[AMBER] Toluene solvent box equilibration

From: Pedro Miguel Reis Figueiredo <pmrfigueiredo.cnc.uc.pt>
Date: Wed, 30 Mar 2022 21:17:02 +0100

Dear Amber Users,

I am trying to equilibrate a cubic box of toluene, but I am facing some
troubles.

I have followed the tutorial "A room-temperature ionic liquid
(ambermd.org [1])", where I have substituted the ACN, BMIM and BF4 units
with toluene molecules only. I have successfully created the mol2 and
frcmod files for toluene and created a box of 150A3 with 19640 molecules
(corresponding to the toluene's density of 0.87g/cm3).

Using the following LEaP input, I have created the topology and
coordinate files:

> _source leaprc.protein.ff14SB_
>
> _loadamberparams tol.frcmod_
>
> _ _
>
> _TOL = loadmol2 toluene.mol2_
>
> _BOX = loadpdb BOX.pdb_
>
> _ _
>
> _setbox BOX centers_
>
> _saveamberparm BOX boxtol.prmtop boxtol.rst7_
>
> _quit_

Then I ran the minimization (min.in [2]) and equilibration (eq.in [3])
simulations, but my system exploded during the equilibration step
without any apparent reason.

min.in [2]:

> _minimization_
>
> _ &cntrl_
>
> _ imin=1, ntpr=100, ntwx=100, maxcyc=10000,_
>
> _ ntb=1,_
>
> _ &end_

eq.in [3]:

> _3ns equilibration box_
>
> _&cntrl_
>
> _ imin=0, ntpr=3000, ntwx=3000,_
>
> _ ntx=1, irest=0,_
>
> _ tempi=293., temp0=293., ntt=3,_
>
> _ gamma_ln=5., ntb=2, ntp=1, taup=1.0, ioutfm=1, nstlim=3000000,_
>
> _ ntwr=1000, dt=.001, ig=-1,_
>
> _/_

Trying to solve this, I have tested the inputs I use for protein MD
simulation:

min1.in:

> _initial minimization of the whole system_
>
> _ &cntrl_
>
> _ imin = 1,_
>
> _ maxcyc = 5000,_
>
> _ ncyc = 2000,_
>
> _ ntb = 1,_
>
> _ ntr = 0,_
>
> _ cut = 10.0_
>
> _ /_

heat.in:

> _500ps MD equilibration_
>
> _ &cntrl_
>
> _ imin = 0,_
>
> _ irest = 0,_
>
> _ ntx = 1,_
>
> _ ntb = 1,_
>
> _ cut = 10.0,_
>
> _ ntr = 0,_
>
> _ ntc = 2,_
>
> _ ntf = 2,_
>
> _ tempi = 0.0,_
>
> _ temp0 = 293.15,_
>
> _ ntt = 3,_
>
> _ gamma_ln = 1.0,_
>
> _ nstlim = 250000, dt = 0.002_
>
> _ ntpr = 5000, ntwx = 5000, ntwr = 1000_
>
> _ /_

In this step, the output print a lot of ***** (not sure if this is
normal, since I have not a solute):

> NSTEP = 95000 TIME(PS) = 190.000 TEMP(K) =********* PRESS = 0.0
>
> Etot = ************** EKtot = ************** EPtot = **************
>
> BOND = -0.0001 ANGLE = 27449525.5430 DIHED = 1805073.3845
>
> 1-4 NB = 0.0000 1-4 EEL = 0.0014 VDWAALS = **************
>
> EELEC = -10309.6839 EHBOND = 0.0000 RESTRAINT = 0.0000

prod.in

> _1ns MD_
>
> _ &cntrl_
>
> _ imin = 0, irest = 1, ntx = 7,_
>
> _ ntb = 2, pres0 = 1.0, ntp = 1,_
>
> _ taup = 2.0,_
>
> _ cut = 10.0, ntr = 0,_
>
> _ ntc = 2, ntf = 2,_
>
> _ tempi = 293.15, temp0 = 293.15,_
>
> _ ntt = 3, gamma_ln = 1.0,_
>
> _ nstlim = 500000, dt = 0.002,_
>
> _ ntpr = 2500, ntwx = 2500, ntwr = 10000_
>
> _ /_

 In the production simulation, the system explode presenting this
message:

> ERROR: Calculation halted. Periodic box dimensions have changed too much from their initial values. Your system density has likely changed by a large amount, probably from starting the simulation from a structure a long way from equilibrium.
>
> [Although this error can also occur if the simulation has blown up for some reason]
>
> The GPU code does not automatically reorganize grid cells and thus you will need to restart the calculation from the previous restart file. This will generate new grid cells and allow the calculation to continue. It may be necessary to repeat this restarting multiple times if your system is a long way from an equilibrated density. Alternatively you can run with the CPU code until the density has converged and then switch back to the GPU code.

Running with CPU or limiting the GPU code with the flag vlimit does not
resolve this. I have observed high shifts in systems' pressure and the
density increase from 0.8 to 1.0:

> NSTEP = 42500 TIME(PS) = 585.000 TEMP(K) =********* PRESS = 31224.2
>
> Etot = ************** EKtot = ************** EPtot = **************
>
> BOND = -0.0001 ANGLE = 28900440.0113 DIHED = 1802556.7630
>
> 1-4 NB = 0.0000 1-4 EEL = 0.0004 VDWAALS = **************
>
> EELEC = -3284.3515 EHBOND = 0.0000 RESTRAINT = 0.0000
>
> EKCMT = -0.0000 VIRIAL = -2527617.1317 VOLUME = 3749235.0471 Density = 0.8014
>
> ------------------------------------------------------------------------------
>
> NSTEP = 45000 TIME(PS) = 590.000 TEMP(K) =********* PRESS = 17159.0
>
> Etot = ************** EKtot = ************** EPtot = **************
>
> BOND = -0.0001 ANGLE = 28883760.4779 DIHED = 1803516.5749
>
> 1-4 NB = 0.0000 1-4 EEL = 0.0004 VDWAALS = **************
>
> EELEC = -13758.8548 EHBOND = 0.0000 RESTRAINT = 0.0000
>
> EKCMT = -0.0000 VIRIAL = -1151566.2531 VOLUME = 3108273.3687 Density = 0.9667
>
> ------------------------------------------------------------------------------
>
> NSTEP = 47500 TIME(PS) = 595.000 TEMP(K) =********* PRESS =-77840.8
>
> Etot = ************** EKtot = ************** EPtot = **************
>
> BOND = -0.0001 ANGLE = 28812866.2413 DIHED = 1802867.7596
>
> 1-4 NB = 0.0000 1-4 EEL = 0.0004 VDWAALS = **************
>
> EELEC = -18520.2184 EHBOND = 0.0000 RESTRAINT = 0.0000
>
> EKCMT = -0.0000 VIRIAL = 4602393.4852 VOLUME = 2738409.5564 Density = 1.0973

Any thoughts on what I am doing wrong?

Thank you in advance!

Kind regards,

Pedro

Links:
------
[1] http://ambermd.org
[2] http://min.in
[3] http://eq.in
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Mar 30 2022 - 13:30:03 PDT
Custom Search