[gmx-users] pressure and potential energy

[abhishek khetan]

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