Relative constraint deviation after LINCS

Hello,

I’m using the eqA.mdp file from the tutorial here:

Despite running several equilibration steps with and without restraints on my ligand I get the error:
Step 3, time 0.006 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.001067, max 0.057983 (between atoms 9731 and 9732)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length

I can see the dummy atoms swinging around like crazy when I visualize the equilibration trajectory. I don’t know why the equilibration doesn’t crash but the run started from the provided mdp file does.

Apparently an mdp setting I need to adjust?

Or is it an issue with my ligand structure? I can’t see why this would be the case since this is a structure I’ve used in many simulations (without the dummy atoms).

Thank you,

Dan

I switched the constraints from “all-bonds” to “h-bonds” and it eliminated the errors.

Switching to h-bond constraints is likely incorrect for the pmx based ligand hybrid topologies. In case, a hydrogen atom is mapped to a dummy in stateA, the constraint will not be generated for the bonds involving this dummy.

I would suggest trying an SD integrator (note that you need to change tau-t in this case).

Okay, thank you. In the case of an SD integrator, I would return constraints = all-bonds?

Thanks, Vytas.

Dan

Yes. It is not guaranteed that the SD integrator will help, as there might be other issues, but it is worth to try it out.

Is this simulation started right after energy minimization? If so, what is the maximal force at the end of the EM?

The first time around I ran nvt and npt equilibration.

This time I went straight from em to the included eq mdp. The maximum force was 2.203x10^3.

Could I go from the eq run with md integrator and my h-bond constraint to the TI runs with the sd integrator all-bond constraint without issue?

I would recommend to retain consistent use of constraints and integrator for both steps: eq and transitions.