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
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
