[gmx-users] pressure coupling

[Justin A. Lemkul]

Shuangxing Dai wrote:
> Thanks for the response.
> Yes, I know that large fluctuations will be expected. But why the 
> average was not even correct? Also, how can I reach some state ( some 
> pressre and temperature) and know/tell that the system is 
> in equilibrium if the average is not correct with large fluctuation? I 
> mean I do not know the criteria of Gromacs. Or maybe my procedure is wrong?
> 

You've only run 100 ps total equilibration.  For an extremely high temperature, 
I might expect it to take quite a bit longer, especially since you're dealing 
with a system that will be subject to large fluctuations (due to its size).

-Justin

> Message: 1
> Date: Mon, 05 Jul 2010 10:39:30 -0400
> From: "Justin A. Lemkul" <jalemkul at vt.edu <mailto:jalemkul at vt.edu>>
> Subject: Re: [gmx-users] pressure coupling
> To: Discussion list for GROMACS users <gmx-users at gromacs.org 
> <mailto:gmx-users at gromacs.org>>
> Message-ID: <4C31EEA2.7070105 at vt.edu <mailto:4C31EEA2.7070105 at vt.edu>>
> Content-Type: text/plain; charset=UTF-8; format=flowed
> 
> 
> 
> 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
>>
> Thanks,
> 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

========================================