# Confusing Density from NPT run **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, mdrun\ **Created:** [August 28, 2024, 7:05am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021 "2024-08-28T07:05:33Z")\ **Posts on this page:** 9\ **Page:** 1
**Author:** ![hehehaha](https://avatars.discourse-cdn.com/v4/letter/h/97f17d/32.png) [@hehehaha](https://gromacs.bioexcel.eu/u/hehehaha)\ **Post date:** [August 28, 2024, 7:05am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/1 "2024-08-28T07:05:33Z")
GROMACS version:2023 (gromacs-amd-gfx90a/2023) GROMACS modification: Yes/No Hi all, I am running MD on a protein-ligand complex. My protein is huge with 6 chains and 277.25 kDa. My forcefield parameterisation was using AmberForce field FF19SB / ambertools. I have decided to use OPC water model to solvate my protein-ligand complex I decided to follow several protocols about my protein (from literature) and decided to do EM (steep descent) → restrained NPT at 20K → restrained NPT at 300K → Unrestrained NPT at 300K → MD [for restraints, i just standardise to 1000 kj/(molnm) (i might improve on this workflow especially with reducing force constants but i have a more important question below). i decided to start low temperature and increase to high temperature because initially, i tried NVT-\>NPT (restrained) at 300K and got densities at ~300kg/m^3 which is far from experimental water or OPC in this case. Since I believe my protein looks stable, i can skip the NVT equilibration and carry out NPT at lower temperatures in hopes that my system would contract and then reach desired density. Not sure if this is the right protocol (EM → Restrained NVT lower temp → Restrained NVT desired temp → Unrestrained NVT desired temp) . Any advice would be appreciated. For all restrained NPT runs, density is kept around 300kg/m^3 and finally at the unrestrained run, density shot to about 1300 kg/m^3. This is far from the water’s density but i would say relatively close to average density reported for proteins. Now i wonder if the density actually reflected the large protein’s density instead of the water or if that is incorrect. I wonder if 1) OPC water model topology was incorrect? 2) Did I not use enough water molcules to solvate the system? (but the water molecules covers the entire protein well and no vacuum bubbles during NPT) Side note: When removing PBC and visualisating the trajectory, the ligand seems stable even in unrestrained runs and seems to be fine so im believe my system is stable but the density is throwing me a little of a little Any advice/help is appreciated. Thanks all OPC water model topology in case anyone’s wants to see it found in [gromacs\_ff/amber14sb\_parmbsc1\_opc\_lmi.ff/opc.itp at master · intbio/gromacs\_ff · GitHub](https://github.com/intbio/gromacs_ff/blob/master/amber14sb_parmbsc1_opc_lmi.ff/opc.itp) [moleculetype] ; molname nrexcl ; OPC model WAT 2 [atoms] ; id at type res nr res name at name cg nr charge mass 1 O 1 SOL O 1 0 16.00000 2 H 1 SOL H1 1 0.67914 1.00800 3 H 1 SOL H2 1 0.67914 1.00800 4 EPW 1 SOL EPW 1 -1.35828 0.00000 #ifndef FLEXIBLE [settles] ; i funct doh dhh 1 1 0.08724 0.13712 #else [bonds] ; i j funct length force.c. 1 2 1 0.08724 502416.0 0.08724 502416.0 1 3 1 0.08724 502416.0 0.08724 502416.0 [angles] ; i j k funct angle force.c. 2 1 3 1 103.6 628.02 103.6 628.02 #endif [virtual\_sites3] ; Vsite from funct a b 4 1 2 3 1 0.147722363 0.147722363 [exclusions] 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 ; The position of the virtual site is computed as follows: ; ; O ; ; V ; ; H H ; ; Ewald OPC: ; const = distance (OV) / [cos (angle(VOH)) \* distance (OH)] ; 0.01594 nm / [cos (51.8 deg) \* 0.0872433 nm] ; then a = b = 0.5 \* const = 0.147722363 ; ; Vsite pos x4 = x1 + a\*(x2-x1) + b\*(x3-x1) npt.mdp : title = Protein-ligand complex NPT equilibration define = -DPOSRES ; position restrain the protein and ligand ; Run parameters integrator = md ; leap-frog integrator nsteps = 500000 ; 2 \* 500000 = 1000 ps = 1ns dt = 0.002 ; 2 fs ; Output control nstenergy = 5000 ; save energies every 10.0 ps nstlog = 5000 ; update log file every 10.0 ps nstxout-compressed = 5000 ; save coordinates every 10.0 ps ; Bond parameters continuation = no ; first NPT run constraint\_algorithm = lincs ; holonomic constraints constraints = h-bonds ; bonds to H are constrained lincs\_iter = 1 ; accuracy of LINCS lincs\_order = 4 ; also related to accuracy ; Neighbor searching and vdW cutoff-scheme = Verlet ns\_type = grid ; search neighboring grid cells nstlist = 20 ; largely irrelevant with Verlet rlist = 1.0 vdwtype = Cut-off vdw-modifier = Potential-shift rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.0 pme\_order = 4 ; cubic interpolation fourierspacing = 0.125 ; grid spacing for FFT ; Temperature coupling tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = protein\_grp non\_protein\_grp ; two coupling groups - more accurate tau\_t = 0.5 0.5 ; time constant, in ps ref\_t = 50 50 ; reference temperature, one for each group, in K ; Pressure coupling pcoupl = C-rescale ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau\_p = 0.5 ; time constant, in ps ref\_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 refcoord\_scaling = com ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction used with OPC water model DispCorr = EnerPres ; Velocity generation gen\_vel = yes ; assign velocities from Maxwell distribution gen\_temp = 50 ; temperature for Maxwell distribution gen\_seed = -1 ; generate a random seed npt2: title = Protein-ligand complex NPT equilibration define = -DPOSRES ; position restrain the protein and ligand ; Run parameters integrator = md ; leap-frog integrator nsteps = 500000 ; 2 \* 500000 = 1000 ps = 1ns dt = 0.002 ; 2 fs ; Output control nstenergy = 5000 ; save energies every 10.0 ps nstlog = 5000 ; update log file every 10.0 ps nstxout-compressed = 5000 ; save coordinates every 10.0 ps ; Bond parameters continuation = yes ; Continue from NPT constraint\_algorithm = lincs ; holonomic constraints constraints = h-bonds ; bonds to H are constrained lincs\_iter = 1 ; accuracy of LINCS lincs\_order = 4 ; also related to accuracy ; Neighbor searching and vdW cutoff-scheme = Verlet ns\_type = grid ; search neighboring grid cells nstlist = 20 ; largely irrelevant with Verlet rlist = 1.0 vdwtype = Cut-off vdw-modifier = Potential-shift rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.0 pme\_order = 4 ; cubic interpolation fourierspacing = 0.125 ; grid spacing for FFT ; Temperature coupling tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = protein\_grp non\_protein\_grp ; two coupling groups - more accurate tau\_t = 1 1 ; time constant, in ps ref\_t = 150 150 ; reference temperature, one for each group, in K ; Pressure coupling pcoupl = C-rescale ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau\_p = 1.0 ; time constant, in ps ref\_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 refcoord\_scaling = com ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction used with OPC water model DispCorr = EnerPres ; Velocity generation gen\_vel = no ; velocity generation off npt3: title = Protein-ligand complex NPT equilibration define = -DPOSRES ; position restrain the protein and ligand ; Run parameters integrator = md ; leap-frog integrator nsteps = 500000 ; 2 \* 500000 = 1000 ps = 1ns dt = 0.002 ; 2 fs ; Output control nstenergy = 5000 ; save energies every 10.0 ps nstlog = 5000 ; update log file every 10.0 ps nstxout-compressed = 5000 ; save coordinates every 10.0 ps ; Bond parameters continuation = yes ; continuing from previous NPT constraint\_algorithm = lincs ; holonomic constraints constraints = h-bonds ; bonds to H are constrained lincs\_iter = 1 ; accuracy of LINCS lincs\_order = 4 ; also related to accuracy ; Neighbor searching and vdW cutoff-scheme = Verlet ns\_type = grid ; search neighboring grid cells nstlist = 20 ; largely irrelevant with Verlet rlist = 1.0 vdwtype = Cut-off vdw-modifier = Potential-shift rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.0 pme\_order = 4 ; cubic interpolation fourierspacing = 0.125 ; grid spacing for FFT ; Temperature coupling tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = protein\_grp non\_protein\_grp ; two coupling groups - more accurate tau\_t = 1 1 ; time constant, in ps ref\_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling pcoupl = C-rescale ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau\_p = 2.0 ; time constant, in ps ref\_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 refcoord\_scaling = com ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction used with OPC water model DispCorr = EnerPres ; Velocity generation gen\_vel = no ; velocity generation off ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [August 28, 2024, 11:16am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/2 "2024-08-28T11:16:33Z")
How did you add water to the system? A condensed-phase system with an initial density of 300 kg/m³ signals a lot of vacuum. In VMD, you can visualize changes in box size if you type “pbc box” in the TCL console. ---
**Author:** ![hehehaha](https://avatars.discourse-cdn.com/v4/letter/h/97f17d/32.png) [@hehehaha](https://gromacs.bioexcel.eu/u/hehehaha)\ **Post date:** [August 28, 2024, 11:29am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/3 "2024-08-28T11:29:55Z")
Two considerations, 1. I solvated using AmberTools tleap but found that the water molecules did not seem to enter the chains to make H-bonding network in my initial EM → NVT → NPT → MD 2. So, i solvated using AmberTools tleap + added the crystallography water molecules pbc box using the original trajectory without centering or after centering? I dont remember seeing any vacuum gaps in any of my runs but ill have a look again ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [August 28, 2024, 11:38am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/4 "2024-08-28T11:38:30Z")
AmberTools has a few options to add solvent, it might be that something failed there. Check how many OPC waters you’ve got in your box, and what the box size is. Centering shouldn’t change anything if you have a vacuum there, so either way ---
**Author:** ![hehehaha](https://avatars.discourse-cdn.com/v4/letter/h/97f17d/32.png) [@hehehaha](https://gromacs.bioexcel.eu/u/hehehaha)\ **Post date:** [August 28, 2024, 11:48am UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/5 "2024-08-28T11:48:26Z")
Just to add more details this was how i solvated water molecules in the tleap solvateoct mol OPCBOX 10.0 addions2 mol Na+ 173 addions2 mol Cl- 76 i added 2393 crystallography WAT and this was tleap output which added 43884 water molecules Scaling up box by a factor of 1.210620 to meet diagonal cut criterion Solute vdw bounding box: 95.644 92.735 131.303 Total bounding box for atom centers: 155.515 155.515 155.515 (box expansion for ‘iso’ is 72.5%) Solvent unit box: 18.865 18.478 19.006 Volume: 1917070.582 A^3 (oct) Total mass 1108755.980 amu, Density 0.960 g/cc Added 43884 residues. The eventual topology seems to have the correct total number of WAT molecules WAT 46277 43884 + 2393 = 46277 if the solvation by AmberTools is causing issues maybe ill just solvate with gromacs then (?) ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [August 28, 2024, 3:02pm UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/6 "2024-08-28T15:02:16Z")
So now the question is, how did you start with a density of 0.96 in Amber and got 0.3 in Gromacs? Perhaps your box dimensions changed, or were converted incorrectly, or you converted a file that wasn’t the final one. ---
**Author:** ![hehehaha](https://avatars.discourse-cdn.com/v4/letter/h/97f17d/32.png) [@hehehaha](https://gromacs.bioexcel.eu/u/hehehaha)\ **Post date:** [August 28, 2024, 3:20pm UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/7 "2024-08-28T15:20:31Z")
This are the commands i ran, based on what you are saying, im getting the idea that i may have messed the very first step which is to define the box? #Define box echo 0 | gmx\_mpi editconf -f topol.gro -o newbox.gro -bt octahedron -d 1.2 #EM gmx\_mpi grompp -f em.mdp -c newbox.gro -p topol.top -o em.tpr -v #Make index gmx\_mpi make\_ndx -f em.gro -o index.ndx ##4 phases of NPT to slowly get desired density## gmx\_mpi grompp -f npt.mdp -c em.gro -t em.trr -r em.gro -p topol.top -n index.ndx -o npt.tpr gmx\_mpi grompp -f npt2.mdp -c npt.gro -t npt.cpt -r npt.gro -p topol.top -n index.ndx -o npt2.tpr gmx\_mpi grompp -f npt3.mdp -c npt2.gro -t npt2.cpt -r npt2.gro -p topol.top -n index.ndx -o npt3.tpr gmx\_mpi grompp -f npt4.mdp -c npt3.gro -t npt3.cpt -r npt3.gro -p topol.top -n index.ndx -o npt4.tpr gmx\_mpi grompp -f md\_10ns.mdp -c npt4.gro -t npt4.cpt -p topol.top -n index.ndx -o md\_10ns.tpr ---
**Author:** ![hehehaha](https://avatars.discourse-cdn.com/v4/letter/h/97f17d/32.png) [@hehehaha](https://gromacs.bioexcel.eu/u/hehehaha)\ **Post date:** [August 28, 2024, 5:17pm UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/8 "2024-08-28T17:17:03Z")
I believe that was the mistake echo 0 | gmx\_mpi editconf -f topol.gro -o newbox.gro -bt octahedron -d 1.2 i shouldnt have done that since my gro was already a octahedron… i proceeded to do gmx\_mpi grompp -f em.mdp -c topol.gro -p topol.top -o em.tpr -v but got insane steric clashes Energy minimization has stopped because the force on at least one atom is not finite. This usually means atoms are overlapping. Modify the input coordinates to remove atom overlap or use soft-core potentials with the free energy code to avoid infinite forces. You could also be lucky that switching to double precision is sufficient to obtain finite forces. Steepest Descents converged to machine precision in 23 steps, but did not reach the requested Fmax \< 1000. Potential Energy = 3.7172511e+21 Maximum force = inf on atom 66275 Norm of force = inf the atom 66275 points to two different water molecules (i have 44k water molecules so there are repeated numbers after 99999 atom number) both of the water molecules (which are outside of protein chains, surrounding the protein) do not seem to be clashing with other molecules in the gro file so how shall i proceed from here…? i could go back to my initial topology without crystallography water molecules (which does energy minimise without clashes interestingly) but the water molecules do not seem to enter into the protein chains so a little concerning? i just tried positional restraining the water molecules but same error seems to happen. Any thoughts? ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [August 28, 2024, 6:06pm UTC](https://gromacs.bioexcel.eu/t/confusing-density-from-npt-run/10021/9 "2024-08-28T18:06:11Z")
Oh yes, that solves it. You added another 1.2 nm margin to your original box (-d). Even if your molecules were not broken across PBC before that, it will make for a strangely behaving system. I remember there’s a tiny difference in the octahedron definition between Gromacs and Amber - some last decimal places in the box angle value - so that if you just transfer them naively, things will not fit exactly. So as you pointed out, if you have the possibility to adjust the box and solvate in Gromacs again, it tends to be the safer option.