I am trying to simulate an organic molecule in Chloroform solvent using Gromos54a7 ff (it is beta sheet mimic-type system so I believe gromos54a7 should be ok). I have created the topology of solvent molecule as well as solute molecule using ATB, and used gmx insert-molecules to insert appropriate number of molecules such that the density of the system is close to the experimental value. Energy minimization ran fine. NVT Equilibration did run without errors while npt equilibration failed with LINCS Warnings, suggestive of system blowing-up. On visualization of nvt.gro in VMD software, I found some bonds being broken (unusual of classical MD). I tried nvt equilibration by restraining first the entire solute molecule and also by restraining only the heavy atoms in solute molecule. The run fails in both the conditions. Even restraining solvent molecules does not help. I have tried with reducing the step size to 1fs. Also, I have tried by using pre-equilibrated chloroform box. But the problem still remains the same.
Here’s how my nvt.mdp looks like:
title = NVT equilibration
define = -DPOSRES -DPOSRES_SOL ; position restrain
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 50000 = 1000 ps
dt = 0.001 ; 1 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.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; 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 = DY4B SJ2H ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 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
Kindly help me resolve the issue and let me know if any further information is required at my end.
Even if it looks like bonds are breaking in VMD it does not mean that they are - in fact they cannot be broken. It’s just a visualization artifact.
At a brief glance, your NVT parameters look OK. I’d suspect that the problem is with the NPT equilibration settings, or that you need to run NVT for longer before continuing to the NPT stage.
Thank you so much MagnusL for your response. I am sharing the npt.mdp file for your reference. I have changed isothermal compressibility value to the literature reported value of chloroform. I am speculative that this value could be the reason for instability of my system. Am I interpreting this parameter correctly or should I change it back to 4.5e-5? Additionally, I have changed the barostat to ‘Berendsen’ and played with different tau values.
My npt.mdp:
title = NPT equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 500000 = 1 ns
dt = 0.001 ; 1 fs
; Output control
nstxout = 5000 ; save coordinates every 10.0 ps
nstvout = 5000 ; save velocities every 10.0 ps
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
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 scheme
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; 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 = DY4B SJ2H ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in 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 = 1.05e-2 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off
First, I would recommend to always use pcoupl = c-rescale.
I think the compressibility of chloroform should be closer to 1e-4. From where did you get your value? High compressibility values should work fine, but I guess (don’t know for sure) that you would might need to raise tau_p.
Thank you so much for your response. I ran a test simulation using the compressibility value of water (4.5 × 10^-5 bar^-1), and the run completed successfully without any errors. It appears that the issue was indeed related to the compressibility value itself.
I obtained this value from the CRC Handbook of Chemistry and Physics, and I am attaching a screenshot of the relevant section for your reference.
The value in the screenshot supports slightly higher than twice as high as water, i.e., ~1e-4 bar^-1. I think you’re just off by a factor 100 in your unit conversions.
Yes, but I am unable to figure out the factor of 100 by unit conversion. Nevertheless, I have set up the system by taking ~1e-4 bar^-1 as the value. Hopefully, it should be ok.
Secondly, in one of the previous posts, I read that if we take c-rescale for pressure coupling during equilibration, it is preferred to use it during production run as well. Since you recommended me to use c-rescale for pressure coupling, I wanted to cross check if its ok to shift to Parrinello Rahman during production run.
According to the post, it does not really matter whether we change isothermal compressibility values or not. What are your thoughts on this?
Different papers reported different values. Also I am unable to figure out the reason behind a factor of 100 as per our previous discussion.
That post does not say that it does not matter if you change the compressibility, but it is not very important. A factor 2 difference will probably not make a noticeable difference (the tau-p recommendations often differ more than that anyhow). But if you know a correct value I think it makes sense to use it. A factor 100 difference can certainly cause problems, as you’ve noticed.
Regarding the factor 100:
Your table says (approximately) 10e-4/MPa, which is 10e-10/Pa, right? This is 10e-5/bar, which is 1e-4/bar.
As a sanity check, do the same for water, based on your table, and you will get:
4.5e-4/MPa, which is 4.5e-10/Pa, which is 4.5e-5/bar.
1e-4/bar is a factor 100 lower than your original setting of 1e-2/bar.
Thank you so much for your response; it was incredibly helpful in clarifying my doubts. Following your advice, I ran the simulations without errors. However, I believe the results may be incorrect. My molecule is a tetramer containing porphyrin rings stacked together, but upon simulation, the pyrrole rings of the porphyrins and the entire porphyrin structure appear distorted and out of plane. This might be a force field issue, though I am not entirely sure. I am also considering whether applying some restraints to the porphyrin rings might be necessary.
When visualizing the molecule, I noticed broken bonds, which could be an artifact, as you suggested earlier. Do you have any recommendations on how to address this?
Here are the steps I followed for my simulation:
I have used Gaussian software to optimize the initial structure of my molecule and then used it to create both the topology of my molecule and solvent using Automated Topology Builder (Gromos54a7 force field). NVT and NPT equilibration using position restraints on all atoms other than hydrogen. Production run in NPT ensemble.
I have attached snapshots of the structures (before and after production run) for your reference. Could you please guide me in identifying the problem and finding a potential solution? If you need any further information, please let me know.
Regarding the broken bonds, you are using the DynamicBonds visualisation option in VMD, right? You can just change the Distance Cutoff to 1.7 or 1.8 Å, that will probably make the bonds show up in the visualisation.
Regarding the structure after simulation, there can be many explanations. First, the force field can certainly be worth considering. Your ligand (monomer) seems quite large, most (semi-)automatic parameter/topology generation tools are far from perfect. There is a high probability that there are moieties that are not perfectly represented. You could try generating parameters for it in GAFF (ACPYPE/Antechember or STaGE) or CGenFF (e.g. Charmm-GUI or STaGE) and compare if there are large discrepancies in the parameters related to, e.g., the pyrrole rings are similar. Also, are you sure that the tetramer would retain its structure in solvent and at the temperature you are using? Another thing to try would of course be to see if it helps running the NPT equilibration with position restraints significantly longer to make sure that the system, at least the solvent, is properly equilibrated in the desired ensemble. I would not recommend using restraints during production simulations, unless that’s the only way you can study the phenomena you are interested in. The question would still be whether your results are be valid or not.
