[gmx-users] Pressure coupling problem

Fabian Casteblanco fabian.casteblanco at gmail.com
Mon Apr 11 22:55:16 CEST 2011


Hi,

I'm still in my first few months of using Gromacs.  I started by creating an
*.itp and *.top file for *Ethanol* 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'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.

I would really appreciate anybody's help!  I'm new to this but I'm eager to
keep getting better.

Thanks.

*NVT SCRIPT  (this works fine and takes me to 298 K)*
File Edit Options Buffers Tools Help
title           =CHARMM ETHANOL  NVT equilibration
;define         =-DPOSRES       ;position restrain the protein
;Run parameters
integrator      =md             ;leap-frog algorithm
nsteps          =50000          ;2 * 50000 = 100 ps
dt              =0.002          ;2fs
;Output control
nstxout         =100            ;save coordinates every 0.2 ps
nstvout         =100            ;save velocities every 0.2 ps
nstenergy       =100            ;save energies every 0.2 ps
nstlog          =100            ;update log file every 0.2 ps
;Bond parameters
continuation    =no             ;first dynamics run
constraint_algorithm=lincs      ;holonomic constraints
constraints     =all-bonds      ;all bonds (even heavy atom-H
bonds)constraind
lincs_iter      =1              ;accuracy of LINCS
lincs_order     =4              ;also related to accuracy
;Neighborhood searching
ns_type         =grid           ;search neighboring grid cells
nstlist         =5              ;10 fs
rlist           =1.0            ;short-range neighborlist cutoff (in nm)
rcoulomb        =1.0            ;short-range electrostatic cutoff (in nm)
rvdw            =1.0            ;short-range van der Waals cutoff (in nm)
;Electrostatics
coulombtype     =PME            ;Particle Mesh Ewald for long-range
electrostat\
;ics
pme_order       =4              ;cubic interpolation
fourierspacing  =0.16           ;grid spacing for FFT
;Temperature coupling is on
tcoupl          =V-rescale      ;modified Berendsen thermostat
tc_grps         =SYSTEM   ;two coupling groups - more accurate
tau_t           =0.1    ;0.1              ;time constant, in ps
ref_t           =298    ;25               ;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
;Dispersion correction
DispCorr        =EnerPres       ;account for cut-off vdW scheme
;Velocity generation
gen_vel         =yes            ;assign velocities from Maxwell distribution
gen_temp        =25             ;temperature for Maxwell distribution
gen_seed        =-1             ;generate a random seed
;END

*NPT SCRIPT*
File Edit Options Buffers Tools Help
title           =Ethanol npt equilibration
;define         =-DPOSRES       ;position restrain the protein
;Run parameters
integrator      =md             ;leap-frog algorithm
nsteps          =50000          ;2 * 50000 = 100 ps
dt              =0.002          ;2fs
;Output control
nstxout         =100            ;save coordinates every 0.2 ps
nstvout         =100            ;save velocities every 0.2 ps
nstenergy       =100            ;save energies every 0.2 ps
nstlog          =100            ;update log file every 0.2 ps
;Bond parameters
continuation    =yes            ;Restarting after NVT
constraint_algorithm=lincs      ;holonomic constraints
constraints     =all-bonds      ;all bonds (even heavy atom-H
bonds)constraind
lincs_iter      =1              ;accuracy of LINCS
lincs_order     =4              ;also related to accuracy
;Neighborhood searching
ns_type         =grid           ;search neighboring grid cells
nstlist         =5              ;10 fs
rlist           =1.0            ;short-range neighborlist cutoff (in nm)
rcoulomb        =1.0            ;short-range electrostatic cutoff (in nm)
rvdw            =1.0            ;short-range van der Waals cutoff (in nm)
;Electrostatics
coulombtype     =PME            ;Particle Mesh Ewald for long-range
electrostat\
;ics
pme_order       =4              ;cubic interpolation
fourierspacing  =0.16           ;grid spacing for FFT
;Temperature coupling is on
tcoupl          =V-rescale      ;modified Berendsen thermostat
tc-grps         =SYSTEM   ;two coupling groups - more accurate
tau_t           =0.1;   0.1               ;time constant, in ps
ref_t           =298;   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 =4.5e-5         ;isothermal compressibility of h2O, 1/bar
;Periodic boundary conditions
pbc             =xyz            ; 3-D PBC
;Dispersion correction
DispCorr        =EnerPres       ;account for cut-off vdW scheme
;Velocity generation
gen_vel         =no             ;Velocity generation is off
;gen_temp       =25             ;temperature for Maxwell distribution
;gen_seed       =-1             ;generate a random seed
;END


-- 
*Best regards,*
**
*Fabian F. Casteblanco*
*Rutgers University -- *
*Chemical Engineering PhD Student*
*C: +908 917 0723*
*E:  **fabian.casteblanco at gmail.com* <fabian.casteblanco at gmail.com>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110411/49df402a/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ethanol.itp
Type: application/octet-stream
Size: 2600 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110411/49df402a/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ethanol.top
Type: application/octet-stream
Size: 713 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110411/49df402a/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: em.gro
Type: application/octet-stream
Size: 449 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110411/49df402a/attachment-0002.obj>


More information about the gromacs.org_gmx-users mailing list