Too high temperature after nvt

User discussions
1

GROMACS version: 2023.1
GROMACS modification: Yes/No
Hello,
After having built a system with 36 1-octanol molecules at the centre surrounded by water solvents (spc216.gro) and at total dimensions of a cube with length 50 A. I run energy minimization and now i tried to run nvt but the average temperature is something like 10000 K (too high). I attach nvt.mdp file, em potential energy file .xvg and energy minimization file (em_2.mdp)
em_2.mdp (1.0 KB)
I also attach the temperature.xvg file
temperature.xvg (36.4 KB)

. Can someone please tell me what i did wrong?
thank you
em_potential.xvg (1.3 KB)

nvt.mdp (2.3 KB)

2

Hi,

Your timestep is just 0.02 fs. Try increasing it to 2 fs.

Petter

3

Hello Peter and thanks for the response. When the timestep was at 2fs i got a segmentation fault because two water molecules were too close

4

*Update: I visualized it and after the energy minimization proccess i see that the octanol molecules’ bonds were disappeared, why is that? I attach the em_2.mdp file.
em_2.mdp (1.0 KB)

5

Then your energy minimization likely did not work or complete properly. Your system is highly unstable which would explain the behavior.

How did you create your system and what force field are you using? Some files that would be useful to see are your configuration (.pdb/.gro), topology (.top) and log files (.log).

It would also be good if you can post your exact commands that were used for both steps (for gmx grompp and mdrun).

6

I have downloaded from automatic topology builder the .itp file for octanol, the .pdb file for octanol and the modified force field gromos54a7 and i also did some adjustments (added to file ffbonded.itp at the last line #include /path/to/octanol.itp and copied at the folder of the force field the octanol.itp file).
I converted the .pdb file with gmx pdb2gmx -f octanol.pdb -o octanol.gro
Then i made a small box
$gmx editconf -f octanol.gro -o octanol_box.gro -c -d 0.6 -bt cubic
I inserted 35 more octanol molecules at this box
$gmx insert-molecules -f octanol_box.gro -ci octanol.gro -o octanol_box_36_octanols.gro -nmol 35 -try 10000 -seed 12345
I made a larger box and solvated it
$gmx solvate -cp larger_boc_36_octanols.gro -cs spc216.gro -o solvated.gro -p topol.top

Then i used the grommacs preproccess
$gmx grompp -f em_2.mdp -c input.gro -p topol.top -o em.tpr -maxwarn 2
$gmx mdrun -s em.tpr -o em.trr -c em.gro

Thats what i get after many tries:
EO:9_6_2023 $ gmx grompp -f em_2.mdp -c octanol_36_box_0.8_solvated.gro -p topol.top -o em_4.tpr -maxwarn 2
:-) GROMACS - gmx grompp, 2023 (-:

Executable: /run/media/EO/data_2TB/Gromacs_2023/Installation_folder/gromacs/bin/gmx
Data prefix: /run/media/EO/data_2TB/Gromacs_2023/Installation_folder/gromacs
Working dir: /mnt/data_2TB/Gromacs_2023/9_6_2023
Command line:
gmx grompp -f em_2.mdp -c octanol_36_box_0.8_solvated.gro -p topol.top -o em_4.tpr -maxwarn 2

Ignoring obsolete mdp entry ‘ns_type’

NOTE 1 [file em_2.mdp]:
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
your simulation.

Setting the LD random seed to -588402949

Generated 210 of the 2556 non-bonded parameter combinations

Excluding 3 bonded neighbours molecule type ‘Other’

turning all bonds into constraints…

Excluding 2 bonded neighbours molecule type ‘SOL’

turning all bonds into constraints…

WARNING 1 [file topol.top, line 1063]:
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. Further information is available at
GROMACS can not reproduce properties with the GROMOS force fields - Redmine #2884 (#2884) · Issues · GROMACS / GROMACS · GitLab, and a longer
explanation of our decision to remove physically incorrect algorithms can
be found at On The Importance of Accurate Algorithms for Reliable Molecular Dynamics Simulations | Theoretical and Computational Chemistry | ChemRxiv | Cambridge Open Engage .

NOTE 2 [file topol.top, line 1063]:
In moleculetype ‘Other’ 972 atoms are not bound by a potential or
constraint to any other atom in the same moleculetype. Although
technically this might not cause issues in a simulation, this often means
that the user forgot to add a bond/potential/constraint or put multiple
molecules in the same moleculetype definition by mistake. Run with -v to
get information for each atom.

Analysing residue names:
There are: 36 Other residues
There are: 3988 Water residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups…
Number of degrees of freedom in T-Coupling group rest is 26841.00
The integrator does not provide a ensemble temperature, there is no system ensemble temperature

The largest distance between excluded atoms is 0.164 nm between atom 5765 and 5766
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 44x44x44, spacing 0.116 0.116 0.116

Estimate for the relative computational load of the PME mesh part: 0.22

This run will generate roughly 5 Mb of data

There were 2 NOTEs

There was 1 WARNING

Back Off! I just backed up em_4.tpr to ./#em_4.tpr.1#

7

Can you upload your .log files somewhere too? Eg. using pastebin.com or whatnot.

8

you mean these files?
em.log (23.9 KB)
em_2.log (19.0 KB)
npt.log (639.8 KB)
nvt.log (639.9 KB)

9

At the end of the em_2.log files it says:
Maximum force = 7.5985640e+07 on atom 784

So there is a very (too) large force acting on atom 784. Check in the structure was is going on there.

10

Thank you for replying!
Maximum force on atom 784 is too large, but the mean energy of all atoms belongs in the same order of magnitude, as a result the whole system is too unstable

11

I don’t understand what you mean with “in the same order of magnitude”. The same as what?

You need to have better initial structure to start MD from.

12

Sorry, i mean in the same order of magnitude as the other atoms. I tried to build my system with the above gromacs commands and i cannot understand why it is so unstable.

13

The forces are extremely large, so there must be overlapping atoms. Have you checked if atom 784 overlaps with another atom in the starting structure for EM?

14

**Update: I visualized it and after the energy minimization proccess i see that the octanol molecules’ bonds were disappeared, why is that? Even when i am not doing energy minimization with steepest descent (i used conjugate gradient).