Thank you Justin and Peter for your responses.  I tried extending the time on the npt equilibration.  It helped but not much.  My final pressure after the MD run was about 1.15 bar compared to ref_p which was set to 1 bar.  Peter, I will try to analyze the potential error using RMSD and Drift.  <br>
<br>Thanks.<br><br><div class="gmail_quote">On Mon, Apr 11, 2011 at 4:55 PM, Fabian Casteblanco <span dir="ltr">&lt;<a href="mailto:fabian.casteblanco@gmail.com">fabian.casteblanco@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<div>Hi,</div>
<div> </div>
<div>I&#39;m still in my first few months of using Gromacs.  I started by creating an *.itp and *.top file for <i>Ethanol</i> using CHARMM force field parameters.  I made the molecule and it looked fine, put 1000 molecules in a box, energy minimized it to a negative potential energy, viewed it on VMD, again looks fine.  When I started running the NVT script, I set it equal to a ref_T of 298 K.  It equilibrated at the temperature.  Then I tried using an NPT script to equilibrate it to a ref_p of 1 bar.  This is where I get the problem.  The output shows the density is close to the actual experimental value of 0.789 g/cm^3.  But for some reason, my pressure never gets an average of 1 bar.  It keeps oscillating, which I understand is normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and 1.45 for 75,000 steps,dt=0.002).  I don&#39;t understand why the ref_p of 1 bar is not working when I run this NPT.mdp script file.  My simple goal is to have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 bar and somewhere near the experimental density.</div>


<div> </div>
<div>I would really appreciate anybody&#39;s help!  I&#39;m new to this but I&#39;m eager to keep getting better.</div>
<div> </div>
<div>Thanks.</div>
<div> </div>
<div><u>NVT SCRIPT  (this works fine and takes me to 298 K)</u></div>
<div>File Edit Options Buffers Tools Help<br>title           =CHARMM ETHANOL  NVT equilibration<br>;define         =-DPOSRES       ;position restrain the protein<br>;Run parameters<br>integrator      =md             ;leap-frog algorithm<br>

nsteps          =50000          ;2 * 50000 = 100 ps<br>dt              =0.002          ;2fs<br>;Output control<br>nstxout         =100            ;save coordinates every 0.2 ps<br>nstvout         =100            ;save velocities every 0.2 ps<br>

nstenergy       =100            ;save energies every 0.2 ps<br>nstlog          =100            ;update log file every 0.2 ps<br>;Bond parameters<br>continuation    =no             ;first dynamics run<br>constraint_algorithm=lincs      ;holonomic constraints<br>

constraints     =all-bonds      ;all bonds (even heavy atom-H bonds)constraind<br>lincs_iter      =1              ;accuracy of LINCS<br>lincs_order     =4              ;also related to accuracy<br>;Neighborhood searching<br>

ns_type         =grid           ;search neighboring grid cells<br>nstlist         =5              ;10 fs<br>rlist           =1.0            ;short-range neighborlist cutoff (in nm)<br>rcoulomb        =1.0            ;short-range electrostatic cutoff (in nm)<br>

rvdw            =1.0            ;short-range van der Waals cutoff (in nm)<br>;Electrostatics<br>coulombtype     =PME            ;Particle Mesh Ewald for long-range electrostat\<br>;ics<br>pme_order       =4              ;cubic interpolation<br>

fourierspacing  =0.16           ;grid spacing for FFT<br>;Temperature coupling is on<br>tcoupl          =V-rescale      ;modified Berendsen thermostat<br>tc_grps         =SYSTEM   ;two coupling groups - more accurate<br>
tau_t           =0.1    ;0.1              ;time constant, in ps<br>
ref_t           =298    ;25               ;reference temperature, one for each \<br>;group, in K<br>;Pressure coupling is off<br>pcoupl          =no             ;no pressure coupling in NVT<br>;Periodic boundary conditions<br>

pbc             =xyz            ; 3-D PBC<br>;Dispersion correction<br>DispCorr        =EnerPres       ;account for cut-off vdW scheme<br>;Velocity generation<br>gen_vel         =yes            ;assign velocities from Maxwell distribution<br>

gen_temp        =25             ;temperature for Maxwell distribution<br>gen_seed        =-1             ;generate a random seed</div>
<div>;END</div>
<div> </div>
<div><u>NPT SCRIPT</u></div>
<div>File Edit Options Buffers Tools Help<br>title           =Ethanol npt equilibration<br>;define         =-DPOSRES       ;position restrain the protein<br>;Run parameters<br>integrator      =md             ;leap-frog algorithm<br>

nsteps          =50000          ;2 * 50000 = 100 ps<br>dt              =0.002          ;2fs<br>;Output control<br>nstxout         =100            ;save coordinates every 0.2 ps<br>nstvout         =100            ;save velocities every 0.2 ps<br>

nstenergy       =100            ;save energies every 0.2 ps<br>nstlog          =100            ;update log file every 0.2 ps<br>;Bond parameters<br>continuation    =yes            ;Restarting after NVT<br>constraint_algorithm=lincs      ;holonomic constraints<br>

constraints     =all-bonds      ;all bonds (even heavy atom-H bonds)constraind<br>lincs_iter      =1              ;accuracy of LINCS<br>lincs_order     =4              ;also related to accuracy<br>;Neighborhood searching<br>

ns_type         =grid           ;search neighboring grid cells<br>nstlist         =5              ;10 fs<br>rlist           =1.0            ;short-range neighborlist cutoff (in nm)<br>rcoulomb        =1.0            ;short-range electrostatic cutoff (in nm)<br>

rvdw            =1.0            ;short-range van der Waals cutoff (in nm)<br>;Electrostatics<br>coulombtype     =PME            ;Particle Mesh Ewald for long-range electrostat\<br>;ics<br>pme_order       =4              ;cubic interpolation<br>

fourierspacing  =0.16           ;grid spacing for FFT<br>;Temperature coupling is on<br>tcoupl          =V-rescale      ;modified Berendsen thermostat<br>tc-grps         =SYSTEM   ;two coupling groups - more accurate<br>
tau_t           =0.1;   0.1               ;time constant, in ps<br>
ref_t           =298;   300               ;reference temperature, one for each \<br>;group, in K<br>;Pressure coupling is on<br>pcoupl          =Parrinello-Rahman      ;Pressure coupling on in NPT<br>pcoupltype      =isotropic              ;uniform scaling of box vectors<br>

tau_p           =2.0                    ;time constant, in ps<br>ref_p           =1.0                    ;reference pressure, in bar<br>compressibility =4.5e-5         ;isothermal compressibility of h2O, 1/bar<br>;Periodic boundary conditions<br>

pbc             =xyz            ; 3-D PBC<br>;Dispersion correction<br>DispCorr        =EnerPres       ;account for cut-off vdW scheme<br>;Velocity generation<br>gen_vel         =no             ;Velocity generation is off<br>

;gen_temp       =25             ;temperature for Maxwell distribution<br>;gen_seed       =-1             ;generate a random seed</div>
<div>;END</div>
<div><br clear="all"><br>-- <br></div>
<div><b><font face="arial,helvetica,sans-serif">Best regards,</font></b></div>
<div><b></b><font face="arial,helvetica,sans-serif"> </font></div>
<div><b><font face="arial,helvetica,sans-serif" size="4">Fabian F. Casteblanco</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">Rutgers University -- </font></b></div>
<div><b><font face="arial,helvetica,sans-serif">Chemical Engineering PhD Student</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">C: +908 917 0723</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">E:  </font></b><a href="mailto:fabian.casteblanco@gmail.com" target="_blank"><b><font face="arial,helvetica,sans-serif">fabian.casteblanco@gmail.com</font></b></a></div>
<br>
</blockquote></div><br><br clear="all"><br>-- <br><div><b><font face="arial,helvetica,sans-serif">Best regards,</font></b></div>
<div><b></b><font face="arial,helvetica,sans-serif"> </font></div>
<div><b><font face="arial,helvetica,sans-serif" size="4">Fabian F. Casteblanco</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">Rutgers University -- </font></b></div>
<div><b><font face="arial,helvetica,sans-serif">Chemical Engineering PhD Student</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">C: +908 917 0723</font></b></div>
<div><b><font face="arial,helvetica,sans-serif">E:  </font></b><a href="mailto:fabian.casteblanco@gmail.com" target="_blank"><b><font face="arial,helvetica,sans-serif">fabian.casteblanco@gmail.com</font></b></a></div><br>