[gmx-users] pressure coupling
Justin A. Lemkul
jalemkul at vt.edu
Mon Jul 5 16:39:30 CEST 2010
Shuangxing Dai wrote:
> Hi, all,
> I am trying to do anisotropic coupling using Parrinello-Rahman. I
> want to reach the equilibrium state at 1000K, 0.1MPa( 1 bar). However,
> I find that the fluctuation is so large ( p=80.22 bar, T=1007.23 K,
> fluctuation is 626 bar for pressure and 11.21 K for temperature.) I run
> energy minimization first, then Berenderson pressure coupling for 50 ps,
> finally Parrinello-Rahman pressure coupling for 50 ps. Anyone has
> experience on anisotropic coupling and know why?
> My system has just 4000 atoms. This is my .mdp file:
>
For a small system, large fluctations in the pressure are to be expected:
http://www.gromacs.org/Documentation/Terminology/Pressure
-Justin
> ;define = -DPOSRES
> ; RUN CONTROL PARAMETERS =
> integrator = sd
> ; start time and timestep in ps =
> tinit = 0
> dt = 0.001
> nsteps = 50000
> ; number of steps for center of mass motion removal =
> nstcomm = 100
> ; OUTPUT CONTROL OPTIONS =
> ; Output frequency for coords (x), velocities (v) and forces (f) =
> nstxout = 0
> nstvout = 0
> nstfout = 0
> ; Output frequency for energies to log file and energy file =
> nstlog = 10
> nstenergy = 10
> ; Output frequency and precision for xtc file =
> nstxtcout = 10
> xtc-precision = 1000
> ; NEIGHBORSEARCHING PARAMETERS =
> ; nblist update frequency =
> nstlist = 20
> ; ns algorithm (simple or grid) =
> ns_type = grid
>
> ;OPTIONS FOR PRESSURE COUPLING
> Pcoupl = Parrinello-Rahman
> pcoupltype = anisotropic
> tau_p = 5
> compressibility = 2.1645e-07 2.1645e-07 2.7322e-07 0 0 0
> ref_p = 1 1 1 0 0 0
> ;OPTIONS FOR TEMPERATURE COUPLING
> tc_grps = system
> tau_t = 1
> ref_t = 1000
> ; OPTIONS FOR BONDS =
> constraints = hbonds
> ; Type of constraint algorithm =
> constraint-algorithm = Lincs
> ; Do not constrain the start configuration =
> unconstrained-start = no
> ; Relative tolerance of shake =
> shake-tol = 0.0001
> ; Highest order in the expansion of the constraint coupling matrix =
> lincs-order = 12
> ; Lincs will write a warning to the stderr if in one step a bond =
> ; rotates over more degrees than =
> lincs-warnangle = 30
> ; Periodic boundary conditions: xyz, no, xy
> pbc = xyz
> periodic_molecules = no
> ; nblist cut-off
> rlist = 1
>
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype = PME
> rcoulomb = 1
> ; Method for doing Van der Waals
> vdw-type = Cut-off
> ; cut-off lengths
> rvdw = 1
>
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> ; EWALD/PME/PPPM parameters
> pme_order = 6
> ewald_rtol = 1e-4
> ewald_geometry = 3d
> epsilon_surface = 0
> optimize_fft = no
>
>
>
>
> I used the same procedure to reach different state ( like 0.01 K, 100K,
> 300K, 500K, 800K, 1000K). Only the 0.01 K got the correct average
> pressure and small fluctuation.
> T p(bar) delta p T (K) delta T
> 0.01 1.03E+00 2.12E+00 9.95E-03 1.21E-04
> 100 5.48359 71.5565 91.198 1.25218
> 300 2.62807 162.926 290.324 3.36594
> 500 16.6356 1048.76 492.616 5.45231
> 800 59.3245 594.773 799.314 9.18876
> 1000 80.2257 626.789 1007.23 11.2159
>
> Thanks in advance,
> Shuangxing Dai
>
--
========================================
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
========================================
More information about the gromacs.org_gmx-users
mailing list