[gmx-users] berendsen P coupling and fluctuation properties

Justin A. Lemkul jalemkul at vt.edu
Wed Aug 17 07:38:16 CEST 2011



Elisabeth wrote:
> Thanks Justin for your thorough description. Is there any specific tool 
> in gromacs to plot histograms or you mean I follow the customary 
> plotting procedure?
> 

You have to write your own code from, e.g., g_energy output.

> I have one more question. I am now trying using PR pressure scheme with 
> v-rescale thermostat for NPT dynamics. Before dynamics run,  can I apply 
> berendsen barostat combined with v-rescale thermostat to equilibrate or 
> since I am using berendsen barostat, T coupling scheme has to be also 
> berendsen?
> 

You can use any combination you like during equilibration.

-Justin

> Thanks,
> 
> 
> On 16 August 2011 13:17, Justin A. Lemkul <jalemkul at vt.edu 
> <mailto:jalemkul at vt.edu>> wrote:
> 
> 
> 
>     Elisabeth wrote:f
> 
>         Hello,
> 
>         Thank you. I am looking at potential energies to calculate
>         vaporization heat. I wanted to know how the fact that berendsen
>         does not lead to correct ensemble is affecting total potential
>         energy of the system. I have an unclear image of "correct
>         ensemble". Does this mean whatever output I am getting from my
>         runs are unreliable?
> 
> 
>     The Berendsen coupling algorithms produce very narrow distributions
>     of temperature and pressure.  These do not correspond to the correct
>     distribution of a true statistical mechanical ensemble; plot a
>     histogram and you'll find that the results are shockingly different
>     between Berendsen and, say, Nose-Hoover for T.  Therefore, the
>     ensemble you're using to measure your properties is not NPT (or NVT,
>     in the case of thermostat only), it's something undefined and likely
>     not real.  Applying better thermostats and barostats that have been
>     shown to produce the correct distributions is more rigorously correct.
> 
>     The effects are usually not terribly noticeable unless you go
>     looking for them.  You'll find that the Berendsen algorithms quite
>     faithfully give you the target T and P that you desire (which does
>     make them very good for initial equilibration) and so you naturally
>     assume that everything is fine.  While that's all well and good in
>     one sense, the fluctuations of T (and thus velocities, and thus
>     kinetic energy) are wrong.  If KE does not fluctuate properly, then
>     neither does PE (since you're then affecting energy flow between PE
>     and KE, while total energy stays fixed, in theory).
> 
>     The larger point is that if you claim to apply an NVT (or NPT)
>     ensemble using Berendsen coupling, you're not correct, strictly
>     speaking.
> 
>     -Justin
> 
>         Appreciate any clarification.
> 
>         Thank you,
>         Best,
> 
>         On 15 August 2011 19:34, Mark Abraham <Mark.Abraham at anu.edu.au
>         <mailto:Mark.Abraham at anu.edu.au>
>         <mailto:Mark.Abraham at anu.edu.__au
>         <mailto:Mark.Abraham at anu.edu.au>>> wrote:
> 
>            On 16/08/2011 7:05 AM, Elisabeth wrote:
> 
>                Dear all,
> 
>                I noticed that applying Parrinello-Rahman (PR) pressure
>             coupling
>                even after equilibration with berendsen does not lead to
>             target
>                value for pressure when
> 
>                ;        Bonds
>                constraints             = none
> 
>                is used.
> 
> 
>            The use of constraints and the integration step size is linked.
>            Roughly speaking, no constraints should accompany a 0.5 fs time
>            step, H-bond constraints with 1fs and all constraints with 2fs.
>            Haphazard changes to .mdp files have all kinds of these
>         "gotchas".
> 
> 
>                I tried Berendsen for to get fixed pressure ( 50 bar) but
>             in the
>                next run PR even for long time is giving 52 bar. This is
>             the case
>                for other target pressures too.
> 
> 
>            You need to be sure to collect statistics only after
>         equilibration,
>            and consider whether the observed variation is consistent with
>            convergence to a given value.
> 
> 
> 
>                So this made me select Berenden which is giving target
>             pressure
>                values but my concern is whether my results are reliable
>             because
>                BR does not give the exact ensemble as PR. I read
>             somewhere on the
>                list that fluctuation properties can not be calculated
>             when BR is
>                used.  What does "fluctuation property" mean?
> 
> 
>            BR does not produce the correct ensemble. I forget the
>         details about
>            why, but there are references in the T-coupling section of the
>            manual you should consider.
> 
> 
>                Does this mean that any property calculated form
>             fluctuations of
>                some other quantity can not be obtained>? like heat
>             capacity which
>                is defined based on enthalpy fluctuations?
> 
> 
>            IIRC, yes.
> 
> 
>                I am interested in potential energy terms (g_energy
>             bonded/non
>                boned terms)
> 
> 
>            I routinely struggle to see why people think they can learn
>         anything
>            from these.
> 
> 
>                and structural properties like rdf for a number of polymer
>                molecules, system size around 3000 atoms.
> 
>                Thank you for your comments.
>                Best,
> 
>                constraints             = none              
>                ;        Run control                       integrator    
>                  =  md                    dt                  =  0.001  
>                             nsteps              =  5000000            
>             nstcomm             =  100          
> 
>                nstenergy           =  100                    nstxout    
>                     =  100                     nstlist             =  10
>                               ns_type             =  grid              
>                coulombtype         =  Shift                      
>              vdw-type            =  Shift                
>             rcoulomb-switch     =  0                       rvdw-switch  
>                   =  0.9 ;0              
>                ;        Cut-offs
>                rlist               =  1.25                     rcoulomb
>                        =  1.0                rvdw                =  1.0
>                          
>                Tcoupl              =  v-rescale                    
>             tc-grps             =  System                 tau_t        
>                   =  0.1                ref_t               =  300    
>                Pcoupl              =  berendsen           Pcoupltype    
>                  =  isotropic                   tau_p               =  1
>                               compressibility     =  3.5e-5            
>              ref_p               =  10                
>                            gen_vel             =  no
>                gen_temp            =  300.0                   gen_seed  
>                      =  173529                
> 
> 
> 
> 
> 
> 
> 
> 
>            --
>            gmx-users mailing list    gmx-users at gromacs.org
>         <mailto:gmx-users at gromacs.org>
>            <mailto:gmx-users at gromacs.org <mailto:gmx-users at gromacs.org>>
> 
>            http://lists.gromacs.org/__mailman/listinfo/gmx-users
>         <http://lists.gromacs.org/mailman/listinfo/gmx-users>
>            Please search the archive at
>            http://www.gromacs.org/__Support/Mailing_Lists/Search
>         <http://www.gromacs.org/Support/Mailing_Lists/Search> before
>         posting!
>            Please don't post (un)subscribe requests to the list. Use the
>            www interface or send it to gmx-users-request at gromacs.org
>         <mailto:gmx-users-request at gromacs.org>
>            <mailto:gmx-users-request at __gromacs.org
>         <mailto:gmx-users-request at gromacs.org>>.
> 
>            Can't post? Read
>         http://www.gromacs.org/__Support/Mailing_Lists
>         <http://www.gromacs.org/Support/Mailing_Lists>
> 
> 
> 
>     -- 
>     ==============================__==========
> 
>     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
>     <tel:%28540%29%20231-9080>
>     http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin
>     <http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
> 
>     ==============================__==========
>     -- 
>     gmx-users mailing list    gmx-users at gromacs.org
>     <mailto:gmx-users at gromacs.org>
>     http://lists.gromacs.org/__mailman/listinfo/gmx-users
>     <http://lists.gromacs.org/mailman/listinfo/gmx-users>
>     Please search the archive at
>     http://www.gromacs.org/__Support/Mailing_Lists/Search
>     <http://www.gromacs.org/Support/Mailing_Lists/Search> before posting!
>     Please don't post (un)subscribe requests to the list. Use the www
>     interface or send it to gmx-users-request at gromacs.org
>     <mailto:gmx-users-request at gromacs.org>.
>     Can't post? Read http://www.gromacs.org/__Support/Mailing_Lists
>     <http://www.gromacs.org/Support/Mailing_Lists>
> 
> 

-- 
========================================

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

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



More information about the gromacs.org_gmx-users mailing list