[gmx-users] NPT simulation with parinello-rahman pressure coupling

[Justin A. Lemkul]

Shyno Mathew wrote:
> Dear gromacs users,
> I am having some difficulty in setting up an NPT simulation, following 
> are the issues
> 
> 1. With Parinello-Rahman coupling, I am getting the error:
> Fatal error:
> The X-size of the box (4.224258) times the triclinic skew factor 
> (1.000000) is smaller than the number of DD cells (5) times the smallest 
> allowed cell size (0.844801)
> As mentioned in the archives, initially I thought my system is too small 
> to split between different processors. When I load the trajectory I can 
> see that the box is shrinking very fast in the x direction.

It may be possible that your system can be split over the chosen number of 
processors, but the margin of error is very small.  If the system is not well 
equilibrated, the box may be changing quite fast.

>  I have performed the same run with Berendsen thermostat and it ran 
> fine. So changing the pcoupling method can cause such a change?

Yes, absolutely.  The Berendsen method leads to faster convergence, while 
Parrinello-Rahman often fluctuates to a larger extent before settling around the 
target pressure values.  It is not suitable for equilibration, generally.  Have 
you done prior equilibration for this system?

> With the berendsen p coupling I have done a benchmark for the system and 
> I am using the optimum number of processors. Will the change in 
> pcoupling affects benchmark as well?
> 

The differences should be minimal.  I doubt there is a significant difference in 
calculation times for the various pressure coupling algorithms, at least when 
considering the other truly time-consuming processes in the MD algorithm.

> some parameters in .mdp file:
> 
> integrator               = sd
> pbc                      = xyz
> rlist                    = 1.2
> coulombtype              = PME
> fourierspacing           = 0.12
> pme_order                = 4
> ewald_rtol               = 1e-05
> optimize_fft             = no
> epsilon_surface          = 0
> ewald_geometry           = 3d
> rcoulomb                 = 1.2
> vdwtype                  = Shift
> rvdw_switch              = 0.9
> rvdw                     = 1.0
> epsilon_r                = 1
> DispCorr                 = EnerPres
> tc-grps                  = system
> tau_t                    = 2.0
> ref_t                    = 310
> pcoupl                   = Parrinello-Rahman
> pcoupltype               = anisotropic
> nstpcouple               = -1
> tau_p                    = 1.0
> 
> 2. I changed some of the above mentioned parameters, for example vdwtype 
> to cut-off and DispCorr to energy as shown below:
> 
> rlist                    = 1.0
> rcoulomb                 = 1.0
> vdwtype                  = Cut-off
> rvdw                     = 1.0
> epsilon_r                = 1
> DispCorr                 = Ener
> 
> Again I get the same error:
> Fatal error:
> The X-size of the box (5.000070) times the triclinic skew factor 
> (1.000000) is smaller than the number of DD cells (5) times the smallest 
> allowed cell size (1.000000)
> 

The fact that this error persists despite other changes suggests to me that 
there is simply inherently wrong with the system, i.e. inadequate minimization 
or equilibration, unstable topology, etc.

Also note that you're not providing several important parameters in your 
pressure coupling section, like reference pressures and compressibilities.  With 
anisotropic pressure coupling, you should make sure that the off-diagonal 
components are not causing compression or distortion of your system, as the 
manual advises.

> Also I thought for obtaining more accurate trajectory, I should apply 
> dispersion corrections for both pressure and energy, is this true?
> 

With a plain cutoff, dispersion correction is often advisable.

> 3. In the mailing list I found that tau_p has to be greater than tau_t 
> to avoid larger fluctuations  and so changed the parameters as shown 
> below, everything remains same as in question2 except the following:
> 
> tau_t                    = 2.0
> ref_t                    = 310
> tau_p                    = 3.0
> 

When citing information found in the archive, providing a link is advisable.  I 
won't categorically say whether that advice is right or wrong, though it is 
generally correct given the nature of barostats vs. thermostats.  Do not 
necessarily believe that tau_p > tau_t will lead to greater stability.  The 
simple fact that pressure oscillates more widely than temperature is the main 
reason for this need, it's simply fundamental.

> 
> Fatal error:
> One of the box vectors has become shorter than twice the cut-off length 
> or box_yy-|box_zy| or box_zz has become smaller than the cut-off.
> 

This is the same problem as before, but by chance a different error was 
triggered based on the algorithm that failed first.  Your system is blowing up. 
  See the explanation and advice here:

http://www.gromacs.org/Documentation/Terminology/Blowing_Up

-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

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