GROMACS version: 2021.5 & 2022.3

Hello,

I am performing MD simulations with different time steps, i.e., 1.0 fs, 0.1 fs, and 0.01 fs, using varying electric fields. I got some surprising results with 0.01fs time step. I observed that box volume gradually increased with time when I performed NPT equilibration with 0.01 fs time step. The initial box dimension was 6 nm x 6 nm x 6 nm. After 1nm NPT equilibration, the box dimension became 82 nm x 82 nm x 82 nm.

Could you please suggest me why I got these kind of observation?

Best,

J

I assume you are running in single precision. This is not sufficiently for taking very small time steps like 0.1 and 0.01 fs, in particular when constraints are present. You should use a double precision build of GROMACS.

Dear Professor Hess,

Thank you very much for your reply. I will try to simulate with a double-precision build of GROMACS and update it here.

Best,

J

UPDATE:

Double precision build of GROMACS produced excellent results with 0.01 fs time step.

Hello,

As mentioned in the previous post, I am trying to perform simulations with 0.01, 0.005, and 0.001 fs using varying electric + magnetic fields with double-precision Gromacs 2022.3. Although the box size was not increasing with double precision Gromacs, the conductivity and cumulative current calculated with 0.01fs were not the same or similar to that calculated with 0.005 and 0.001 fs under an extremely high applied magnetic field. The literature and numerical calculations suggest that the 0.01 fs time step is sufficient to simulate an ion-water system in the presence of 5e7 Tesla and 1V applied fields.

Since I use the same system and conditions, I should expect the same conductivity and cumulative current with a 0.01fs, 0.005fs, and 0.001 fs time step.

Could you please suggest to me what could be the problem?

Is the double-precision version of Gromacs perfect to simulate with 0.005 fs and 0.001 fs time steps?

Best,

J

I think you should get the same results. Did you run equally long simulation times?

And do you have constraints that are not using settle?

Dear Professor Hess,

Thank you very much for your reply. Since 0.001 fs is very small and takes longer to simulate, I compared 1.0 ns data of 0.001 fs with 5.0 ns data of 0.01 fs and 0.005 fs. I will compare it with 5.0 ns data. However, there were no fluctuations in the cumulative current of 1.0 ns with 0.001 fs, which was perfectly linear.

I used LINCS as the constraint_algorithm. I will use the SETTLE algorithm.

Best regards,

J

Settle only works for water molecules. If you can use that then that is best. For other bond constraints you should use LINCS. But the results might be sensitive to the LINCS parameters. What LINCS parameters are you using?

Dear Professor Hess,

Thank you very much for the suggestion. I will do it as you suggested.

Please see the details of the LINCS parameter:

; OPTIONS FOR BONDS

constraints = h-bonds

; Type of constraint algorithm

constraint_algorithm = LINCS

; Do not constrain the start configuration

continuation = yes

; Use successive overrelaxation to reduce the number of shake iterations

Shake-SOR = no

; Relative tolerance of shake

shake-tol = 0.0001

; Highest order in the expansion of the constraint coupling matrix

lincs-order = 4

; Number of iterations in the final step of LINCS. 1 is fine for

; normal simulations, but use 2 to conserve energy in NVE runs.

; For energy minimization with constraints it should be 4 to 8.

lincs-iter = 1

; Lincs will write a warning to the stderr if in one step a bond

; rotates over more degrees than

lincs-warnangle = 30

; Convert harmonic bonds to morse potentials

morse = no

Best Regards,

J

What is the energy drift with the different time steps?

Dear Professor,

Thank you so much for your reply and suggestion. Please see all the details.

Looking forward to your reply.

For a system of 1M KCl solution in a box dimension of 4nm x 4 nm x 4nm, all the production MD simulations were performed in an NVT ensemble using a V-rescale thermostat. Systems were equilibrated initially with 1.0 fs time step, then with 0.1 fs time step, and finally with 0.01 fs time step in an NVT ensemble. Thereafter, NPT equilibration with 0.01fs time step. NPT equilibrated system was considered for the final MD simulations. The energy drifts calculated using “gmx energy” are as follows:

## dt=0.01 fs

Energy Average Err.Est. RMSD Tot-Drift

Total Energy -80459.5 – 1100.65 59.2548 (kJ/mol)

## dt=0.005 fs

Energy Average Err.Est. RMSD Tot-Drift

Total Energy -84839.7 27 5168.36 40.4147 (kJ/mol)

## dt=0.001 fs

Energy Average Err.Est. RMSD Tot-Drift

Total Energy -87421.4 71 129193 337.759 (kJ/mol)

The total energy of the system started with an exceptionally very high value, but it reached a value of -85000 kJ/mol to -91000 kJ/mol after 1-3 ps, which is probably the reason for this very drift. Or something wrong?

Please see the starting energy values:

dt=0.01fs

@ s0 legend “Total Energy”

0.000000 2134080.486576

0.001000 -101724.000466

0.002000 -101708.973963

0.003000 -101697.222019

0.004000 -101680.479681

dt=0.005fs

@ s0 legend “Total Energy”

0.000000 8841561.863444

0.000500 -101724.095149

0.001000 -101716.830464

0.001500 -101710.787462

0.002000 -101704.581309

0.002500 -101697.289736

0.003000 -101689.809808

dt=0.001fs

@ s0 legend “Total Energy”

0.000000 223480966.012610

0.000100 -101724.130147

0.000200 -101722.777638

0.000300 -101721.465892

0.000400 -101720.100881

0.000500 -101718.838124

0.000600 -101717.382092

0.000700 -101715.623553

FYI,

nstcalcenergy = 100

nstenergy = 100

Best regards,

J

I meant of the conserved energy quantity. The drift is reported somewhere close to the end to the log file.