[gmx-users] Pressure coupling problem
Peter C. Lai
pcl at uab.edu
Mon Apr 11 23:41:14 CEST 2011
So your density graph looks stabilized? I also tend to look for changes in
box x, y, z as well since the scale of their changes is easier to track.
Sometimes it helps to look at the error vs. rmsd vs total drift statistics as
well for such parameters that are easier to track - again if density shows
stability for 0-150ps then check your boxes, it might be shrinking or
growing due to the pressure perturbation, and you can use the average rate
of change of those to find your equilibration point instead of trying to
do something with a -300 to 300 bar pressure smear or whatnot.
On 2011-04-11 04:17:59PM -0500, Justin A. Lemkul wrote:
>
>
> Fabian Casteblanco wrote:
> > 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.
> >
>
> Your equilibration period (100-150 ps) is rather short, and the systematic
> increase suggests that you're simply not equilibrated yet.
>
> Also bear in mind that a quantity that is prone to pressure fluctuations in the
> hundreds to thousands can only be so accurate. There was a very thorough
> discussion about the statistical significance of pressure values that are not
> equal to ref_p just some time ago. You may want to look through the archives to
> find this discussion.
>
> -Justin
>
> > 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* <mailto:fabian.casteblanco at gmail.com>
> >
>
> --
> ========================================
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> ========================================
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
===============================================================
Peter C. Lai | University of Alabama-Birmingham
Programmer/Analyst | BEC 257
Genetics, Div. of Research | 1150 10th Avenue South
pcl at uab.edu | Birmingham AL 35294-4461
(205) 690-0808 |
===============================================================
More information about the gromacs.org_gmx-users
mailing list