Leonard Jones Potential, Leonard Jones Force minimization problem

Question

I am working on a minimization problem.

this is formula for Leonard Jones Potential

this is formula for Leonard Jones Force

r is the distance between two atoms

for our use case, assume value of epsilon and sigma to be 1.

I have to produce the following output

But this is the output I am getting currently

I have gotten the energy (first part) right but I have got minimization (second part) and minimization with jacob = force function (third) part wrong.

I have been provided with the list of 38 points

pos = np.array(
[[9.1912408,8.85347993,2.52084286],
[1.01912764,3.87522545,4.62970322],
[5.24551702,3.08277634,7.05478086],
[7.45460644,1.82528748,9.1244322],
[2.23631807,6.50856483,8.17087062],
[2.38413586,2.45351112,9.09679809],
[3.89190792,7.0473434,8.64435811],
[8.18056541,5.82798785,7.06721024],
[7.77931597,1.1913546,8.189075],
[6.01469145,7.25783933,2.85902436],
[3.20436434,6.12880544,8.40349576],
[1.17308183,5.01219025,5.18666853],
[9.56018292,8.35520658,4.53423182],
[8.96446002,2.54045198,1.18540885],
[6.71330128,1.47359293,8.48820804],
[7.76428251,4.82036586,6.64556474],
[3.51015797,6.52003299,9.52943662],
[6.36423462,8.27639182,3.20621038],
[2.86670404,7.02224667,8.93137555],
[9.66500032,7.80058222,6.18588381],
[0.86644517,6.0966345,4.93413621],
[0.51555606,4.29384818,3.60742217],
[3.86087918,2.44117332,0.94168493],
[9.53757102,4.99027756,9.32440563],
[0.7141156,7.10720338,4.38772771],
[6.92954254,0.05736179,4.72862505],
[3.03575535,7.41135827,7.91044752],
[4.54490467,6.15200622,9.18212839],
[8.17499498,1.06876349,6.98800152],
[8.80925458,1.61273032,1.76078201],
[2.54396165,1.3325331,6.15136613],
[6.22496233,7.9443845,2.08881114],
[1.69322599,6.85444596,9.02678238],
[6.53529894,5.56049505,0.42834861],
[5.41818267,8.17326036,2.80320575],
[6.37606463,3.60233175,3.17719802],
[1.37927686,3.51987394,3.55278354],
[1.32779946,2.79836489,4.6765608],]
)

There are x, y and z component for each point.

I am computing every possible unique distance among these 38 points and send them as an initial guess for minimization.

following is the code that I have implemented so far.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# defined to handle numpy arrays
def LJ_Vectorized(r, sigma=1, epsilon=1) :

  r6 = r**6
  r12 = r6*r6
  s6 = sigma ** 6
  s12 = s6 * s6

  return np.sum(4*epsilon*(s12 / r12 - s6/r6 ))

# defined to handle numpy arrays
def LJ_Force_Vectorized(r, sigma=1, epsilon=1):
  r13 = r ** 13
  r7 = r ** 7
  s6 = sigma ** 6
  s12 = s6 * s6

  return np.sum(48 * epsilon * (s12/r13 - 0.5 * (s6/r7)))

pos = np.array(
[[9.1912408,8.85347993,2.52084286],
[1.01912764,3.87522545,4.62970322],
[5.24551702,3.08277634,7.05478086],
[7.45460644,1.82528748,9.1244322],
[2.23631807,6.50856483,8.17087062],
[2.38413586,2.45351112,9.09679809],
[3.89190792,7.0473434,8.64435811],
[8.18056541,5.82798785,7.06721024],
[7.77931597,1.1913546,8.189075],
[6.01469145,7.25783933,2.85902436],
[3.20436434,6.12880544,8.40349576],
[1.17308183,5.01219025,5.18666853],
[9.56018292,8.35520658,4.53423182],
[8.96446002,2.54045198,1.18540885],
[6.71330128,1.47359293,8.48820804],
[7.76428251,4.82036586,6.64556474],
[3.51015797,6.52003299,9.52943662],
[6.36423462,8.27639182,3.20621038],
[2.86670404,7.02224667,8.93137555],
[9.66500032,7.80058222,6.18588381],
[0.86644517,6.0966345,4.93413621],
[0.51555606,4.29384818,3.60742217],
[3.86087918,2.44117332,0.94168493],
[9.53757102,4.99027756,9.32440563],
[0.7141156,7.10720338,4.38772771],
[6.92954254,0.05736179,4.72862505],
[3.03575535,7.41135827,7.91044752],
[4.54490467,6.15200622,9.18212839],
[8.17499498,1.06876349,6.98800152],
[8.80925458,1.61273032,1.76078201],
[2.54396165,1.3325331,6.15136613],
[6.22496233,7.9443845,2.08881114],
[1.69322599,6.85444596,9.02678238],
[6.53529894,5.56049505,0.42834861],
[5.41818267,8.17326036,2.80320575],
[6.37606463,3.60233175,3.17719802],
[1.37927686,3.51987394,3.55278354],
[1.32779946,2.79836489,4.6765608],]
)

# compute distances between 38 atoms and store unique distances in a numpy array

import itertools

distances = []
for pair in itertools.combinations(pos, 2):
    p1, p2 = pair
    distance = np.linalg.norm(p1 - p2)
    distances.append(distance)

# Remove duplicates to get unique distances
distances = np.array(list(set(distances)))

print('Energy ', LJ_Vectorized(distances))
res = minimize(LJ_Vectorized, distances , method='CG' , tol =1e-3)
print( 'After minimization' , res.fun )
res = minimize(LJ_Vectorized , distances , jac=LJ_Force_Vectorized , method='CG' , tol =1e-3)
print('After minimization with force' , res.fun )