<div dir="ltr">Hi, all,<div> 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? </div>
<div> My system has just 4000 atoms. This is my .mdp file:</div><div><br></div><div><div>;define = -DPOSRES</div><div>; RUN CONTROL PARAMETERS =</div><div>integrator = sd</div><div>; start time and timestep in ps =</div>
<div>tinit = 0</div><div>dt = 0.001</div><div>nsteps = 50000</div><div>; number of steps for center of mass motion removal =</div><div>nstcomm = 100</div>
<div>; OUTPUT CONTROL OPTIONS =</div><div>; Output frequency for coords (x), velocities (v) and forces (f) =</div><div>nstxout = 0</div><div>nstvout = 0</div><div>nstfout = 0</div>
<div>; Output frequency for energies to log file and energy file =</div><div>nstlog = 10</div><div>nstenergy = 10</div><div>; Output frequency and precision for xtc file =</div><div>nstxtcout = 10</div>
<div>xtc-precision = 1000</div><div>; NEIGHBORSEARCHING PARAMETERS =</div><div>; nblist update frequency =</div><div>nstlist = 20</div><div>; ns algorithm (simple or grid) =</div><div>ns_type = grid</div>
<div><br></div><div>;OPTIONS FOR PRESSURE COUPLING</div><div>Pcoupl = Parrinello-Rahman</div><div>pcoupltype = anisotropic</div><div>tau_p = 5</div><div>compressibility = 2.1645e-07 2.1645e-07 2.7322e-07 0 0 0</div>
<div>ref_p = 1 1 1 0 0 0</div><div>;OPTIONS FOR TEMPERATURE COUPLING</div><div>tc_grps = system</div><div>tau_t = 1</div><div>ref_t = 1000</div><div>
; OPTIONS FOR BONDS =</div><div>constraints = hbonds</div><div>; Type of constraint algorithm =</div><div>constraint-algorithm = Lincs</div><div>; Do not constrain the start configuration =</div><div>
unconstrained-start = no</div><div>; Relative tolerance of shake =</div><div>shake-tol = 0.0001</div><div>; Highest order in the expansion of the constraint coupling matrix =</div><div>lincs-order = 12</div>
<div>; Lincs will write a warning to the stderr if in one step a bond =</div><div>; rotates over more degrees than =</div><div>lincs-warnangle = 30</div><div>; Periodic boundary conditions: xyz, no, xy</div><div>
pbc = xyz</div><div>periodic_molecules = no</div><div>; nblist cut-off </div><div>rlist = 1</div><div><br></div><div>; OPTIONS FOR ELECTROSTATICS AND VDW</div><div>; Method for doing electrostatics</div>
<div>coulombtype = PME</div><div>rcoulomb = 1</div><div>; Method for doing Van der Waals</div><div>vdw-type = Cut-off</div><div>; cut-off lengths </div><div>rvdw = 1</div>
<div><br></div><div>; Spacing for the PME/PPPM FFT grid</div><div>fourierspacing = 0.12</div><div>; FFT grid size, when a value is 0 fourierspacing will be used</div><div>fourier_nx = 0</div><div>
fourier_ny = 0</div>
<div>fourier_nz = 0</div><div>; EWALD/PME/PPPM parameters</div><div>pme_order = 6</div><div>ewald_rtol = 1e-4</div><div>ewald_geometry = 3d</div><div>epsilon_surface = 0</div>
<div>optimize_fft = no</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div>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.</div>
<div><table border="0" cellpadding="0" cellspacing="0" width="375" style="border-collapse:
collapse">
<col width="75" span="5">
<tbody><tr height="13">
<td height="13" width="75">T</td>
<td width="75">p(bar)</td>
<td width="75">delta p</td>
<td width="75">T (K)</td>
<td width="75">delta T</td>
</tr>
<tr height="13">
<td height="13" align="right">0.01</td>
<td class="xl65" align="right">1.03E+00</td>
<td class="xl65" align="right">2.12E+00</td>
<td class="xl65" align="right">9.95E-03</td>
<td class="xl65" align="right">1.21E-04</td>
</tr>
<tr height="13">
<td height="13" align="right">100</td>
<td align="right">5.48359</td>
<td align="right">71.5565</td>
<td align="right">91.198</td>
<td align="right">1.25218</td>
</tr>
<tr height="13">
<td height="13" align="right">300</td>
<td align="right">2.62807</td>
<td align="right">162.926</td>
<td align="right">290.324</td>
<td align="right">3.36594</td>
</tr>
<tr height="13">
<td height="13" align="right">500</td>
<td align="right">16.6356</td>
<td align="right">1048.76</td>
<td align="right">492.616</td>
<td align="right">5.45231</td>
</tr>
<tr height="13">
<td height="13" align="right">800</td>
<td align="right">59.3245</td>
<td align="right">594.773</td>
<td align="right">799.314</td>
<td align="right">9.18876</td>
</tr>
<tr height="13">
<td height="13" align="right">1000</td>
<td align="right">80.2257</td>
<td align="right">626.789</td>
<td align="right">1007.23</td>
<td align="right">11.2159</td>
</tr>
</tbody></table>
</div><div>Thanks in advance,<br>Shuangxing Dai<br>
</div></div>