Hi, I am currently performing NVT equilibration of i-pma in a water-CCl4 box. Four moles of pre-equilibrated pma with a length of 30 units were randomly inserted into the pre-equilibrated water-CCl4 box. The system then underwent NVT equilibration for 500 ps, resulting in a temperature of 299.9 K (reference temperature = 300 K). Subsequently, I ran the NVT equilibration for an additional 1 ns, but encountered the following error:
Fatal error:
2 particles communicated to PME rank 1 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x. This usually indicates that your system is not well equilibrated.
How can I resolve this issue? Your assistance would be greatly appreciated. The nvt.mdp file is given below:
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = LBER M486 SOL ; two coupling groups - more accurate
tau_t = 0.1 0.1 0.1 ; time constant, in ps
ref_t = 300 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
Please help. Thanks in advance
You are using three temperature coupling groups. Are they very large, all of them? Small temperature coupling groups can lead to instabilities. In most cases it’s better to just use one temperature coupling group. You should probably also use comm-mode linear.
In case you still get crashes, I would suggest reducing dt to 0.001 or 0.0005 and see if it helps, i.e. in the NVT equilibration only.
Not related to the crashes:
- What force field are you using? rcoulomb = 1.4 and rvdw = 1.4 is higher than what is usually used.
- Why are you writing coordinates, velocities and energies so often? And why are you writing uncompressed coordinates?
Thank you for your help @MagnusL . The temperature coupling groups are CCl4 (M486), H2O (SOL) and PMA(LBER). here, PMA is the polymer and CCl4 and H2O are solvent molecules. As you mentioned I have used PMA and WATER as temperature coupling groups but getting error :
Fatal error:
6060 atoms are not part of any of the T-Coupling groups
which indicates the CCl4 group is missing. How can I solve this issue?
@MagnusL I am using Gromos 54 a7 force field. Currently i am reproducing a research paper. There they have mentioned the values for rcoulomb and rvdw as 1.4. Thats why i have chosen the same value. I didn’t get the next question. Could you explain once more?
Regarding temperature coupling groups, I often recommend using one (system), but if there is not enough temperature/velocity exchange between separate parts of the system you could use two groups (e.g. water_solutes_and_ions protein_and_lipids). You will have to define the selection/index groups yourself (see gmx make_ndx). If you need three or more temperature coupling groups, I think they will risk getting too small for temperature coupling to work well, but it all depends on your system of course.
OK. If those cutoffs are recommended for the force field then you should use them.
What I meant with my second (off-topic) comment was about:
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
Why are you using nstxout and nstvout? There’s nothing wrong with it at all, but you should ask yourself why you need it. In many cases you will just get very large files that you don’t have any use for. Why are you writing energies so often? Do you need as many as 1000 observation points/frames in your analyses? If you need to analyse why a simulation is crashing it can be good to write coordinates often. But in many cases, for 500,000 steps I’d say something similar to:
nstxout-compressed = 25000 ; save compressed coordinates every 50 ps
nstenergy = 10000 ; save energies every 20 ps (gives 50 output frames)
nstlog = 10000 ; update log file every 20 ps
would be OK. If you want to make movies, or if states you want to observe are very short-lived, you might want more coordinate frames. If you think plots in your analyses don’t look good you might want to have more energy output frames.
Dear sir,
I have tried NVT equilibration for 1ns with
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
But again I am getting the error:
Fatal error:
1 particles communicated to PME rank 0 are more than 2/3 times the cut-off out
of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.
Could you please help with my issue?
The first thing I would do it to set tau-t to 1 or 2 ps and see if that helps. If that is not enough to solve the problem, I would try to set dt to 0.001. I.e., if you did not already to that. You could see if it helps with two separate comm-grps (e.g., one with the i-pma (assuming that you mean isopropyl myristate) and one with the rest).
Thank you sir, I will try
@MagnusL , thank you for your help. I have changed dt=0.002 to dt=0.001. Its working now.
Thank you
Once the system is stable you should be able to switch back to dt=0.002.
Dear Sir,
As per your instructions, I have conducted the NPT equilibration . The simulation was started from 100 ps. At 1600 ps, I observed a pressure value of 0.9156 bar, with a reference pressure set at 1 bar. Based on this pressure value, I am uncertain whether the system has stabilized completely. Would it be advisable to run an additional 100 ps simulation to further assess stability?
Regarding the density, the average value obtained was 1252.77 kg/m³. To compare this value meaningfully, I need guidance on which standard density value to use for comparison. Would it be appropriate to compare this density value with that of water only, or are there other reference values I should consider?
Your guidance in this matter would be greatly appreciated.
Thank you, sir.
When you say that you “observed a pressure value of 0.9156 bar”, do you mean a single point observation or do you mean the average pressure over the NPT simulation (preferably discarding the first 100 or 200 ps using the -b option in gmx energy)? If it is the average, what is the uncertainty? 0.92+/-??? sounds very reasonable, especially if the system is fairly small and since the simulation time was quite short. See, e.g., Pressure is not going to 1 bar - #2 by cblau.
With a system consisting mostly of water I’d expect the density to be close to the density of the water model that you are using. I’m not sure how well, e.g., proteins and lipids in the commonly used force fields reproduce empirical densities. What is the empirical density of your polymer?
In the log file, I got the value for pressure=9.15914e-01 from averages :
<====== ############### ==>
<==== A V E R A G E S ====>
<== ############### ======>
Pressure (bar)
9.15914e-01
Here, number of molecules: water= 8483, CCl4=1212 and pma(poly methacrylic acid )=4 chains (each chain of poly methacrylic acid consisting of 30 repeating units).
The density of pma remains within the range around 1,100 to 1,300 kilograms per cubic meter (kg/m³).
Which density should I consider?
please help sir.
Thank you
After running the NPT equilibration for 1600 ps (average pressure= 0.9156 bar), an additional 100 ps NPT equilibration was conducted, totaling 1700 ps. Despite this, the average pressure remains at -1.11193e+00 bar. The pressure continues to vary over time. How can I confirm that the system has stabilized
Since you’ve got the densities of the components I’m sure you can calculate an approximate theoretical density as well as I can.
Use gmx energy to get the statistics of the pressure. Just one average value does not say very much.
Only you can decide when you think your system is stable enough and you should be able to motivate it. I would recommend gmx energy to check that, e.g., the volume, temperature and potential energy are stable (there is no clear rule what is stable) in the last half of the equilibration simulation (or for a period you think is appropriate to make you confident that you are sampling what you are interested in). The values will fluctuate over time, but they should fluctuate around the value you expect.
The time required for equilibrating a system depends on the complexity of the system and the starting structure. It can take anything from a few hundred ps (for simple solvents) to at least a few hundred ns.
Dear sir,
Thank you for your assistance. I have completed the NPT equilibration with an average pressure of 8.01649e-01 at 28 ns. Now, I am conducting a 1 ns MD simulation. During the production run, the system’s temperature was coupled with a Nose-Hoover thermostat to 300 K, and the pressure coupling of the system to 1 bar was performed semi-isotropically using the Parrinello-Rahman method while keeping the x- and y-dimensions constant.
md.mdp file is given:
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 500000 = 1000 ps (1 ns)
dt = 0.001 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = nose-hoover ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; Apply pressure coupling separately in X-Y and Z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
But I am getting error :
Fatal error:
The X-size of the box (5.044856) times the triclinic skew factor (1.000000) is
smaller than the number of DD cells (5) times the smallest allowed cell size
(1.008875)
So, I have relaxed the barostat coupling time to allow for a
more gentle equilibration by applying tau_p=5. still it shows the same error. How to solve the issue without changing the box size of the system?
Thanks in advance
what is your mdrun flags? how many cores you are running? or gpu?
@scinikhil thank you for your reply. I am suing the mdrun command:
mdrun -deffnm md
40 cores were used for he job running.
The script which I am using is given below:
#SBATCH --job-name=packmol_test # Job name
#SBATCH --ntasks=40 # Number of MPI tasks (i.e. processes)
#SBATCH --cpus-per-task=1 # Number of cores per MPI task
#SBATCH --ntasks-per-node=20
#SBATCH --nodes=2 # Maximum number of nodes to be allocated
#SBATCH --partition=mediumq
#SBATCH --output=md-15ns-run.log
module load apps/gromacs/2020.7/gromacs-2020.7-gpu
#name of the executable
exe=“gmx_mpi”
export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
#run the application
mpirun -bootstrap slurm -n $SLURM_NTASKS $exe mdrun -deffnm md-1ns
How did you paramterise these molecules? and for which force field? i suspect the issue is in the paramters. Can you share the itp files for the parameters?
I would reparameterise carefully and minimise the system a couple of times and equilibrate for longer time.
I have parametrised using ATB. gromos 54a7 force field is used as force field. I have used original geometry all-atom pdb file of polymer downloaded from ATB. I am trying to upload the itp files here. But since .itp files can not be uploaded here, I am uploading them as .log files. The extenstion is changed from .itp to .log just for the sake of uploading. LBER corresponds to the polymer itp file and M486 corresponds to the CCl4 itp file
LBER.log (87.6 KB)
M486.log (5.3 KB)
Thank you