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: 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?
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.
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:
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?