Simulation unstable with external potential

GROMACS version: 2021.3
GROMACS modification: Yes

Hi,

I added an extenal potential to the force field but it seems that the simulation is not stable. Since the h-bonds are constrained but the external potential have additional forces on hydrogens, there are some LINCS warnings:

Step 140957, time 281.914 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 3.462555, max 16.844141 (between atoms 18 and 19)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
9 10 90.0 0.1647 0.2452 0.1080
18 19 35.4 0.8695 1.9825 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
18 20 48.7 0.1051 0.1155 0.1111
26 27 90.0 1.0481 0.2800 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111
26 28 90.0 0.1168 0.1443 0.1111

Then the simluation is stopped:

step 140957: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
/var/slurm/slurmd.spool/job1834770/slurm_script: line 12: 18320 Segmentation fault

I tried smaller integration step (1 fs) but it still has the same problem. Does anyone know how to run the stable simulation with h-bond constraints and external potentials?

Thanks in advance,
Pan

What do you mean with “external potential”? An electric field?

My guess would be that magnitude of the external potential is one or more orders of magnitude too large.

The external potential is a machine learning potential. The magnitude is usually between -1 kJ/mol and -10^3 kJ/mol, which should be 2 or more orders of magnitude smaller than the potential energy.

So you changed the GROMACS code to add you own potential? If that is the case then very likely there is an issue in your code.

Yes, the GROMACS code was changed. I didn’t find the bug in my code and the force calculation is also okay according to the numerical differentiation test.

However, I have a question about calculating MM single point energy. Since my machine learning model is training on the energy difference between QM and MM energies, I calculated MM energy using the mdp file below:

integrator = md
nsteps = 0
dt = 0.001
nstxout = 1
nstvout = 1
nstfout = 1
nstenergy = 1
nstlog = 1
nstxout-compressed = 1
compressed-x-grps = System
continuation = no
constraints = none
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rcoulomb = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rlist = 1.2
rvdw = 1.2
rvdw-switch = 1.0
DispCorr = no
coulombtype = cutoff
coulomb-modifier = none
tcoupl = V-rescale
tc-grps = Protein
tau_t = 0.1
ref_t = 300
pbc = xyz
gen_vel = yes
gen_temp = 300

Is this okay to calculate MM energy of a molecule in vaccum? I didn’t use constraint because I found the energy will be different if the h-bond constraint is used, and the MM energy without constraint seems to be the same with the MM energy calculated from other program like NAMD after comparison. However, the production simulation is running with h-bond constraint, and machine learning energy (difference between QM energy and MM energy without constraint) is added. Will this cause the LINCS warnings of the bonds rotating more than 30 degrees and affect the stability of simulations?