[gmx-users] Pressure coupling problem
Fabian Casteblanco
fabian.casteblanco at gmail.com
Fri Apr 15 21:15:11 CEST 2011
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.
Thanks.
On Mon, Apr 11, 2011 at 4:55 PM, Fabian Casteblanco <
fabian.casteblanco at gmail.com> 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.
>
> 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>
>
>
--
*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/20110415/45815a16/attachment.html>
More information about the gromacs.org_gmx-users
mailing list