Mdrun Module of nvt Wrote Unusual pdb Files

User discussions
1

GROMACS version: 2023.1
GROMACS modification: Yes/No

Hello.

I am performing MD simulation of a pure 500 molecules of ionic liquid. The pdb files of the cation and anion are written separately. Packmol is used to packed 500 cations and anions together in a single pdb file. The packmol pdb file was changed to .gro and then inserted into a cubic box of 1nm. Initial pdb files of cation and anion were uploaded to ATB server to obtained their all-atoms itp files which was included in the topology file, together with gromos 54a7_atb.ff.
Energy minimization was achieved but with unusual nature of potential graph. Md run of nvt do wrote unusual pdb files.The promped statements from the end of energy minimization to nvt equilibriation, together with my files are provided below:

At the end of energy minimization step
Steepest Descents converged to Fmax < 1000 in 1729 steps
Potential Energy = -4.7715419e+05
Maximum force = 9.7373108e+02 on atom 19615
Norm of force = 7.8541580e+01

Files
topol.top (338 Bytes)
nvt.mdp (2.2 KB)
em.mdp (1.0 KB)

mdrun of nvt prompt this

gmx mdrun -v -deffnm nvt

Reading file nvt.tpr, VERSION 2023.1 (single precision)
Changing nstlist from 10 to 100, rlist from 1 to 1.029

Using 1 MPI thread
Using 4 OpenMP threads

starting mdrun ‘pure ionic liquid’
50000 steps, 100.0 ps.
step 0
Step 2, time 0.004 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000110, max 0.010312 (between atoms 10446 and 10447)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
10446 10447 38.0 0.1121 0.1132 0.1120
10446 10448 39.3 0.1121 0.1131 0.1120

Step 3, time 0.006 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000028, max 0.002499 (between atoms 10446 and 10447)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
10446 10447 61.2 0.1132 0.1123 0.1120
10446 10448 65.7 0.1131 0.1123 0.1120

Step 4, time 0.008 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.038587, max 4.566490 (between atoms 10446 and 10448)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
25127 25128 57.7 0.1121 0.1106 0.1120
25127 25129 51.9 0.1121 0.1104 0.1120
29654 29656 32.6 0.1090 0.1078 0.1090
29654 29657 90.0 0.1090 0.1454 0.1090
10446 10447 90.0 0.1123 0.4301 0.1120
10446 10448 90.0 0.1123 0.6234 0.1120
Wrote pdb files with previous and current coordinates

Step 5, time 0.01 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1.818926, max 252.313431 (between atoms 10446 and 10447)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
25127 25128 90.0 0.1106 2.0233 0.1120
25127 25129 90.0 0.1104 0.2656 0.1120
29654 29656 30.1 0.1078 0.1117 0.1090
29654 29657 81.5 0.1454 0.1087 0.1090
10446 10447 90.0 0.4301 28.3711 0.1120
10446 10448 90.0 0.6234 1.1008 0.1120
10485 10486 90.0 0.1120 0.8523 0.1120
10485 10487 90.0 0.1120 2.3489 0.1120
10488 10489 90.0 0.1120 0.7347 0.1120
10488 10490 90.0 0.1120 0.4482 0.1120
Wrote pdb files with previous and current coordinates

Step 6, time 0.012 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 2995.472168, max 405171.218750 (between atoms 29199 and 29202)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
25127 25128 90.0 2.0233 33.0502 0.1120
25127 25129 90.0 0.2656 2232.0037 0.1120
25150 25152 90.0 0.1130 0.9243 0.1130
29199 29200 90.0 0.1090 1482.6108 0.1090
29199 29201 90.0 0.1090 1487.4843 0.1090
29199 29202 90.0 0.1090 44163.7695 0.1090
29654 29657 90.0 0.1087 0.1187 0.1090
10446 10447 90.0 28.3711 79.4997 0.1120
10446 10448 90.0 1.1008 8137.3271 0.1120
10443 10444 90.0 0.1120 1.3141 0.1120
10443 10445 90.0 0.1121 1.7986 0.1120
10442 10443 90.0 0.1120 1.5135 0.1120
8985 8986 90.0 0.1130 7793.2314 0.1130
8985 8987 90.0 0.1130 192.4110 0.1130
10449 10450 90.0 0.1120 2.6645 0.1120
10449 10451 90.0 0.1120 1.8408 0.1120
10482 10483 90.0 0.1130 1.7628 0.1130
10482 10484 90.0 0.1130 21.9267 0.1130
10485 10486 90.0 0.8523 42.7724 0.1120
10485 10487 90.0 2.3489 7.0602 0.1120
10472 10473 41.4 0.1120 0.1120 0.1120
10491 10492 90.0 0.1120 0.5516 0.1120
10491 10493 90.0 0.1120 0.3977 0.1120
10491 10494 90.0 0.1120 0.4314 0.1120
10488 10489 90.0 0.7347 29.3122 0.1120
10488 10490 90.0 0.4482 26.6822 0.1120
Wrote pdb files with previous and current coordinates
Segmentation fault (core dumped)

The mdrun of nvt module wrote step 4b and 4c pdb files up to step 6b and 6c pdb files only.

Please, how can I rectify the system?
Thank You.

2

Hi,

you got LINCS warning from the very early stage of the simulations. That might indicate that your starting structure is not properly energy minimised for that force field or something else is wrong is your starting structure.

You can visualize the structure to see what is the issue with atom id 10446 and 10447.
Once an LINCS occurs, the application mdrun writes a pdb file label step_ID.pdb
I hope it helps
Alessandra

3

Hello.
Thanks a lot for the guidance.

4

Hi @avilla @jalemkul
I could finish the MD simulation successfully. I see no discrepancy in the MD movie. Whole structure remains intact through out the simulation.
But, I have also got two pdbs “step11c_n4.pdb” and “step11b_n4.pdb” generated after Energy Minimization step.

Note: Those pdbs were not generated when I used dodecahedron type box. (same structure, same mdp file, same system configuration). but I was facing periodic boundary condition violation, so I switched box type to cubic.

So should I consider my simulation successful or it has failed in the very beginning stage.
Please find the terminal output (terminal_output.dat) as attached. (open in notepad)
terminal_output.dat (531.6 KB)
Also let me know if you need any additional info to look into my problem.

Thanks in advance for considering my query.

5

LINCS warning (and associated step*.pdb files) are okay during the energy minimization (EM). It means your initial structure was far from ideal, but resolving this is why we do EM. Out of an abundance of caution, you can check around the involved atoms the final geometry, but your simulation should be fine, barring other possible issues of course.

On the other hand: one, or worst several, LINCS warnings during the simulation (whether equilibration or production) are concerning - see above messages for explanation.

6

@ebriand Thankyou so much for replying.

I didn’t get any LINCS warning during Minimization step. At Step 11 of minimization, It just shows

"step 11: 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
Step=   11, Dmax= 6.2e-02 nm, Epot=          inf Fmax= 1.08415e+05, atom= 916"

and wrote two Pdb files named “step11c_n4.pdb ” and “step11b_n4.pdb ”
What does above error mean, how to solve that as it asked to check for bad contacts/reduce time step?

Also, as I have mentioned, with the same Pdb and MDPs but the dodecahedron type of box, it didn’t produce any additional files (step.pdb). Please consider this also.

Thanks again

7

@jalemkul @avilla
I tried minimization with smaller timestep dt=0.001 (previously dt=0.002), it didn’t produce any additional pdbs. Though it is not reproducible all the time (50:50).
Upon investigating a little, I have found that the smaller dt provides a comparative accurate and controlled simulation but with the increased time of simulation.

So it is valid to have different timesteps for different process in a single MD simulation? i.e. smaller time step (dt=0.001) for minimization and the equilibration as they take little time to finish, and large time step (dt=0.002) for production?

Thanks

8

One or more water molecules can not be settled means that is is also a constraint error - while not a LINCS warning, because a slightly different constraint solver is used, it is the same kind of error.

It is important to note that energy minimization and equilibration are qualitatively different steps in the protocol. In energy minimization (where integrator = steep), the gradient descent algorithm optimizes the particle position to get a lower energy conformation. There is no concept of time, nor timestep, since this is not MD proper.

On the other hand, equilibration and production are actual simulations, where particles moves according to the law of mechanics (integrator = md), but, in the case of equilibration, typically with restraints on their movement.

This means that during energy minimization (where integrator = steep), the dt parameter does not affect anything. The relevant parameters are as follow:

integrator  = steep     ; Algorithm (steep = steepest descent minimization)
emtol       = 300.0     ; Stop minimization when the maximum force < 300.0 kJ/mol/nm
emstep      = 0.01      ; Energy step size
nsteps      = 50000     ; max steps

You might want to decrease your emstep.

9

@ebriand Thanks for your kind explanation.

Still I am wondering that why with dodecahedron type box it didn’t produce any additional step.pdb while cubic type of box produces. Is it related to box type?

Thanks

10

It is not directly related to the box type, however the positions of the water created by gmx solvate depends on the shape or size of the box (its is copy/pasting a preset small water box until the box is filled - simple but reliable).

Some of the water might end up placed in relatively bad position, leading to a few warnings at the beginning of the energy minimization (which, again, are fine - whereas warnings during actual simulation would generally not be).

11

@ebriand Thank you so much. I really appreciate your assistance in addressing my inquiry.