Issue of High Potential Energy on Running Minimization

User discussions
, , , ,
21

Dear milosz,

I sent the files to your email yesterday.

Thanks

22

Got them, I’m just working on some minor details regarding the formatting of OPLS topologies.

23

So I did the following test to make sure all your periodic bonds are in place:

p = gml.Pdb('FER.pdb')

cutoff = 1.8 # in angstroms
pairs = [] # all pairs going through the PBC

for a1 in p.atoms:
  for a2 in p.atoms[a1.serial:]:
    if p._atoms_dist(a1, a2) > min(p.box[:3]) and p._atoms_dist_pbc(a1, a2) < cutoff:
      pairs.append((a1.serial, a2.serial))

t = gml.Top('FER.top', ignore_ifdef=True)
zeo = t.molecules[0]
for atom_pair in pairs:
  print(atom_pair in zeo.bonds)

and all return True, which means your topology is formed correctly.

Running a quick minimization with FER.top and FER.pdb also goes well. I do

gml.calc_gmx_energy('minimized.gro', 'FER.top', terms='all')

to check the energy contributions, and they come out as

{'bond': [15956.757812], 'angle': [24865.548828], 'proper-dih.': [3996.352051], 'lj-(sr)': [0.0], 'coulomb-(sr)': [7405572.5], 'potential': [7450391.0], 'box-x': [5.7054], 'box-y': [4.2909], 'box-z': [4.5246], 'volume': [110.768105], 'density': [1751.036743]}

In short, it all looks more or less fine. Since the topology doesn’t use Coulomb attenuation of 1-4 pairs (fudgeQQ set to 0 and no [ pairs ], the magnitude of Coulombic repulsion is large, but you should check that against the original parameterization