# System blowing up during equilibration step **URL:** \ **Category:** User discussions\ **Created:** [October 15, 2024, 11:42am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344 "2024-10-15T11:42:41Z")\ **Posts on this page:** 12\ **Page:** 1
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 15, 2024, 11:42am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/1 "2024-10-15T11:42:41Z")
GROMACS version: 2020.4-MODIFIED Dear Gromacs Users, 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. ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 15, 2024, 1:05pm UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/2 "2024-10-15T13:05:20Z")
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. ---
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 15, 2024, 5:43pm UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/3 "2024-10-15T17:43:38Z")
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 ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 16, 2024, 6:26am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/4 "2024-10-16T06:26:44Z")
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`. ---
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 17, 2024, 5:52am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/5 "2024-10-17T05:52:51Z")
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. Thanks > **(attachments)** > > ![CHCL3_Isothermal_Compressibility.JPG](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/6/61321fe947511910bfab7163c84d6361613f8977.jpeg) ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 17, 2024, 6:33am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/6 "2024-10-17T06:33:51Z")
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. ---
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 17, 2024, 6:51am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/7 "2024-10-17T06:51:34Z")
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. Thanks ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 17, 2024, 6:53am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/8 "2024-10-17T06:53:32Z")
Yes, `c-rescale` is recommended both for equilibration and production. ---
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 17, 2024, 6:57pm UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/9 "2024-10-17T18:57:31Z")
Thank you so much MagnusL for your response. I was going through some of the previous posts and found this: [https://gromacs.bioexcel.eu/t/on-compressibility-under-pressure-coupling-section-of-mdp/5421](https://gromacs.bioexcel.eu/t/on-compressibility-under-pressure-coupling-section-of-mdp/5421) 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. Thanks ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 18, 2024, 7:14am UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/10 "2024-10-18T07:14:25Z")
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. ---
**Author:** ![NidhiM](https://avatars.discourse-cdn.com/v4/letter/n/3bc359/32.png) [@NidhiM](https://gromacs.bioexcel.eu/u/NidhiM)\ **Post date:** [October 19, 2024, 4:11pm UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/11 "2024-10-19T16:11:22Z")
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. Looking forward to your response. Thanks > **(attachments)** > > ![Initial.jpg](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/f/fd78ef26015e4fedc80dc66b407062c23398e45e.jpeg) > ![After200ns.jpg](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/c/c6c965262669cc68409ea4f1fbe853c44d531abe.jpeg) ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 20, 2024, 7:34pm UTC](https://gromacs.bioexcel.eu/t/system-blowing-up-during-equilibration-step/10344/12 "2024-10-20T19:34:42Z")
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.