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