Re: [gmx-users] Lambda Weights from Expanded Ensemble Code

[Michael Shirts]

The simulated tempering options haven't been as well tested as the
hamiltonian expanded ensemble version.  The weights SHOULD be showing up in
the column that says -nan, but clearly they aren't.  If you file a redmine
issue, I may be able to take a look, but it might take a while to address.

On Wed, Jan 15, 2020 at 8:52 PM Daniel Kozuch  wrote:

> Hello,
>
> I am interested in using simulated tempering in GROMACS (2019.5) under the
> expanded ensemble options. Is there a way to monitor the ensemble weights
> as the simulation progresses? I think in theory they are supposed to be
> printed out in the log file, but it is only printing 0, -nan, and inf:
>
> MC-lambda information
>   N  Temp.(K)Count   G(in kT)  dG(in kT)
> ...
>  36  359.105  118   -nan   -nan <<
>  37  366.880   96   -nan   -nan
>  38  374.852  107   -nan   -nan
>  39  383.026  129   -nan   -nan
>  40  391.407  166   -nan   -nan
>  41  400.000  199   -nan0.0
>
> Here are my relevant mdp settings:
> simulated-tempering = yes
> nstexpanded = 500
> simulated-tempering-scaling = exponential
> lmc-stats   = metropolis-transition
> lmc-move= metropolis
>
> Any suggestions?
>
> Best,
> Dan Kozuch
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Implementing Hamiltonian replica exchange simulations

[Michael Shirts]

Hamilton replica exchange is implemented in GROMACS (unless it's been
broken in some way).  See
http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy
for
an older tutorial, but it should be relatively straightforward to adjust to
current call patterns.

On Tue, Oct 15, 2019 at 8:18 AM Searle Duay  wrote:

> Hi Matthew,
>
> Thank you very much! I appreciate it.
>
> Best,
> Searle
>
> On Tue, Oct 15, 2019 at 8:44 AM Matthew Fisher <
> matthew.fis...@stcatz.ox.ac.uk> wrote:
>
> > It's possible to implement HREX-MD with the Plumed patch. A tutorial for
> > it is here - https://www.plumed.org/doc-v2.5/user-doc/html/hrex.html
> >
> > If you're interested in the scientific basis behind the methodology then
> > take a look here
> > https://www.tandfonline.com/doi/full/10.1080/00268976.2013.824126
> >
> > Best,
> > Matthew
> >
> > Matthew Fisher
> > DPhil Candidate in Inorganic Chemistry (L.L.Wong Group)
> >
> > Inorganic Chemistry Laboratory
> > South Parks Road
> > Oxford
> > OX1 3QR
> > UNITED KINGDOM
> >
> > Tel (Office): +44 (0)1865 272679
> >
> > 
> > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> > gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Searle
> > Duay 
> > Sent: 15 October 2019 12:02
> > To: gmx-us...@gromacs.org 
> > Subject: Re: [gmx-users] Implementing Hamiltonian replica exchange
> > simulations
> >
> > Thank you Prof. Lemkul. I haven't seen any discussion either on doing
> HREX.
> >
> > -Searle
> >
> > On Tue, Oct 15, 2019 at 5:18 AM Justin Lemkul  wrote:
> >
> > >
> > >
> > > On 10/14/19 3:48 PM, Searle Duay wrote:
> > > > Hi everyone,
> > > >
> > > > I want to do Hamiltonian replica exchange simulations on my helical
> > > peptide
> > > > transitioning from transmembrane state to surface-bound state. My
> > > > collective variable is the angle between the z-axis and the peptide
> > > helical
> > > > axis. I have 30 windows with different initial configurations. Here
> is
> > an
> > > > example of mdp file for a window:
> > > >
> > > > integrator  = md
> > > > dt  = 0.002
> > > > nsteps  = 5000
> > > > nstlog  = 50
> > > > nstxout = 1000
> > > > nstvout = 1000
> > > > nstfout = 1000
> > > > nstcalcenergy   = 100
> > > > nstenergy   = 1000
> > > > ;
> > > > cutoff-scheme   = Verlet
> > > > nstlist = 20
> > > > rlist   = 1.2
> > > > coulombtype = pme
> > > > rcoulomb= 1.2
> > > > vdwtype = Cut-off
> > > > vdw-modifier= Force-switch
> > > > rvdw_switch = 1.0
> > > > rvdw= 1.2
> > > > ;
> > > > tcoupl  = Nose-Hoover
> > > > tc_grps = PROT   MEMB   SOL_ION
> > > > tau_t   = 1.01.01.0
> > > > ref_t   = 300 300 300
> > > > ;
> > > > pcoupl  = Parrinello-Rahman
> > > > pcoupltype  = semiisotropic
> > > > tau_p   = 5.0
> > > > compressibility = 4.5e-5  4.5e-5
> > > > ref_p   = 1.0 1.0
> > > > ;
> > > > constraints = h-bonds
> > > > constraint_algorithm= LINCS
> > > > continuation= yes
> > > > ;
> > > > nstcomm = 100
> > > > comm_mode   = linear
> > > > comm_grps   = PROT   MEMB   SOL_ION
> > > > ;
> > > > refcoord_scaling= com
> > > > ;
> > > > pull= yes
> > > > pull_ncoords= 1
> > > > pull_ngroups= 2
> > > > pull-group1-name= n_com
> > > > pull-group2-name= c_com
> > > > pull-coord1-type= umbrella
> > > > pull-coord1-geometry= angle-axis
> > > > pull-coord1-groups  = 1 2
> > > > pull-coord1-vec = 0 0 -1
> > > > pull-coord1-start   = no
> > > > pull-coord1-init= 3.0
> > > > pull-coord1-k   = 1000.0
> > > >
> > > > The mdp file for each window varies only in the pull-coord1-init
> > > parameter
> > > > (ranging from 3 to 90 degrees, in 3-degree intervals). Then, I grompp
> > to
> > > > get 30 different tpr files.
> > > >
> > > > This is the command that I use for running the simulations:
> > > >
> > > > srun --mpi=pmi2 gmx_mpi mdrun -s window.tpr -replex 1000 -multidir
> > angle0
> > > > angle1 angle2 angle3 angle4 angle5 angle6 angle7 angle8 angle9
> angle10
> > > > angle11 angle12 angle13 angle14 angle15 angle16 angle17 angle18
> angle19
> > > > angle20 angle21 angle22 angle23 angle24 angle25 angle26 angle27
> angle28
> > > > angle29 -maxh 119.5
> > > >
> > > > And I'm getting the following error:
> > > > The properties of the 30 systems are all the same, there is nothing
> to
> > > > exchange
> > > >
> > > > I'm wondering where the error comes from. My systems should be
> > different
> > > > such that the target angle is 

Re: [gmx-users] Pressure coupling in expanded ensemble simulations

[Michael Shirts]

Ah, sorry, I thought there was more information coming.

Have you considered just using temperature replica exchange?  It's not that
much less efficient, and is easier to deal with.  Replica exchange should
be working with NPT (as long as you use Parrinello-Rahman and a reasonable
temperature control algorithm).  The size scaling is about the same; i.e.
if you need a lot of replicas, you will also need a lot of expanded
ensemble intermediates.




On Fri, May 17, 2019 at 8:06 AM Gregory Man Kai Poon  wrote:

> Hi Michael,
>
> I am just following up on your thoughts on how carrying out expanded
> ensemble at NVT and converting back to NPT on the mailing list.  Again I
> appreciate your advice in this area.
>
> Best wishes,
>
> Gregory
>
>
> On 5/8/2019 12:01 PM, Michael Shirts wrote:
>
> Yeah, this is an unfortunately place in the code where not all combinations
> work - very long story.  Hopefully this will be working better in 2020.
>
> What I would recommend is, if possible, performing the expanded ensemble
> simulation at NVT.  Everything should work fine there (paper coming out
> hopefully soon comparing a bunch of free energy methods).  Once can always
> correct the free energy at lambda=1 from NVT to NPT.  I Can fill in the
> details.
>
> You do NOT want to do Berendsen for NPT when running expanded ensemble.
> The results will be incorrect (as I have learned by sad experience_
>
> On Wed, May 8, 2019 at 8:14 AM Gregory Man Kai Poon  ><mailto:gp...@gsu.edu> wrote:
>
>
>
> Hi all:
>
> We are interested to do expanded ensemble simulations (such as simulated
> tempering) on GROMACS.  Extensive fiddling with the settings and
> googling on other people's experience suggests that these simulations
> must use the md-vv integrator, which in turn is compatible with
> Berendsen or MTTK coupling for pressure.  However, MTTK does not work
> with constraints, which are needed for the forcefields.  Berendsen can
> handle constraints but is not recommended for preserving thermodynamic
> ensembles.  Any ideas on how one should proceed?
>
> Many thanks for your thoughts.
>
> Gregory
>
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
>
> https://nam03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.gromacs.org%2FSupport%2FMailing_Lists%2FGMX-Users_Listdata=02%7C01%7Cgpoon%40gsu.edu%7C0691e8847a2b444b3a8208d6d3ce834d%7C515ad73d8d5e4169895c9789dc742a70%7C0%7C1%7C636929281266198778sdata=6gqdkFvDH2oX6U17pgPfcUvs34mmPP6OCyPnUQbgwKE%3Dreserved=0
> before
> posting!
>
> * Can't post? Read
> https://nam03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.gromacs.org%2FSupport%2FMailing_Listsdata=02%7C01%7Cgpoon%40gsu.edu%7C0691e8847a2b444b3a8208d6d3ce834d%7C515ad73d8d5e4169895c9789dc742a70%7C0%7C1%7C636929281266198778sdata=xJsYvSXLGH4aeubUlxygOvzszCraVOEywg%2BwFcnAIQA%3Dreserved=0
>
> * For (un)subscribe requests visit
>
> https://nam03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmaillist.sys.kth.se%2Fmailman%2Flistinfo%2Fgromacs.org_gmx-usersdata=02%7C01%7Cgpoon%40gsu.edu%7C0691e8847a2b444b3a8208d6d3ce834d%7C515ad73d8d5e4169895c9789dc742a70%7C0%7C1%7C636929281266198778sdata=1zXWxg16PNlz5%2FFgFsmfO79eangf1sRB%2BkzpcuT3OFU%3Dreserved=0
> or
> send a mail to gmx-users-requ...@gromacs.org gmx-users-requ...@gromacs.org>.
>
>
> --
>
> Gregory M. K. Poon, PhD, RPh
> Associate Professor
> Departments of Chemistry and Nutrition | Georgia State University
> NSC 414/415/416 | 50 Decatur St. SE, Atlanta, GA 30302
> P.O. Box 3965 | Atlanta, GA 30303
> Ph (404) 413-5491 | gp...@gsu.edu<mailto:gp...@gsu.edu>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] TIP4P molecules stuck together

[Michael Shirts]

If you have "unprotected" electrostatic sites (i.e. with nonzero repulsive
terms directly on top of the charge), then there will always be some
configurations with essentially infinitely negative energy.  For vanilla
MD, these are essentially impossible to reach kinetically at any reasonable
temperature or reasonable timestep, because the R^-12 of the neighboring
atoms creates such an enormously large barrier.  But with certain
accelerated sampling approaches, you can skip over the barrier, and access
these sites, which will either blow up the simulation, or get stuck
forever.  If you cap the forces, then weirder things will happen.  Is your
cap smoothly varying?  If not, then your dynamics on hitting the cap will
be unphysical.  How are the forces propagated into the energies (if grad U
=/= F, then weird non-newtonian physics will also happen).

> the forces on each atom are massive (on the order of 10^7).

What are the energies? Are they lower or higher than zero?

On Thu, May 9, 2019 at 8:43 AM John Whittaker <
johnwhitt...@zedat.fu-berlin.de> wrote:

> Hi all,
>
> I have a rather strange question that I hope someone can shed some light
> on.
>
> Before I begin, I want to note that I am pioneering some new developments
> of the Adaptive Resolution Simulation technique
> (https://doi.org/10.1002/adts.201900014), so the simulations/techniques I
> am performing/implementing are fairly non-standard with respect to normal
> atomistic simulations.
>
> With that in mind, I am simulating a box of TIP4P water and calculating
> structural/static properties. My simulations utilize a force-cap of 2000
> kJ/(mol nm) at each time step - i.e., when the force on an atom is larger
> than +/- 2000, the force is automatically normalized to +/- 2000 to
> prevent explosive forces due to atomic overlaps.
>
> For the most part, this works for the purposes of my simulations but I
> have observed some water molecules "sticking" together in the
> configuration shown here:
>
> https://www.dropbox.com/s/p5rkximspp25flf/tip4pDimer.jpg?dl=0
>
> with a corresponding O-H radial distribution function (unnormalized) shown
> here:
>
> https://www.dropbox.com/s/ez56db4qggv1iii/rdf_OH_long.jpg?dl=0
>
> where there is a clear (albeit, small) probability of finding a hydrogen
> atom an extremely short distance from an oxygen.
>
> The molecules travel together like this for several ps and then, for
> seemingly no reason, split apart and carry on perfectly fine for the rest
> of the simulation.
>
> I have performed a single-point energy calculation on this configuration
> in vacuum and have found, as one would expect, the forces on each atom are
> massive (on the order of 10^7). Yet, the molecules do not repel and seem
> to prefer this configuration for a short time.
>
> I have a feeling that this configuration is allowed when the forces are
> normalized to 2000 and the molecules become trapped there.
>
> I am wondering if anyone may have some experience with TIP4P water
> molecules taking on unphysical configurations for non-negligible times. I
> have not tried this same simulation using TIP3P yet, so I'm unsure if it
> has something to do with electrostatic interactions with the virtual site,
> but I will test this tomorrow.
>
> Thank you for any information/speculation/guesses as to why this is
> happening.
>
> - John
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure coupling in expanded ensemble simulations

[Michael Shirts]

Yeah, this is an unfortunately place in the code where not all combinations
work - very long story.  Hopefully this will be working better in 2020.

What I would recommend is, if possible, performing the expanded ensemble
simulation at NVT.  Everything should work fine there (paper coming out
hopefully soon comparing a bunch of free energy methods).  Once can always
correct the free energy at lambda=1 from NVT to NPT.  I Can fill in the
details.

You do NOT want to do Berendsen for NPT when running expanded ensemble.
The results will be incorrect (as I have learned by sad experience_

On Wed, May 8, 2019 at 8:14 AM Gregory Man Kai Poon  wrote:

> Hi all:
>
> We are interested to do expanded ensemble simulations (such as simulated
> tempering) on GROMACS.  Extensive fiddling with the settings and
> googling on other people's experience suggests that these simulations
> must use the md-vv integrator, which in turn is compatible with
> Berendsen or MTTK coupling for pressure.  However, MTTK does not work
> with constraints, which are needed for the forcefields.  Berendsen can
> handle constraints but is not recommended for preserving thermodynamic
> ensembles.  Any ideas on how one should proceed?
>
> Many thanks for your thoughts.
>
> Gregory
>
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Performance of GROMACS on GPU's on ORNL Titan?

[Michael Shirts]

Does anyone have experience running GROMACS on GPU's on Oak Ridge National
Labs Titan or Summit machines, especially parallelization over multiple
GPUs? I'm looking at applying for allocations there, and am interested in
experiences that people have had. We're probably mostly looking at systems
in the 100-200K atoms range, but we need to get to long timescales
(multiple microseconds, at least) for some of the phenomena we are looking
at.

Thanks!
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Bennet error in FEP-calculations for charged ligands.

[Michael Shirts]

As free energies get larger, then the error gets less accurate. So if it is
reporting 18.51 +/-  2.95 for one of the intervals, that likely suggest
there is very little overlap in that area.
A total free energy difference of  in -6.12 +/-  3.19 for the
transformation indicates that the result is not very certain; you're within
two standard deviations, and again, the Bennett error formula is inaccurate
for larger error. I would suggest adding 1-2 more intermediates between
those points.

https://github.com/MobleyLab/alchemical-analysis provides some of the tools
to check overlap.

I would also suggest doing a multi stage transformation to increase overlap
- for atoms that are disappearing or appearing that are charged, turn off
charges off, change vdw, then turn charges back on. In some cases,
especially when the atoms involved have different sizes, sc-coul can lead
to some semi-pathological cases when vdw softcore and coul softcore result
in some very low/high potentials at intermediates.

On Wed, Jan 23, 2019 at 8:21 AM Artem Shekhovtsov 
wrote:

> Hi all!
> I encountered an error in my relative free energy calculations and do not
> know how to fix it.
> Molecules for which I want to carry out calculations contain a carboxyl
> group.
> To validate the protocol, I tried to run the fep-calculations of symmetric
> molecules for which the change in free energy will be zero.
> During validation, I was faced with the fact that the convergence error for
> charged ligands greatly exceeds that for neutral ones.
> For example, for the case of m-methyl benzoic acid [O-]C(= O)c1cc(C)ccc1:
> Mutation (-CH3 + H) for one meta position (-H + CH3) for another meta
> position relative to the carboxyl group.
>
> Ionized acid (solvent leg, 5 ns):
> point  0 -  1,   DG 28.74 +/-  0.02
> point  1 -  2,   DG 53.79 +/-  0.41
> point  2 -  3,   DG 27.76 +/-  2.45
> point  3 -  4,   DG 18.51 +/-  2.95
> point  4 -  5,   DG  8.80 +/-  0.79
> point  5 -  6,   DG  4.04 +/-  0.07
> point  6 -  7,   DG -0.06 +/-  0.04
> point  7 -  8,   DG -2.62 +/-  0.33
> point  8 -  9,   DG -8.95 +/-  0.27
> point  9 - 10,   DG -23.89 +/-  1.61
> point 10 - 11,   DG -29.32 +/-  1.29
> point 11 - 12,   DG -54.15 +/-  0.08
> point 12 - 13,   DG -28.79 +/-  0.03
>
> total  0 - 13,   DG -6.12 +/-  3.19
>
> Unionized acid (solvent leg, 5 ns):
> point  0 -  1,   DG  0.08 +/-  0.01
> point  1 -  2,   DG -6.88 +/-  0.10
> point  2 -  3,   DG -14.69 +/-  0.05
> point  3 -  4,   DG -21.07 +/-  0.03
> point  4 -  5,   DG -13.21 +/-  0.02
> point  5 -  6,   DG -7.46 +/-  0.02
> point  6 -  7,   DG -0.09 +/-  0.03
> point  7 -  8,   DG  7.37 +/-  0.02
> point  8 -  9,   DG 13.24 +/-  0.03
> point  9 - 10,   DG 21.15 +/-  0.06
> point 10 - 11,   DG 14.72 +/-  0.03
> point 11 - 12,   DG  6.98 +/-  0.07
> point 12 - 13,   DG -0.05 +/-  0.01
>
> total  0 - 13,   DG  0.09 +/-  0.22
>
> For a neutral molecule containing a charged carboxyl and amino groups
> ([O-]C(=O)c1cc(C)c([NH3+])cc1) the error is small:
> point  0 -  1,   DG -4.65 +/-  0.01
> point  1 -  2,   DG -19.95 +/-  0.04
> point  2 -  3,   DG -28.67 +/-  0.13
> point  3 -  4,   DG -57.84 +/-  0.19
> point  4 -  5,   DG -44.18 +/-  0.08
> point  5 -  6,   DG -26.84 +/-  0.05
> point  6 -  7,   DG  0.10 +/-  0.20
> point  7 -  8,   DG 26.92 +/-  0.06
> point  8 -  9,   DG 44.28 +/-  0.13
> point  9 - 10,   DG 57.79 +/-  0.13
> point 10 - 11,   DG 28.68 +/-  0.09
> point 11 - 12,   DG 19.98 +/-  0.12
> point 12 - 13,   DG  4.66 +/-  0.03
>
> total  0 - 13,   DG  0.28 +/-  0.21
>
> Adding Na and Cl ions to ([O-]C(=O)c1cc(C)c([NH3+])cc1) does not cause an
> increase in error.
> point  0 -  1,   DG -4.61 +/-  0.01
> point  1 -  2,   DG -19.96 +/-  0.07
> point  2 -  3,   DG -28.68 +/-  0.04
> point  3 -  4,   DG -57.67 +/-  0.04
> point  4 -  5,   DG -44.24 +/-  0.01
> point  5 -  6,   DG -26.91 +/-  0.04
> point  6 -  7,   DG  0.03 +/-  0.03
> point  7 -  8,   DG 26.89 +/-  0.03
> point  8 -  9,   DG 44.20 +/-  0.08
> point  9 - 10,   DG 57.75 +/-  0.09
> point 10 - 11,   DG 28.65 +/-  0.06
> point 11 - 12,   DG 20.02 +/-  0.05
> point 12 - 13,   DG  4.67 +/-  0.01
>
> total  0 - 13,   DG  0.12 +/-  0.18
>
> FEP-part of *.mdp:
> free-energy   = yes
> sc-power  = 1
> sc-alpha  = 0.5
> sc-coul   = yes
> fep-lambdas   = 0.0 0.01  0.05  0.1  0.2  0.3   0.4  0.6   0.7
> 0.8  0.9  0.95 0.99  1.0
> coul-lambdas  = 0.0 0.01  0.05  0.1  0.2  0.3   0.4  0.6   0.7
> 0.8  0.9  0.95 0.99  1.0
> 

Re: [gmx-users] Pressure control with constraints incompatible with expanded ensemble

[Michael Shirts]

Hi, Jacob-
Unfortunately, that is correct.  We had though we'd have functionality that
could support this combination in 2018 (a Monte Carlo barostat and/or
expanded ensemble in alchemical space working with leapfrog), but those did
not make it in.  The support for MTTK with constraints was removed in order
to simplify the inner loop of MD for other improvements.

The good thing is that in virtually all situations at P near 1 atm, \Delta
A and \Delta G for solvation calculations are pretty much identical, since
P \Delta V is so smal.  We can chat more if you think that you would need
additional help on this.  We are working with some improvements for
expanded ensemble, so contact me off the list if you want to discuss those
and see if they will help.

We DO know for sure that expanded ensemble or replica exchange with
Berendsen pressure control has some nasty issues (will be mentioned in a
paper hopefully being submitted soon, so it's better documented). So don't
do that.



On Wed, Sep 5, 2018 at 12:54 PM Jacob Monroe 
wrote:

> Hi all,
>
> I would like to use Gromacs (2018 versions on GPUs) to compute solvation
> free energies using expanded ensemble simulations.  I have no problems in
> the NVT ensemble, but when I try and implement pressure control I run into
> errors.
>
> Specifically, if I want to use expanded ensemble, I must use the md-vv
> integrator.  If I use this integrator, I must select MTTK pressure control
> (I’d rather not use Berendsen for production runs), but if I want to
> constrain hydrogen bonds, I cannot use MTTK (grompp throws an error).
> Previous versions of Gromacs appear to have allowed MTTK as long as SHAKE
> was used for constraints (see the mdp files here:
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Ethanol_solvation_with_expanded_ensemble
> ).
>
> So is the only way to run NPT, expanded ensemble simulations with current
> versions of Gromacs to use the Berendsen barostat?
>
> Thanks,
> Jacob
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] MBAR bootstrap

[Michael Shirts]

This is not the right place to ask about MBAR.  You could pose the question
on https://github.com/choderalab/pymbar if you have specific questions
about how to implement something there.
https://github.com/MobleyLab/alchemical-analysis is a wrapper for MBAR that
can process gromacs files.

See:
http://www.alchemistry.org/wiki/Analyzing_Simulation_Results#Bootstrap_Sampling
for a discussion of bootstrap sampling in the context of free energies.

On Tue, Jan 9, 2018 at 8:01 PM, Thanh Le  wrote:

> Hi everyone,
> I have a question regarding using bootstrap for MBAR.
> After running 30 windows for the complex and 20 windows for the ligand,
> 100 ns each, I have a problem trying to do bootstrap (based on what I have
> read, bootstrap sampling is recommended to analyze results from MBAR)
> Can anyone point me to the correct direction to implement bootstrap? Or
> any advice on how to analyze MBAR?
> Thanks,
> Thanh Le
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free energy calculation

[Michael Shirts]

You don't have a well-posed scientific question yet.  Come back when you've
thought about the hypotheses you want to test, and have written out what
you are planning to do.  Try to do the calculations yourself, and come back
if they are not working.  This is a help discussion forum for people who
are having trouble running Gromacs.  You haven't gotten to the point you
are running gromacs yet!

On Wed, Sep 20, 2017 at 4:13 PM, Sheikh Imamul Hossain <
s.imamul...@gmail.com> wrote:

> Hi All,
> I am interested to calculate the free energy of coarse-grained lipid
> monolayer which contains CG lipid, cholesterol, water and nanoparticles.
> Any suggestions regarding Metadynamics or  Umbrella sampling? Any tutorial
> that will be helpful for learning? I am trying Metadynamics learning with
> Belfast tutorial: plumed and umbrella with Bevan Lab.
> Thanks in advanced.
> Regards
> Imamul
>
> --
>
> Sincerely Your’s
> *Sheikh Imamul Hossain*
> PhD Research Student
>
> Mechanical Engineering | School of Chemistry, Physics and Mechanical
> Engineering | Faculty of Science and Engineering | Queensland University of
> Technology| 2 George St Brisbane | GPO Box 2434 Brisbane, Queensland,
>  Australia 4001
> phone: +61407258933 Office Email: s2.hoss...@.qut.edu.au
>
> Office:   P701B, P Block, Gardens Point
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Energy decomposition from FEP

[Michael Shirts]

More precisely, the van der Waals energy depends on whether you turn off
the charges at the same time or after.  It's a path-dependent quantity, not
a state function.

On Wed, Sep 20, 2017 at 12:33 PM, David van der Spoel 
wrote:

> On 20/09/17 19:51, Nikhil Maroli wrote:
>
>> Dear all,
>>
>> How to compute van der waals and electrostatic contribution to free energy
>> from fep simulation. I have used slow growth procedure using the
>> appropriate mdp files and followed fep tutorial.
>> Thanks in advance.
>>
>> This is physically meaningless. You want to decompose into enthalpy and
> entropy instead.
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Head of Department, Cell & Molecular Biology, Uppsala University.
> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
> http://www.icm.uu.se
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support
> /Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] FEP calculations on multiple nodes

[Michael Shirts]

Significantly more information is be needed to understand what happened.

* What does the logfile say that was output?
* What command are you using to run on multiple nodes?
* What is the .mdp file?
* How many nodes are you running on?
* What version of the program?

And so forth.

On Mon, Aug 21, 2017 at 4:49 AM, Vikas Dubey 
wrote:

> Hi everyone,
>
> I am trying runa  FEP calculation with a system of ~25 particles. I
> have 20 windows and I am currently running my simulations on 1 node each.
> Since, my system is big, I just get 2.5ns in day. So, I thought to run each
> of my window on multiple nodes but for some reason, it crashes immediately
> after starting with an error.
>
>
> *Segmentation fault (core dumped)*
>
> Simulations run smoothly on one node. No error there. I tried to see file
> but there was nothing written there. Any help would be very much
> appreciated.
>
>
> Thanks,
> Vikas
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] How to orientationally restrain a molecule?

[Michael Shirts]

I'm interested carrying out a long simulation of a thin, narrow protein.
Clearly, the volume-minimizing way to do this would be to put it in a thin,
narrow box of solvent.  The issue with this is that it could rotate around
and interact with itself. if the orientation of the protein changed
relative to the orientation of the box.

Is there a good way to restrain the orientation in GROMACS currently?  I
looked a little at the orientation restraint for NMR, but it didn't seem
like that would really work well (or at least, it was not obvious it would
work well), for restraining the overall orientation.  Am I missing
something?

Any other thoughts, suggestions for doing this?  Has anyone else done this?

Thanks!
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: Error bar in alchemical free energy calculation, alchemical_analysis

[Michael Shirts]

I'd take a look at the paper describing the tool:
http://dx.doi.org/10.1007/s10822-015-9840-9.


On Tue, Apr 18, 2017 at 12:07 PM, Alex  wrote:

> Any idea please !
> Thanks
> Alex
>
> -- Forwarded message --
> From: Alex 
> Date: Sun, Apr 16, 2017 at 11:29 AM
> Subject: Error bar in alchemical free energy calculation,
> alchemical_analysis
> To: gmx-us...@gromacs.org
>
>
> Dear gromacs user,
>
> I was wondering how error bar in relative binding free energy calculation
> is estimated? I mean the thermodynamic integration (TI) method using the
> alchemical_analysis python tools.
> Any reference also would be highly appreciated.
>
> Regards,
> Alex
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Potential energy in Zwanzig relationship

[Michael Shirts]

Take a look at how the python code does EXP.  The information you want is
in the dhdl.xvg. Look at the documentation for this file, and read the
header.  Come back with specific questions about the file if the header and
documentation are not enough.

On Fri, Feb 10, 2017 at 6:47 AM, gozde ergin <gozdeeer...@gmail.com> wrote:

> Dear Michael.
>
> Thanks for the reply. I have already used this python code however I would
> like to calculate myself by using the potential energies.
> Because I need to reweigh the free energy by using exponential
> re-weighting technique.
> > On 10 Feb 2017, at 14:44, Michael Shirts <mrshi...@gmail.com> wrote:
> >
> > https://github.com/MobleyLab/alchemical-analysis
> >
> > Takes gromacs dhdl.xvg output and calculate free energies by many
> different
> > methods, including BAR, MBAR and Zwanzig.  See
> > http://www.alchemistry.org/wiki/Main_Page for more information.
> >
> > On Fri, Feb 10, 2017 at 5:26 AM, gozde ergin <gozdeeer...@gmail.com>
> wrote:
> >
> >> Dear all,
> >>
> >> I run thermodynamic integration simulation with gromacs and got the free
> >> energy by g_bar command.
> >> I also would like to estimate this free energy by using Zwanzig
> >> relationship of \DeltaG = -RT ln (<exp{-\Delta U_(ij)/RT}>_i
> >> Here U is the potential energy, right?
> >> However the results that I am getting with g_bar (82 kcal/mol) and with
> >> Zwanzig equation (406 kcal/mol) are totally different.
> >> Do anyone has an idea why I am getting that different results?
> >>
> >> Bests
> >> --
> >> Gromacs Users mailing list
> >>
> >> * Please search the archive at http://www.gromacs.org/
> >> Support/Mailing_Lists/GMX-Users_List before posting!
> >>
> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>
> >> * For (un)subscribe requests visit
> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >> send a mail to gmx-users-requ...@gromacs.org.
> >>
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Potential energy in Zwanzig relationship

[Michael Shirts]

https://github.com/MobleyLab/alchemical-analysis

Takes gromacs dhdl.xvg output and calculate free energies by many different
methods, including BAR, MBAR and Zwanzig.  See
http://www.alchemistry.org/wiki/Main_Page for more information.

On Fri, Feb 10, 2017 at 5:26 AM, gozde ergin  wrote:

> Dear all,
>
> I run thermodynamic integration simulation with gromacs and got the free
> energy by g_bar command.
> I also would like to estimate this free energy by using Zwanzig
> relationship of \DeltaG = -RT ln (_i
> Here U is the potential energy, right?
> However the results that I am getting with g_bar (82 kcal/mol) and with
> Zwanzig equation (406 kcal/mol) are totally different.
> Do anyone has an idea why I am getting that different results?
>
> Bests
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Diffusion coefficient water

[Michael Shirts]

See the following paper
(http://pubs.acs.org/doi/abs/10.1021/ct400109a) for some examples of
diffusion calculations in GROMACS.

The finite size correction for diffusion coefficients is something
different than the long range energy corrections. A discussion of this
(I can't vouch for all the details, I just googled) is here:
http://pubs.acs.org/doi/abs/10.1021/jp0756247

On Thu, Aug 18, 2016 at 6:37 AM, Mark Abraham  wrote:
> Hi,
>
> I would also consider using a pressure-coupling algorithm that isn't known
> to create the wrong velocity distribution.
>
> Mark
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] quick clarification on free energy and domain decomposition

[Michael Shirts]

Hi, all-

I wanted to check quickly.  I noticed a difference in domain
decomposition when I ran with free energy on and off (using the
coupl-mol functionality), with a relatively large molecule (8 residue
peptide) being decoupled.  Domain decomposition failed for the free
energy on 4 cores, and looking at the log file I noticed:

No free energy:
< two-body bonded interactions: 0.434 nm, LJ-14, atoms 340 348

---
Free energy:
> two-body bonded interactions: 2.586 nm, LJC Pairs NB, atoms 1727 1809

I'm sure people have noticed this before, but I'm wondering if there
is some workaround?  If I mutate the atoms directly rather than using
coupl-mol, would I be able to get around it?

Thanks!
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] PSCP

[Michael Shirts]

See http://pubs.acs.org/doi/abs/10.1021/acs.jctc.6b00397 for an
example of a PSCP in GROMACS (although it was from crystal to
crystal).  It takes series of different .mdps to do it.

On Fri, Jul 22, 2016 at 5:34 AM, Justin Lemkul  wrote:
>
>
> On 7/21/16 11:20 PM, kar...@shell.com wrote:
>>
>> Hello Justin!
>>
>> In recent literature, a pseudo-supercritical path has been used to model
>> the melting point of crystals. To do so, the molecule starts in the liquid
>> state, (1) has the intermolecular forces scaled down (which I did using
>> lambda and free energy techniques), (2) the volume of the box is reduced to
>> that of the solid state, (3) an external potential field is introduced, and
>> then (4) the field is removed and the IMF are restored.
>>
>> Does GROMACS have the capacity to introduce an FCC lattice of 3D Gaussian
>> Wells (step 3)? I know this is also supposed to use lambda-integration as
>> well, but I do not know what the syntax is for the mdp file.
>>
>
> If an option isn't here:
> http://manual.gromacs.org/documentation/5.1.3/user-guide/mdp-options.html or
> in the printed manual, it doesn't exist.
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 629
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] md-vv and md

[Michael Shirts]

Could you file an issue on redmine with sample input files you use? It
may be some particular combination of md-vv and some other options
that is causing the problem, and it is very hard to debug anything
without a precise copy of the inputs that cause the problem.

On Thu, Jul 21, 2016 at 9:28 AM, gozde ergin <gozdeeer...@gmail.com> wrote:
> Dear all,
>
> I am using md-vv integrator in 5.1 and still get this all bulk moving problem.
> Before I was using 4.6 and informed that this dribble is fixed 5.1 however I 
> still got the same problem.
> My bulk is moving in z direction however if I use md integrator I do not face 
> with this problem.
> Any idea?
>> On 18 Jul 2015, at 04:13, Michael Shirts <mrshi...@gmail.com> wrote:
>>
>>> Has problem been solved about molecules drifting in md-vv integrator?
>>
>> We think so, but some additional tests are running to address any
>> other lingering integrator issues.  The fix is in 5.1 beta, but not in
>> the 4.6 branch yet (it's a pull request, but not yet merged).
>>
>>> Also can you please check my potential energy second pic in
>>> http://imgur.com/6aJkRoQ#0
>>> . <http://imgur.com/6aJkRoQ#0>
>>> What could be the reason this big jump in potential energy (integrator :
>>> md-vv)?
>>> And why I do not see the same jump by using md integrator?
>>
>> I don't think there's sufficient information here to answer that question.
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
>> mail to gmx-users-requ...@gromacs.org.
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: Conversion from GROMACS .top .itp to CHARMM .psf .par combo

[Michael Shirts]

Try ParmEd.

http://parmed.github.io/ParmEd/html/index.html

On Sat, May 28, 2016 at 7:57 AM, Sourav Ray  wrote:
> Hello
>
> Currently I have a combination of:
>
> 1. .pdb file (four copies of a monomer)
> 2. .itp file defining the topology of one monomer (applied on all four with
> the help of the .top file described below).
> 3. .top file listing the constituents of the entire pdb:
>
> 
>
> #include "martini_v2.2P.itp"
>
> #include "Protein_A.itp"
>
> [ system ]
> ; name
> Martini system from fa_qpr.pdb
>
> [ molecules ]
> ; namenumber
> Protein_A  1
> Protein_A  1
> Protein_A  1
> Protein_A  1
>
> -
>
> It would be really helpful if someone could refer to a software or script
> maybe that takes the .top and .itp files as input and converts them to .psf
> and .par/.prm files.
>
> Thanks and regards
> Sourav
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Expanded ensemble error

[Michael Shirts]

If you are doing solvation in water, I would suggest not bothering
with expanded ensemble. I would just use standard molecular dynamics
at a range of lambda values.  We have some information that it can be
useful with binding problems, but with solvation alone, you are
probably just fine running a series of alchemical simulations in
parallel.  The expanded ensemble simulations don't really help
sampling that much.

We have recently done some expanded ensemble simulations with drug
binding, and we will try to put a tutorial in the next couple of
months.  There are some subtleties that can cause problems.

On Sat, Apr 30, 2016 at 7:37 AM, Marzieh Saeedi Masineh
 wrote:
> Dear Gromacs users,
> I wanted to use expanded ensemble for calculating solvation free energy of a 
> drug molecule in the water based on "Ethanol solvation with expanded ensemble 
> in gromacs". in this example the integrator, t coupl, pcoupl, and 
> constraint-algorithm were set to md-vv, Nose-Hoover, MTTK and shake, 
> respectively. When I used these parameters in Gormacs 4.6, I faced an error 
> related to the shake algorithm:(SHAKE is not supported with domain 
> decomposition and constraint that cross charge group boundaries, use LINCS)
> I tried different options for these parameters, but I my runs were crashed. 
> Can any person help about the exact combination of these parameters ( the 
> integrator, t coupl, pcoupl, and constraint-algorithm) that works properly. 
> it is noted that, I searched in the mailing list, but i didn't find any 
> answer about it. Thanks in advance.
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Looking for the source of the H-REMD slowdown when decoupling a lot of atoms

[Michael Shirts]

'd expect ~2x in PME and less in
> constraints), but it could be that this too is simply because most
> interaction in the system are perturbed which trigger non-optimized
> code-paths.
>
> Cheers,
> --
> Szilárd
>
> On Wed, Mar 30, 2016 at 3:05 PM, Mark Abraham <mark.j.abra...@gmail.com>
> wrote:
>
>> Hi,
>>
>> Yeah, unfortunately Michael's pretty much right - the free-energy kernel is
>> currently that cousin nobody talks about (or to). It's essentially
>> unchanged since GROMACS 4.0 days, except that the Verlet scheme has some
>> kludge so that it can call the same kernel that the group scheme used to
>> call. But it has none of the optimizations and also some simplifying
>> pessimizations. The result is highly unsuitable for Chris's use case, but
>> not horrible for the more normal case of perturbing a small part of a
>> system. For now, I can only suggest trying a run with more than one OpenMP
>> thread per rank, but there's nothing in the log file snippets that fills me
>> with any hope that it would be noticeably faster.
>>
>> We have a pile of infrastructure built and in the works that will lead to
>> being able to offer similarly optimized free-energy kernels, but they won't
>> see the light of day until next year, I'm afraid. A set of sample .tprs at
>> Redmine 742 would be most welcome, however - it's very good for us to know
>> when/that we're optimizing a workflow someone actually wants to run, but
>> currently has reason to avoid.
>>
>> Mark
>>
>> On Wed, Mar 30, 2016 at 8:05 AM Michael Shirts <mrshi...@gmail.com> wrote:
>>
>> > Hi, Chris: I'm pretty sure that it's because the nonbonded free
>> > energies are much slower than the standard free energies.  You state:
>> >
>> > > I took a look at gmxlib/nonbonded/nb_free_energy.c in v.5.1.2, but I
>> was
>> > unable to find a function called "gmx_waste_time_here()" and beyond that
>> I
>> > was out of my depth.
>> >
>> > But it's much more the fact that the non-free energy nonbondeds are
>> > SUPER optimized.
>> >
>> > I don't see a particularly viable way around this for now. The only
>> > thing I can think of splitting the neighborlists into two force calls,
>> > and scaling the forces and energies coming out of those.  That would
>> > be a huge pain.
>> >
>> > On Tue, Mar 29, 2016 at 9:50 PM, Christopher Neale
>> > <chris.ne...@alum.utoronto.ca> wrote:
>> > > Dear Users:
>> > >
>> > > I am trying to do some Hamiltonian replica exchange (H-REMD) in gromacs
>> > 5.1.2 and am running up against really large slowdowns when decoupling a
>> > large number of atoms. I am decoupling 5360 atoms out of the 15520 atoms
>> in
>> > my system. The goal is not to get a PMF, but to enhance sampling using
>> the
>> > REST approach to partially decouple lipids in a bilayer. This approach
>> > enhances lipid relaxation times (
>> > http://pubs.acs.org/doi/pdf/10.1021/ct500305u ) though the authors of
>> > that paper modified the gromacs code to do their own H-REMD in order to
>> > avoid the really slow speed they also got when decoupling lots of atoms
>> via
>> > the free energy code.
>> > >
>> > > I have already posted details here
>> http://redmine.gromacs.org/issues/742
>> > , which includes .mdp options and some timing output. I compare the
>> timing
>> > output to a standard temperature REMD (T-REMD) run. For my usage, the
>> > slowdown is about 12x for H-REMD vs. T-REMD.
>> > >
>> > > I am motivated to find a solution within gromacs because the
>> alternative
>> > is to use gromacs 4.6.7 with plumed (or with the aforementioned modified
>> > code, which is also gromacs v4). Normally that would be a viable option,
>> > but I am using the charmm force field and the charmm TIP3P water and I
>> > would rather not give up the speed boost that I see in gromacs v5.1.2,
>> > which allows the use of the verlet cutoff scheme and has been tested and
>> > shown to give the correct reproduction of charmm forces (vs. the forces
>> one
>> > would get using the charmm simulation software).
>> > >
>> > > I took a look at gmxlib/nonbonded/nb_free_energy.c in v.5.1.2, but I
>> was
>> > unable to find a function called "gmx_waste_time_here()" and beyond that
>> I
>> > was out of my depth.
>> > >
>> > > Thank you for any pointers,

Re: [gmx-users] Looking for the source of the H-REMD slowdown when decoupling a lot of atoms

[Michael Shirts]

Hi, Chris: I'm pretty sure that it's because the nonbonded free
energies are much slower than the standard free energies.  You state:

> I took a look at gmxlib/nonbonded/nb_free_energy.c in v.5.1.2, but I was 
> unable to find a function called "gmx_waste_time_here()" and beyond that I 
> was out of my depth.

But it's much more the fact that the non-free energy nonbondeds are
SUPER optimized.

I don't see a particularly viable way around this for now. The only
thing I can think of splitting the neighborlists into two force calls,
and scaling the forces and energies coming out of those.  That would
be a huge pain.

On Tue, Mar 29, 2016 at 9:50 PM, Christopher Neale
 wrote:
> Dear Users:
>
> I am trying to do some Hamiltonian replica exchange (H-REMD) in gromacs 5.1.2 
> and am running up against really large slowdowns when decoupling a large 
> number of atoms. I am decoupling 5360 atoms out of the 15520 atoms in my 
> system. The goal is not to get a PMF, but to enhance sampling using the REST 
> approach to partially decouple lipids in a bilayer. This approach enhances 
> lipid relaxation times ( http://pubs.acs.org/doi/pdf/10.1021/ct500305u ) 
> though the authors of that paper modified the gromacs code to do their own 
> H-REMD in order to avoid the really slow speed they also got when decoupling 
> lots of atoms via the free energy code.
>
> I have already posted details here http://redmine.gromacs.org/issues/742 , 
> which includes .mdp options and some timing output. I compare the timing 
> output to a standard temperature REMD (T-REMD) run. For my usage, the 
> slowdown is about 12x for H-REMD vs. T-REMD.
>
> I am motivated to find a solution within gromacs because the alternative is 
> to use gromacs 4.6.7 with plumed (or with the aforementioned modified code, 
> which is also gromacs v4). Normally that would be a viable option, but I am 
> using the charmm force field and the charmm TIP3P water and I would rather 
> not give up the speed boost that I see in gromacs v5.1.2, which allows the 
> use of the verlet cutoff scheme and has been tested and shown to give the 
> correct reproduction of charmm forces (vs. the forces one would get using the 
> charmm simulation software).
>
> I took a look at gmxlib/nonbonded/nb_free_energy.c in v.5.1.2, but I was 
> unable to find a function called "gmx_waste_time_here()" and beyond that I 
> was out of my depth.
>
> Thank you for any pointers,
> Chris.
>
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Separate lambdas for repulsive and dispersive free energies?

[Michael Shirts]

Correct, currently a WCA decomposition has not been implemented.

I'm not sure if I understand this comment, though:

>  For this reason, I only plan to calculate E_vdw_dispersive and E_elec for my 
> system.

E_vdw_repulsive is a gigantic contribution to the solvation free
energy of the protein. Why would you not calculate it?  Without it,
you're calculating something else entirely.

On Sun, Mar 27, 2016 at 8:44 PM, Nuo Wang  wrote:
> Hi all,
>
> I am trying to calculate the hydration/solvation free energy for proteins
> using free energy perturbation:
> E_hydration = E_vdw + E_elec
> E_vdw (WCA decomposition) = E_vdw_repulsive + E_vdw_dispersive
>
> However, unlike small molecules, if I turn off both elec and vdw for a
> normal-sized protein, a huge cavity will be created and it will take
> unfeasible amount of time to equilibrate the system. For this reason, I
> only plan to calculate E_vdw_dispersive and E_elec for my system.
>
> I have read the mailing list posts of implementing the repulsive WCA
> potential in Gromacs, but I have not found any discussion that try to
> implement the WCA-decomposed potential, i.e. repulsive + dispersive, and
> further scaling them separately in free energy calculations.
>
> According the the manual, it seems that it is impossible to implement the
> decomposition with the Gromacs function form g(x) + f(x). On top of this,
> vdw is scaled together by one lambda. I wonder if I am right on this and if
> there is an alternative way that I am not aware of?
>
> FYI, the WCA decomposed scaling protocol has been implemented in CHARMM:
> http://pubs.acs.org/doi/abs/10.1021/jp048502c
>
> Thanks,
> Nuo
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] μVT simulations and genbox

[Michael Shirts]

In theory, https://github.com/GromPy/GromPy can do it
(http://www.ncbi.nlm.nih.gov/pubmed/22370965), but it's old and I've
never tried it.  Not in the main code currently.

On Wed, Feb 10, 2016 at 11:38 AM, Irem Altan  wrote:
> Hi,
>
> I have two questions:
>
> - Is it possible to do μVT simulations?
> - When I copy vdwradii.dat to my working folder and modify it, does it affect 
> the MD simulations? I used it to modify the number of water molecules genbox 
> adds, and then ran an MD simulation without removing the vdwradii.dat file. 
> Would that affect anything?
>
> Best,
> Irem
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: FW: PME settings in free energy calculations

[Michael Shirts]

Hi, Dries-

Can you print out what the corresponding section of mdout.mdp look like?

Have you  tried using the cutoffs above with group cutoffs, and seen
what the difference is?

On Thu, Feb 4, 2016 at 4:59 AM, Dries Van Rompaey
<dries.vanromp...@gmail.com> wrote:
> Hi,
>
> I'm using PME with the default Potential-Shift-Verlet modifier. If I'm
> interpreting the manual and comments on the mailinglist correctly, the
> settings for rcoulomb-switch shouldn't affect this and PME-switching
> shouldn't be necessary (although there seems to have been some discussion
> about this, https://gerrit.gromacs.org/#/c/2220/). Has a consensus been
> reached elsewhere?
>
> Current nonbonded settings in my .mdp files are:
> ; Electrostatics
> coulombtype  = PME
> rcoulomb= 1.0
>
> ; van der Waals
> vdwtype= cutoff
> vdw-modifier= Potential-switch
> rvdw-switch  = 0.9
> rvdw = 1.0
>
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr  = EnerPres
>
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing   = 0.08
>
> ; EWALD/PME/PPPM parameters
> pme_order= 6
> ewald_rtol = 1e-06
> epsilon_surface= 0
>
> Thanks in advance,
>
> Dries
>
> On 4 February 2016 at 00:56, Michael Shirts <mrshi...@gmail.com> wrote:
>
>> Hi, Dries-
>>
>> Questions like this are probably best answered on the gmx-users list.
>> I can't say too much for the Verlet scheme -- I know that it was
>> relatively recently adapted for free energies, and there may be some
>> combinations of settings that could give unanticipated results.
>>
>> Pretty much all of our experience about nonbonded calculations and
>> free energies is collected in the following paper:
>> http://pubs.acs.org/doi/abs/10.1021/ct4005068. We used the group
>> cutoff scheme since that was the only one that was supported in
>> GROMACS at the time.
>>
>> Since you haven't sent any files, it's hard to tell what is actually
>> going on. The one thing that has a tendency to happen with the more
>> recent update schemes is if you set a potential modifier that is a
>> switch, but don't set the distance the switch starts at, then it is
>> automatically set to zero.  Check the mdout.mdp to see if this is
>> happening.
>>
>> > From: Van Rompaey Dries <dries.vanromp...@uantwerpen.be>
>> Date: Wednesday, February 3, 2016 at 12:34 PM
>> To: Michael Shirts <michael.shi...@colorado.edu>
>> Subject: PME settings in free energy calculations
>>
>> Dear prof. Shirts,
>>
>> I'm currently trying to figure out the PME settings for a relative
>> free binding energy simulation I'm working on. I took the parameters
>> from Matteo Aldeghi's
>> paper(
>> http://pubs.rsc.org/en/content/articlelanding/2015/sc/c5sc02678d#!divAbstract
>> )
>> as a starting point (adapted for the Verlet scheme by scaling the
>> fourierspacing to 0.08 and setting coulomb = rvdw = 1.0). I then tried
>> to verify these settings by comparing a single point energy
>> calculation with these settings and one with very long coulomb cutoffs
>> as recommended on alchemistry.org.
>>
>> Unfortunately, I can't seem to get this quite right. I'm getting
>> differences in the hundreds of kj/mol, leading me to suspect I'm doing
>> something wrong. I'm calculating the energy values by extracting
>> CoulombSR and Coul-recip from the energy.xvg files. I've tried
>> calculating the energy with coulomb cutoffs at 3 nm and 10 nm, but
>> agreement with the PME results remains rather poor. David Mobley
>> mentioned you performed extensive research on this topic, and I'm
>> hoping you could point me in the right direction.
>>
>> Thanks in advance,
>>
>> Dries
>>
>>
>>
>> 
>> Michael Shirts
>> Associate Professor
>> michael.shi...@colorado.edu
>> Phone: (303) 735-7860
>> Office: JSCBB D317
>> Department of Chemical and Biological Engineering
>> University of Colorado Boulder
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
> --
> Gromacs Users m

Re: [gmx-users] Fwd: FW: PME settings in free energy calculations

[Michael Shirts]

Also, can you be more precise about the inconsistency in the outputs
you are seeing (with numerical output from GROMACS and analysis code).

On Thu, Feb 4, 2016 at 6:17 AM, Michael Shirts <mrshi...@gmail.com> wrote:
> Hi, Dries-
>
> Can you print out what the corresponding section of mdout.mdp look like?
>
> Have you  tried using the cutoffs above with group cutoffs, and seen
> what the difference is?
>
> On Thu, Feb 4, 2016 at 4:59 AM, Dries Van Rompaey
> <dries.vanromp...@gmail.com> wrote:
>> Hi,
>>
>> I'm using PME with the default Potential-Shift-Verlet modifier. If I'm
>> interpreting the manual and comments on the mailinglist correctly, the
>> settings for rcoulomb-switch shouldn't affect this and PME-switching
>> shouldn't be necessary (although there seems to have been some discussion
>> about this, https://gerrit.gromacs.org/#/c/2220/). Has a consensus been
>> reached elsewhere?
>>
>> Current nonbonded settings in my .mdp files are:
>> ; Electrostatics
>> coulombtype  = PME
>> rcoulomb= 1.0
>>
>> ; van der Waals
>> vdwtype= cutoff
>> vdw-modifier= Potential-switch
>> rvdw-switch  = 0.9
>> rvdw = 1.0
>>
>> ; Apply long range dispersion corrections for Energy and Pressure
>> DispCorr  = EnerPres
>>
>> ; Spacing for the PME/PPPM FFT grid
>> fourierspacing   = 0.08
>>
>> ; EWALD/PME/PPPM parameters
>> pme_order= 6
>> ewald_rtol = 1e-06
>> epsilon_surface= 0
>>
>> Thanks in advance,
>>
>> Dries
>>
>> On 4 February 2016 at 00:56, Michael Shirts <mrshi...@gmail.com> wrote:
>>
>>> Hi, Dries-
>>>
>>> Questions like this are probably best answered on the gmx-users list.
>>> I can't say too much for the Verlet scheme -- I know that it was
>>> relatively recently adapted for free energies, and there may be some
>>> combinations of settings that could give unanticipated results.
>>>
>>> Pretty much all of our experience about nonbonded calculations and
>>> free energies is collected in the following paper:
>>> http://pubs.acs.org/doi/abs/10.1021/ct4005068. We used the group
>>> cutoff scheme since that was the only one that was supported in
>>> GROMACS at the time.
>>>
>>> Since you haven't sent any files, it's hard to tell what is actually
>>> going on. The one thing that has a tendency to happen with the more
>>> recent update schemes is if you set a potential modifier that is a
>>> switch, but don't set the distance the switch starts at, then it is
>>> automatically set to zero.  Check the mdout.mdp to see if this is
>>> happening.
>>>
>>> > From: Van Rompaey Dries <dries.vanromp...@uantwerpen.be>
>>> Date: Wednesday, February 3, 2016 at 12:34 PM
>>> To: Michael Shirts <michael.shi...@colorado.edu>
>>> Subject: PME settings in free energy calculations
>>>
>>> Dear prof. Shirts,
>>>
>>> I'm currently trying to figure out the PME settings for a relative
>>> free binding energy simulation I'm working on. I took the parameters
>>> from Matteo Aldeghi's
>>> paper(
>>> http://pubs.rsc.org/en/content/articlelanding/2015/sc/c5sc02678d#!divAbstract
>>> )
>>> as a starting point (adapted for the Verlet scheme by scaling the
>>> fourierspacing to 0.08 and setting coulomb = rvdw = 1.0). I then tried
>>> to verify these settings by comparing a single point energy
>>> calculation with these settings and one with very long coulomb cutoffs
>>> as recommended on alchemistry.org.
>>>
>>> Unfortunately, I can't seem to get this quite right. I'm getting
>>> differences in the hundreds of kj/mol, leading me to suspect I'm doing
>>> something wrong. I'm calculating the energy values by extracting
>>> CoulombSR and Coul-recip from the energy.xvg files. I've tried
>>> calculating the energy with coulomb cutoffs at 3 nm and 10 nm, but
>>> agreement with the PME results remains rather poor. David Mobley
>>> mentioned you performed extensive research on this topic, and I'm
>>> hoping you could point me in the right direction.
>>>
>>> Thanks in advance,
>>>
>>> Dries
>>>
>>>
>>>
>>> 
>>> Michael Shirts
>>> Associate Professor
>>> michael.shi...@colorado.edu
>>&g

[gmx-users] Fwd: FW: PME settings in free energy calculations

[Michael Shirts]

Hi, Dries-

Questions like this are probably best answered on the gmx-users list.
I can't say too much for the Verlet scheme -- I know that it was
relatively recently adapted for free energies, and there may be some
combinations of settings that could give unanticipated results.

Pretty much all of our experience about nonbonded calculations and
free energies is collected in the following paper:
http://pubs.acs.org/doi/abs/10.1021/ct4005068. We used the group
cutoff scheme since that was the only one that was supported in
GROMACS at the time.

Since you haven't sent any files, it's hard to tell what is actually
going on. The one thing that has a tendency to happen with the more
recent update schemes is if you set a potential modifier that is a
switch, but don't set the distance the switch starts at, then it is
automatically set to zero.  Check the mdout.mdp to see if this is
happening.

> From: Van Rompaey Dries <dries.vanromp...@uantwerpen.be>
Date: Wednesday, February 3, 2016 at 12:34 PM
To: Michael Shirts <michael.shi...@colorado.edu>
Subject: PME settings in free energy calculations

Dear prof. Shirts,

I'm currently trying to figure out the PME settings for a relative
free binding energy simulation I'm working on. I took the parameters
from Matteo Aldeghi's
paper(http://pubs.rsc.org/en/content/articlelanding/2015/sc/c5sc02678d#!divAbstract)
as a starting point (adapted for the Verlet scheme by scaling the
fourierspacing to 0.08 and setting coulomb = rvdw = 1.0). I then tried
to verify these settings by comparing a single point energy
calculation with these settings and one with very long coulomb cutoffs
as recommended on alchemistry.org.

Unfortunately, I can't seem to get this quite right. I'm getting
differences in the hundreds of kj/mol, leading me to suspect I'm doing
something wrong. I'm calculating the energy values by extracting
CoulombSR and Coul-recip from the energy.xvg files. I've tried
calculating the energy with coulomb cutoffs at 3 nm and 10 nm, but
agreement with the PME results remains rather poor. David Mobley
mentioned you performed extensive research on this topic, and I'm
hoping you could point me in the right direction.

Thanks in advance,

Dries



~~~~~~~~
Michael Shirts
Associate Professor
michael.shi...@colorado.edu
Phone: (303) 735-7860
Office: JSCBB D317
Department of Chemical and Biological Engineering
University of Colorado Boulder
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Contact_Angles_Calculation

[Michael Shirts]

See:

http://pubs.acs.org/doi/abs/10.1021/acs.langmuir.5b00842

For a study using GROMACS.  They code can be downloaded (though it
might not be trivial to use).

On Thu, Jan 21, 2016 at 6:32 PM, HZJ <1047334...@qq.com> wrote:
> Dear GMX users,
>
>
> Recently, I want to analyze the contact angle of a gas bubble adsorbed on a 
> surface, during the simulation,
> the bubbule can move freely on the surface. I aim to calculate the contact 
> angle for each frame. Did anyone
> do such analysis and suggest a simple way to do it? I read some paper, one 
> method is to calculate the 2D density of the bubble, and then to fit the 
> edges of the curve to a function type and then get the angle. It seems this 
> method is simple only for one frame. Thanks!
>
>
>
> Best wishes!
>
>
> Dr. He
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Contact_Angles_Calculation

[Michael Shirts]

I believe the last line of the paper is: "We have made the code used
for this analysis available on the GitHub site:
https://github.com/shirtsgroup/droplet-analysis.;

On Thu, Jan 21, 2016 at 10:02 PM, HZJ <1047334...@qq.com> wrote:
> Thanks. I have come across this paper. It introduce how to calculate this 
> angle. No code is provided with this
> paper. Can they share this code? I use a smooth surface.
>
>
> Best wishes!
>
>
>
>
>
> ------ Original --
> From:  "Michael Shirts";<mrshi...@gmail.com>;
> Date:  Fri, Jan 22, 2016 11:40 AM
> To:  "Discussion list for GROMACS users"<gmx-us...@gromacs.org>;
> Cc:  "gromacs.org_gmx-users"<gromacs.org_gmx-users@maillist.sys.kth.se>;
> Subject:  Re: [gmx-users] Contact_Angles_Calculation
>
>
>
> See:
>
> http://pubs.acs.org/doi/abs/10.1021/acs.langmuir.5b00842
>
> For a study using GROMACS.  They code can be downloaded (though it
> might not be trivial to use).
>
> On Thu, Jan 21, 2016 at 6:32 PM, HZJ <1047334...@qq.com> wrote:
>> Dear GMX users,
>>
>>
>> Recently, I want to analyze the contact angle of a gas bubble adsorbed on a 
>> surface, during the simulation,
>> the bubbule can move freely on the surface. I aim to calculate the contact 
>> angle for each frame. Did anyone
>> do such analysis and suggest a simple way to do it? I read some paper, one 
>> method is to calculate the 2D density of the bubble, and then to fit the 
>> edges of the curve to a function type and then get the angle. It seems this 
>> method is simple only for one frame. Thanks!
>>
>>
>>
>> Best wishes!
>>
>>
>> Dr. He
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
>> mail to gmx-users-requ...@gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] amber force field ff12SB and ff14SB for gromacs uploaded

[Michael Shirts]

> Yes, I just try to show to Alexey Shvetsov or some gromas users who want
> to do check the conversion.

The issue is that the hardest and most time consuming part is not the
conversion itself, but the testing.  Without the testing, the
conversion isn't that useful.  I appreciate the work you've done, but
the hard part is still ahead.

> It is too much thing to post for the detail of
> residue comparing. I agree with your suggestion about testing. I will try
> to do test and report it as soon as possible.

The information here gives a good example of good testing.

http://ffamber.cnsm.csulb.edu/ffamber.php

In fact, you could offer your files to Prof. Sorin and see if he wants
to collaborate on testing.

> Here I really want other
> gromacs users do check my conversion and take apart the testing to make
> the conversion completely successful. It will give us more freedom in
> doing with GROMACS.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Errors from alchemical_analysis.py from analyzing .xvg files

[Michael Shirts]

You should file an issue on the alchemical-analysis issue tracker.

https://github.com/MobleyLab/alchemical-analysis

On Tue, Dec 15, 2015 at 8:35 PM, Nathan K Houtz  wrote:
> Hello,
>
>
> This request is related to python programs used for analyzing gromacs 
> outputs, not gromacs itself. If there is a more appropriate forum for this 
> question, let me know and I apologize for posting here.
>
>
> I have performed a bunch of thermodynamic-integration simulations for liquid 
> water in Gromacs, and am trying to analyze the .xvg files with the python 
> script mentioned in the ethanol solvation tutorial: 
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy
>  (I think the tutorial is out of date - the python script is now called 
> alchemical_analysis.py and I think it is called with slightly different 
> flags). It may be relevant that I'm using Windows 10 and installed Anaconda2 
> (64 bit) with the gui installer and then pymbar through Anaconda2 before 
> downloading alchemical-analysis-master.
>
>
> In a Windows command prompt, I entered:
>
>>python alchemical_analysis.py -d C:\\ -t 200 -p 1 -v  > results
>
>
> This is the result:
>
> Traceback (most recent call last):
>   File "alchemical_analysis.py", line 1224, in 
> main()
>   File "alchemical_analysis.py", line 1166, in main
> nsnapshots, lv, dhdlt, u_klt = parser_gromacs.readDataGromacs(P)
>   File 
> "C:\<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", 
> line 187, in readDataGromacs
> fs = sorted(fs, key=F.sortedHelper)
>   File 
> "C:\<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", 
> line 48, in sortedHelper
> self.state = l[0] = int(l[0]) # Will be of use for selective MBAR 
> analysis.
> IndexError: list index out of range
>
>
> where <...> is just the directory to which I downloaded 
> alchemical-analysis-master. I don't understand the error and cannot find help 
> for it online. In the directory where I have the .xvg files, all the files 
> are simply named integers: 0.xvg, 1.xvg, 2.xvg,  20.xvg, and there are no 
> other files in that folder. I don't want to tinker with the python code and 
> accidentally mess something up. Am I calling the script correctly? Could it 
> be the operating system (Windows vs. Linux)? I might be able to get on a 
> Linux machine if it would help.
>
>
> Thanks very much,
>
> N. H.
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free energy of ice crystals (Thermodynamic Integration)

[Michael Shirts]

IF it's NVT, then you should just be able to add harmonic potentials
using the position restraints (not constraints) -- the harmonic
restraints can be controlled by lambda.  You should only need to alter
the charge and VDW interactions -- no need to adjust the bonds.

We have been doing some work on crystals of small organics-- email me
off the list and we might be able to give you some example scripts.

On Fri, Dec 11, 2015 at 8:46 PM, Nathan K Houtz  wrote:
> Hello,
>
>
> I would like to perform some free energy calculations of ice Ic (cubic ice) 
> at various temperatures. However, I'm having difficulty finding information 
> on how Gromacs can be used to do this. The goal is to gradually turn off VDW 
> and Coulomb interactions, just as one would do in a liquid, but also turn on 
> virtual spring-like potentials to create an Einstein crystal, which is where 
> I'm having trouble. I know you can fix atom locations, and considering bond 
> potentials are spring-like, my best idea is to create a virtual copy of each 
> molecule (not virtual sites, but just atoms that have no charge and don't 
> interact with other atoms/molecules), fix it in the crystal lattice, and 
> create bonds between corresponding atoms with an equilibrium distance of 0. 
> But this means I have to simulate twice as many atoms, and more importantly I 
> don't know how to couple only some of a molecule's bonds with lambda. Am I 
> even on the right track, or is there a more straightforward way to do this? 
 Ar
>  e there any tutorials for free energies of solids that I could take a look 
> at? (I couldn't find any)
>
>
> Thanks for your help!
>
> N. H.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] AMBER force fields, ff12SB and ff14SB for GROMACS

[Michael Shirts]

Great!  It's important to have all data validated carefully. Have you
put together a set of tests, with coverage of all the atom types and
interaction types, showing that the molecules give the same energies
in both GROMACS and AMBER?  Validating conversions is quite important
so that people can trust the results.

On Thu, Dec 10, 2015 at 4:01 PM, Man Hoang Viet  wrote:
> Dear GROMACS Admins,
>
> I have implemented AMBER force field ff12SB and ff14SB into GROMACS and
> want to share them with all gromacs users. I have tried to upload it via
> below link
> http://www.gromacs.org/Downloads/User_contributions/Force_fields
> However, I can not do register to upload. Could you please tell me how I
> can share the amber force fields to gromacs users.
>
> Thank you very much.
>
> Yours sincerely,
> --
> Viet Man,
>
> Postdoctoral Research Associate,
> Physics Department, NCSU
>
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Example of H REMD

[Michael Shirts]

Was this not useful?

> Most of this should still be applicable for 5.0.4.
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy

On Mon, Nov 30, 2015 at 8:51 AM, Sanja Zivanovic
 wrote:
> Dear users,
>
> Could you please send me some example of inputs for Hamiltonian Replica
> exchange? It would be great if you can write me command line and attach mdp
> file. I would like to make H REMD for following ligand.
>
> Thanks in advance a lot
>
> Best Regards, Sanja
>
>
> --
> *Sanja Zivanovic*
> PhD Student
> https://www.linkedin.com/in/sanjazivanovic
>
> Molecular Modelling and Bioinformatics Laboratory
> Institute for Research in Biomedicine (IRB Barcelona)
> Parc Científic de Barcelona
> Baldiri Reixac 10, 08028 Barcelona, SPAIN
> Tel. (+34) 9340 37155
> http://www.irbbarcelona.org
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Hamiltonian REMD

[Michael Shirts]

Most of this should still be applicable for 5.0.4.

http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy

On Wed, Nov 25, 2015 at 6:23 AM, Sanja Zivanovic
 wrote:
> Could you please tell me how I can perform Hamiltonian replica exchange
> simulations in Gromacs version 5.0.4 without using PLUMED?
>
> Thanks
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free Energy of Liquid Water

[Michael Shirts]

If ALL the particles are changing with the free energy coupling
parameter, then GROMACS will slow down quite a bit.  If only one
molecule is changing, then it shouldn't be that much slower (20-30%
slower?)  But without .mdp and.top files, it's hard to say exactly
what is happening.


On Wed, Oct 7, 2015 at 11:08 PM, Nathan K Houtz <nho...@purdue.edu> wrote:
> Thanks for explaining, (and Professor Farais de Moura as well) that makes 
> sense.
>
> On a different note, I've noticed a huge difference in performance between my 
> equillibration run and my actual simulation, and I'm wondering if this is 
> normal. For comparison, I ran a couple small (meaningless) simulations where 
> the only difference was the collection of free energy information and Gromacs 
> slows down by a factor of 39. I had been under the impression that I would be 
> able to do each simulation in just under an hour but now it looks like it 
> will take the better part of two days, which is too much. I'll just have to 
> run shorter or fewer simulations if I can't improve that, but I'm hoping I'm 
> just doing something wrong. I'm essentially using modified versions of the 
> .mdp file I found for this tutorial: 
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy.
>  I changed things like turning off the pressure coupling, changing the 
> temperature, etc. The settings for free energy calculations are unchanged. Is 
> the calculation of dh/dl 
 re
>  ally that expensive, or is this abnormal?
>
> Thank you for the help as always. Regards,
> Nathan H.
>
> - Original Message -
> From: "Michael Shirts" <mrshi...@gmail.com>
> To: "Discussion list for GROMACS users" <gmx-us...@gromacs.org>
> Cc: "gromacs org gmx-users" <gromacs.org_gmx-users@maillist.sys.kth.se>
> Sent: Tuesday, October 6, 2015 3:05:38 PM
> Subject: Re: [gmx-users] Free Energy of Liquid Water
>
> For a pure fluid, G  = N \mu.  And \mu = (dG/dN)_(T,P).  So you only
> need to change one molecule to ideal gas to get the change in free
> energy. The free energy of transfer of water from liquid to gas is
> indeed the free energy of solvation of one water molecule in bath of
> water.  So there's a reason why you're just finding tutorials of
> solvation free energies!
>
> On Thu, Oct 1, 2015 at 10:44 PM, Nathan K Houtz <nho...@purdue.edu> wrote:
>> Hi everyone,
>>
>> I would like to use Gromacs to do Thermodynamic Integration (TI) from liquid 
>> water (TIP4P model) to an ideal gas, to find the relative free energy. To do 
>> this, I believe  one generally integrates above the critical point by 
>> increasing the temperature above the critical temperature and then relaxing 
>> the pressure until the system is a diffuse gas. The mdp options 
>> documentation is helpful, and I went through an ethanol solvation tutorial, 
>> but there doesn't appear to be a "pressure-lambda" or a "volume-lambda" 
>> option that I could use to do the second part. How can I get Gromacs to 
>> calculate the dh/dl derivative while relaxing the pressure?
>>
>> In addition, all of the tutorials I found for thermodynamic integration were 
>> for finding solvation free energies. The coulomb and VDW forces are 
>> essentially changed from "completely on" to "completely off". But in my 
>> case, I'd like to change the temperature and pressure between two nonzero 
>> values. I don't want to begin my simulation at 0K and 0atm, but lambda 
>> *must* go from 0 to 1. How can I define both starting and ending points for 
>> the temperature and pressure (or volume, or density)?
>>
>> Thanks for your help!
>> Nathan H.
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
>> mail to gmx-users-requ...@gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>

Re: [gmx-users] Free Energy of Liquid Water

[Michael Shirts]

For a pure fluid, G  = N \mu.  And \mu = (dG/dN)_(T,P).  So you only
need to change one molecule to ideal gas to get the change in free
energy. The free energy of transfer of water from liquid to gas is
indeed the free energy of solvation of one water molecule in bath of
water.  So there's a reason why you're just finding tutorials of
solvation free energies!

On Thu, Oct 1, 2015 at 10:44 PM, Nathan K Houtz  wrote:
> Hi everyone,
>
> I would like to use Gromacs to do Thermodynamic Integration (TI) from liquid 
> water (TIP4P model) to an ideal gas, to find the relative free energy. To do 
> this, I believe  one generally integrates above the critical point by 
> increasing the temperature above the critical temperature and then relaxing 
> the pressure until the system is a diffuse gas. The mdp options documentation 
> is helpful, and I went through an ethanol solvation tutorial, but there 
> doesn't appear to be a "pressure-lambda" or a "volume-lambda" option that I 
> could use to do the second part. How can I get Gromacs to calculate the dh/dl 
> derivative while relaxing the pressure?
>
> In addition, all of the tutorials I found for thermodynamic integration were 
> for finding solvation free energies. The coulomb and VDW forces are 
> essentially changed from "completely on" to "completely off". But in my case, 
> I'd like to change the temperature and pressure between two nonzero values. I 
> don't want to begin my simulation at 0K and 0atm, but lambda *must* go from 0 
> to 1. How can I define both starting and ending points for the temperature 
> and pressure (or volume, or density)?
>
> Thanks for your help!
> Nathan H.
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Converting Amber frcmod and mol2 files to Gromacs

[Michael Shirts]

ParmEd (http://parmed.github.io/ParmEd/html/index.html) is an
interesting new tool that can help automate such conversions so they
don't need to be done as manually in the future.

On Wed, Sep 30, 2015 at 1:37 PM, Simon Dürr  wrote:
> Hi,
>
> we had the pleasure to convert the *.mol2 and *.frcmod for the HEME
> and CYP prosthetic groups from the following paper  for use with
> GROMACS:
> Shahrokh, K., Orendt, A., Yost, G. S. and Cheatham, T. E. (2012),
> Quantum mechanically derived AMBER-compatible heme parameters for
> various states of the cytochrome P450 catalytic cycle. J. Comput.
> Chem., 33: 119–133. doi: 10.1002/jcc.21922
>
> As the protocol for doing so is rather lengthy we want to share our
> protocol with you guys.
> If you can spare some time we would be grateful If some of you could
> check our procedure.
> If you find any errors let us know. However so far the trajectories
> look good (angles and bond distances are good).
>
> After further testing of the trajectories we want to release the
> rtp/hdb/itp - files.
> Where would be the best way to upload them?
>
> The Protocol is as follows:
>
>
> #
> ##  Table of Contents##
> #
> Part 1: RTP-file
> Part 2: ffbondend-file
> Part 3: ffnonbonded
> Part 4: residuetypes.dat
> Part 5: specbond.dat
> Part 6: hdb-file
> Part 7: atomtypes.atp
>
> ##
> ##  PART 1 ##
> ##
>
> #  BEGIN RTP-File ###
>
> [ HEME ]
>  [ atoms]
>  NC   ncr   0.00660 1
>C1C   ccr  -0.01920 2
>C4C   ccr  -0.01930 3
>C2C   ccr   0.01140 4
>C3C   ccr  -0.03570 5
>…
> [ bonds ]
> NC   C1C
> NC   C4C
> NCFE
>[…]
> ## END RTP-FILE #
>
> The information to fill this file is found in the HEM.mol2 file.
>
> ## HEM.mol2 ##
> @ATOM
>   1 NC  2.946000   -0.972000   -0.732000 nc 1 HEM  0.0066 
>   2 C1C 4.282000   -0.756000   -0.518000 cc 1 HEM -0.0192 
>   4 C2C 5.025000   -1.989000   -0.70 cc 1 HEM  0.0114 
>   5 C3C 4.097000   -2.957000   -1.046000 cc 1 HEM -0.0357 
>   […]
> @BOND
> 1 1 2 1
> 2 1 3 1
> 3 128 1
> ## END HEM.mol2 #
>
> To convert to Gromacs
> 1. take the uppercase name in column 2 of the mol2 which is the name
> found in the PDB input file and put it in column 1 of the [atoms]
> section of the .rtp
> 2. Add the lowercase identifier for GAFF (column 6 in the mol2 ) to
> the second column of the atoms section.
> 3. As we don't want to overwrite any of the normal definitions from
> AMBER99SB-ILDN with GAFF definitions we use a unique name for the GAFF
> atoms. Here we just added "r" at the end of each atom name.
> 4. Add the column 6 of the mol2 with the charge as third column to the
> atoms section
> 5 Group  the atoms to charge groups in the fourth column of the
> [atoms]-section. We put every atom in a new group.
>
> Look at the bonds section in the mol2:
> The first column is the index of the bond, the second column is the id
> of the first atom of the bond, the third column the id of the second
> atom of the bond. The last is the bond type (single or double bond)
>
> 7.Take the id of the first atom and replace it with the uppercase
> identifier from the  ATOM section of the mol2. Do the same for the id
> of the second atom.
>
> The columns 1,3,4,5,7 and 8 in the ATOM section are not relevant for Gromacs.
>
>
> ##
> ##  PART 2 ##
> ##
>
>
> # BEGIN FFBONDED.ITP #
> [ bondtypes ]
> ; ij  func   b0  kb
>  fer ncr10.2   83680 ; IC6.frcmod
> [..]
> [ angletypes ]
>  ;  ijk  func   th0   cth
>  ncr fer ndr   1 85.521 1104.576 ; IC6.frcmod
> [..]
>  [ dihedraltypes ]
>  ;i  j   k  lfunc  phasekd  pn
>  X   SH  fer X   9 180.037.656002 ; IC6.frcmod
> [..]
> [ dihedraltypes ]
>  X   orr crr orr 4 180.04.60240 2 ; GAFF.dat improper
>
>  END FFBONDED.ITP #
>
> The information for this file can be found in the IC6.frcmod file OR
> in the gaff.dat
>
>  IC6.frcmod #
> BOND
> fe-nc 100.000   2.000
> [..]
> ANGLE
> nc-fe-nd  132.000  85.521 # average angle
> [..]
> DIHEDRAL
> X -SH-fe-X   1 9.000 180.000   2.000
> [..]
> IMPROPER
>
> NONBON
> fe 1.8 0.01 1.0
>  END IC6.frcmod 
>
>  BEGIN GAFF.dat ###
> […]
> X -c -c -X41.200   180.000   2.000
> X -o -c -o  1.1   180. 2. JCC,7,(1986),230
>  END GAFF.dat #
>
> To convert to Gromacs
> 1. Add the bonds section:
>  a) Copy the two names of the atoms of the bond. Rename to new
> identifier and separate by a single space.
>  b) Set func to 1 

Re: [gmx-users] HREX simulation

[Michael Shirts]

Apologies: I missed the part where it was stated that the simulation
was restarted.  Perhaps I shouldn't have tried to steal a few min to
answer mailing lis questions .  . .

Averaging would give the correct answer.

On Thu, Sep 24, 2015 at 4:49 PM, Mark Abraham <mark.j.abra...@gmail.com> wrote:
> Hi,
>
> It can't be an accumulated matrix. It's not written to a checkpoint, and
> nothing in GROMACS can read one.
>
> At some point, someone needs to roll their sleeves up, write such data to
> fields in TNG format, and stop dumping potentially useful data to a log
> file that'll get truncated upon restart-with-append.
>
> Mark
>
> On Thu, Sep 24, 2015 at 10:42 PM Michael Shirts <mrshi...@gmail.com> wrote:
>
>> > Now, to determine the acceptance rate for the whole 10ns run I should
>> > average these numbers as shown below?
>>
>> No, the matrix at 10ns is the accumulated matrix.  You should use the
>> last matrix.
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] HREX simulation

[Michael Shirts]

Ugh, sorry for missing that.  Yes, the replica exchange transition
matrix is not stored.  Some sort of unification should be possible . .
.

On Thu, Sep 24, 2015 at 5:33 PM, Mark Abraham <mark.j.abra...@gmail.com> wrote:
> Hi,
>
> Expanded ensemble transition matrix and REMD transition matrix are
> different things. Shyno wasn't very specific, but his samples of output
> include the string "Repl" which suggests the latter, whose data is not
> checkpointed.
>
> Mark
>
> On Thu, Sep 24, 2015 at 11:19 PM Michael Shirts <mrshi...@gmail.com> wrote:
>
>> static int do_cpt_df_hist(XDR *xd, int fflags, df_history_t *dfhist, FILE
>> *list)
>>...
>>
>> case edfhTIJ:  ret = do_cpte_nmatrix(xd, cptpEDFH, i,
>> fflags, nlambda, dfhist->Tij, list); break;
>> case edfhTIJEMP:   ret = do_cpte_nmatrix(xd, cptpEDFH, i,
>> fflags, nlambda, dfhist->Tij_empirical, list); break;
>>
>> On Thu, Sep 24, 2015 at 5:06 PM, Mark Abraham <mark.j.abra...@gmail.com>
>> wrote:
>> > Hi,
>> >
>> > Some things are accumulated there, yes. But not the transition history.
>> >
>> > Mark
>> >
>> > On Thu, Sep 24, 2015 at 11:04 PM Michael Shirts <mrshi...@gmail.com>
>> wrote:
>> >
>> >> I take it back, I think it is actually accumulated.  See (in
>> >> checkpoint.cpp):
>> >>
>> >> static int do_cpt_df_hist(XDR *xd, int fflags, df_history_t *dfhist,
>> FILE
>> >> *list)
>> >>
>> >>
>> >> On Thu, Sep 24, 2015 at 4:59 PM, Michael Shirts <mrshi...@gmail.com>
>> >> wrote:
>> >> > Apologies: I missed the part where it was stated that the simulation
>> >> > was restarted.  Perhaps I shouldn't have tried to steal a few min to
>> >> > answer mailing lis questions .  . .
>> >> >
>> >> > Averaging would give the correct answer.
>> >> >
>> >> > On Thu, Sep 24, 2015 at 4:49 PM, Mark Abraham <
>> mark.j.abra...@gmail.com>
>> >> wrote:
>> >> >> Hi,
>> >> >>
>> >> >> It can't be an accumulated matrix. It's not written to a checkpoint,
>> and
>> >> >> nothing in GROMACS can read one.
>> >> >>
>> >> >> At some point, someone needs to roll their sleeves up, write such
>> data
>> >> to
>> >> >> fields in TNG format, and stop dumping potentially useful data to a
>> log
>> >> >> file that'll get truncated upon restart-with-append.
>> >> >>
>> >> >> Mark
>> >> >>
>> >> >> On Thu, Sep 24, 2015 at 10:42 PM Michael Shirts <mrshi...@gmail.com>
>> >> wrote:
>> >> >>
>> >> >>> > Now, to determine the acceptance rate for the whole 10ns run I
>> should
>> >> >>> > average these numbers as shown below?
>> >> >>>
>> >> >>> No, the matrix at 10ns is the accumulated matrix.  You should use
>> the
>> >> >>> last matrix.
>> >> >>> --
>> >> >>> Gromacs Users mailing list
>> >> >>>
>> >> >>> * Please search the archive at
>> >> >>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >> >>> posting!
>> >> >>>
>> >> >>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >> >>>
>> >> >>> * For (un)subscribe requests visit
>> >> >>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
>> or
>> >> >>> send a mail to gmx-users-requ...@gromacs.org.
>> >> >>>
>> >> >> --
>> >> >> Gromacs Users mailing list
>> >> >>
>> >> >> * Please search the archive at
>> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >> posting!
>> >> >>
>> >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >> >>
>> >> >> * For (un)subscribe requests visit
>> >> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
>> or
>> >> send a mail to gmx-users-requ...@gromacs.org.
>> >> --
>> >> Gromacs Users mailing list
>> >>
>

Re: [gmx-users] HREX simulation

[Michael Shirts]

static int do_cpt_df_hist(XDR *xd, int fflags, df_history_t *dfhist, FILE *list)
   ...

case edfhTIJ:  ret = do_cpte_nmatrix(xd, cptpEDFH, i,
fflags, nlambda, dfhist->Tij, list); break;
case edfhTIJEMP:   ret = do_cpte_nmatrix(xd, cptpEDFH, i,
fflags, nlambda, dfhist->Tij_empirical, list); break;

On Thu, Sep 24, 2015 at 5:06 PM, Mark Abraham <mark.j.abra...@gmail.com> wrote:
> Hi,
>
> Some things are accumulated there, yes. But not the transition history.
>
> Mark
>
> On Thu, Sep 24, 2015 at 11:04 PM Michael Shirts <mrshi...@gmail.com> wrote:
>
>> I take it back, I think it is actually accumulated.  See (in
>> checkpoint.cpp):
>>
>> static int do_cpt_df_hist(XDR *xd, int fflags, df_history_t *dfhist, FILE
>> *list)
>>
>>
>> On Thu, Sep 24, 2015 at 4:59 PM, Michael Shirts <mrshi...@gmail.com>
>> wrote:
>> > Apologies: I missed the part where it was stated that the simulation
>> > was restarted.  Perhaps I shouldn't have tried to steal a few min to
>> > answer mailing lis questions .  . .
>> >
>> > Averaging would give the correct answer.
>> >
>> > On Thu, Sep 24, 2015 at 4:49 PM, Mark Abraham <mark.j.abra...@gmail.com>
>> wrote:
>> >> Hi,
>> >>
>> >> It can't be an accumulated matrix. It's not written to a checkpoint, and
>> >> nothing in GROMACS can read one.
>> >>
>> >> At some point, someone needs to roll their sleeves up, write such data
>> to
>> >> fields in TNG format, and stop dumping potentially useful data to a log
>> >> file that'll get truncated upon restart-with-append.
>> >>
>> >> Mark
>> >>
>> >> On Thu, Sep 24, 2015 at 10:42 PM Michael Shirts <mrshi...@gmail.com>
>> wrote:
>> >>
>> >>> > Now, to determine the acceptance rate for the whole 10ns run I should
>> >>> > average these numbers as shown below?
>> >>>
>> >>> No, the matrix at 10ns is the accumulated matrix.  You should use the
>> >>> last matrix.
>> >>> --
>> >>> Gromacs Users mailing list
>> >>>
>> >>> * Please search the archive at
>> >>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >>> posting!
>> >>>
>> >>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >>>
>> >>> * For (un)subscribe requests visit
>> >>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> >>> send a mail to gmx-users-requ...@gromacs.org.
>> >>>
>> >> --
>> >> Gromacs Users mailing list
>> >>
>> >> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>> >>
>> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >>
>> >> * For (un)subscribe requests visit
>> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Amber to GROMACS structure/topology

[Michael Shirts]

http://parmed.github.io/ParmEd/html/index.html

https://code.google.com/p/acpype/

On Thu, Sep 24, 2015 at 5:14 PM, Gustavo Avelar Molina
 wrote:
> Hi,
>
> I would like to know if there is a way to convert Amber structures and
> topologies to the GROMACS format.
>
> Thank you very much for your time.
>
> Best regards,
>
> Gustavo
>
> ==
> Gustavo Avelar Molina, B.Sc. Chem.
> M.Sc. Chem. Student
>
> Department of Chemistry
> Faculty of Philosophy, Sciences and Literature of Ribeirão Preto
> Protein Biochemistry and Biophysics Laboratory
> University of São Paulo, Ribeirão Preto, São Paulo, Brazil
>
> +55 16 994311221 | +55 11 949874141
>
> avelarmolinagust...@gmail.com | gustavoavelarmol...@usp.br
>
> https://lbbpusp.wordpress.com/ |
>
> 
> 
> 
>
> ==
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] HREX simulation

[Michael Shirts]

> Now, to determine the acceptance rate for the whole 10ns run I should
> average these numbers as shown below?

No, the matrix at 10ns is the accumulated matrix.  You should use the
last matrix.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] HREX simulation

[Michael Shirts]

I take it back, I think it is actually accumulated.  See (in checkpoint.cpp):

static int do_cpt_df_hist(XDR *xd, int fflags, df_history_t *dfhist, FILE *list)


On Thu, Sep 24, 2015 at 4:59 PM, Michael Shirts <mrshi...@gmail.com> wrote:
> Apologies: I missed the part where it was stated that the simulation
> was restarted.  Perhaps I shouldn't have tried to steal a few min to
> answer mailing lis questions .  . .
>
> Averaging would give the correct answer.
>
> On Thu, Sep 24, 2015 at 4:49 PM, Mark Abraham <mark.j.abra...@gmail.com> 
> wrote:
>> Hi,
>>
>> It can't be an accumulated matrix. It's not written to a checkpoint, and
>> nothing in GROMACS can read one.
>>
>> At some point, someone needs to roll their sleeves up, write such data to
>> fields in TNG format, and stop dumping potentially useful data to a log
>> file that'll get truncated upon restart-with-append.
>>
>> Mark
>>
>> On Thu, Sep 24, 2015 at 10:42 PM Michael Shirts <mrshi...@gmail.com> wrote:
>>
>>> > Now, to determine the acceptance rate for the whole 10ns run I should
>>> > average these numbers as shown below?
>>>
>>> No, the matrix at 10ns is the accumulated matrix.  You should use the
>>> last matrix.
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>> posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-requ...@gromacs.org.
>>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
>> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fixing periodicity effects on trajectory file

[Michael Shirts]

Maybe the ligand actually is leaving the binding site and moving
around, and that's what your simulation is telling you!  Hard to say
with the posted information.

On Tue, Sep 15, 2015 at 11:24 AM, Homa rooz  wrote:
> Dear Justin
> I have followed your advice but it didn't solve my problem.
> Ligand is still out of binding site.
>
> --
> Homa Ahmadi
> Shahed University
> Faculty of Basic Sciences
> Alternative email> h.ahm...@shahed.ac.ir 
> --
> Gromacs Users mailing list
>
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] HREX simulation

[Michael Shirts]

 The exchange attempt forces replicas to synchronize, and this forces the
 whole simulation system to pay the cost of any localized slow-down (e.g.
 network traffic congestion, internal simulation load imbalance from
 diffusion, even throttling of core clocks from heating effects) at each
 synchronization point. Longer periods between exchange attempts use
 resources more efficiently if the observed variations are approx random.

Excellent point.  I usually do low-core number replica exchange, but
when is running on lots of cores,  I can see how frequent exchange
would slow things down.

 That ought to be a problem, but in fact the replica-exchange implementation
 already synchronizes much more frequently than it should (
 https://gerrit.gromacs.org/#/c/4312/ needs more testing before this is
 solved).

Thanks for pointing that out, I had missed that proposed change.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] HREX simulation

[Michael Shirts]

 I am new to doing replica exchange simulations.
 I am trying to understand the empirical transition matrix printed in the
 log files. I tried exchange intervals 2 ps and 0.2 ps. However, later I
 read gromacs recommends exchange interval not to be under 1ps.

I don't know that there are firm rules for this.  It seems like faster
exchange is generally better, up to a point (if there are specific
code reasons one might thing this is not the case, anyone else on the
list, speak up).

 the system is protein in water and I am only showing parts of the matrix:

 Exchange interval: 0.2 ps
 Repl1234 5

 Repl  0.9528  0.0472  0.  0.  0.
 Repl  0.0472  0.9017  0.0511  0.  0.
 Repl  0.  0.0511  0.8946  0.0543  0.
 Repl  0.  0.  0.0543  0.8898  0.0558

 1. If I look at line 3 Repl  0.0472  0.9017  0.0511 does this mean
 exchange between replica 1, 2 is accepted 4.72 %  and exchange between 2,3
 is accepted 5.11 % ?

Yes.

 Repl1234 5


 Repl  0.9504  0.0496   0.  0.  0.

 Repl  0.0496  0.8982  0.0522   0.  0.

 Repl  0.  0.0522  0.8942  0.0536   0.

 Repl  0.  0.  0.0536  0.  0.0576

 2. With exchange interval 2 ps, I am getting slightly better acceptance
 rate. For example, again looking at line 3 Repl  0.0496  0.8982  0.0522
 the acceptance rate for exchange 1 to 2 is 4.96 %, and he acceptance rate
 for exchange 2 to 3 is 5.22 %

Those are statistically insignificant results.  The success rate
exchange should not change as the interval changes, which is what you
are seeing, which is good.

You probably want to make your spacing narrower in order to improve
the exchange rate; otherwise, replicas move too slowly in state space.

The optimal rate depends on a few factors, but generally you want to
keep the diagonal around 0.60 to 0.70. not 0.9.  One minus the second
smallest eigenvalue of the transition matrix expressed the rate at
which the system moves in state space (see
http://scitation.aip.org/content/aip/journal/jcp/135/19/10.1063/1.3660669).
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] chemical potential differences between two component liquids and crystal structures in gromacs

[Michael Shirts]

Chemical potential in the fluid would be almost exactly the sanme to a
solvation calculation -- as long as you can count on it to be
metastable enough to not crystalize / phase separate during the time
of the simulation.

Chemical potential in the crystal is . . .  harder.  Check out the
pseudo-supercritical path approach of Maginn as currently decent way
to do this.  We are writing a paper applying it in Gromacs, but we
don't have the files and setup ready to share yet.

On Tue, Aug 25, 2015 at 12:25 PM, Nathan K Houtz nho...@purdue.edu wrote:
 Hello everyone,

 I am interested in calculating the change in chemical potential between a 
 metastable liquid, consisting or solvent and solute, and the crystal 
 structure of the pure solute. For instance, one system I am looking at is 
 tetrolic acid in CCL4 as my metastable liquid, and one of the polymorphs of 
 tetrolic acid as the final crystal structure. I have looked through a couple 
 tutorials on free energy calculations using gromacs, but I get the impression 
 that this type of calculation may be less routine than say a solvation free 
 energy which I came across. Can anyone tell me if this type of calculation is 
 possible in gromacs, and where else I could be looking for directions as to 
 how to proceed?

 Thank you,
 Nathan
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Coulomb barriers and Coulomb Softcore Potential

[Michael Shirts]

Adjusting the soft core is a gigantic pain.  I wouldn't recommend it,
and it's likely not necesssary. What settings are you using now?  Note
that the manual describes how it is defined.

Can you post the alchemical-analysis output for N/N_k?  There's an
analysis quirk we are looking at where the correlation times are
measured to be longer than they actually are, though that is when LJ
is changing and coul is turned off.

Look at the time correlation for the variable that is changing -- does
it look stationary, or does it look like it's changing slowly, or has
just a few discrete steps?  If the autocorrelation time is slow, then
N/N_k actually is that long, and the code is doing it's job.

On Mon, Aug 24, 2015 at 9:55 AM, Andreas Mecklenfeld
a.mecklenf...@tu-braunschweig.de wrote:
 Dear Gromacs-Users,

 I want to calculate the solvation free energy of water in an aqueous ionic
 solution and I'm using the Python tool alchemical-analysis.py by
 Klimovich, Shirts and Mobley for evaluation. This tool demonstrates a very
 high N/N_k ratio (up to 7000) while decreasing the electrostatic potential
 (Lennard Jones fully active).

 Considering energetic barriers, I would like to adjust the Coulomb Softcore
 Potential. Does this seem plausible and if so, how is the Coulomb Softcore
 Potential defined in Gromacs?
 Naden and Shirts provide a concept by equation (A.2) in Linear Basis
 Function Approach to Efficient Alchemical Free Energy Calculations. 2.
 Inserting and Deleting Particles with Coulombic Interactions
 (http://pubs.acs.org/doi/abs/10.1021/ct501047e) - is this the formula
 intended?

 Best regards,
 Andreas

 --
 M. Sc. Andreas Mecklenfeld
 Stipendiat

 Technische Universität Braunschweig
 Institut für Thermodynamik
 Hans-Sommer-Straße 5
 38106 Braunschweig
 Deutschland / Germany

 Tel: +49 (0)531 391-2634
 Fax: +49 (0)531 391-7814

 www.ift-bs.de

 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Is dVremain/dl related to fep-lambdas?

[Michael Shirts]

That could be documented better.  It is all of the changes in the
molecule that are not explicitly occurring in their own lambda
variable.  For example, in setup 2, then all the changes have separate
dhdl terms, so there is nothing left to be in dVremain/dl.  Note that
in setup 1) dVremain/dl is just the sum of the other three terms in
setup 2).

On Mon, Aug 24, 2015 at 3:55 AM, Igor Leontyev ileont...@ucdavis.edu wrote:
 Dear gmx-users,
 When I compute with gromacs-5.1 relative free energies using mixed
 topologies (i.e. by setting typeA, typeB, bond kA, kB ant etc) I see that
 some undocumented value dVremain/dl is printed in log- and xvg-files.

 Please, clarify the exact meaning of the dVremain/dl quantity. How it is
 used in dG computation?



 == Some Thoughts ===
 Playing with different setup by setting parameters fep-lambdas,
 coul-lambdas, vdw-lambdas and bonded-lambdas, I figured out that
 dVremain/dl of Setup-1 (when set only fep-lambdas) is equal to sum of all
 dV/dl components of Setup-2 (when all coul-lambdas, vdw-lambdas and
 bonded-lambdas are set explicitly):
 dVremain/dl(1) = Vcoul/dl(2) + dVvdw/dl(2) + dVbonded/dl(2)

 Setup-1)
 # if in MDP-file
; coul-lambdas = 0.0 0.10
; vdw-lambdas = 0.0 0.10
; bonded-lambdas = 0.0 0.10
fep-lambdas = 0.0 0.10
 # then
dVremain/dl
1.52124e+03
 #endif

 Setup-2)
 # if in MDP-file
coul-lambdas = 0.0 0.10
vdw-lambdas = 0.0 0.10
bonded-lambdas = 0.0 0.10
fep-lambdas = 0.0 0.10
 # then
dVremain/dl dVcoul/dl dVvdw/dl dVbonded/dl
0.0e+00 1.56939e+03 -6.97050e+01 2.15596e+01
 #endif
 ===




 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Add new integrator

[Michael Shirts]

Adding an integrator for your own work is hard but doable, especially
if you don't need it parallelized.  Look at the do_update_X routines
for examples.   It will be harder for 5th/6th order since you probably
need values of the forces for times in the past, which are not
automatically stored in gromacs -- you'll have to add arrays to store
them yourself (which will get rather hard with parallelization).

Adding something that other people could use as a general release --
probably not something you would want to try now. Some plans in the
works to make it easier, but it will be a few months at the soonest.

Just wondering about adding predictor/corrector integrators.  They
tend to be quite bad over longer time scales.  Is there a particular
reason?

On Tue, Aug 11, 2015 at 1:24 PM, Eric Smoll ericsm...@gmail.com wrote:
 Hello GROMACS users,

 I am interested in adding a fifth/sixth order predictor/corrector
 integrator to GROMACS. Is this possible? If so, any guidance on which files
 need to be updated is much appreciated.

 Best,
 Eric
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Question about forcefield files

[Michael Shirts]

The identifiers are just strings, so as long as they are unique, there
shouldn't be a problem (I can't remember offhand if capitalization is
respected or not, though).

I'd name them something other than opls_XXX, though, if you are adding
your own, just to avoid confusion of you share the files with others.

On Tue, Jul 28, 2015 at 3:43 PM, Eric Smoll ericsm...@gmail.com wrote:
 Hello GROMACS users,

 In the OPLS atomtypes.atp file, is it problematic to define atom type names
 (e.g., opls_XXX) with nonconsecutive integers? For example,

 opls_1000 ...
 opls_1100 ...
 opls_1050 ...

 Best,
 Eric
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] md-vv and md

[Michael Shirts]

 Has problem been solved about molecules drifting in md-vv integrator?

We think so, but some additional tests are running to address any
other lingering integrator issues.  The fix is in 5.1 beta, but not in
the 4.6 branch yet (it's a pull request, but not yet merged).

 Also can you please check my potential energy second pic in
 http://imgur.com/6aJkRoQfjXGmXu#0
 . http://imgur.com/6aJkRoQfjXGmXu#0
 What could be the reason this big jump in potential energy (integrator :
 md-vv)?
 And why I do not see the same jump by using md integrator?

I don't think there's sufficient information here to answer that question.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free energy calculations (FEP) and soft core potential 1-1-48

[Michael Shirts]

I put in 1-1-48, and although somewhat more efficient (20-30%), it has
a tendency to have floating point errors in single precision. I'll
probably take it out sometime in the next few versions, as the benefit
is not really worth the hassle.  If you are off by a lot, it's almost
certainly not going to help.

Whether something matches experiment is independent of whether the
calculation is performed correctly. If it's off by 2-3 kcal/mol, that
could just be the nature of the problem.  If it's off by 5-10
kcal/mol, especially for a relative calculation, there is probably
something else wrong.  Possibly it has to do with your top file, but
of course hard to say.

100 lambda states is almost certainly more than needed.  When I
disappear 60 heavy atoms, I use only 40 lambda points!  If you are
changing 1-5 heavy atoms, no more than 20 lambda would be needed,
10-12 is probably enough.
If there is a lot of protein motion, 1 ns is likely not long enough.

 I am using umbrella pulling to maintain
an interaction between one residue and the ligand (especially important
during SA) but I remove it for MD because I am afraid it would influence my
results (even though I don't know it for sure).

It probably would affect something during MD, but what the effect
would be depends on many factors.

On Wed, Jul 15, 2015 at 10:20 AM, Julian Zachmann
frankjulian.zachm...@uab.cat wrote:
 Dear Gromacs Users,

 I am running free energy calculations (FEP) to estimate relative binding
 affinities. So far, my results don't match the experimental results (not
 even close), so I must be doing something wrong.

 My protocol is the following:

 - Energy minimisation
 - Leap frog minimisation
 - 0.1ns  NVT
 - 2nsSimulated Annealing up to 350K
 - 0.1ns  NPT
 - 1ns MD production run

 During the whole equilibration time I am using umbrella pulling to maintain
 an interaction between one residue and the ligand (especially important
 during SA) but I remove it for MD because I am afraid it would influence my
 results (even though I don't know it for sure).

 The simulation time for 1ns is quite short of course but I have 'nstdhdl'
 put to 10 to get a lot of output.

 The other settings are pasted below. So far I was using 1-1-6 for the soft
 core potential but because the results were not good, I was thinking to try
 1-1-48. However in this case the simulations always crash. Has anybody run
 successfully 1-1-48 soft core potentials in Gromacs and how did you chose
 your settings? Like me pasted below?

 Any other suggestions about what I am doing wrong?
 Could I maintain the pull distance constraint during the MD simulation?

 Best regards,
 Julian

 free_energy  = yes
 init_lambda_state= $i
 ;   0  1  2  3  4  5  6
  7  8  9 10 11 12 13 14 15 16 17
   18 19 20 21 22 23 24 25 26 27 28
 29 30 31 32 33 34 35 36 37 38
 39 40 41 42 43 44 45 46 47 48 49
   50 51 52 53 54 55 56 57 58 59 60
 61 62 63 64 65 66 67 68 69 70
 71 72 73 74 75 76 77 78 79 80 81
   82 83 84 85 86 87 88 89 90 91 92
 93 94 95 96 97 98 99
 fep_lambdas = 0. 0.0300 0.0600 0.0900 0.1200 0.1500 0.1800
 0.2100 0.2400 0.2700 0.3000 0.3300 0.3600 0.3900 0.4200 0.4500 0.4800
 0.5100 0.5400 0.5700 0.6000 0.6300 0.6600 0.6900 0.7200 0.7500 0.7800
 0.8100 0.8400 0.8700 0.9000 0.9300 0.9600 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1.
 vdw_lambdas = 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0.0500 0.1000 0.1500
 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500
 0.7000 0.7500 0.8000 0.8200 0.8500 0.8700 0.8800 0.8900 0.9000 0.9100
 0.9200 0.9300 0.9350 0.9400 0.9450 0.9500 0.9530 0.9570 0.9600 0.9630
 0.9670 0.9700 0.9730 0.9770 0.9800 0.9820 0.9840 0.9860 0.9880 0.9900
 0.9910 0.9920 0.9930 0.9940 0.9950 0.9955 0.9960 0.9965 0.9970 0.9975
 0.9980 0.9983 0.9987 0.9990 0.9992 0.9993 0.9994 0.9995 0.9996 0.9997
 0.9998 0. 1.
 ; 1-1-6
 sc-coul= no
 sc-alpha = 0.5
 

Re: [gmx-users] Number of Compute Nodes in Gromacs Expanded Ensemble Simulations

[Michael Shirts]

Expanded ensemble doesn't require any special number of nodes (states are
visited sequentially, not in parallel) so the rules are the same as for any
other free energy simulation.

On Thu, Jun 18, 2015 at 4:46 AM, Andreas Mecklenfeld 
a.mecklenf...@tu-braunschweig.de wrote:

 Hi,

 I'm new to expanded ensemble simulations in Gromacs. Is there any rule of
 thumb for the most efficient number of compute nodes?

 Thanks in advance
 Andreas

 --
 M. Sc. Andreas Mecklenfeld
 Stipendiat

 Technische Universität Braunschweig
 Institut für Thermodynamik
 Hans-Sommer-Straße 5
 38106 Braunschweig
 Deutschland / Germany

 Tel: +49 (0)531 391-2634
 Fax: +49 (0)531 391-7814

 http://www.ift-bs.de

 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] compressibility

[Michael Shirts]

The 'compressibility' in the pressure coupling is not actually the
compressibility of the fluid; it's another way to write the mass of the
fictitious piston.  It affects the inertia of the piston and the period of
the oscillations of the box, but won't affect the magnitude of fluctuations
of the box.  Having a value near the compressibility of the fluid, I
believe, minimizes the chances of unstable simulations and weird
resonances, but is not required.

On Fri, Jun 12, 2015 at 3:12 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Hi,

 I simulate a molecule in NPT ensemble, which I can not find any
 experimental data about it's compressibility, what should I put in .mdp
 file for it's compressibility in pressure coupling?

 Best regards
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Different Cv and Cp

[Michael Shirts]

I can't diagnose what is wrong in every case; I can just say that
statistical mechanics says that both methods for calculating Cv and both
methods for calculating Cp should agree within noise.

On Wed, Jun 10, 2015 at 8:54 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Dear Michael,

 Could you please comment on my last question, thanks a lot.
 I have noticed that when I run both NVT and NPT simulations from a same
 .gro file (obtained from energy minimization) I obtain same Cv and Cp for
 two ensembles. However, so far I was running NVT after NPT (meaning I used
 .gro file obtained from NPT as initial configuration for NVT, since I
 wanted to have correct volume). The second way gives that big difference in
 Cv and Cp.
 I run a long simulations in all cases, but why both ways produce such
 difference?



 Best regards


 On Tue, Jun 9, 2015 at 5:45 PM, Michael Shirts mrshi...@gmail.com wrote:

  The variance formula is derived from the derivative formula + the
  assumption the distribution in Boltzmann, so they must agree if the
  distribution is Boltzmann.
 
  On Tue, Jun 9, 2015 at 11:43 AM, Faezeh Pousaneh fpoosa...@gmail.com
  wrote:
 
   Thank you so much for the reply.
  
   Yes, I use contact applied pressure and I am careful about units. I
  checked
   and average enthalpy and U are close, meaning that PV is negligible.
 But
   the point is variance of enthalpy in NPT differs from variance of
 energy
  in
   NVT and that causes the difference.
  
   You had given me an article showing that for example for Benzene the
   difference in Cv and Cp is 25%, and here Ii get similar ( my molecule
 is
   also carbon ring). But still I can answer why two both ways does not
 give
   same Cv-Cp.
   I will follow your suggestion.
  
  
  
   Best regards
  
  
   On Tue, Jun 9, 2015 at 5:29 PM, Michael Shirts mrshi...@gmail.com
  wrote:
  
If the simulation are generating configurations with the Boltzmann
probability distribution, the results should the same up to error.
   
Cv and Cp should not be exactly the same, though for liquids at room
temperature, they are pretty close (look up the precise numbers for
 the
fluid you are interested in).
   
Are you calculating enthalpy as U + PV, where P is the constant
 APPLIED
pressure, not the instantaneous pressure, and PV is in same units?
  Since
the PV term should be quite low for liquids, the two heat capacities
   should
be relatively close (within noise -- fluctuation based calculations I
   think
are noisier).
   
Else you need to check if the Boltzmann distributions are being
  correctly
generated: See the code and paper describing it linked here:
https://github.com/shirtsgroup/checkensemble
   
   
   
   
On Tue, Jun 9, 2015 at 11:23 AM, Faezeh Pousaneh 
 fpoosa...@gmail.com
wrote:
   
 Dear Michael,

 Can I ask a question concerning your previous email,
 I followed

 Cv = (dU/dT)_V

 Cp = (dH/dT)_P

 for my lutidine molecule, and I get same values for Cv and Cp. But
   when I
 test with

 kB T^2 c_P = Var (Enthalpy)
 kB T^2 c_V = Var (Energy)

 I get 40 J/mol.K difference in Cv and Cp.

 Mean that fluctuation play big role. Which way of checking I can
  rely?


 Best regards


 On Tue, May 26, 2015 at 3:38 PM, Michael Shirts 
 mrshi...@gmail.com
 wrote:

  By definition (more fundamental that fluctuation formulas)
 
  Cv = (dU/dT)_V
 
  Cp = (dH/dT)_P
 
  Run two simulations at different T and estimate the derivatives.
 
  On Tue, May 26, 2015 at 5:12 AM, Faezeh Pousaneh 
   fpoosa...@gmail.com
  wrote:
 
   Dear Michael,
  
   I still would like to know what was your method you mentioned
 on
   last
   paragraph, just for learning:
  
   ''Also, to be sure, you should double check by calculating both
   heat
   capacities by finite difference formulas as well with two
   simulations
 at
   T+dt/2 and T-dt/2 -- if the fluctuation and finite difference
   resutls
  don't
   agree within propagated error, then something is off.''
  
   ?
   thanks
  
  
   Best regards
  
  
   On Mon, May 25, 2015 at 5:10 PM, Faezeh Pousaneh 
fpoosa...@gmail.com
   wrote:
  
Dear Andre,
   
thank you for the link, you are probably right, It seems that
  my
  molecule
has the difference Cp-Cv in the same range as benzene (since
 it
   has
  also
ring structure).
   
   
Best regards
   
   
On Mon, May 25, 2015 at 4:44 PM, Faezeh Pousaneh 
 fpoosa...@gmail.com
wrote:
   
Dear Michael,
   
I use  Parrinello-Rahman for barostat and v-rescale for
thermostat.
   
Sorry, could you explain more the second paragraph please? I
  did
not
  get

Re: [gmx-users] Different Cv and Cp

[Michael Shirts]

The variance formula is derived from the derivative formula + the
assumption the distribution in Boltzmann, so they must agree if the
distribution is Boltzmann.

On Tue, Jun 9, 2015 at 11:43 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Thank you so much for the reply.

 Yes, I use contact applied pressure and I am careful about units. I checked
 and average enthalpy and U are close, meaning that PV is negligible. But
 the point is variance of enthalpy in NPT differs from variance of energy in
 NVT and that causes the difference.

 You had given me an article showing that for example for Benzene the
 difference in Cv and Cp is 25%, and here Ii get similar ( my molecule is
 also carbon ring). But still I can answer why two both ways does not give
 same Cv-Cp.
 I will follow your suggestion.



 Best regards


 On Tue, Jun 9, 2015 at 5:29 PM, Michael Shirts mrshi...@gmail.com wrote:

  If the simulation are generating configurations with the Boltzmann
  probability distribution, the results should the same up to error.
 
  Cv and Cp should not be exactly the same, though for liquids at room
  temperature, they are pretty close (look up the precise numbers for the
  fluid you are interested in).
 
  Are you calculating enthalpy as U + PV, where P is the constant APPLIED
  pressure, not the instantaneous pressure, and PV is in same units?  Since
  the PV term should be quite low for liquids, the two heat capacities
 should
  be relatively close (within noise -- fluctuation based calculations I
 think
  are noisier).
 
  Else you need to check if the Boltzmann distributions are being correctly
  generated: See the code and paper describing it linked here:
  https://github.com/shirtsgroup/checkensemble
 
 
 
 
  On Tue, Jun 9, 2015 at 11:23 AM, Faezeh Pousaneh fpoosa...@gmail.com
  wrote:
 
   Dear Michael,
  
   Can I ask a question concerning your previous email,
   I followed
  
   Cv = (dU/dT)_V
  
   Cp = (dH/dT)_P
  
   for my lutidine molecule, and I get same values for Cv and Cp. But
 when I
   test with
  
   kB T^2 c_P = Var (Enthalpy)
   kB T^2 c_V = Var (Energy)
  
   I get 40 J/mol.K difference in Cv and Cp.
  
   Mean that fluctuation play big role. Which way of checking I can rely?
  
  
   Best regards
  
  
   On Tue, May 26, 2015 at 3:38 PM, Michael Shirts mrshi...@gmail.com
   wrote:
  
By definition (more fundamental that fluctuation formulas)
   
Cv = (dU/dT)_V
   
Cp = (dH/dT)_P
   
Run two simulations at different T and estimate the derivatives.
   
On Tue, May 26, 2015 at 5:12 AM, Faezeh Pousaneh 
 fpoosa...@gmail.com
wrote:
   
 Dear Michael,

 I still would like to know what was your method you mentioned on
 last
 paragraph, just for learning:

 ''Also, to be sure, you should double check by calculating both
 heat
 capacities by finite difference formulas as well with two
 simulations
   at
 T+dt/2 and T-dt/2 -- if the fluctuation and finite difference
 resutls
don't
 agree within propagated error, then something is off.''

 ?
 thanks


 Best regards


 On Mon, May 25, 2015 at 5:10 PM, Faezeh Pousaneh 
  fpoosa...@gmail.com
 wrote:

  Dear Andre,
 
  thank you for the link, you are probably right, It seems that my
molecule
  has the difference Cp-Cv in the same range as benzene (since it
 has
also
  ring structure).
 
 
  Best regards
 
 
  On Mon, May 25, 2015 at 4:44 PM, Faezeh Pousaneh 
   fpoosa...@gmail.com
  wrote:
 
  Dear Michael,
 
  I use  Parrinello-Rahman for barostat and v-rescale for
  thermostat.
 
  Sorry, could you explain more the second paragraph please? I did
  not
get
  the method. What I checked so far is checking if gromacs
 correctly
gives
  Cv,Cp= Var(Energy or Enthalpy)/kBT^2 , and I find that it gives.
 
 
 
 
  Best regards
 
 
  On Mon, May 25, 2015 at 4:11 PM, Michael Shirts 
  mrshi...@gmail.com
   
  wrote:
 
  Are you running with the Berendsen thermostat or barostat?  The
gromacs
  g_energy functions for heat capacity use the fluctuation
 formula,
   and
 the
  fluctuations with both of these algorithms are wrong (as should
  be
  printed
  in the log file warning message). Make sure you use
ensemble-preserving
  thermostats if you want fluctuation properties.
 
  Also, to be sure, you should double check by calculating both
  heat
  capacities by finite difference formulas as well with two
   simulations
 at
  T+dt/2 and T-dt/2 -- if the fluctuation and finite difference
   resutls
  don't
  agree within propagated error, then something is off.
 
 
  On Mon, May 25, 2015 at 5:59 AM, Faezeh Pousaneh 
fpoosa...@gmail.com
  wrote:
 
   Hi,
  
   I do not know why I obtain two difference cp and cv from NVT

Re: [gmx-users] Different Cv and Cp

[Michael Shirts]

If the simulation are generating configurations with the Boltzmann
probability distribution, the results should the same up to error.

Cv and Cp should not be exactly the same, though for liquids at room
temperature, they are pretty close (look up the precise numbers for the
fluid you are interested in).

Are you calculating enthalpy as U + PV, where P is the constant APPLIED
pressure, not the instantaneous pressure, and PV is in same units?  Since
the PV term should be quite low for liquids, the two heat capacities should
be relatively close (within noise -- fluctuation based calculations I think
are noisier).

Else you need to check if the Boltzmann distributions are being correctly
generated: See the code and paper describing it linked here:
https://github.com/shirtsgroup/checkensemble




On Tue, Jun 9, 2015 at 11:23 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Dear Michael,

 Can I ask a question concerning your previous email,
 I followed

 Cv = (dU/dT)_V

 Cp = (dH/dT)_P

 for my lutidine molecule, and I get same values for Cv and Cp. But when I
 test with

 kB T^2 c_P = Var (Enthalpy)
 kB T^2 c_V = Var (Energy)

 I get 40 J/mol.K difference in Cv and Cp.

 Mean that fluctuation play big role. Which way of checking I can rely?


 Best regards


 On Tue, May 26, 2015 at 3:38 PM, Michael Shirts mrshi...@gmail.com
 wrote:

  By definition (more fundamental that fluctuation formulas)
 
  Cv = (dU/dT)_V
 
  Cp = (dH/dT)_P
 
  Run two simulations at different T and estimate the derivatives.
 
  On Tue, May 26, 2015 at 5:12 AM, Faezeh Pousaneh fpoosa...@gmail.com
  wrote:
 
   Dear Michael,
  
   I still would like to know what was your method you mentioned on last
   paragraph, just for learning:
  
   ''Also, to be sure, you should double check by calculating both heat
   capacities by finite difference formulas as well with two simulations
 at
   T+dt/2 and T-dt/2 -- if the fluctuation and finite difference resutls
  don't
   agree within propagated error, then something is off.''
  
   ?
   thanks
  
  
   Best regards
  
  
   On Mon, May 25, 2015 at 5:10 PM, Faezeh Pousaneh fpoosa...@gmail.com
   wrote:
  
Dear Andre,
   
thank you for the link, you are probably right, It seems that my
  molecule
has the difference Cp-Cv in the same range as benzene (since it has
  also
ring structure).
   
   
Best regards
   
   
On Mon, May 25, 2015 at 4:44 PM, Faezeh Pousaneh 
 fpoosa...@gmail.com
wrote:
   
Dear Michael,
   
I use  Parrinello-Rahman for barostat and v-rescale for thermostat.
   
Sorry, could you explain more the second paragraph please? I did not
  get
the method. What I checked so far is checking if gromacs correctly
  gives
Cv,Cp= Var(Energy or Enthalpy)/kBT^2 , and I find that it gives.
   
   
   
   
Best regards
   
   
On Mon, May 25, 2015 at 4:11 PM, Michael Shirts mrshi...@gmail.com
 
wrote:
   
Are you running with the Berendsen thermostat or barostat?  The
  gromacs
g_energy functions for heat capacity use the fluctuation formula,
 and
   the
fluctuations with both of these algorithms are wrong (as should be
printed
in the log file warning message). Make sure you use
  ensemble-preserving
thermostats if you want fluctuation properties.
   
Also, to be sure, you should double check by calculating both heat
capacities by finite difference formulas as well with two
 simulations
   at
T+dt/2 and T-dt/2 -- if the fluctuation and finite difference
 resutls
don't
agree within propagated error, then something is off.
   
   
On Mon, May 25, 2015 at 5:59 AM, Faezeh Pousaneh 
  fpoosa...@gmail.com
wrote:
   
 Hi,

 I do not know why I obtain two difference cp and cv from NVT and
  NPT
 simulations.
 What I do is, I take 1000 lutidne molecules, and I do firstly an
   energy
 minimization with steep integrator, then NPT simulation at T=300
  and
P=1
 atm for 10ns, (I obtain Cp= 230), then I run NVT for 10 ns with
  same
mdp
 file except no pressure coupling, and with initial .gro file
  obtained
from
 NPT run, (I obtain Cv=180).
 Does some one know where is my mistake? (In both runs, I obtain
 Cv
   and
Cp
 from g_energy and in different time intervals and after
  equilibrited
time)


 Best regards
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List
 before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit

 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
   or
 send a mail to gmx-users-requ...@gromacs.org.

--
Gromacs Users mailing list
   
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!
   
* Can't post? Read http

Re: [gmx-users] Different Cv and Cp

[Michael Shirts]

By definition (more fundamental that fluctuation formulas)

Cv = (dU/dT)_V

Cp = (dH/dT)_P

Run two simulations at different T and estimate the derivatives.

On Tue, May 26, 2015 at 5:12 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Dear Michael,

 I still would like to know what was your method you mentioned on last
 paragraph, just for learning:

 ''Also, to be sure, you should double check by calculating both heat
 capacities by finite difference formulas as well with two simulations at
 T+dt/2 and T-dt/2 -- if the fluctuation and finite difference resutls don't
 agree within propagated error, then something is off.''

 ?
 thanks


 Best regards


 On Mon, May 25, 2015 at 5:10 PM, Faezeh Pousaneh fpoosa...@gmail.com
 wrote:

  Dear Andre,
 
  thank you for the link, you are probably right, It seems that my molecule
  has the difference Cp-Cv in the same range as benzene (since it has also
  ring structure).
 
 
  Best regards
 
 
  On Mon, May 25, 2015 at 4:44 PM, Faezeh Pousaneh fpoosa...@gmail.com
  wrote:
 
  Dear Michael,
 
  I use  Parrinello-Rahman for barostat and v-rescale for thermostat.
 
  Sorry, could you explain more the second paragraph please? I did not get
  the method. What I checked so far is checking if gromacs correctly gives
  Cv,Cp= Var(Energy or Enthalpy)/kBT^2 , and I find that it gives.
 
 
 
 
  Best regards
 
 
  On Mon, May 25, 2015 at 4:11 PM, Michael Shirts mrshi...@gmail.com
  wrote:
 
  Are you running with the Berendsen thermostat or barostat?  The gromacs
  g_energy functions for heat capacity use the fluctuation formula, and
 the
  fluctuations with both of these algorithms are wrong (as should be
  printed
  in the log file warning message). Make sure you use ensemble-preserving
  thermostats if you want fluctuation properties.
 
  Also, to be sure, you should double check by calculating both heat
  capacities by finite difference formulas as well with two simulations
 at
  T+dt/2 and T-dt/2 -- if the fluctuation and finite difference resutls
  don't
  agree within propagated error, then something is off.
 
 
  On Mon, May 25, 2015 at 5:59 AM, Faezeh Pousaneh fpoosa...@gmail.com
  wrote:
 
   Hi,
  
   I do not know why I obtain two difference cp and cv from NVT and NPT
   simulations.
   What I do is, I take 1000 lutidne molecules, and I do firstly an
 energy
   minimization with steep integrator, then NPT simulation at T=300 and
  P=1
   atm for 10ns, (I obtain Cp= 230), then I run NVT for 10 ns with same
  mdp
   file except no pressure coupling, and with initial .gro file obtained
  from
   NPT run, (I obtain Cv=180).
   Does some one know where is my mistake? (In both runs, I obtain Cv
 and
  Cp
   from g_energy and in different time intervals and after equilibrited
  time)
  
  
   Best regards
   --
   Gromacs Users mailing list
  
   * Please search the archive at
   http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
   posting!
  
   * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  
   * For (un)subscribe requests visit
   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
 or
   send a mail to gmx-users-requ...@gromacs.org.
  
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 
 
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Different Cv and Cp

[Michael Shirts]

Are you running with the Berendsen thermostat or barostat?  The gromacs
g_energy functions for heat capacity use the fluctuation formula, and the
fluctuations with both of these algorithms are wrong (as should be printed
in the log file warning message). Make sure you use ensemble-preserving
thermostats if you want fluctuation properties.

Also, to be sure, you should double check by calculating both heat
capacities by finite difference formulas as well with two simulations at
T+dt/2 and T-dt/2 -- if the fluctuation and finite difference resutls don't
agree within propagated error, then something is off.


On Mon, May 25, 2015 at 5:59 AM, Faezeh Pousaneh fpoosa...@gmail.com
wrote:

 Hi,

 I do not know why I obtain two difference cp and cv from NVT and NPT
 simulations.
 What I do is, I take 1000 lutidne molecules, and I do firstly an energy
 minimization with steep integrator, then NPT simulation at T=300 and P=1
 atm for 10ns, (I obtain Cp= 230), then I run NVT for 10 ns with same mdp
 file except no pressure coupling, and with initial .gro file obtained from
 NPT run, (I obtain Cv=180).
 Does some one know where is my mistake? (In both runs, I obtain Cv and Cp
 from g_energy and in different time intervals and after equilibrited time)


 Best regards
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: Bad performance in free energy calulations

[Michael Shirts]

 somewhat off-topic but I wonder why in your free energy protocol you
 only vary the vdW and electrostatic lambdas.  What about the others?
 Your mutation also transforms bonded terms and masses.



Minor point - if you are taking the difference between two mutations (say,
with and without a ligand), then it's better to leave the masses untouched
-- any contribution will cancel.  There are some weird issues with changing
masses in the constraint terms and their contribution to the free energies
that need to be accounted for post-hoc.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: Bad performance in free energy calulations

[Michael Shirts]

 As Michael said it is better to leave most things untouched and just change
 VDW and LJ - at least this is what I have been reading so far.


Well, I didn't quite say this.  Changes in the bonded terms do not entirely
cancel.  Changes in masses do 100% cancel because of the seperability of
momenta and position, and the masses only contributed to the total energy
through the momenta.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Fwd: Bad performance in free energy calulations

[Michael Shirts]

Yep, that's what I generally do.  Almost all alchemical changes involve a
difference in two calculations (since the alchemical change itself is
unphysical).  Even one-calculation solvation free energy calculations are
actually calculating the difference in free energy from liquid to vapor
state.

On Tue, May 19, 2015 at 11:58 AM, Hannes Loeffler 
hannes.loeff...@stfc.ac.uk wrote:

 On Tue, 19 May 2015 11:37:11 -0400
 Michael Shirts mrshi...@gmail.com wrote:

  
   somewhat off-topic but I wonder why in your free energy protocol you
   only vary the vdW and electrostatic lambdas.  What about the others?
   Your mutation also transforms bonded terms and masses.
  
 
 
  Minor point - if you are taking the difference between two mutations
  (say, with and without a ligand), then it's better to leave the
  masses untouched -- any contribution will cancel.  There are some
  weird issues with changing masses in the constraint terms and their
  contribution to the free energies that need to be accounted for
  post-hoc.

 So you suggest to leave the mass lambdas simply at their initial (or
 final) value to be safe in case constraints are involved (and I am too
 lazy to check)?
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Thermostats and MSD

[Michael Shirts]

For #1, see: http://pubs.acs.org/doi/abs/10.1021/ct400109a

For #2, I think you'll have to figure that out for yourself.

On Wed, Mar 11, 2015 at 5:42 AM, sujithkakkat . sujithk...@gmail.com
wrote:

 Dear all,

   I am simulating a system of small hydrophobic solutes ( CH4)  dissolved
 in water. Simulations are performed in the NPT simulation, and I would like
 to determine the diffusion coefficients from MSD plots. I have two related
 questions;

 (1)  Is the thermostat going to seriously affect the diffusivity of the
 molecules. I am using Nose-Hoover thermostat with a time constant of 0.2
 ps. Is the choice of time constant very crucial while studying diffusivity.
 Should I avoid thermostat and go for NVE simulations for studying
 diffusivity.

 (2) From radial distribution functions, it is found that there is
 association between the hydrophobic solvents. I would like to know, if in
 case a small cluster of the hydrophobic solutes are formed, then the
 whether the rotational motion of these clusters contribute to dispalcement
 in the MSD calculation.

 Please comment.

 Regards,
 Sujith.
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] bennet acceptance ratio

[Michael Shirts]

I didn't write the g_bar script, so I can't answer all of the questions
about it.

Equation 3 in that paper is a general equation true for any weight function
alpha.  It is not Bennett's acceptance ratio. Equation 4 is the one that's
actually Bennett's equations, which chooses an alpha that gives a minimum
variance answer.

So you'll have to ask others about the exact usage of g_mbar. The tutorials
on the Alchemistry.org web page show how to use the tool I wrote for
analyzing free energy calculations with gromacs.  Gromacs has plenty of
functionality for printing out the energy differences between states.

On Wed, Oct 29, 2014 at 8:03 AM, Ahmet yıldırım ahmedo...@gmail.com wrote:

 Dear Prof.Shirts,

 I have a question about your paper titled An Introduction to Best
 Practices in Free Energy
 Calculations.There are two equation (equation 3 and equation 4) of Bennet
 Acceptance Ratio in your paper.


 Which equation does Gromacs g_bar tool use from your paper? If it uses
 equation 3, then g_bar computes potential energy difference between two
 states. This potential energy difference means free energy difference
 between two states?
 Also I saw the same equation on alchemistry.org

 Another interesting thing, the dhdl.xvg files includes dh/dl values. And
 when g_bar calculates the free energy, it uses these files (dHdl.xvg).
 Whereas g_bar needs the potential energy differences between states/window?

 H=U+K, where H:Hamiltonian, U:Potential energy, T:Kinetic energy. You know
 g_bar gives it(Delta-G) on output screen. T=3/2kT, T=constant temperature,
 then H=U right. Then dH equals to dU?

 Thanks in advance

 --
 Ahmet Yıldırım
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Conserved quantity in NPT simulations

[Michael Shirts]

Those conserved quantities are indeed included in the reported
information.  For thermostats/barostats where there is such a quantity,
it's included in the .edr file.

On Fri, Oct 24, 2014 at 4:23 AM, Paolo Franz paolo.fr...@gmail.com wrote:

 Hi,
 I am wondering how to test the conserved quantity in an NPT simulation fron
 the ener.edr file. As you well know for npt with a Nose-Hoover thermostat
 and a Parrinello-Rahman barostat, the conserved quantity is not just K+P,
 but contains additional terms function of the additional npt variables. Are
 those energies placed somewhere? I am using a somewhat new potential for a
 system and I want to know how well it conserves the energy at constant
 pressure and temperature.
 Thnaks in advance!

 Paolo
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Velocity Verlet integrator

[Michael Shirts]

Because the 'start' of the vv integrator step is halfway through the loop.
This is a byproduct of 1) putting leapfrog and velocity verlet in the same
loop and 2) minimizing communication and output.  It is not as elegant as
it should be.  There are efforts to clean this up, but it's a lot of
reorganization, and has gone slowly.

On Wed, Oct 15, 2014 at 8:46 AM, Mario Fernández Pendás mariof...@gmail.com
 wrote:

 Yes, I understand that. But my question is more about why the two velocity
 updates are implemented before the position update and not the other way
 round?

 From the theoretical point of view I would think more in one of the
 following schemes:


1. Calculate: [image: \vec{v}\left(t + \tfrac12\,\Delta t\right) =
\vec{v}(t) + \tfrac12\,\vec{a}(t)\,\Delta t\,]
2. Calculate: [image: \vec{x}(t + \Delta t) = \vec{x}(t) +
\vec{v}\left(t + \tfrac12\,\Delta t\right)\, \Delta t\,]
3. Derive [image: \vec{a}(t + \Delta t)] from the interaction potential
using [image: \vec{x}(t + \Delta t)]
4. Calculate: [image: \vec{v}(t + \Delta t) = \vec{v}\left(t +
\tfrac12\,\Delta t\right) + \tfrac12\,\vec{a}(t + \Delta t)\Delta t,]



1. Calculate: [image: \vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}(t)\,
\Delta t+\tfrac12 \,\vec{a}(t)\,\Delta t^2]
2. Derive [image: \vec{a}(t + \Delta t)] from the interaction potential
using [image: \vec{x}(t + \Delta t)]
3. Calculate: [image: \vec{v}(t + \Delta t) = \vec{v}(t) +
\tfrac12\,\left(\vec{a}(t)+\vec{a}(t + \Delta t)\right)\Delta t\,]


 This is why my confusion arises.


 2014-10-15 14:05 GMT+02:00 Justin Lemkul jalem...@vt.edu:

 
 
  On 10/15/14 7:30 AM, Mario Fernández Pendás wrote:
 
  Dear all,
 
  I am still interested in some integrator related issues.
 
  I understand that the easiest way to implement velocity Verlet was to
  split
  the updates in two updates. But I don't understad the order of those
  updates.
  I mean why there are two updates for velocities and then the update for
  positions?
  My intuitive idea would be to update first one half for velocities,
 then a
  full step for positions and, finally, using these new positions the
 second
  half for velocities.
 
 
  Yes, there are two separate half-step updates for velocities.  The
  comments in the md.cpp code are quite verbose if you want to trace
 through.
 
  -Justin
 
  --
  ==
 
  Justin A. Lemkul, Ph.D.
  Ruth L. Kirschstein NRSA Postdoctoral Fellow
 
  Department of Pharmaceutical Sciences
  School of Pharmacy
  Health Sciences Facility II, Room 629
  University of Maryland, Baltimore
  20 Penn St.
  Baltimore, MD 21201
 
  jalem...@outerbanks.umaryland.edu | (410) 706-7441
  http://mackerell.umaryland.edu/~jalemkul
 
  ==
 
  --
  Gromacs Users mailing list
 
  * Please search the archive at http://www.gromacs.org/
  Support/Mailing_Lists/GMX-Users_List before posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Velocity Verlet integrator

[Michael Shirts]

Yes, I've been using the theory here:

http://arxiv.org/abs/1301.3800

Which describes how to concantenate integrator steps in a formal way.

I can say that the time savings you get by concatenating integrators is
VERY small.  The only time it is nonnegligible is when there is a LOT of
communication, and even then, there are better ways to make simulations
faster.  It's not an area where there is a lot of improvement that can be
made.

When doing temperature and pressure control, there are many cases that you
cannot really join the steps as well.

On Wed, Oct 15, 2014 at 10:05 AM, Mario Fernández Pendás 
mariof...@gmail.com wrote:

 Thank you very much Professor Shirts.

 I have these doubts because I am trying to implement new integrators based
 in the concatenation of two VV steps to make a single step. The idea
 follows the integrators suggested in
 http://web.mit.edu/~ripper/www/research/efficient_md_integrators.pdf

 This is why it is important for me to know where each step starts and
 finishes.

 2014-10-15 15:52 GMT+02:00 Michael Shirts mrshi...@gmail.com:

  Because the 'start' of the vv integrator step is halfway through the
 loop.
  This is a byproduct of 1) putting leapfrog and velocity verlet in the
 same
  loop and 2) minimizing communication and output.  It is not as elegant as
  it should be.  There are efforts to clean this up, but it's a lot of
  reorganization, and has gone slowly.
 
  On Wed, Oct 15, 2014 at 8:46 AM, Mario Fernández Pendás 
  mariof...@gmail.com
   wrote:
 
   Yes, I understand that. But my question is more about why the two
  velocity
   updates are implemented before the position update and not the other
 way
   round?
  
   From the theoretical point of view I would think more in one of the
   following schemes:
  
  
  1. Calculate: [image: \vec{v}\left(t + \tfrac12\,\Delta t\right) =
  \vec{v}(t) + \tfrac12\,\vec{a}(t)\,\Delta t\,]
  2. Calculate: [image: \vec{x}(t + \Delta t) = \vec{x}(t) +
  \vec{v}\left(t + \tfrac12\,\Delta t\right)\, \Delta t\,]
  3. Derive [image: \vec{a}(t + \Delta t)] from the interaction
  potential
  using [image: \vec{x}(t + \Delta t)]
  4. Calculate: [image: \vec{v}(t + \Delta t) = \vec{v}\left(t +
  \tfrac12\,\Delta t\right) + \tfrac12\,\vec{a}(t + \Delta t)\Delta
 t,]
  
  
  
  1. Calculate: [image: \vec{x}(t + \Delta t) = \vec{x}(t) +
  \vec{v}(t)\,
  \Delta t+\tfrac12 \,\vec{a}(t)\,\Delta t^2]
  2. Derive [image: \vec{a}(t + \Delta t)] from the interaction
  potential
  using [image: \vec{x}(t + \Delta t)]
  3. Calculate: [image: \vec{v}(t + \Delta t) = \vec{v}(t) +
  \tfrac12\,\left(\vec{a}(t)+\vec{a}(t + \Delta t)\right)\Delta t\,]
  
  
   This is why my confusion arises.
  
  
   2014-10-15 14:05 GMT+02:00 Justin Lemkul jalem...@vt.edu:
  
   
   
On 10/15/14 7:30 AM, Mario Fernández Pendás wrote:
   
Dear all,
   
I am still interested in some integrator related issues.
   
I understand that the easiest way to implement velocity Verlet was
 to
split
the updates in two updates. But I don't understad the order of those
updates.
I mean why there are two updates for velocities and then the update
  for
positions?
My intuitive idea would be to update first one half for velocities,
   then a
full step for positions and, finally, using these new positions the
   second
half for velocities.
   
   
Yes, there are two separate half-step updates for velocities.  The
comments in the md.cpp code are quite verbose if you want to trace
   through.
   
-Justin
   
--
==
   
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
   
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
   
jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
   
==
   
--
Gromacs Users mailing list
   
* Please search the archive at http://www.gromacs.org/
Support/Mailing_Lists/GMX-Users_List before posting!
   
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
   
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
 or
send a mail to gmx-users-requ...@gromacs.org.
   
   --
   Gromacs Users mailing list
  
   * Please search the archive at
   http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
   posting!
  
   * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  
   * For (un)subscribe requests visit
   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
   send a mail to gmx-users-requ...@gromacs.org.
  
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org

Re: [gmx-users] Use NVT to mimic NVE

[Michael Shirts]

 I guess, if I pick a coupling constant that is just small enough to keep the
energy conserved, I would get a NVT simulation that is as close as a NVE
simulation as possible.

 Is this correct?

Yes, but then at that point the thermostat isn't actually thermostatting.
The Bussi comment is merely to show that his thermostat correctly reduces
to Newton's law in the limit, not that it would be useful to run it in that
limit.

On Mon, Oct 13, 2014 at 10:28 AM, Johnny Lu johnny.lu...@gmail.com wrote:

 On page 014101-3, the Bussi paper (http://dx.doi.org/10.1063/1.2408420)
 mentioned: On the other hand, for coupling constant approaching
 infinity,the Hamiltonian dynamics is recovered.
 Does that means that for a large enough coupling constant, the velocities
 are nearly not rescaled, and the dynamics (like rate of motion) would be
 same as that of NVE?

 A larger coupling constant, means a smaller diffusion coefficient in the
 axillary dynamics by equation 6.

 While the effects of the velocity rescaling at each step will accumulate, a
 larger coupling constant means the thermostat perturb less of the dynamics,
 and the resulting dynamics is closer to a NVE simulation.
 There is no worry that the thermostat would suddenly rescale the dynamics
 every x step, because in the procedure of the thermostat, the velocities
 are rescaled every step, regardless of the coupling constant.

 I guess, if I pick a coupling constant that is just small enough to keep
 the energy conserved, I would get a NVT simulation that is as close as a
 NVE simulation as possible.

 Is this correct?
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Use NVT to mimic NVE

[Michael Shirts]

 The point was that if the energy leaves the system, then the system
is never in equilibrium. Not vice versa, please.

This is incorrect logic.  f there are errors in the integrator, then
F=gradient of the potential is not true. If that is not true, energy is not
conserved, and energy can leave or enter the system.


On Sun, Oct 12, 2014 at 1:58 PM, Dr. Vitaly Chaban vvcha...@gmail.com
wrote:

 The point was that if the energy leaves the system, then the system is
 never in equilibrium. Not vice versa, please.



 Dr. Vitaly V. Chaban

 Виталий Витальевич ЧАБАН


 On Sun, Oct 12, 2014 at 7:20 PM, Johnny Lu johnny.lu...@gmail.com wrote:
  Why non-equilibrium can cause lost of total energy in NVE simulation?
  (baring the case that the non-equilibrium causes very large values that
  causes over flow / getting out of the stable range of algorithms of
  gromacs).
 
  A closed system should have constant energy, regardless of whether it is
 at
  equilibrium.
 
  I guess people use the word converge instead of equilibrium in
 molecular
  dynamics simulation of protein, possibly because the simulation may never
  reach equilibrium with current computing power.
  http://pubs.rsc.org/en/content/articlepdf/2012/cp/c2cp23961b
 
  I may try older version of gromacs and/or not using constraint, after
  checking the bug list and performance. I had never thought about that
 route.
 
 
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Use NVT to mimic NVE

[Michael Shirts]

 If we use thermostats, the thermostats will correct the energy drop/raise
caused by numerical error (which acts as a fluctuating source/sink that may
not even follow the normal distribution), without causing additional
artifacts.

Thermostats correct SOME of the artifacts.  It is not always clear which
ones.

 Regardless, thermostats causes artifacts, especially in dynamics and
rates.

In SOME cases it may be negligible.  In others, not.  It is very hard to
say.

See our paper linked before http://pubs.acs.org/doi/abs/10.1021/ct400109a
Effects of Temperature Control Algorithms on Transport Properties and
Kinetics in Molecular Dynamics Simulations  which covers most short-time
scale effects.

Whether the thermostats that do not affect short time scale properties to
any noticable degree still affect long-time scale properties through some
sort of error accumulation is not known, nor which properties they might
affect.


On Sun, Oct 12, 2014 at 3:07 PM, Johnny Lu johnny.lu...@gmail.com wrote:

 oh, i still mix up things. the numerical error might act as dissipation.

 On Sun, Oct 12, 2014 at 3:01 PM, Johnny Lu johnny.lu...@gmail.com wrote:

  Sorry that I misunderstood your post.
 
  Can compensating energy lost caused by numerical error with thermostat
  cause error in the probability distribution of thermodynamic variables
 and
  distribution of conformations?
 
  Below is just a guess. I don't understand the fluctuation dissipation
  theorem enough at this point.
 
  I would guess so.. unless the fluctuation of the velocity of the atom
  doesn't follow maxwell-Boltzmann distribution (I guess that can happen
 in a
  crystal made of harmonic oscillator of a single frequency, or if some
 part
  of the protein has low degree of freedom. This is a flaw found in
  nose-hover thermostats.)
 
  Otherwise, some thermostats already try to provide the correct
  fluctuation, so that the distribution of thermodynamic variables is
 correct
  (as in fluctuation-dissipation theorem). If we use thermostats, the
  thermostats will correct the energy drop/raise caused by numerical error
  (which acts as a fluctuating source/sink that may not even follow the
  normal distribution), without causing additional artifacts.
 
  Regardless, thermostats causes artifacts, especially in dynamics and
 rates.
 
  On Sun, Oct 12, 2014 at 2:07 PM, Michael Shirts mrshi...@gmail.com
  wrote:
 
   The point was that if the energy leaves the system, then the system
  is never in equilibrium. Not vice versa, please.
 
  This is incorrect logic.  f there are errors in the integrator, then
  F=gradient of the potential is not true. If that is not true, energy is
  not
  conserved, and energy can leave or enter the system.
 
 
  On Sun, Oct 12, 2014 at 1:58 PM, Dr. Vitaly Chaban vvcha...@gmail.com
  wrote:
 
   The point was that if the energy leaves the system, then the system is
   never in equilibrium. Not vice versa, please.
  
  
  
   Dr. Vitaly V. Chaban
  
   Виталий Витальевич ЧАБАН
  
  
   On Sun, Oct 12, 2014 at 7:20 PM, Johnny Lu johnny.lu...@gmail.com
  wrote:
Why non-equilibrium can cause lost of total energy in NVE
 simulation?
(baring the case that the non-equilibrium causes very large values
  that
causes over flow / getting out of the stable range of algorithms of
gromacs).
   
A closed system should have constant energy, regardless of whether
 it
  is
   at
equilibrium.
   
I guess people use the word converge instead of equilibrium in
   molecular
dynamics simulation of protein, possibly because the simulation may
  never
reach equilibrium with current computing power.
http://pubs.rsc.org/en/content/articlepdf/2012/cp/c2cp23961b
   
I may try older version of gromacs and/or not using constraint,
 after
checking the bug list and performance. I had never thought about
 that
   route.
   
   
   
   --
   Gromacs Users mailing list
  
   * Please search the archive at
   http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
   posting!
  
   * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  
   * For (un)subscribe requests visit
   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
   send a mail to gmx-users-requ...@gromacs.org.
  
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail

Re: [gmx-users] Use NVT to mimic NVE

[Michael Shirts]

Following up on fluctuations and ensembles:

The following paper discusses tests of whether an ensemble has the correct
distribution of energies, including fluctuations of all scale:

Simple quantitative tests to validate sampling from thermodynamic
ensembles
http://pubs.acs.org/doi/abs/10.1021/ct300688p

Turns out that v-rescale behaves very well, at least for fluctuations of
the total energy.

In the limit of large scale, for NVE and NVT (if carried out at the
temperature where E_NVT = E_NVE, then all system averages (including
pressure, etc) will be equal.  But the fluctuations (which include 2nd
derivative properties like heat capacity, etc). will not be equal.  So it
depends on what you are interested in.

If you are interested in thermodynamics, just go ahead and run the ensemble
you are interested in.  If you are interested in dynamics, you can check
out this paper which partly answers the questions of what the artifacts are
created by thermostats.  Again, vrescale is quite good --  Langevin or
anything with stochastic noise is not so good.

Effects of Temperature Control Algorithms on Transport Properties and
Kinetics in Molecular Dynamics Simulations
http://pubs.acs.org/doi/abs/10.1021/ct400109a


Best,

Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu

(434)-243-1821
On Fri, Oct 10, 2014 at 7:40 PM, Antonio Baptista bapti...@itqb.unl.pt
wrote:

 On Fri, 10 Oct 2014, Mark Abraham wrote:

  On Fri, Oct 10, 2014 at 7:41 PM, Johnny Lu johnny.lu...@gmail.com
 wrote:

  Hi.

 Is it a good idea to mimic NVE by a NVT simulation with a large
 temperature
 coupling time constant, to reduce the effect of the thermostat ?


 Not if your observable is the total energy. But for most interesting
 observables, this is more than fine.


  If I use V-rescale thermostat, what artifacts will the simulation get if
 I
 use a large coupling time (say, 500 ps) and single precision gromacs ?


 Probably nobody knows :-) Thermodynamic ensemble differences disappear as
 the system becomes large. But quantifying that from converged ensembles
 for
 a given observable is work.

 Mark


 Note also that, unlike averages, fluctuations are ensemble-dependent.
 Although v-rescaling (as most other thermostats) does indeed sample from
 the NVE ensemble in the mathematical limit of infinite coupling constant,
 this does not imply that, in practice, fluctuations converge to those in
 the NVE ensemble for large but finite coupling constants.

 As you can check in the v-rescaling paper, fluctuations were found to be
 very robust to changes of the coupling parameter (at least for the kinetic
 and potential energies). This is generally a good thing, but maybe not what
 you want to run a fake NVE simulation as you described. So, if you are
 interested in the fluctuations of some properties, you will have to check
 what is in practice their asymptotic behavior.

 Best,
 Antonio





  Thank you.
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

  --
 Gromacs Users mailing list

 * Please search the archive at http://www.gromacs.org/
 Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

  --
 Gromacs Users mailing list

 * Please search the archive at http://www.gromacs.org/
 Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Reminder about Gromacs user survey

[Michael Shirts]

Dear Gromacs users/developers,

I'd like to remind you about the currently running Gromacs feature and user
survey Mark Abraham emailed out last week. The survey can be found at:

https://www.surveymonkey.com/s/BD63LRP

We've had a lot of great feedback so far, but would love to hear your
thoughts and opinions as well!

Those filling out the survey will be eligible for the previously detailed
random drawing of two $50 Amazon gift certificates.  Due to Labor Day
weekend in the US cutting the beginning of this week short, we will extend
the original eligibility deadline from 11:59 p.m. GMT tonight (September
2nd) until 11:59 p.m. Pacific Daylight Time tomorrow, Wednesday, September
3rd.

Best,

Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Gromacs feature and usage survey

[Michael Shirts]

Survey is now with 100% more African representation! Thanks for catching
that smaller but important category.  I'm not sure it's possible to edit
old entries, but it's there for the rest of the week.


On Tue, Aug 26, 2014 at 8:24 AM, MUSYOKA THOMMAS 
mutemibiochemis...@gmail.com wrote:

 Hey Mark,
 In ref to part 3 of the survey, I am surprised that it does not capture
 Africa.


 On Tue, Aug 26, 2014 at 1:39 PM, Mark Abraham mark.j.abra...@gmail.com
 wrote:

  Dear Gromacs users/developers,
 
  Gromacs developers Profs. Michael Shirts and Peter Kasson at the
 University
  of Virginia are conducting a survey on Gromacs usage and to gauge
 interest
  in potential additional features to Gromacs.  The survey can be found at:
 
  https://www.surveymonkey.com/s/BD63LRP
 
  Survey participants submitted by 11:59 p.m. GMT on Tuesday, September
 2nd,
  will be eligible for a random drawing of two $50 Amazon gift certificates
  (one response per person, excluding participants in any location where
 such
  a drawing runs afoul of any local law, and those in the Shirts and Kasson
  labs, Amazon country site of your choice).
 
  Thanks!
 
  Mark Abraham
  Gromacs development manager
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] How To Turn Off Electrostatics Linearly in TI Simulations

[Michael Shirts]

 However, if I were to linearly decouple the Coulombic interactions, does
it mean I slowly reduce the cut-off for  Coulombic interaction to 0
angstroms?

No, it is reducing the contribution of q_i q_j energy terms in the atoms
being coupled to zero.


On Sun, Jul 13, 2014 at 10:59 PM, Yip Yew Mun yipy0...@gmail.com wrote:

 I read the mdp options for sc-coul, and I understand that by setting to
 the default option, I’m not applying the soft core energy to the Coulomb
 interaction. However, if I were to linearly decouple the Coulombic
 interactions, does it mean I slowly reduce the cut-off for Coulombic
 interaction to 0 angstroms? If not, what is the technical action to execute
 in the mdp file?

 Thanks.

 On 7/13/14, 10:10 AM, Yip Yew Mun wrote:
  Hi all, I?m planning to run thermodynamic integration (TI) simulations
 with a few lambda values. I read articles and found out that the one of the
 better ways to deal with these lambda simulations is:
 
  1. Disable all electrostatics linearly.
  2. Disable LJ terms by soft core.
 
  Therefore, what I wish to ask is how do I disable electrostatics
 ?linearly? in the mdp files of these lambda values?
 

 Linear decoupling of Coulombic interactions is the default behavior;
 please read
 the description of the options in the manual (sc-coul being the most
 directly
 relevant to this issue).

 -Justin

 --
 ==

 Justin A. Lemkul, Ph.D.
 Ruth L. Kirschstein NRSA Postdoctoral Fellow

 Department of Pharmaceutical Sciences
 School of Pharmacy
 Health Sciences Facility II, Room 601
 University of Maryland, Baltimore
 20 Penn St.
 Baltimore, MD 21201

 jalem...@outerbanks.umaryland.edu | (410) 706-7441
 http://mackerell.umaryland.edu/~jalemkul

 ==

 Yip Yew Mun (Mr) | PhD Research Scholar | Division of Chemistry 
 Biological Chemistry
 School of Physical  Mathematical Sciences | Nanyang Technological
 University | Singapore 639798
 Tel: (+65) 97967803 | Email: yipy0...@e.ntu.edu.sg | GMT+8h

 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] NPT MD and Berendsen

[Michael Shirts]

A nonisotropic system will not necessarily end up with exactly the applied
pressure -- things like surface tension and stresses/strains end up getting
mixed in and throwing things off.


On Mon, Jun 23, 2014 at 3:32 PM, Chetan Mahajan chetanv...@gmail.com
wrote:

 Thanks.

 Average was taken over 1000 ps. Error estimate is 2.9 which is smaller than
 average pressure of 19.87 bar. Here is the complete statistics:

 Energy  Average Err.Est.  RMSD  Tot-Drift

 ---
 Temperature  299.968  0.0131.98625   0.062225  (K)
 Pressure19.86992.9   726.92413.5906
  (bar)

 Interestingly, average pressure is close to one (0.86 or 0.8 bar) when
 refcoord scaling of -all is used for position restraints. Above simulation
 where average pressure is so large was carried out with refcoord-scaling of
 -com option. I am simulating TiO2 crystal solvated by water, formate and
 sodium ion. Now, one important observation is that with -com option,
 crystal slab moves a distance during change of box dimensions in NPT
 equilibration, whereas slab DOES NOT move with -all option. However, -com
 option is more favorable since very less water manage to get on the side
 surfaces of slab compared to that in -all option, where lot of water get on
 the side which is not desirable.

 What is happening here? Is it because of slab movement that average
 pressure is annoyingly large in simulations with -com option? what can be
 done to rectify or is it okay to proceed anyways?

 Thanks a lot!
 regards
 Chetan



 On Mon, Jun 23, 2014 at 1:47 AM, Francis Jing francij...@gmail.com
 wrote:

  How long was the average taken over? If the error is larger than the
  pressure itself, you cannot make judgment.
  Generally Berendsen is only inaccurate about the fluctuations, not the
 mean
  value.
 
 
  Francis
 
 
  On Mon, Jun 23, 2014 at 1:50 PM, Chetan Mahajan chetanv...@gmail.com
  wrote:
 
   Dear All,
  
   I am using Berendsen barostat and thermostat for NPT equilibration MD
 run
   at 1 atm and 300 K. While it runs fine, average pressure at the end of
  the
   run is 19 bar. I am confused whether this is way off 1 atm target, or
 is
   justified recalling gromacs manual that says Berendsen barostat does
 not
   yield correct thermodynamic ensemble.
  
   IN addition, initial box elongates to have space at both ends in
   z-direction as can be seen in following link:
  
   https://www.dropbox.com/sh/dahnyy3zy6vb1xh/AAD5ppxrm18rB3DHZ0Qji0qva
  
   Again, I am not sure if this is a case for concern or just artifact of
   using Berendsen.
  
   Is Parrinello-Rahman better than Berendsen, even for initial NPT
   equilibration?
  
   Thanks
   Chetan
   --
   Gromacs Users mailing list
  
   * Please search the archive at
   http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
   posting!
  
   * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  
   * For (un)subscribe requests visit
   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
   send a mail to gmx-users-requ...@gromacs.org.
  
 
 
 
  --
  Zhifeng (Francis) Jing
  Graduate Student in Physical Chemistry
  School of Chemistry and Chemical Engineering
  Shanghai Jiao Tong University
  http://sun.sjtu.edu.cn
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] large free energy difference of mutating ATP to GTP in solution using FEP method in Gromacs-4.6.5

[Michael Shirts]

There are many ways for free energy calculations to go wrong. Look at
http://www.alchemistry.org for some discussions.

One thing I noticed (I didn't have a chance to look at everything) -
you only talk about one simulation. When doing classical simulations,
then you can't just mutate atoms, because classical simulations
neglect quantum mechanical terms. You can only look at differences in
alchemical simulations.  For example, you can look at the differences
between solvation free energies by performing the mutation in both
water and vacuum.  You have to think VERY carefully about the
thermodynamic end states that you are using.

On Sat, May 24, 2014 at 1:18 AM, dbaogen dbao...@gmail.com wrote:
 Dear all,

I want to calculate the free energy difference of mutating ATP to 
 GTP in solution using FEP method. Firstly, the hybrid topology and structure 
 files for A (ATP) and B (GTP) state using dummy atoms were constructed. 
 Secondly,  the system is running for 10 ns to reach an equilibrium state.  
 And then the structure at 10 ns is as the starting structure to carry out FEP 
 calculation. In the course of FEP, the coulomb interaction was firstly 
 changed, and then the VDW interactions. Total 32 lambda points are set in the 
 mdp file shown in the following:
 integrator   = sd
 nsteps   = 10
 dt  = 0.002
 nstenergy= 1000
 nstlog   = 5000
 nstcalcenergy= 100
 nstcomm  = 1
 cutoff-scheme= group
 rlist= 1.2
 dispcorr = EnerPres
 vdw-type = switch
 ;cut-off lengths
 rvdw = 1.1
 rvdw-switch  = 1
 ; Coulomb interactions
 coulombtype  = pme
 rcoulomb = 1.2
 fourierspacing   = 0.12
 ; Constraints
 constraints  = all-bonds
 ; set temperature to 310K
 tcoupl   = v-rescale
 tc-grps  = system
 tau-t= 1.0
 ref-t= 310
 ; pressure control
 pcoupl  = Parrinello-Rahman
 ref-p= 1
 compressibility= 4.5e-5
 tau-p= 0.5

 ; and set the free energy parameters
 free-energy  = yes
 sc-power = 1
 sc-sigma = 0.3
 sc-alpha = 0.5
 sc-coul  = no
 sc-r-power   = 6
 ; we still want the molecule to interact with itself at lambda=0
 couple-intramol  = no
 couple-lambda0   = vdw-q
 couple-lambda1   = vdw-q
 ; $LAMBDA$ changed from 0 to 32
 init-lambda-state= $LAMBDA$
 nstdhdl  = 100
 calc-lambda-neighbors= 1
 fep-lambdas =  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.03 0.05 0.1 
 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 
 1.0
 ;change electrostatic and then LJ interaction
 coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.00 1.00 1.0 
 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 
 1.0
 vdw-lambdas  = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.03 0.05 0.1 
 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 
 1.0

 atom part of hybrid topology file is :
 [ atoms ]
 ;   nr   type  resnr residue  atom   cgnr charge   mass  typeB
 chargeB
  1 O3  1RAGO1G  1   -0.95260  16.00
  2  P  1RAG PG  21.26500  31.00
  3 O3  1RAGO2G  3   -0.95260  16.00
  4 O3  1RAGO3G  4   -0.95260  16.00
  5 OS  1RAGO3B  5   -0.53220  16.00
  6  P  1RAG PB  61.38520  31.00
  7 O2  1RAGO1B  7   -0.88940  16.00
  8 O2  1RAGO2B  8   -0.88940  16.00
  9 OS  1RAGO3A  9   -0.56890  16.00
 10  P  1RAG PA 101.25320  31.00
 11 O2  1RAGO1A 11   -0.87990  16.00
 12 O2  1RAGO2A 12   -0.87990  16.00
 13 OS  1RAGO5* 13   -0.59870  16.00
 14 CT  1RAGC5* 140.05580  12.00
 15 H1  1RAGH50 150.06790   1.008000
 16 H1  1RAGH51 160.06790   1.008000
 17 CT  1RAGC4* 170.10650  12.00
 18 H1  1RAGH40 180.11740   1.008000
 19 OS  1RAGO4* 19   -0.35480  16.00
 20 CT  1RAGC1* 200.03940  12.00   CT 
 0.01910  12.00
 21 H2  1RAGH10 210.20070   1.008000   H2 
 0.20060   1.008000
 22 N*  1RAG N9 22   -0.02510  

Re: [gmx-users] binding sites with MD

[Michael Shirts]

Look for papers by D. E. Shaw and Gianni di Fabritiis.  They have done
this, but it generally takes at least microseconds or milliseconds to
converge.

On Sat, May 24, 2014 at 6:24 AM, Nidhi Katyal nidhikatyal1...@gmail.com wrote:
 Hi all,

 I would like to ask if unbiased MD in nanoseconds time scale be used to
 find the potential binding sites of ligand with protein?

 I have simulated for 50ns, 1:14 and 1:24 protein:ligand simultaneously with
 random placement of ligand initially. In the time interval between 40 to
 50ns, movement of ligand molecules can be seen around certain sites only
 for both the runs. Can these sites be considered as binding sites? Also,
 ligand molecules are involved in both hydrophobic and hydrogen bonding
 interactions with these sites.

 Thanks in advance.
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] binding sites with MD

[Michael Shirts]

Vmd can convert trajectories, as can other tools. Google is your friend!

Sent from my iPhone

 On May 24, 2014, at 15:13, Albert mailmd2...@gmail.com wrote:
 
 this plug-in looks interesting.
 
 However, as far as I saw from the tutorial that even we run simulation from 
 Gromacs, this plugin can only import NAMD psf and DCD to analysis. That's  a 
 headache issue. Can you improve this vmd plugin so that it can also read 
 Gromacs format files directly?
 
 best
 Albert
 
 
 On 05/24/2014 07:02 PM, Tom Dodson wrote:
 Nidhi,
 
 You might be interested in this tool:http://prody.csb.pitt.edu/drugui/
 Here are slide from a talk about the method:
 http://mmbios.org/images/workshops/HandsOn2014/Lecture3.pdf
 And here is a tutorial:
 http://prody.csb.pitt.edu/tutorials/drugui_tutorial/drugui_tutorial.pdf
 
 I am sure that you could use GROMACS to do the simulated annealing step
 instead of NAMD.
 
 Tom
 
 -- 
 Gromacs Users mailing list
 
 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
 
 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Reproducibility in free energy calculations.

[Michael Shirts]

Mean of the three calculations is 4.35, standard error is 0.19.  So this
seems statistically consistent.  If you want to get more consistent
results, you'll need to run longer until the estimated statistical
uncertainty is lower.  Also, the estimated uncertainty is often an
underestimate (BAR often underestimates by 20-30%).


On Sat, Apr 19, 2014 at 9:05 AM, Justin Lemkul jalem...@vt.edu wrote:



 On 4/19/14, 1:16 AM, sujithkakkat . wrote:

 Hello Justin,

 I tried the free energy calculations in your tutorial at different
 sets
 of conditions. Also I repeated a simulation at the same conditions with
 same parameters, thrice on the same computer on same number of processors.
 The results in the three cases where different (4.16 +/- 0.19 , 4.36
 +/-0.14 , 4.54 +/-0.15 kJ/mol ) . I am not surprised, since I thought that
 this might happen to a small system (241 water + 1 methane). However, I
 want to be sure that this is OK. I guess in cases like this I should go
 for
 an average value of many runs.


 Offhand, I would have expected slightly better agreement between the runs,
 but there are tons of factors at play, so the outcome seems very plausible
 and reasonably consistent.

 -Justin

 --
 ==

 Justin A. Lemkul, Ph.D.
 Ruth L. Kirschstein NRSA Postdoctoral Fellow

 Department of Pharmaceutical Sciences
 School of Pharmacy
 Health Sciences Facility II, Room 601
 University of Maryland, Baltimore
 20 Penn St.
 Baltimore, MD 21201

 jalem...@outerbanks.umaryland.edu | (410) 706-7441
 http://mackerell.umaryland.edu/~jalemkul

 ==

 --
 Gromacs Users mailing list

 * Please search the archive at http://www.gromacs.org/
 Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Velocity Verlet integrator

[Michael Shirts]

Putting both velocity Verlet and leapfrog Verlet both in Gromacs turns
out to be non-trivial for the bookkeeping.  The easiest way to do this
was split the velocity Verlet updates.

Also, the additional computational cost of two half steps for
velocities is trivial compared to the cost of the forces for almost
all systems.

On Sun, Mar 23, 2014 at 8:39 AM, Mario Fernández Pendás
mariof...@gmail.com wrote:
 Dear all,

 In terms of computational efficiency, why the velocity Verlet integrator is
 implemented in GROMACS in one full step for positions and two half steps
 for velocities?

 Would it be more efficient to merge the second halft step for velocities
 with the first half step of the following scheme, ie, integrating in one
 full step for velocities?

 Thank you very much,
 Mario Fernánez-Pendás
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Rigid water model

[Michael Shirts]

Setting the compressibility in gromacs actually is just as way adjust the 
stiffness of the barostat. 

Berendsen thermostat guarantees a Cp that is too low. Try v-rescale or 
nose-Hoover.

See http://pubs.acs.org/doi/abs/10.1021/ct300688p

For more info.

Sent from my iPhone

 On Mar 16, 2014, at 8:33, Mark Abraham mark.j.abra...@gmail.com wrote:
 
 On Mar 14, 2014 11:06 PM, Michelle Aranha mara...@utk.edu wrote:
 
 Hi!
 I tried to reproduce results for rigid water with 895 molecules~1000 kg/m3
 density, 10 ns NPT simulation (298 K  1 bar)with Parinello Rahman
 Barostat
 and Berendsen thermostat.
 
 Using a thermostat that produces the wrong velocity distribution is asking
 for trouble. No idea if that's your main problem, though.
 
 Here is my result for SPC water
 
 Kappa (J/m3)5.689E-10
 Cp  22.8054
 Cp- Cv  5.94102
 
 The density, diffusion coefficient, the average configuration energy agree
 with previous results. However the heat capacity value is quite low as
 compared to 75 J/mol K of real water or even that has been previously
 found
 with water models. I compared my values with a 2 ns and 10 ns simulations,
 but the heat capacity value shows no real improvement. I do not know much
 about quantum corrections hence I have not tried it. Will applying quantum
 corrections give the correct heat capacity?
 The isothermal compressibility  reported by GROMACS is in J/m3 and the
 adiabatic bulk modulus is reported in m^3/J, could that be a mistake and
 should it be rather in  (J/m^3)^-1 and J/m3 respectively.
 And, why does the isothermal compressibility change during the course of
 simulation if it has been fixed in the mdp file.
 
 A parameter for an algorithm and an observed property are two different
 things!
 
 Mark
 
 
 Best wishes,
 Michelle
 
 
 --
 View this message in context:
 http://gromacs.5086.x6.nabble.com/Rigid-water-model-tp5015158.html
 Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
 --
 Gromacs Users mailing list
 
 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
 
 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.
 -- 
 Gromacs Users mailing list
 
 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
 
 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Dubious results with NPT

[Michael Shirts]

No full solutions to this problem -- I'll write a few notes.

* With large systems run for relatively small amounts of time, the
ensemble could be statistically indistinguishble from the true
distribution even if the averages don't line up correctly. See:

http://pubs.acs.org/doi/abs/10.1021/ct300688p

For ways to examine the full distribution of volumes.

* Parrinello-Rahman isn't quite correct in GROMACS because virial and
kinetic energy are offset by 1/2 a step, though getting something more
like 1.2 atm when running at 1 atm is more common.

* Note that with systems that are heterogeneous in direction (not
uniform in x y and z), any errors in the pressure may be magnified.

* If you have a property that changes when the pressure changes from 1
to 5 atm, odds are that molecular simulation is not the right way to
study it -- the force fields are not that accurate.


On Tue, Feb 25, 2014 at 6:48 AM, Dr. Vitaly Chaban vvcha...@gmail.com wrote:
 in that case the average is reasonable, but you should state in the
 manuscript that all the reported results correspond to 5.15 bar.


 Dr. Vitaly V. Chaban


 On Tue, Feb 25, 2014 at 10:42 AM, sujithkakkat . sujithk...@gmail.com wrote:
 Hi,

   In fact the average pressure 5.15 bar was obtained  for the case where the
 reference pressure was 5 bar.  I had mentioned it in one of the previous
 mails. But still the average value is not that good.

I had run an earlier simulation with the same parameters for a pure water
 system for 5ns and reference pressure 5bar, and things worked fine there
 with the average pressure at 5.15bar.


 Sujith.


 On Tue, Feb 25, 2014 at 2:29 PM, Dr. Vitaly Chaban vvcha...@gmail.com
 wrote:

 your average pressure is the pressure you should report in the
 publication. if you got 5 bars instead of 1 bar, you should write

 I simulated the system at 5 bars


 Dr. Vitaly V. Chaban


 On Tue, Feb 25, 2014 at 6:01 AM, sujithkakkat . sujithk...@gmail.com
 wrote:
  On Mon, Feb 24, 2014 at 5:37 AM, sujithkakkat . sujithk...@gmail.com
  wrote:
  Hello,
 
Thank you both for  the comments. I am using gromos96 forcefield . I
  read
  a little bit  and as you said the nonbonded cutoff has to be higher.
  The tau_p=5ps was chosen , since the manual mentions that the value has
  to
  be raised by 4-5 times on going from berendsen to parrinello-rahman
  barostat, though I did not completely follow the reasons behind it. I
  will
  try with lower values.
 
 I had run an earlier simulation with the same parameters for a pure
  water
  system for 5ns and reference pressure 5bar, and things worked fine
  there
  with the average pressure at 5.15bar.
 
 
  No, this average pressure is not fine.
 
 
 
   Why isn't it fine? I admit it is not very good. I guess due to huge
  fluctuations the average pressure is tend to vary a bit to either side
  of
  the reference value, unless I continue the equilibration for a very long
  time. Also the average was calculated with g_energy, which I guess
  considered the entire simulation, where the first few time steps, when
  the
  system is away from equilibrium,  would have poor pressure values with
  respect to reference , which would reflect in the average pressure
  calculated.
  Also using g_analyze to calculate average pressure between different
  times
  during the simulation did not give the same value.  So I am not sure if
  I
  should aim at an average pressure which is equal to the reference
  pressure.
 
  Thanks,
  Sujith.
 
 
 
  On Mon, Feb 24, 2014 at 4:26 PM, Justin Lemkul jalem...@vt.edu wrote:
 
 
 
  On 2/23/14, 11:37 PM, sujithkakkat . wrote:
 
  Hello,
 
 Thank you both for  the comments. I am using gromos96 forcefield .
  I
  read
  a little bit  and as you said the nonbonded cutoff has to be higher.
  The tau_p=5ps was chosen , since the manual mentions that the value
  has
  to
  be raised by 4-5 times on going from berendsen to parrinello-rahman
  barostat, though I did not completely follow the reasons behind it. I
  will
  try with lower values.
 
  I had run an earlier simulation with the same parameters for a
  pure
  water system for 5ns and reference pressure 5bar, and things worked
  fine
  there with the average pressure at 5.15bar. I guess the sigma  for the
  case
  of water was low and therefore the small cut-off of 0.8nm did not
  matter.
  However the case of cyclohexane alone remains to be tried.
 
  I guess Dr Vitaly was saying about using a switch/shift function.
  I will try the simulation with the new settings and see.
 
 
  The Gromos force fields were parametrized with a plain cutoff for van
  der
  Waals.  They do not require switching or shifting.  It's true that
  plain
  cutoffs are a bit harsh, but switching and shifting have their own
  issues,
  and one should always stick to the parametrization methods unless there
  is
  demonstrable evidence that the new method is better.  I have never
  known
  Gromos96 force fields to need anything other than 

Re: [gmx-users] aMD for Gromacs?

[Michael Shirts]

Accelerated MD does not give correct dynamics either.  Systems are not
trapped in minima the appropriate amount of time; you can't directly
map the faster movement through space back to the movement in the
original space unless the cutoff energy is so low as to not really
provide any advantage.

I would say that both metadynamics and accelerated MD are not
necessarily the most precise of names for methods out there . . .

On Tue, Feb 25, 2014 at 8:28 PM, Francis Jing francij...@gmail.com wrote:
 In fact there is one major difference: metadynamics does not have dynamics!
 Although a recent version have dealt with it, its implementation may
 require some time.
 Anyway, it's fine if you do not care about dynamics.

 Francis


 On Wed, Feb 26, 2014 at 6:39 AM, Michael Shirts mrshi...@gmail.com wrote:

 Your best bet is probably something like the PLUMED plugin to do
 metadynamics.  aMD raises up the potential if it is below a given
 cutoff; metadynamics puts in bumps in the potential where you have
 visited.  There are definitely subtle differences, but the general
 effects are similar.

 On Tue, Feb 25, 2014 at 1:31 PM, Justin Lemkul jalem...@vt.edu wrote:
 
 
  On 2/25/14, 1:03 PM, Albert wrote:
 
  Hello:
 
  Does anybody knows can we do Accelerated molecular dynamics (aMD,
  developed in
  Mccammon group) in Gromacs? I seldom
  see someone did this in Gromacs
 
 
  Gromacs has flooding, which I think is similar.  But as far as aMD, no.
 
  -Justin
 
  --
  ==
 
  Justin A. Lemkul, Ph.D.
  Ruth L. Kirschstein NRSA Postdoctoral Fellow
 
  Department of Pharmaceutical Sciences
  School of Pharmacy
  Health Sciences Facility II, Room 601
  University of Maryland, Baltimore
  20 Penn St.
  Baltimore, MD 21201
 
  jalem...@outerbanks.umaryland.edu | (410) 706-7441
  http://mackerell.umaryland.edu/~jalemkul
 
  ==
 
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a
  mail to gmx-users-requ...@gromacs.org.
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.




 --
 Zhifeng (Francis) Jing
 Graduate Student in Physical Chemistry
 School of Chemistry and Chemical Engineering
 Shanghai Jiao Tong University
 http://sun.sjtu.edu.cn
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Temperature coupling

[Michael Shirts]

See the following for more information:

Effects of Temperature Control Algorithms on Transport Properties and
Kinetics in Molecular Dynamics Simulations

http://pubs.acs.org/doi/abs/10.1021/ct400109a

When discussing velocity profile, I believe this is referring to any
external driving force that results in net movement of the system.

On Sat, Feb 22, 2014 at 4:21 PM, Dr. Vitaly Chaban vvcha...@gmail.com wrote:
 I do not know what the manual says, but if you are talking about
 viscosity computation through cos-like acceleration, then g_energy is
 doing all the work for you.


 Dr. Vitaly V. Chaban


 On Fri, Feb 21, 2014 at 10:38 PM, Marcelo Vanean vanea...@gmail.com wrote:
 I am trying to calculate the viscosity of water (spc model) by periodic
 perturbation method. The manual says (manual 4.6.5, pages 169 - 170): To
 obtain the correct value for the viscosity the generated velocity profile
 should not be coupled to the heat bath, also the velocity profile should be
 excluded from the kinetic energy and The heat generated by the viscous
 friction is removed by coupling to a heat bath. I don't understand.
 How is possible
 to couple the system to the thermal bath, but not the velocity profile?

 Thank you in advance.
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Temperature coupling

[Michael Shirts]

I'm pretty sure that viscous friction is not what is referred to when
it's talking about a velocity profile; it's the COM velocities.

On Sat, Feb 22, 2014 at 7:46 PM, Marcelo Vanean vanea...@gmail.com wrote:
 Vitaly V. Chaban, this is the question. How to remove the heat
 generated by viscous friction by heat bath without coupling the
 velocity profile to heat bath? Also, why to remove the velocity
 profile of the kinetic energy?

 On Sat, Feb 22, 2014 at 9:08 PM, Marcelo Vanean vanea...@gmail.com wrote:
 Michael Shirts, in this article is not discussed about the periodic
 perturbation method to calculate viscosity. Anyway, thank you.

 On Sat, Feb 22, 2014 at 7:16 PM, Marcelo Vanean vanea...@gmail.com wrote:
 Dr. Vitaly Chaban, I  know  that g_energy is doing all the work for me.
 However I want to understand the concepts.


 On Fri, Feb 21, 2014 at 6:38 PM, Marcelo Vanean vanea...@gmail.com wrote:

 I am trying to calculate the viscosity of water (spc model) by periodic
 perturbation method. The manual says (manual 4.6.5, pages 169 - 170): To
 obtain the correct value for the viscosity the generated velocity profile
 should not be coupled to the heat bath, also the velocity profile should be
 excluded from the kinetic energy and The heat generated by the viscous
 friction is removed by coupling to a heat bath. I don't understand. How is
 possible to couple the system to the thermal bath, but not the velocity
 profile?

 Thank you in advance.


 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free energy calculation problems (bug?)

[Michael Shirts]

Two thoughts that probably don't explain the problem, since they
should have similar behavior between versions:

* vdw-type = Cut-off

is almost never a good idea.  If you are doing this to be
reproducible, it's fine, but I'd avoid it in the future.

* You almost certainly want DispCorr = EnerPres, or else the density
will change with cutoff, and dispersion interactions beyond the cutoff
will be neglected (also changing the energy with cutoff).

Can you file a redmine (http://redmine.gromacs.org) with .top  and
.gro files as well?  One never knows where a problem might be coming
from (maybe only with certain top files?), so having those helps as
well.  It would be better to have the softcore version, since that is
likely to be better behaved anyway, and it has the same difference.


On Sat, Feb 8, 2014 at 12:00 AM, Олег Титов tito...@qsar.chem.msu.ru wrote:
 Dear GROMACS users!

 I've encountered a problem in thermodynamic integration with GROMACS. I'm
 interested in calculating of solvation free energy differences between
 holbenzenes and benzene. To calculate this values I use thermodynamic
 integration approach (5 point formula).

 Half a year ago I did my calculation with GROMACS 4.6.1. The results looked
 good and were consistent with literature data. Recently, I've discovered
 that versions 4.6.0 and 4.6.1 had a bug in free energy code, so I decided
 to recalculate everything with 4.6.5 version of GROMACS. New results look
 not so good and, moreover, are not consistent with literature.

 I use GAFF forcefield with TIP3P water models. To obtain GROMACS topology
 and coordinate files I convert AMBER topologies (generated with tleap) with
 amb2gmx.pl script and this topologies look fine.

  For phenylbromide molecule calculated solvation free difference should be
 near -0.1 - -0.2 kcal/mol (Experimental data is 0.59 kcal/mol).
 4.6.1 version gave me -0.1+-0.13 kcal/mol - Fine!
 I've performed the same calculation with AMBER11  - -0.27 +-0.16 kcal/mol -
 also fine.
 4.6.5 version for the same system gives me 1.1 - 1.2 kcal/mol

 I've played with some options (tried adding soft-core potentials and turned
 on/off couple-intramol option). The results are:
 no_sc no_intramol 1.27 +- 0.18
 no_sc intramol 1.20 +- 0.20
 sc no_intramol 1.35 +- 0.21
 sc intramol 1.22 +- 0.19

 So now I don't know what to do. Old buggy version produced good results
 consistent with literature and other MD packages and new fixed version
 produces different results. The dH/dl curve for two versions of GROMACS
 looks similar for lambdas 0-0.7 (mutation is PhBr - PhH ) but becomes
 different at lambda = 0.9 So something changed there... Moreover different
 elecrostatic approaches for PhBr molecule are consistent with each other
 (dddG are fine), but new version of GROMACS gives me overestimated ddG as
 shown above.

 Can anyone give me some advice what to do with it? May be there is some bug
 in new version? May be I'm doing something wrong?

 Attaching sample .mdp (for solvated and vacuum system) and .top files.

 Thanks for the help,
 Oleg

 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure coupling in NPT equilibration step

[Michael Shirts]

Note that Parrinello-Rahman contains some very minor approximations
(virial and kinetic energies off by 1/2 a step) that might affect
things by 0.3 atm. But that is almost certainly not the limiting
factor in the accuracy of the simulation.

On Mon, Dec 16, 2013 at 3:09 PM, Justin Lemkul jalem...@vt.edu wrote:
 On Mon, Dec 16, 2013 at 2:42 PM, Mahboobeh Eslami 
 mahboobeh.esl...@yahoo.com wrote:

 hi GMX user
 please help me

 i want to study protein ligand complex by gromacs4.6.3. i repaired my
 protein by MODELER and docked my ligand by AUTODOCK software.

 i use amber99sb-ildn force field.  after Energy Minimization i use
 following option forTemperature coupling in NVT equilibration step:

 ; Temperature coupling
 tcoupl= v-rescale   ; Couple temperature to external heat
 bath according to velocity re-scale method
 tc-grps   = Protein_UNK NA_SOL  ; two coupling groups - more accurate
 tau_t = 0.10.1  ; Coupling time constant, controlling
 strength of coupling
 ref_t = 300 300 ; Temperature of heat bath

  then i use following option for Pressure coupling  in NPT equilibration
 step

 ; Pressure coupling
 pcoupl   = Parrinello-Rahman   ; pressure coupling is on for
 NPT
 Pcoupltype   = Isotropic  ; uniform scaling of box
 vectors
 tau_p= 0.5  ;  time constant, in ps
 compressibility  = 4.5e-5; isothermal compressibility
 of water, bar^-1
 ref_p= 1.0 ; reference pressure,
 in bar
 refcoord_scaling = com

  i obtain temperature average 299.63 for NVT step bu i obtain pressure
 average 1.35 for NPT stepI tried Various Thermostats and Barostatsbut i
 don’t get good pressure averageClose to 1bar.


 You don't consider 1.35 bar close to 1 bar?  For pressure, I'd call that
 exceptionally good, especially if you pay attention to the fluctuations.

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

 -Justin

 --

 ==

 Justin A. Lemkul, Ph.D.
 Postdoctoral Fellow

 Department of Pharmaceutical Sciences
 School of Pharmacy
 Health Sciences Facility II, Room 601
 University of Maryland, Baltimore
 20 Penn St.
 Baltimore, MD 21201

 jalem...@outerbanks.umaryland.edu | (410) 706-7441


 ==
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure coupling constants

[Michael Shirts]

The most physical thing would be to prepare a system in a small NVT
box, then make the box bigger along one dimension, then do an NVT
simulation of that system and see what happens along that dimension
over time. Barostats are simply not designed to do nonequilibrium
right.  They magically create volume at all points in the box
simultaneously. An equilibrium simulation doesn't care where in the
box the volume appears; a nonequilibrium simulation very much does.

On Thu, Dec 12, 2013 at 4:01 PM, kpsanto santo...@gmail.com wrote:
 Hi
 Thanks a lot. But if one wants to simulate the sudden expansion (in a crude
 way ,of course) smaller coupling constant would be more appropriate. Isn't
 it?

 santo

 -
 K. P. Santo
 Post doctoral fellow
 Department of Chemistry
 University of North Carolina at Chapel Hill
 Chapel Hill, NC, USA

 --
 View this message in context: 
 http://gromacs.5086.x6.nabble.com/Pressure-coupling-constants-tp5013292p5013300.html
 Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure coupling constants

[Michael Shirts]

You're right, NVE would probably be better!  Otherwise, you would not
get expansionary cooling.

On Thu, Dec 12, 2013 at 7:08 PM, kpsanto santo...@gmail.com wrote:
 guess that right! but why do you think NVT scaling of velocities is physical
 which should also affect the non-equilibrium expansion process?

 thanks in advance
 Santo

 -
 K. P. Santo
 Post doctoral fellow
 Department of Chemistry
 University of North Carolina at Chapel Hill
 Chapel Hill, NC, USA

 --
 View this message in context: 
 http://gromacs.5086.x6.nabble.com/Pressure-coupling-constants-tp5013292p5013303.html
 Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Compilation issue with F77_FUNC functions?

[Michael Shirts]

And looks like the solution was to make sure that mpic++ was used as
the compiler instead of mpicc.



On Mon, Dec 9, 2013 at 5:22 AM, Mark Abraham mark.j.abra...@gmail.com wrote:
 On Mon, Dec 9, 2013 at 1:44 AM, Michael Shirts mrshi...@gmail.com wrote:

 So, I'm trying to compile with MPI using mpich3.  Previous
 installations worked, and installations without MPI worked. I'm
 getting errors like:

 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 warning: parameter names (without types) in function declaration


 Apparently F77_FUNC is getting #defined by someone else after config.h,
 because vec.h transitively #includes the entire world, including mpi.h.
 That could be proved by inspecting the preprocessed source, e.g. by hacking
 the command you can get from make VERBOSE=1. Assigning blame there would
 tax Solomon. Possible fixes include trying a more recent compiler, hacking
 out the #include vec.h (which I think is even superfluous), or getting
 rid of the F77_FUNC stuff there (which is probably a relic from when some
 fortran kernel needed the data).

 Mark

 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 error: function ‘F77_FUNC’ is initialized like a variable
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: braces around scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: invalid initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)

 And it keeps on like that for a long while.

 Any suggestions?  Perhaps something wrong in the way mpicc is handling
 Fortran code?
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] pressure based replica exchange

[Michael Shirts]

Pressure replica exchange is not in currently. It will be added
eventually (it's not that hard, just takes time, and hasn't been a
priority.).

Note that pressure temperature control with berendsen pressure control
will be 100% wrong.  For pressure replica exchange to work, the
barostat has to give the correct volume fluctuations, and berendsen
does not give correct volume fluctuations.

Of course, that's irrelevant since pressure replica exchange isn't in yet.

The exchange criteria is always of the form min(1,exp(-X), where X =
\beta(Delta (beta U) + Delta (beta P V)).  Where P, beta, and U can
vary with state.

What is in is temperature replica exchange with NPT, where beta
changes, but P stays the same.

On Thu, Dec 12, 2013 at 1:40 AM, Jianguo Li ljg...@yahoo.com.sg wrote:
 Thanks very much. As you suggested, making T with small increments such as 
 300.0, 300.0001, 300.0002 ... leads to complaints :-)

 I am trying to do replica exchange at different surface tensions for my 
 membrane system.
 To achieve different surface tensions for different replica, I only changed 
 Pxy, but keep Pz the same for all replica.
 pcoupl  = berendsen
 pcoupltype  = semi-isotropic
 tau_p   = 0.5 0.5
 ref_p   = 1.0 1.0  ; Pxy and Pz
 compressibility = 4.5e-5  4.5e-5

 I noticed that in the manual, the eq (3.141) shows the exchange probability. 
 Can you explain what does the P1 and P2 represent? Pxy or Pz?

 Also, the manual mentioned about using Hamiltonian and temperature replica 
 exchange simultaneously (eq 3.143). But there is no pressure term in the 
 acceptance criteria (eq 3.143). I am wondering if 
 pressure+temperature+Hamiltonian replica exchange has been implemented in 
 Gromacs. If yes, what is the general equation of acceptance criteria?

 Cheers,
 Jianguo







 On Tuesday, 10 December 2013, 17:49, Mark Abraham mark.j.abra...@gmail.com 
 wrote:






 On Tue, Dec 10, 2013 at 1:47 PM, Jianguo Li ljg...@yahoo.com.sg wrote:

 Dear All,

Is it possible to do pressure based replica exchange simulations in gromacs? 
Basically I want to do replica exchange simulations for my membrane system at 
different surface tensions. If I just set different pressures in the mpd 
file, then the mdrun will complain systems are all the same, there is 
nothing to exchange. Is there any other parameters I need to set?


 Yes, the implementation is a bit broken, and requires that one of ref-t or 
 lambda varies before it even considers the pressure. If you want T constant 
 and P as the replica control variable, you can can only do that by fooling 
 the implementation by having T vary in tiny increments.

 Mark


Thanks a lot!

Jianguo
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Compilation issue with F77_FUNC functions?

[Michael Shirts]

Apologies! I certainly didn't post enough information.  mpicc is using:

gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-3)

Which is exactly what the non-mpi versions used, which did not have
this problem.

I ran cmake with:

cmake ../gromacs -DGMX_MPI=ON -DGMX_GPU=OFF -DGMX_DOUBLE=ON
-DGMX_CPU_ACCELERATION=AVX_256
-DCMAKE_INSTALL_PREFIX=/h3/n1/shirtsgroup/gromacs_46/install
-DFFTW_INCLUDE_DIR=/h3/n1/shirtsgroup/software/fft3w/include
-DFFTW_LIBRARY='/h3/n1/shirtsgroup/software/fft3w/lib/libfftw3.a;/h3/n1/shirtsgroup/software/fft3w/lib/libfftw3.so;/usr/lib64/libm.so;/share/apps/mpich3/gnu/lib/libmpich.so;'

All the extra libraries for FFTWF appear to be necessary for some
reason, but I don't think that's it . . .

Commentinf out the line F77_FUNC int src/config.h.cmakein and
rerunning cmake did not change anything.

On Sun, Dec 8, 2013 at 10:29 PM, Roland Schulz rol...@utk.edu wrote:
 Hi,

 what compiler is used by mpicc? What does mpicc -showme and mpicc
 --version show? Does it help to uncomment the line containing F77_FUNC int
 src/config.h.cmakein?

 Roland


 On Sun, Dec 8, 2013 at 7:44 PM, Michael Shirts mrshi...@gmail.com wrote:

 So, I'm trying to compile with MPI using mpich3.  Previous
 installations worked, and installations without MPI worked. I'm
 getting errors like:

 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 warning: parameter names (without types) in function declaration
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 error: function ‘F77_FUNC’ is initialized like a variable
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: braces around scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: invalid initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)

 And it keeps on like that for a long while.

 Any suggestions?  Perhaps something wrong in the way mpicc is handling
 Fortran code?
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.







 --
 ORNL/UT Center for Molecular Biophysics cmb.ornl.gov
 865-241-1537, ORNL PO BOX 2008 MS6309
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Compilation issue with F77_FUNC functions?

[Michael Shirts]

And FWIW, it's being compiled on CentOS 6.4.

On Sun, Dec 8, 2013 at 11:16 PM, Michael Shirts mrshi...@gmail.com wrote:
 Apologies! I certainly didn't post enough information.  mpicc is using:

 gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-3)

 Which is exactly what the non-mpi versions used, which did not have
 this problem.

 I ran cmake with:

 cmake ../gromacs -DGMX_MPI=ON -DGMX_GPU=OFF -DGMX_DOUBLE=ON
 -DGMX_CPU_ACCELERATION=AVX_256
 -DCMAKE_INSTALL_PREFIX=/h3/n1/shirtsgroup/gromacs_46/install
 -DFFTW_INCLUDE_DIR=/h3/n1/shirtsgroup/software/fft3w/include
 -DFFTW_LIBRARY='/h3/n1/shirtsgroup/software/fft3w/lib/libfftw3.a;/h3/n1/shirtsgroup/software/fft3w/lib/libfftw3.so;/usr/lib64/libm.so;/share/apps/mpich3/gnu/lib/libmpich.so;'

 All the extra libraries for FFTWF appear to be necessary for some
 reason, but I don't think that's it . . .

 Commentinf out the line F77_FUNC int src/config.h.cmakein and
 rerunning cmake did not change anything.

 On Sun, Dec 8, 2013 at 10:29 PM, Roland Schulz rol...@utk.edu wrote:
 Hi,

 what compiler is used by mpicc? What does mpicc -showme and mpicc
 --version show? Does it help to uncomment the line containing F77_FUNC int
 src/config.h.cmakein?

 Roland


 On Sun, Dec 8, 2013 at 7:44 PM, Michael Shirts mrshi...@gmail.com wrote:

 So, I'm trying to compile with MPI using mpich3.  Previous
 installations worked, and installations without MPI worked. I'm
 getting errors like:

 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 warning: parameter names (without types) in function declaration
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:53:
 error: function ‘F77_FUNC’ is initialized like a variable
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: braces around scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:56:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: invalid initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 error: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: excess elements in scalar initializer
 /h3/n1/shirtsgroup/gromacs_46/gromacs/src/gmxlib/cinvsqrtdata.c:57:
 warning: (near initialization for ‘F77_FUNC’)

 And it keeps on like that for a long while.

 Any suggestions?  Perhaps something wrong in the way mpicc is handling
 Fortran code?
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.







 --
 ORNL/UT Center for Molecular Biophysics cmb.ornl.gov
 865-241-1537, ORNL PO BOX 2008 MS6309
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Regarding soft core parameter: Relative solvation free energy calculation

[Michael Shirts]

The fluctuation is often of the order of magnitude of dh/dl itself.
It's dh/dl that matters.



On Thu, Nov 28, 2013 at 1:26 PM, bipin singh bipinel...@gmail.com wrote:
 Thanks for you reply Prof. Shirts.

 I have plotted the dh/dl values (link mentioned below), but I am not sure
 whether the fluctuation is in acceptable range at lambda=0.0. Please have a
 look at the plot and let me know your thoughts.

 http://researchweb.iiit.ac.in/~bipin.singh/dhdl.png


 On Thu, Nov 28, 2013 at 8:44 PM, Michael Shirts mrshi...@gmail.com wrote:

 It's very possible that this is entirely physical. dh/dl fluctuates a
 lot.  You may want to look at graphs of dh/dl after a few 10's of ps
 to see if it looks reasonable.

 On Thu, Nov 28, 2013 at 10:04 AM, bipin singh bipinel...@gmail.com
 wrote:
  I am mentioning below the average and standard deviation of dh/dl values
 at
  lambda=0.0 and lambda=1.0, for sc-alpha=0.5
 
  Avg Std.
   lambda0.0-2.055812e+01   2.730571e+01
 
   lambda1.07.086960e+017.670135e+00
 
 
  On Thu, Nov 28, 2013 at 7:34 PM, bipin singh bipinel...@gmail.com
 wrote:
 
  Thanks for the reply.
 
  By large fluctuations, I mean the standard deviation of dh/dl values
  during the simulation at lambda 0.0 to 0.5 in comparison to lambda
 close to
  1, when using sc-alpha=0.5
 
 
  On Thu, Nov 28, 2013 at 7:27 PM, Michael Shirts mrshi...@gmail.com
 wrote:
 
  Define large fluctuations.  They might be physical large fluctuations!
 
 
 
  On Thu, Nov 28, 2013 at 1:19 AM, bipin singh bipinel...@gmail.com
  wrote:
   Hello All,
  
   I am trying to calculate relative solvation free energy for p-Cresol
 and
   p-Chlorophenol using Gromacs 4.6.3.
  
   I am using equispaced lamda values viz. 0.1, 0.05, 0.1,, 1.0
  
   During one transformation (lambda 0.0 to 1.0) I am switching off the
 vdW
   terms (Changed to DUM atoms) for the -CH3 and -Cl group in p-Cresol
 amd
   p-Chlorophenol respectively.
  
   Thus I am using a constant sc-alpha=0.5 at each lambda values (0.0 to
  1.0)
   for both the molecules.
  
   But I am getting a large fluctuation in dh/dl values only at initial
  lambda
   values i.e. from 0.0 to 0.5, and only in the case of p-Cresol (-CH3
 is
   being changed to DUM).
  
   If I change the sc-alpha=0 for initial lamda ( i.e. from 0.0 to
 0.5), I
  get
   consistent values of dh/dl (i.e. less deviation during simulation).
  
  
   Can anyone suggest the reason for this behaviour and how to select
   appropriate values of sc-alpha based on transformation. And can we
 use
   different sets of sc-alpha values for two molecules to compute their
   relative solvation free energies.
  
  
   
   I have used the following mdp settings:
  
   sc-power = 1
   sc-alpha = 0.5
   sc-r-power   = 6
   sc-coul  = no
  
  
   and
  
  
   sc-power = 1
   sc-alpha = 0.0
   sc-r-power   = 6
   sc-coul  = no
  
   ###
  
  
  
  
  
   *Thanks and Regards,Bipin Singh*
   --
   Gromacs Users mailing list
  
   * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
  
   * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  
   * For (un)subscribe requests visit
   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-usersor
  send a mail to gmx-users-requ...@gromacs.org.
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 
 
 
  --
 
 
 
  *Thanks and Regards, Bipin Singh*
 
 
 
 
  --
 
 
 
  *Thanks and Regards,Bipin Singh*
  --
  Gromacs Users mailing list
 
  * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.




 --



 *Thanks and Regards,Bipin Singh*
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org

Re: [gmx-users] Regarding soft core parameter: Relative solvation free energy calculation

[Michael Shirts]

Define large fluctuations.  They might be physical large fluctuations!



On Thu, Nov 28, 2013 at 1:19 AM, bipin singh bipinel...@gmail.com wrote:
 Hello All,

 I am trying to calculate relative solvation free energy for p-Cresol and
 p-Chlorophenol using Gromacs 4.6.3.

 I am using equispaced lamda values viz. 0.1, 0.05, 0.1,, 1.0

 During one transformation (lambda 0.0 to 1.0) I am switching off the vdW
 terms (Changed to DUM atoms) for the -CH3 and -Cl group in p-Cresol amd
 p-Chlorophenol respectively.

 Thus I am using a constant sc-alpha=0.5 at each lambda values (0.0 to 1.0)
 for both the molecules.

 But I am getting a large fluctuation in dh/dl values only at initial lambda
 values i.e. from 0.0 to 0.5, and only in the case of p-Cresol (-CH3 is
 being changed to DUM).

 If I change the sc-alpha=0 for initial lamda ( i.e. from 0.0 to 0.5), I get
 consistent values of dh/dl (i.e. less deviation during simulation).


 Can anyone suggest the reason for this behaviour and how to select
 appropriate values of sc-alpha based on transformation. And can we use
 different sets of sc-alpha values for two molecules to compute their
 relative solvation free energies.


 
 I have used the following mdp settings:

 sc-power = 1
 sc-alpha = 0.5
 sc-r-power   = 6
 sc-coul  = no


 and


 sc-power = 1
 sc-alpha = 0.0
 sc-r-power   = 6
 sc-coul  = no

 ###





 *Thanks and Regards,Bipin Singh*
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Regarding soft core parameter: Relative solvation free energy calculation

[Michael Shirts]

It's very possible that this is entirely physical. dh/dl fluctuates a
lot.  You may want to look at graphs of dh/dl after a few 10's of ps
to see if it looks reasonable.

On Thu, Nov 28, 2013 at 10:04 AM, bipin singh bipinel...@gmail.com wrote:
 I am mentioning below the average and standard deviation of dh/dl values at
 lambda=0.0 and lambda=1.0, for sc-alpha=0.5

 Avg Std.
  lambda0.0-2.055812e+01   2.730571e+01

  lambda1.07.086960e+017.670135e+00


 On Thu, Nov 28, 2013 at 7:34 PM, bipin singh bipinel...@gmail.com wrote:

 Thanks for the reply.

 By large fluctuations, I mean the standard deviation of dh/dl values
 during the simulation at lambda 0.0 to 0.5 in comparison to lambda close to
 1, when using sc-alpha=0.5


 On Thu, Nov 28, 2013 at 7:27 PM, Michael Shirts mrshi...@gmail.comwrote:

 Define large fluctuations.  They might be physical large fluctuations!



 On Thu, Nov 28, 2013 at 1:19 AM, bipin singh bipinel...@gmail.com
 wrote:
  Hello All,
 
  I am trying to calculate relative solvation free energy for p-Cresol and
  p-Chlorophenol using Gromacs 4.6.3.
 
  I am using equispaced lamda values viz. 0.1, 0.05, 0.1,, 1.0
 
  During one transformation (lambda 0.0 to 1.0) I am switching off the vdW
  terms (Changed to DUM atoms) for the -CH3 and -Cl group in p-Cresol amd
  p-Chlorophenol respectively.
 
  Thus I am using a constant sc-alpha=0.5 at each lambda values (0.0 to
 1.0)
  for both the molecules.
 
  But I am getting a large fluctuation in dh/dl values only at initial
 lambda
  values i.e. from 0.0 to 0.5, and only in the case of p-Cresol (-CH3 is
  being changed to DUM).
 
  If I change the sc-alpha=0 for initial lamda ( i.e. from 0.0 to 0.5), I
 get
  consistent values of dh/dl (i.e. less deviation during simulation).
 
 
  Can anyone suggest the reason for this behaviour and how to select
  appropriate values of sc-alpha based on transformation. And can we use
  different sets of sc-alpha values for two molecules to compute their
  relative solvation free energies.
 
 
  
  I have used the following mdp settings:
 
  sc-power = 1
  sc-alpha = 0.5
  sc-r-power   = 6
  sc-coul  = no
 
 
  and
 
 
  sc-power = 1
  sc-alpha = 0.0
  sc-r-power   = 6
  sc-coul  = no
 
  ###
 
 
 
 
 
  *Thanks and Regards,Bipin Singh*
  --
  Gromacs Users mailing list
 
  * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.




 --



 *Thanks and Regards, Bipin Singh*




 --



 *Thanks and Regards,Bipin Singh*
 --
 Gromacs Users mailing list

 * Please search the archive at 
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
 mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.