[gmx-users] Pressure coupling problem
Fabian Casteblanco
fabian.casteblanco at gmail.com
Mon Apr 11 23:20:31 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
--
Best regards,
Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
C: +908 917 0723
E: fabian.casteblanco at gmail.com
-------------- next part --------------
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
-------------- 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/5c87edcd/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/5c87edcd/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/5c87edcd/attachment-0002.obj>
More information about the gromacs.org_gmx-users
mailing list