[gmx-users] pressure and potential energy

[Mark Abraham]

Hi,

Did grompp warn you that your time step was too large for a system without
constraints? What did you do about that?

Mark

On Thu, 24 Mar 2016 00:28 abhishek khetan <askhetan at gmail.com> wrote:

> Dear gmxers,
>
> I've been tryin to simulate ion-pairing in non-aqueous solvents. I first
> did the energy of the pure solvent box without any dynamics and
> minimization to emtol= 10 and noticed that the last details of the system
> were:
>
>    Energies (kJ/mol)
>            Bond          Angle Ryckaert-Bell.          LJ-14     Coulomb-14
>     2.52698e+01    1.44838e+02    2.55442e+01    3.67978e+02    2.19911e+03
>         LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      LJ recip.
>    -1.11939e+03   -1.08246e-06   -1.28234e+03    5.66549e+01   -1.33420e+03
>       Potential Pres. DC (bar) Pressure (bar)
>    -9.16532e+02    0.00000e+00   -1.63119e-06
>
> Where the potential energy is clearly negative. I then use this an an input
> for the NVT MD run for 50 ns (in 5 successive runs of 10ns, each time
> generating a new *tpr by using the *gro and *cpt files, using a time
> step=0.002 ps and ref-t=300K). The output that i get at the end is:
>
>     <======  ###############  ==>
>     <====  A V E R A G E S  ====>
>     <==  ###############  ======>
>
>     Statistics over 5000001 steps using 50001 frames
>
>     1.43175e+03    1.99400e+03    4.10265e+02    5.00595e+02    2.70866e+03
>         LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      LJ recip.
>    -7.35653e+02   -1.08246e-06   -1.89995e+03    9.78779e+01   -1.32252e+03
>       Potential    Kinetic En.   Total Energy  Conserved En.    Temperature
>     3.18503e+03    3.82951e+03    7.01454e+03    7.09544e+03    3.00153e+02
>  Pres. DC (bar) Pressure (bar)
>    -1.63119e-06    3.13881e+00
>
>    Total Virial (kJ/mol)
>     1.27257e+03   -1.98606e-01    3.62448e+00
>    -1.98665e-01    1.27190e+03   -4.08938e-01
>     3.62448e+00   -4.08919e-01    1.28192e+03
>
>    Pressure (bar)
>     1.23205e+01   -4.20579e-01   -9.97975e+00
>    -4.20402e-01    1.42463e+01    9.48858e-01
>    -9.97974e+00    9.48803e-01   -1.71505e+01
>
> Where, as you can clearly see, the potential energy is positive, which i've
> read at many places should never be positive, because that means that the
> system is not cohesive/stable.I want to understand what could be going
> wrong here ??
>
> Next, I go to NPT and try to equilibrate everything at 1 bar. When I take
> timestep =0.001 ps, then box , whose original size is 2.22 nm , starts to
> blow up and become 2.7 nm. However, when i use timestep=0.002 ps, the size
> of the box comes back at 2.3 nm. I would have thought that decreasing the
> timestep would increase accuracy, but here it is quite the opposite. The
> current output pressure at the end of log file, after having already done
> 20 ns of NPT looks like
>
> Energies (kJ/mol)
>            Bond          Angle Ryckaert-Bell.          LJ-14     Coulomb-14
>     1.07452e+03    1.97490e+03    6.76611e+02    4.87055e+02    2.63633e+03
>         LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      LJ recip.
>    -3.63376e+02   -8.84997e-07   -1.75096e+03    1.14315e+02   -1.12029e+03
>       Potential    Kinetic En.   Total Energy  Conserved En.    Temperature
>     3.72912e+03    3.81619e+03    7.54531e+03    5.13130e+06    2.99109e+02
>  Pres. DC (bar) Pressure (bar)
>    -1.09082e-06    2.34149e+03
>
> Where, as you can see, the potential energy is again positive and the
> pressure is in the order of 1000 bars. I know that what matters is the
> averaged pressure, but the average pressure is also of the order of 100
> bars. I just try to keep extending the simulation, so that I can hit an
> equilibrated stage but it doesn't seem to arrive. I understand that the
> pressure fluctuations can be huge but the mean pressure is so so far away
> from 1 bar. Could you please help me diagnose this problem? I use
> nose-hoover and MTTK couplings with md-vv integrator.
>
> Thanks in advance
>
> --
> MfG,
> abhishek
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>