[gmx-users] pressure coupling

[Shuangxing Dai]

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:

;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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100705/78d8aef7/attachment.html>