Constraint problems in GROMACS with CP2K QMMM interface

Hi all,

I have encountered problems with the interface between GROMACS and CP2K. I setup my calculation standardwise for GROMACS and did the whole standard procedure, minimization, NVT equilibration and NPT equilibration. After that I want to do QM/MM dynamics with NPT-ensemble for my system. The force fiield I used was OPLS, and the water model was tip4p. The time step was 2 fs.

The first problem is related to the generation of tpr files. As you can see, if you put in the bond parameters,

; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

Then you will get the following warning:

“WARNING 1 [file md.mdp]:
Your QM subsystem has a lot of constrained bonds. They probably have been
generated automatically. That could produce an artifacts in the
simulation. Consider constraints = none in the mdp file.”

If you select none instead of hbonds in constraints, then the calculation will fail after several step. In this case after 6 steps.

“step 6: 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”

When you make the time-step smaller, then the calculation will run a bit longer but also fail after some steps.

“step 156: 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”

I tried to redo the simulation with AMBER force field. The same error occured after 345 steps.

My question is:

  1. will the constraints be applied on the QM subsystem, so that I can ignore the warning?

I try to modify the LJ-Parameter as the CP2K QMMM tutorial suggested.

With OPLSFF the calculation just failed because the potential energy is just too high.
“Step 0: The total potential energy is 1.33998e+17, which is extremely high.
The LJ and electrostatic contributions to the energy are 1.33998e+17 and
-107853, respectively. A very high potential energy can be caused by
overlapping interactions in bonded interactions or very large coordinate
values. Usually this is caused by a badly- or non-equilibrated initial
configuration, incorrect interactions or parameters in the topology.”

With AMBERFF, the calculation fails at the first step, and gave the same mistake “step 0: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.”

So basically I think the problem with water settling is not related with the parameters itself but rather with the constraint.

Unfortunately, as a new user, I cannot attach my input files here.

Best regards

Hi thanks for the interest in our new interface,

  1. I found a major problem in your index file:

[ QMAtoms ]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
46 47
46 47

  1. 2 fs time-step is too big for QM/MM in general (I suggest to use 1 fs)

  2. Also you are trying xTB for the QM system, I’m not sure that it works well now in CP2K. At least what I see that it starts failing to converge at some point. Also it could not use periodicity and GEEP electrostatics.

I have checked you system with DFT(PBE functional) after fixing index.ndx file and reducing time-step to 1 fs and seems like it works fine for me.

I try to modify the LJ-Parameter as the CP2K QMMM tutorial suggested.

LJ is managed by GROMACS in our interface, so that modifications are not needed

  1. will the constraints be applied on the QM subsystem, so that I can ignore the warning?\

Regarding this question - yes they will be applied, that’s why one should not use automatic constraints generation with QMMM (which is also stated in the warning you see). This is not the case however for water molecules, as they typically constrained with SETTLE. If you put water molecule into QM part, then SETTLE constraint will be removed from that water by the interface.

Also be careful of using Pressure coupling with QMMM as it currently accounts only for pressure due to MM-MM atoms interactions.

Hi dmorozov,

thanks for your reply. I accidentally gave you the wrong index.ndx file, sorry for that. The correct one should not have duplicate atomic number 46 47. In any case, I guess GROMMP will fail if the index file is wrong.

Okay, I did not know that xTB cannot use PBC, and yes it fails to converge after some points. However if you use constraints, then the calculation ran just fine, even with 2fs. So I thought this is more related to the constraints rather than periodicity. Then I guess, I need to modify the topology and use [settles] for the water molecules. About GEEP, yes I knew only Coulomb can be used.

That is good to hear that the LJ is already managed in the interface.

So if I understood correctly the pressure coupling is not working really well with QM subunit. Is this the reason, why the QMMM tutorial for CP2K also suggested NVT instead of NPT?

Many thanks

Okay, I did not know that xTB cannot use PBC, and yes it fails to converge after some points. However if you use constraints, then the calculation ran just fine, even with 2fs. So I thought this is more related to the constraints rather than periodicity. Then I guess, I need to modify the topology and use [settles] for the water molecules. About GEEP, yes I knew only Coulomb can be used.

I’m not really an expert in xTB so I could not explain why it fails without constraints, sorry.
SETTLE for water molecules are used automatically, regardless of constraints option (water molecules are always constrained by default in GROMACS).

So if I understood correctly the pressure coupling is not working really well with QM subunit. Is this the reason, why the QMMM tutorial for CP2K also suggested NVT instead of NPT?

Yes, NPT formulation for QMMM is very ad-hoc, to calculate pressure one need to evaluate force virial, which is essentially depends on pair-wise interactions between particles in your system. However there is now way to decompose QMMM forces into individual contributions between pairs of atoms. So basically in QMMM pressure is coupled and calculated only for MM part.