GROMACS version: 2023.3
GROMACS modification: No
Here post your question
My system is composed by ionic surfactants and their counterions soaked in water to form wormlike micelles. I’ve been looking for the appropriate interaction parameters, specifically the cutoff distances for the relevant type of interactions within the system (repulsion, dispersion and electrostatics). For repulsion and dispersion I formerly used a plain cutoff scheme in direct space, whilst PME for electrostatics. I noticed that for VdW interactions (repulsion and dispersion) increasing rvdw in a plain cutoff scheme affected the energy results, specifically the average pressure throughout the NVT simulation. Thus, I found that for non-homogeneous systems outside the cutoff distance, the GROMACS manual suggests to apply PME for VdW interactions as well and divide short and long parts in direct and reciprocal spaces respectively. This resulted in a different average pressure output in terms of labeling an data.
The former simulation’s relevant input and output is:
vdw-type = Cut-off
rvdw = 1.4
DispCorr = EnerPres
Energy Average Err.Est. RMSD Tot-Drift
Pressure 8.5272 1.5 207.754 7.70906 (bar)
While the later simulation’s relevant input and output is:
vdw-type = PME
rvdw = 1.4
DispCorr = EnerPres
fourierspacing = 0.12
pme-order = 4
ewald-rtol-lj = 0.001
lj-pme-comb-rule = geometric
Energy Average Err.Est. RMSD Tot-Drift
Pres. DC 1.53899e-06 9.2e-10 2.77028e-08 -5.49326e-09 (bar)
The denomination of “Pressure” changed to “Pres. DC” and the results are completely different. Does anyone know why is this happening?
It might be important to add that GROMACS delivered the following note for the mismatch between PME and the Force Field used (GROMOS from ATB server) and for the combination rule (geometric), while generating the .tpr executable:
“You are using geometric combination rules in LJ-PME, but your non-bonded
C6 parameters do not follow these rules. This will introduce very small
errors in the forces and energies in your simulations. Dispersion
correction will correct total energy and/or pressure for isotropic
systems, but not forces or surface tensions.”
You are right, I just typed the same entry (17) as I used to. But now with reciprocal VdW the entry “11 LJ-recip” was added. Therefore “Pressure” is just simply the number 18. I thought pressure was being miscalculated due to an improper Force-Field. Thanks!
Now the only issue arising would be whether or no the dispersion correction in the pressure is too large, what might turn unrealistic, deriving in the need of Force-Field parametrization rather than the use of corrections, and thus need to rely back in cutoff distances. Also if it is feasible to use the geometric combination rule rather than the Lorentz-Berthelot for performance improvement. Few computer-consuming run stages missing to deduce both.
If you or anyone have experience on this, any suggestions are welcome.
With LJ-PME the dispersion correction only corrects for the difference between the LB and geometric rule when using the geometric rule. This difference will be very small. I don’t think there are reasons not to use the geometric rule for the grid part.
Several notes coming from GROMACS executions and instructions made me cautious to not blindly rely on the GROMOS54A7 Force-Field at issue. I am aware that depending on the FF and the source (in this case the ATb server), small or large deviations in thermodynamic, structural and dynamic output may arise.
Although the .itp files of the molecules were recently generated (they used quantum mechanical optimization), and parameter files were updated and replaced for those available in the ATb server; the need of validation tests for specific applications was coined by themselves. Particularly, the GROMACS execution note I quoted at the end of the discussion opening mentions that the C6 (V) parameter of the geometric rule (that would change to σ for the LB rule) do not follow the geometric combination. Regarding deviations due to the dispersion correction, I quote the following GROMACS execution note:
“The GROMOS force fields have been parametrized with a physically
incorrect multiple-time-stepping scheme for a twin-range cut-off. When
used with a single-range cut-off (or a correct Trotter
multiple-time-stepping scheme), physical properties, such as the density,
might differ from the intended values. Since there are researchers
actively working on validating GROMOS with modern integrators we have not
yet removed the GROMOS force fields, but you should be aware of these
issues and check if molecules in your system are affected before
proceeding.”
So far I have obtained the behaviour I was expecting, but I will try to validate the accuracy of the output once establishing the most suitable input to reproduce the physical trajectory that best fits to the parameters provided by the ATb server, but which might not be necesarily correct.
The GROMOS force fields are somewhat special in that many pairs do not follow combination rules; they use a scaling factor. The dispersion correction will correct for this when using the geometric rule for LJ-PME. Note that using LB rules will only make things worse, as GROMOS follows geometric rules for many pairs. I don’t know how large the effects of geometric rule assumption in LJ-PME are with the GROMOS force field.
The note on the force field is mainy in reaction to papers published by some members of the GROMOS community that claim that GROMACS implements the GROMOS force field incorrectly. The GROMOS force fiels is, or at least has been, used extensively with GROMACS and I have never seen issues in results.