Shuangxing Dai wrote:
Well, the final average pressure is not correct:
* **<====== ############### ==>*
* **<==== A V E R A G E S ====>*
* **<== ############### ======>*
*
*
* Energies (kJ/mol)*
* Buck.ham (SR) Coulomb (SR) Coul. recip. Potential Kinetic
En.*
* 1.60633e+06 -5.69243e+06 -3.45003e+06 -7.53612e+06
1.80719e+04*
* Total Energy Temperature Pressure (bar)*
* -7.51805e+06 3.62347e+02 1.11758e+04*
*
*
*
*
* Pressure (bar)*
* 1.21051e+04 -2.94017e-01 1.30456e+00*
* -2.94017e-01 9.03597e+03 -8.10850e-01*
* 1.30456e+00 -8.10850e-01 1.23863e+04*
*
*
*
*
* **<====== ############################### ==>*
* **<==== R M S - F L U C T U A T I O N S ====>*
* **<== ############################### ======>*
*
*
* Energies (kJ/mol)*
* Buck.ham (SR) Coulomb (SR) Coul. recip. Potential Kinetic
En.*
* 3.92376e+05 4.00326e+05 4.71365e+01 3.02978e+04
2.35039e+03*
* Total Energy Temperature Pressure (bar)*
* 3.01091e+04 4.71259e+01 1.71117e+05*
*
*
*
*
* Pressure (bar)*
* 1.82323e+05 4.91513e+02 3.07107e+02*
* 4.91513e+02 2.02402e+05 3.13857e+02*
* 3.07107e+02 3.13857e+02 1.30179e+05*
*
*
*Also the fluctuation is very large. I think 200ps is long enough. But
still I cannot find the equilibrium.*
Then 200 ps is clearly not enough, or your entire procedure is not correct.
Note too that your temperature (362 K) is not at all close to what you want (300
K). You may be setting tau_t incorrectly; check the archives. This value, when
using the sd integrator, is not the same as the normal interpretation of
relaxation time, so instead of 0.1, you may want 1.0 or larger.
*I do not know whether the parameters are wrong or I ran the anisotropic
pressure coupling wrong. Is *Parrinello-Rahman pressure coupling
sensitive to the initial pressure? I mean if the initial pressure ( in
order of GPa), object pressure is 1 bar ( 1e5 Pa).
If your system is far from equilibrium, it may be beneficial to relax it with a
weak coupling method (i.e., Berendsen) first. Otherwise, a system that is not
close to equilibrium may oscillate very wildly and not converge to the desired
value.
-Justin
*
*
Thanks,
Shuangxing Dai
On Mon, Jun 14, 2010 at 8:50 PM, <gmx-users-requ...@gromacs.org
<mailto:gmx-users-requ...@gromacs.org>> wrote:
Message: 3
Date: Mon, 14 Jun 2010 19:18:14 -0400
From: "Justin A. Lemkul" <jalem...@vt.edu <mailto:jalem...@vt.edu>>
Subject: Re: [gmx-users] anisotropic pressure coupling
To: Discussion list for GROMACS users <gmx-users@gromacs.org
<mailto:gmx-users@gromacs.org>>
Message-ID: <4c16b8b6.4010...@vt.edu <mailto:4c16b8b6.4010...@vt.edu>>
Content-Type: text/plain; charset=UTF-8; format=flowed
Shuangxing Dai wrote:
> Hi all:
> I was using Parrinello-Rahman for pressure coupling and here
is the
> mdp file:
>
> ; RUN CONTROL PARAMETERS =
> integrator = sd
> ; start time and timestep in ps =
> tinit = 0
> dt = 0.001
> nsteps = 200000
> ; number of steps for center of mass motion removal =
> nstcomm = 100
> ; NEIGHBORSEARCHING PARAMETERS =
> ; nblist update frequency =
> nstlist = 50
> ; ns algorithm (simple or grid) =
> ns_type = grid
>
> ;OPTIONS FOR PRESSURE COUPLING
> Pcoupl = Parrinello-Rahman
> pcoupltype = anisotropic
> tau_p = 1
> compressibility = 2.1645e-09 2.1645e-09 2.7322e-09 0 0 0
> ref_p = 1 1 1 0 0 0
> ;OPTIONS FOR TEMPERATURE COUPLING
> tc_grps = system
> tau_t = 0.1
> ref_t = 300
> ; 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
>
> It seems that the temperature coupling is fast and fluctuation
is less
> that 1%. But the pressure coupling fails. The pressure fluctuates
from
> +1e3 bar to -1e3 bar. Anyone can help?
>
I don't call that failure. Depending on the size of your system,
you can
probably expect very wide fluctuations in the pressure. Is the
average correct?
http://www.gromacs.org/Documentation/Terminology/Pressure
-Justin
>
> 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 <http://vt.edu> | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
========================================
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