Re: [gmx-users] Fwd: Relative free energy perturbation

[Hannes Loeffler]

On Thu, 6 Jul 2017 09:17:18 +0200
Davide Bonanni <davide.bona...@unito.it> wrote:

> I am apologize if I resend the same message but I forgot to change the
> object.
> 
> Dear Hannes,
> 
> Thank you very much for your reply, really appreciated.
> 
> 
> 
> > Date: Mon, 26 Jun 2017 12:09:45 +0100
> > From: Hannes Loeffler <hannes.loeff...@stfc.ac.uk>
> > To: <gromacs.org_gmx-users@maillist.sys.kth.se>
> > Cc: gmx-us...@gromacs.org
> > Subject: Re: [gmx-users] Fwd: Relative free energy perturbation
> > Message-ID: <20170626120945.75215...@stfc.ac.uk>
> > Content-Type: text/plain; charset="US-ASCII"
> >
> > On Mon, 26 Jun 2017 12:25:19 +0200
> > Davide Bonanni <davide.bona...@unito.it> wrote:
> >  
> > > 1) Can I perform the calculation in a single step with soft core
> > > potential enabled? I mean, is it correct to transform directly the
> > > hydrogen into a chlorine instead of using 2 topologys and 2
> > > complexes, where in the first step I transform the hydrogen into
> > > dummy atom, and in the second I transform the dummy atom into
> > > chlorine.  
> >
> > Technically speaking you can perfectly do that but in practice it
> > can be much more efficient to directly and linearly transform one
> > atom type into another (single topology approach).  There is no
> > need for a softcore potential in this case.  Those would only be
> > activated for atoms that either appear or disappear i.e. atoms with
> > zero vdW parameters.  The input and topology files from FEsetup
> > should be all you need.  
> 
> 
> Can I have problems if I keep active the softcore potential whether
> is not needed?

I don't think vdW softcores would be activated as said above. In any
case you can just sc-alpha=0 to switch softcores off all together.
Even if they were on, all that would happen is that you would choose an
alternative path (so need to be consistently used in the other leg(s)
of the thermodynamic cycle).  So you should get the same overall ddG
albeit in a possible less efficient fashion because you completely
decouple on atom while recoupling its counter-part.  It just doesn't
make sense to me to do it that way.

Cheers,
Hannes.
-- 
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: Relative free energy perturbation

[Hannes Loeffler]

On Mon, 26 Jun 2017 12:25:19 +0200
Davide Bonanni  wrote:

> 1) Can I perform the calculation in a single step with soft core
> potential enabled? I mean, is it correct to transform directly the
> hydrogen into a chlorine instead of using 2 topologys and 2
> complexes, where in the first step I transform the hydrogen into
> dummy atom, and in the second I transform the dummy atom into
> chlorine.

Technically speaking you can perfectly do that but in practice it can
be much more efficient to directly and linearly transform one atom type
into another (single topology approach).  There is no need for a
softcore potential in this case.  Those would only be activated for
atoms that either appear or disappear i.e. atoms with zero vdW
parameters.  The input and topology files from FEsetup should be all
you need.


> 2) Referring to BevanLab Tutorial 6: Free Energy Calculation (
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gm
> x-tutorials/free_energy/index.html), I have to perform every step of
> molecular dynamics at every Lambda value, is that right?

Yes, you will need to do a simulation for every and each lambda step.


> 3) I run a test minimization step (mdp file attached) of my complex
> at the last "init-lambda-state", 15 in my case. Looking at the .trr
> output file I can see that the bond between the carbon and the
> hydrogen which should be trasformed is longer than a normal C-H bond,
> but the atom is still recognized as hydrogen (picture
> "http://tinypic.com/r/2wp11dh/9": purple -> init-lambda-state = 0 ;
> blue -> init-lambda-state = 15). I was wondering if this is what I am
> supposed to see, so if gromacs is considering the state B of my
> system where I have chlorine bound to carbon instead hydrogen.

What does "recognized as hydrogen" mean?  I suspect that what you are
referring to is the output of some visualisation program because you
instructed it to interpret that particular atom to be a hydrogen.

What you need to expect to see is that a C-H bond is transformed into a
C-Cl bond and accordingly the bond length increases.
-- 
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] Perturbation Thermodynamic Integration

[Hannes Loeffler]

 = 5000
> 
> ;Output Control
> nstxout = 0
> nstvout = 0
> 
> ;Cutoff Schemes
> cutoff-scheme= group
> rlist= 1.0
> vdw-type = cut-off
> rvdw = 2.0
> 
> ;Coulomb interactions
> coulombtype  = pme
> rcoulomb = 1.0
> fourierspacing   = 0.4
> 
> ;Constraints
> constraints  = none
> 
> ;Temperature coupling
> tcoupl   = v-rescale
> tc-grps  = system
> tau-t= 0.1
> ref-t= 300
> 
> ;Pressure coupling
> pcoupl = parrinello-rahman
> ref-p = 1.01325
> compressibility = 4.5e-5
> tau-p = 5
> 
> ;Free energy calculation
> free-energy  = yes
> init-lambda-state= 0
> delta-lambda = 0
> fep-lambdas  =
> calc-lambda-neighbors= 1
> vdw_lambdas  = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> coul_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> bonded_lambdas   = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> restraint_lambdas= 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> mass_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> nstdhdl  = 10
> 
> --Results--
> A(aq)->B(aq)
> Means
> ; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda |
> restraint-lambda | dH/dlambda to before | dH/dlambda to this |
> dH/dlambda to next | pV | lambda
> 585.258 29.695 138.389 -106.111 0.000 0.000 0.000 78.180 14.434 0.0
> 194.905 27.926 109.373 -113.567 0.000 -8.514 0.000 35.428 14.433 0.1
> 201.785 23.966 90.522 -123.135 0.000 -5.856 0.000 32.987 14.437 0.2
> 101.067 21.244 76.286 -130.204 0.000 6.727 0.000 20.621 14.435 0.3
> 67.782 19.301 64.038 -136.090 0.000 12.170 0.000 15.391 14.436 0.4
> 54.102 17.403 54.579 -140.360 0.000 15.207 0.000 12.567 14.438 0.5
> 45.855 14.376 46.120 -146.966 0.000 17.949 0.000 10.041 14.437 0.6
> 39.769 11.953 39.144 -151.300 0.000 20.040 0.000 8.168 14.438 0.7
> 35.384 11.090 33.911 -151.862 0.000 21.252 0.000 7.172 14.439 0.8
> 31.411 9.077 28.113 -154.575 0.000 22.809 0.000 5.829 14.440 0.9
> 28.398 9.974 24.842 -151.347 0.000 23.126 0.000 0.000 14.439 1.0
> Variances
> ; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda |
> restraint-lambda | dH/dlambda to before | dH/dlambda to this |
> dH/dlambda to next | pV | lambda
> 17764.592 28.101 889.212 7624.759 0.000 0.000 0.000 265.656 0.002 0.0
> 1580.485 25.381 492.489 7085.318 0.000 93.157 0.000 95.406 0.002 0.1
> 4269.905 26.694 347.400 7025.386 0.000 118.370 0.000 120.552 0.002 0.2
> 1329.740 26.648 259.096 7083.151 0.000 85.371 0.000 87.603 0.002 0.3
> 249.921 24.406 195.233 7149.494 0.000 77.239 0.000 79.441 0.002 0.4
> 121.269 17.111 156.460 7039.120 0.000 74.030 0.000 76.206 0.002 0.5
> 87.000 15.901 128.267 7054.185 0.000 73.556 0.000 75.723 0.002 0.6
> 65.818 14.577 104.742 7096.497 0.000 73.225 0.000 75.402 0.002 0.7
> 55.143 14.713 91.124 7081.810 0.000 72.768 0.000 74.919 0.002 0.8
> 41.118 13.472 76.214 7198.105 0.000 73.496 0.000 75.675 0.002 0.9
> 33.836 15.572 65.528 7162.912 0.000 72.932 0.000 0.000 0.002 1.0
> 
> A(gas)->B(gas)
> Means
> ; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda |
> restraint-lambda | dH/dlambda to before | dH/dlambda to this |
> dH/dlambda to next | pV | lambda
> 1110.166 32.509 114.451 -115.748 0.000 0.000 0.000 127.586 4.377 0.0
> 398.592 26.708 101.772 -133.086 0.000 -26.055 0.000 52.957 3.912 0.1
> 237.629 27.141 84.628 -137.244 0.000 -7.766 0.000 34.879 5.057 0.2
> 84.897 20.358 76.144 -133.012 0.000 8.728 0.000 18.621 14.436 0.3
> 131.475 20.205 64.964 -158.091 0.000 7.811 0.000 19.736 3.927 0.4
> 107.720 20.655 57.763 -156.171 0.000 10.776 0.000 16.984 4.815 0.5
> 91.054 25.769 51.617 -146.044 0.000 11.632 0.000 16.326 4.020 0.6
> 78.003 21.697 45.907 -154.841 0.000 14.905 0.000 13.273 4.114 0.7
> 69.346 15.392 41.194 -167.694 0.000 18.277 0.000 10.139 6.244 0.8
> 62.263 17.129 37.528 -164.429 0.000 18.954 0.000 9.667 6.133 0.9
> 56.582 18.052 34.094 -162.977 0.000 19.727 0.000 0.000 4.456 1.0
> Variances
> ; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda |
> restraint-lambda | dH/dlambda to before | dH/dlambda to this |
> dH/dlambda to next | pV | lambda
> 40348.430 3.943 514.703 7085.063 0.000 0.000 0.000 481.126 0.030 0.0
> 6978.993 1.390 352.654 7972.572 0.000 151.629 0.000 154.128 0.000 0.1
> 2160.691 4.261 219.874 7020.561 0.000 94.026 0.000 96.219 0.263 0.2
> 301.074 32.826 249.228 7137.750 0.000 78.224 0.000 80.446 0.002 0.3
> 694.265 1.765 123.119 7225.511 0.000 80.863 0.000 83.117 0.000 0.4
> 452.408 7.903 97.154 7151.103 0.000 78.148 0.000 80.364 0.099 0.5
> 347.707 1.1

Re: [gmx-users] Perturbation Thermodynamic Integration

[Hannes Loeffler]

On Tue, 16 May 2017 15:13:10 -0400
Dan Gil  wrote:


> If you do this via decoupling ("absolute" transformation) you do that
> > once for molecule A and once for molecule B.  
> 
> 
> I believe you are referring to the BAR method? I am trying to see if
> I get the same results as another group. They used thermodynamic
> integration from A to B, so I decided to spend some more time with
> this. I will try later if I get consistent results for both methods.

No, I am not referring to a specific analysis methods.  I commented on
an alternative pathway.


> transforming the masses may interact badly
> > with bond constraints that are applied to alchemically transformed
> > bonds  
> 
> 
> Thank you for this information! Do you know if there are there
> publications that refer to this issue?

Not that I know of.
-- 
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] Perturbation Thermodynamic Integration

[Hannes Loeffler]

On Tue, 16 May 2017 10:28:08 -0400
Dan Gil  wrote:

> Thank you for the advice on the cut-off schemes and PME methods.
> 
> What is the physical meaning of a non-interacting final state
> > that has different masses from the initial state?  
> 
> 
> These free energy options was just from me trying to figure out why
> mass has any contributions at all. I am going from molecule A to B
> described in the topology. In the actual simulations, I am planning
> on using this: vdw_lambdas  = 0 0.1 0.2 0.3 0.4 0.5 0.6
> 0.7 0.8 0.9 1 coul_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6
> 0.7 0.8 0.9 1 bonded_lambdas   = 0 0.1 0.2 0.3 0.4 0.5 0.6
> 0.7 0.8 0.9 1 restraint_lambdas= 0 0.1 0.2 0.3 0.4 0.5 0.6
> 0.7 0.8 0.9 1 mass_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6
> 0.7 0.8 0.9 1
> 
> To get the solvation free energy difference between molecule A and B.

If you do this via decoupling ("absolute" transformation) you do that
once for molecule A and once for molecule B.  I don't see why you would
want to change masses in this case.

If you want to do this via a relative transformation you would not use
couple-moltype at all but fill in the B columns in your topology file
for all force field parameters that change.  In this case you _can_
change the masses but you would have to do the same transformation in
vacuum. Then the mass contributions should cancel perfectly (at least in
the limit of infinite sampling, I guess).  But Michael Shirts has
commented a while back that transforming the masses may interact badly
with bond constraints that are applied to alchemically transformed
bonds (we have seen problems with this too). So it is just simplest to
not use mass-lambdas precisely because of the aforementioned argument.


Cheers,
Hannes.
-- 
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] Perturbation Thermodynamic Integration

[Hannes Loeffler]

I have not really followed the previous email exchange but from this
mdp file I wonder what you are trying to achieve.  You seem to want to
decouple all atoms of your HEPT molecule (couple-moltype,
couple-intramol) from its environment but then you also change the
masses.  What is the physical meaning of a non-interacting final state
that has different masses from the initial state?  I believe the mass
contributions are supposed to cancel in a closed thermodynamic cycle
but what is the cycle you are simulating?


On Tue, 16 May 2017 09:30:08 -0400
Dan Gil  wrote:

> Sorry, here is the mdp file:
> 
> ;Integration Method and Parameters
> integrator   = sd
> nsteps   = 10
> dt = 0.002
> nstenergy= 1000
> nstlog   = 5000
> 
> ;Output Control
> nstxout = 0
> nstvout = 0
> 
> ;Cutoff Schemes
> cutoff-scheme= group
> rlist= 1.0
> vdw-type = cut-off
> rvdw = 2.0
> 
> ;Coulomb interactions
> coulombtype  = pme
> rcoulomb = 1.0
> fourierspacing   = 0.4
> 
> ;Constraints
> constraints  = all-bonds
> 
> ;Temperature coupling
> tcoupl   = v-rescale
> tc-grps  = system
> tau-t= 0.1
> ref-t= 300
> 
> ;Pressure coupling
> pcoupl = parrinello-rahman
> ref-p = 1.01325
> compressibility = 4.5e-5
> tau-p = 5
> 
> ;Free energy calculation
> free-energy  = yes
> init-lambda-state= 8
> delta-lambda = 0
> fep-lambdas  =
> calc-lambda-neighbors= 1
> vdw_lambdas  = 0 0   0   0   0   0   0   0   0   0   0
> coul_lambdas = 0 0   0   0   0   0   0   0   0   0   0
> bonded_lambdas   = 0 0   0   0   0   0   0   0   0   0   0
> restraint_lambdas= 0 0   0   0   0   0   0   0   0   0   0
> mass_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
> couple-moltype   = HEPT
> couple-lambda0   = vdwq
> couple-lambda1   = none
> couple-intramol  = no
> nstdhdl  = 10
> 
> 
> On Tue, May 16, 2017 at 1:02 AM, Mark Abraham
>  wrote:
> 
> > Hi,
> >
> > What use are you making of constraints? Justin suggested sharing a
> > full mdp file, which I think may help. We discovered last year that
> > you can get equipartition failure for (IIRC) all-bonds constraints
> > for moieties like -CH2Cl, and latest grompp now detects this.
> >
> > Mark
> >
> > On Tue, 16 May 2017 01:16 Dan Gil  wrote:
> >  
> > > Hello,
> > >
> > > The last thread was getting too big, and the conversation evolved
> > > to a topic different from my original question, so I decided to
> > > start a new thread.
> > >
> > > We were discussing thermodynamic integration, and why the
> > > mass_lambdas would have any contribution to the derivative of the
> > > Hamiltonian.
> > >
> > > I found a source (link below) which derives the Gibbs free energy
> > > change  
> > as  
> > > a function of lambda. I learned that the mass contribution is
> > > often  
> > assumed  
> > > to be small and negligible, given that the mass difference
> > > between the  
> > two  
> > > lambda states are small.
> > > http://www.tandfonline.com/doi/abs/10.1080/00268970600893060
> > >
> > > I think that the mass of the two lambda states that equation (14)
> > > is referring to is the total mass (mass of solvent plus solute).
> > > My system  
> > is  
> > > 1 solute (~40 atoms) infinitely diluted in solvent (23500). I
> > > wonder if I am getting nonzero mass contributions (in my dhdl.xvg
> > > output) because of finite-size effects? Would completely
> > > neglecting the mass contributions  
> > be  
> > > acceptable? Does doing this technically change the system to one
> > > that is  
> > 1  
> > > solute and an infinite number of solvent molecules where the mass
> > > contributions limit is zero?
> > >
> > > Best Regards,
> > >
> > > Dan
> > > --
> > > 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 

Re: [gmx-users] Increasing sigma

[Hannes Loeffler]

On Fri, 5 May 2017 15:09:09 +0200
Pallavi Banerjee  wrote:

> I intend to perform thermodynamic integration of my system, but not to
> calculate the free energy. I want to gradually increase the radius of
> the carbon atoms of the alkyl chains of my lipid molecules. This is
> because I could see a hole forming in my lipid bilayer along the
> edges of the simulation box and the water molecules pass through it.
> 
> Once the system with the bigger radius equilibrates, I would then
> gradually beat down the radius to the original value.

I don't think that would solve your problem.  You would be pretending
that your lipids (and thus the box size) are larger then they are but
what the bilayer really wants to do is to contract.  And water is very
quick in filling the voids.  So if you decrease the radii again you may
just be back where you have started because your lipids still need to
attain the desired density.
-- 
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] Alchemical free energy

[Hannes Loeffler]

On Tue, 4 Apr 2017 18:57:50 +0200
Alex  wrote:

> In relative binding free energy calculation via alchemical free
> energy, I noticed that either the word "Decouple" or "Annihilate" is
> used. What is the difference between them?

Actually, this is typically used in the context of absolute free
energies. The meaning depends on the author, see e.g.
http://www.alchemistry.org/wiki/Decoupling_and_annihilation


> When should we use which?

See above.  In any case, make sure that you explain what your intended
meaning is.
-- 
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] query about structural consistency

[Hannes Loeffler]

On Tue, 21 Mar 2017 13:09:46 +0530
abhisek Mondal  wrote:

> Hi,
> 
> I just have a basic query.
> 
> I'm working with a protein-ligand system with a goal of performing
> Umbrella sampling in gromacs. So first, after I build the complex, I'm
> going for NVT and followed by NPT equilibration of 20ns each. After
> equilibration, is it mandatory that my equilibrated protein-ligand
> complex matches exactly with the crystal structure used for
> equlibration (I mean before and after) ?
> 
> Please be comprehensive. I just need to understand something.

It's easy to ask for a "comprehensive" answer on a mailing list but
what you need to realise is that you are really asking about is the
basics of molecular simulation.  There is no quick answer to this and
whole text books have been written on the topic.

The way "equilibration" is typically used is to really mean to reach a
nearby local minimum and expect to see a convergent behaviour with
respect to a chosen target property.  This paper,
dx.doi.org/10.1021/acs.jctc.5b00784 may be an interesting starting
point.

You may also want to abandon words like "exact" or "correct" in this
context.


Cheers,
Hannes.
-- 
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] Error in free energy calculation

[Hannes Loeffler]

On Thu, 16 Mar 2017 10:09:58 +
Vytautas Rakeviius  wrote:

> You should try to add PLUMED (plumed.org) to GROMACS for such
> things.Bonomi, M., Branduardi, D., Bussi, G., Camilloni, C., Provasi,
> D., Raiteri, P., ... & Parrinello, M. (2009). PLUMED: A portable
> plugin for free-energy calculations with molecular dynamics. Computer
> Physics Communications, 180(10), 1961-1972. 

I'm not sure how this is helping the original request.  Does Plumed
handle the method being used here?  If not, what would the alternative
be and what would be the advantage of that?


> On Thursday, March 16, 2017 11:50 AM, Tasneem Kausar
>  wrote: 
> 
>  Dear All
> 
> I am trying to perform free energy calculations of a protein drug
> system. I am using Gromacs-5.1.4 for this purpose. decoulpe molecule
> type i.e. drug molecule is defined in separate drug.itp file. From
> the discussions of mailing list I have found that this is not a
> problem. To obtain the equilibrium values for protein-ligand
> restraint 15 ns simulation was performed. Here is the last section of
> topology file. ; Include Position restraint file
> #ifdef POSRES
> #include "posre.itp"
> #endif
> #include "drug.itp"
> [ intermolecular_interactions ]
> [ bonds ]
> ;    i    j  type    r0A    r1A    r2A    fcA    r0B    r1B    r2B
> fcB
>   1444  3118    10    0.591  0.591  10.0  0.0    0.591  0.591  10.0
> 4184.000
> 
> [ angle_restraints ]
> ;  ai    aj    ak    al  type    thA      fcA    multA  thB      fcB
> multB
>   1431  1444  3118  1444    1    52.35    0.0    1    52.35
> 41.840    1 1444  3118  3120  3118    1  118.14    0.0    1  118.14
>   41.840    1
> 
> [ dihedral_restraints ]
> ;  ai    aj    ak    al  type    phiA    dphiA  fcA    phiB      dphiB
> fcB
>   1444  1431  1444  3118    1    -99.41  0.0    0.0    -99.41    0.0
> 41.840
>   1431  1444  3118  3120    1    17.49    0.0    0.0    17.49    0.0
> 41.840
>   1444  3118  3120  3119    1    -6.49  0.0    0.0    -6.49    0.0
> 41.840
> ; Include water topology
> #include "gromos54a7.ff/spc.itp"
> 
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct      fcx        fcy        fcz
>   1    1      1000      1000      1000
> #endif
> 
> ; Include topology for ions
> #include "gromos54a7.ff/ions.itp"
> 
> [ system ]
> ; Name
> DUMM in water t= 1.0
> [ molecules ]
> ; Compound        #mols
> Protein_chain_A    1
> DRG                1
> SOL            17685
> 
> Is there any default position for [ intermolecular_interaction ] in
> topology file.
> Since I am using a drug.itp file for drug molecule. Numbering in itp
> file starts from 1. Is atom number in itp file needed to be modified?
> 
> 
> grompp gives many warning  and 1 error
> Back Off! I just backed up mdout.mdp to ./#mdout.mdp.9#
> 
> WARNING 1 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_geometry' in parameter file
> 
> 
> 
> WARNING 2 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_dim' in parameter file
> 
> 
> 
> WARNING 3 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_start' in parameter file
> 
> 
> 
> WARNING 4 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_init1' in parameter file
> 
> 
> 
> WARNING 5 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_ngroups' in parameter file
> 
> 
> 
> WARNING 6 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_group0' in parameter file
> 
> 
> 
> WARNING 7 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_group1' in parameter file
> 
> 
> 
> WARNING 8 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_k1' in parameter file
> 
> 
> 
> WARNING 9 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_kB1' in parameter file
> 
> 
> 
> NOTE 1 [file em_steep_0.mdp]:
>   With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20.
> Note that with the Verlet scheme, nstlist has no effect on the
> accuracy of your simulation.
> 
> Setting the LD random seed to 1259182625
> Generated 226 of the 1711 non-bonded parameter combinations
> 
> ---
> Program gmx grompp, VERSION 5.1.4
> Source code file:
> /home/shahid/Software/gromacs-5.1.4/src/gromacs/gmxpreprocess/topio.c,
> line: 755
> 
> Fatal error:
> Syntax error - File complex.top, line 19956
> Last line read:
> '[ intermolecular_interactions ]'
> Invalid order for directive intermolecular_interactions
> For more information and tips for troubleshooting, please check the
> GROMACS website at http://www.gromacs.org/Documentation/Errors
> ---
> 
> Here are mdp parameters
> 
> ; Options for the decoupling
> sc-alpha                = 0.5
> sc-coul                  = no
> sc-power                = 1
> sc-sigma                = 0.3
> couple-moltype          = DRG
> couple-lambda0          = none
> couple-lambda1          = vdw-q
> 

Re: [gmx-users] Error in free energy calculation

[Hannes Loeffler]

I can't help you much because I have not really paid attention how
restraints in alchemical free energy simulations are handled in newer
versions of Gromacs.  I understand that [intermolecular_interactions]
is a new section that allows using restraints based on global indices.

The error message suggest that you have put this section at the wrong
place in the topology file.  I guess the manual will have information
on where grompp expects this to be.

The warnings come from the fact to you have set pull=no but other pull
control parameters are found (and ignored as not needed).

I guess from your input that you are trying to use something like
Stefan Boresch' VBA method.  The bond restraint seems excessively
high.  I also note that you attempt to recouple the DRG molecule in one
step (couple-lambda1 = vdw-q).  This will probably mean that you should
also use a softcore potential for the Coulomb interactions.
Alternatively, you may want to recouple electrostatic and vdW
interactions separately.  There is literature that suggests the latter
approach to be more efficient.


On Thu, 16 Mar 2017 15:19:48 +0530
Tasneem Kausar  wrote:

> Dear All
> 
> I am trying to perform free energy calculations of a protein drug
> system. I am using Gromacs-5.1.4 for this purpose. decoulpe molecule
> type i.e. drug molecule is defined in separate drug.itp file. From
> the discussions of mailing list I have found that this is not a
> problem. To obtain the equilibrium values for protein-ligand
> restraint 15 ns simulation was performed. Here is the last section of
> topology file. ; Include Position restraint file
> #ifdef POSRES
> #include "posre.itp"
> #endif
> #include "drug.itp"
> [ intermolecular_interactions ]
> [ bonds ]
> ;i j  type r0A r1A r2AfcAr0B r1B
> r2B fcB
>   1444  311810 0.591   0.591   10.0   0.00.591   0.591
> 10.0 4184.000
> 
> [ angle_restraints ]
> ;   aiajakal  typethA  fcAmultA  thB  fcB
> multB
>   1431  1444  3118  1444 152.350.0152.35
> 41.8401 1444  3118  3120  3118 1   118.140.01
> 118.1441.8401
> 
> [ dihedral_restraints ]
> ;   aiajakal  typephiA dphiA  fcAphiB
> dphiB fcB
>   1444  1431  1444  3118 1-99.41   0.00.0-99.410.0
> 41.840
>   1431  1444  3118  3120 1 17.490.00.017.490.0
> 41.840
>   1444  3118  3120  3119 1 -6.49   0.00.0 -6.490.0
> 41.840
> ; Include water topology
> #include "gromos54a7.ff/spc.itp"
> 
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct   fcxfcyfcz
>11   1000   1000   1000
> #endif
> 
> ; Include topology for ions
> #include "gromos54a7.ff/ions.itp"
> 
> [ system ]
> ; Name
> DUMM in water t= 1.0
> [ molecules ]
> ; Compound#mols
> Protein_chain_A 1
> DRG 1
> SOL 17685
> 
> Is there any default position for [ intermolecular_interaction ] in
> topology file.
> Since I am using a drug.itp file for drug molecule. Numbering in itp
> file starts from 1. Is atom number in itp file needed to be modified?
> 
> 
> grompp gives many warning  and 1 error
> Back Off! I just backed up mdout.mdp to ./#mdout.mdp.9#
> 
> WARNING 1 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_geometry' in parameter file
> 
> 
> 
> WARNING 2 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_dim' in parameter file
> 
> 
> 
> WARNING 3 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_start' in parameter file
> 
> 
> 
> WARNING 4 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_init1' in parameter file
> 
> 
> 
> WARNING 5 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_ngroups' in parameter file
> 
> 
> 
> WARNING 6 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_group0' in parameter file
> 
> 
> 
> WARNING 7 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_group1' in parameter file
> 
> 
> 
> WARNING 8 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_k1' in parameter file
> 
> 
> 
> WARNING 9 [file em_steep_0.mdp, line 78]:
>   Unknown left-hand 'pull_kB1' in parameter file
> 
> 
> 
> NOTE 1 [file em_steep_0.mdp]:
>   With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20.
> Note that with the Verlet scheme, nstlist has no effect on the
> accuracy of your simulation.
> 
> Setting the LD random seed to 1259182625
> Generated 226 of the 1711 non-bonded parameter combinations
> 
> ---
> Program gmx grompp, VERSION 5.1.4
> Source code file:
> /home/shahid/Software/gromacs-5.1.4/src/gromacs/gmxpreprocess/topio.c,
> line: 755
> 
> Fatal error:
> Syntax error - File complex.top, line 19956
> Last line read:
> '[ intermolecular_interactions ]'
> Invalid order for directive 

Re: [gmx-users] Adding FE2+ parameters

[Hannes Loeffler]

On Tue, 10 Jan 2017 12:19:15 -0500
Justin Lemkul  wrote:

> On 1/10/17 12:09 PM, CROUZY Serge 119222 wrote:
> > Hello Hannes
> >
> > I'm perfectly aware how you need to be careful in using metal
> > parameters - checking for which solvent and which coordination they
> > have been created for. In my case structural metal coordinated to
> > protein amino acids... I did try to see in Li/Merz parameters which
> > line could resemble the parameters in Gromacs for Zn2+ but with no
> > success - This brings me back to the original question, where do Zn
> > parameters in Gromacs 54a7 come from ?  Knowing this I could try
> > deduce parameters for Fe2+ 
> 
> This came up recently, and the best answer that I found was: these
> parameters have simply always existed in GROMOS.  The only reference
> I could find was the GROMOS manual itself, which requires payment of
> a license fee to even read.  If that's not correct, someone please
> speak up!

I can only correct you about the availability of GROMOS.  Since the mid
last year the source code and the manual are available from their web
page.  There only seems the need for a license when one wants to use it
in connection with CPMD.


Cheers,
Hannes.
-- 
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] Adding FE2+ or FE3+ into the system

[Hannes Loeffler]

What are you planning to do with those parameters?

You could have a look into the Li/Merz parameters (and papers!)
available with the AmberTools (may not have been converted yet to
Gromacs formats and the 12-6-4 sets would need support in the code).
Generally, you should be wary when using simple ion force fields and
check carefully how they have been parameterised and what for.
Distinguishing different valencies of the same transition metal atom
with only vdW/Coulomb parameters will most likely not capture their
complex chemistry.


On Tue, 10 Jan 2017 16:24:47 +
CROUZY Serge 119222  wrote:

> I'm interested as well in knowing how to get decent parameters for
> Fe2+
> 
> From gromacs-5.1.2/share/top/gromos54a7.ff/ffnonbonded.itp
> We have Zn2+  (What is a reference for these parameters ?)
> 
> name   bondtype   mass   charge   ptype CA
> ;name  at.num   mass  charge  ptype   c6   c12
> ZN2+   30  0.000  0.000 A  0.0004182025  9.4400656e-09
> 
> From which I deduced
> if LJ = v(rij)=C(12)/r12 - C(6)/r6 = A/r^12 -C/r^6
> to be compared to the standard V=4*eps*[(sig/r)^12 - (sig/r)^6]
> A =4*eps*sig^12
> C=4*eps*sig^6
> sig = (A/C)^1/6
> eps=C^2 / (4 A)
> For Zn sig=0.168112  eps=4.631677 kJ/mol
> 
> But no Fe2+
> 
> Regards
> 
> Serge Crouzy PhD HDR
> Groupe de Modélisation et Chimie Théorique
> Laboratoire de Chimie et Biologie des Métaux 
> Institut de Biosciences et Biotechnologies de Grenoble
> CEA Grenoble  UMR  CEA/CNRS/UJF 5249
> 17, rue des martyrs
> 38054 Grenoble Cedex 9
> Bat. K  pièce 110   
> Tel (33)0438782963
> Fax (33)0438785487
> http://big.cea.fr/drf/big/CBM/GMCT
> 
> 
> 
> -Message d'origine-
> De : gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> [mailto:gromacs.org_gmx-users-boun...@maillist.sys.kth.se] De la part
> de liming_52 Envoyé : mardi 10 janvier 2017 16:54 À :
> gmx-us...@gromacs.org Objet : [gmx-users] Adding FE2+ or FE3+ into
> the system
> 
> Dear gromacs users,
> 
> I want to compare the protein in the buffers with FE2+ and FE3+,
> respectively. 
> 
> How can I add FE2+ or FE3+ into the system? What is the command?
> 
> Thanks in advance.
> 
> 
> 
> 
> 
> --
> 
> With my best wishes,
> Ming Li, PhD
> Chinese Academy of Agricultural Sciences, Beijing, China

-- 
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] forward and backward states for the binding free energy energy

[Hannes Loeffler]

On Mon, 2 Jan 2017 14:17:56 +0300
Qasim Pars  wrote:

> 1-) The mean of the forward state is ~15 kcal/mol. That is too big,
> right? Maybe the ligand doesn't get decoupled in the forward state?

It will be hard to guess whether this is meaningful or not...  You
switch with 100 ps and carry out 50 replicates.  Is that maybe to fast
or will you need more repetitions?  Experimentation (with monitoring of
the results an convergence), comparison to equilibrium runs, and
literature search should help.

There are many reasons why your current results may not match
expectations. You don't quite say where the "literature" results come
from.  You also don't say if you have closed the thermodynamical cycle
properly or if you take any measurements regarding the vanishing end
state.


> 2-) Is the free energy control stuff used in the em.mdp, nvt.mdp,
> npt.mdp and md.mdp correct for the forward and backward/reversible
> state simulations? I must use the same free energy control parameters
> in both states except init-lambda?

What you need to do is to prepare a number (equal to the number of
trajectories you plan to run in the fast-switching runs) of equilibrated
starting points for each end state. You have various ways to do that
with Gromacs. If you want to use the same control parameters as in the
free energy runs, as you plan to do it, you can use init-lambda=0.0 with
a positive delta-lambda increment and init-lambda=1.0 with a negative
increment.

 
> 3-) Is the free energy control stuff used for 100 ps simulations of
> the forward and backward/reversible state simulations correct? The
> only difference between both states must be delta-lambda?

Plus init-lambda as in answer to 2.


> 5-) I added the state B parameters to the
> [ intermolecular_interactions ] parts. That should be enough to
> define the state B of a molecule in GROMACS2016?

I think the state B topology will be created automatically with
couple-moltype.  Why do you feel like modifying this by hand?
-- 
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] Continuing a FEP alchemical

[Hannes Loeffler]

On Wed, 28 Dec 2016 17:24:08 +0100
Alex  wrote:

> Thank for your response.
> 
> You mean for the TI analysis, now problem if one .xvg file (e.g.
> case.3.xvg) has just 15 columns while another .xvg file (e.g.
> case.18.xvg) has 20 columns? and still the "alchemical analysis
> package" could be used to calculate the free energy change via TI
> method?

I suppose you mean David Mobley's tool.  Yes, it can do that.  You can
easily check that yourself though.
-- 
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] Continuing a FEP alchemical

[Hannes Loeffler]

On Wed, 28 Dec 2016 16:40:49 +0100
Alex  wrote:

> Hello Gromacs user,
> 
> I have a free energy simulation converged results harvested  by
> alchemcial analysis in 15 lambada windows, now, I would like to
> increase the number of lambada windows to 20. Would you please let me
> know how I can continue my old simulation with 15 windows to 20
> windows?

Just add the additional lambda points to your lambda vector and run
the simulations for those.  For TI analysis this is all you need to do.
For FEP style analysis you would need to make sure that you also have
the energy projections for the neighbouring lambda(s).  You should be
able to use the -rerun feature to compute those (but think about
nstxout != nstdhdl).
-- 
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 issue

[Hannes Loeffler]

On Fri, 11 Nov 2016 09:04:35 +0100
gozde ergin  wrote:

> Dear all,
> 
> I follow James Barnett tutorial of “Methane Free Energy of Solvation”
> however I use betaine molecule other than methane. In order to
> analyse the data I use alchemical_analysis.py code however I get an
> error :

You should use the issue tracker of that software
https://github.com/MobleyLab/alchemical-analysis/issues

It's not a Gromacs problem.


> pymbar.utils.ParameterError: Warning: Should have \sum_k N_k W_nk =
> 1.  Actual row sum for sample 0 was 1.307134. 663638 other rows have
> similar problems
> 
> Do you have any idea about the error or anyone here face with this
> error? Thanks in advance.
> Bests 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.

Re: [gmx-users] Free energy - FEP

[Hannes Loeffler]

On Mon, 17 Oct 2016 18:06:01 +0200
Alex  wrote:

> What I have already done in Free energy simulation was to use the
> output configuration of the last lambada windows as input
> configuration (lambada_n.gro) for the new lambada windows and so on,
> which this make simulation so time demanding as one lambada should
> start after the other one. However, I was wondering if it is
> possible(Meaningful) to simulate each windows associated with each
> lambda in FEP, completely independent to each other, so that all the
> windows could start simultaneously and then much more time would be
> saved?

Principally, each window can be considered independent of all others.
The practical problem will be that you can make sure that each window is
equilibrated adequately.  When you know both end states you may be
able to make an educated guess if the transformation will go smoothly
and you can run from the same starting state. In case you do not know
the end state, you may have to consider that that state is very
different from the initial state so "dragging" the transformation along
lambda does make some sense (but which may also fail if the end state is
totally different). In any case, try to avoid "growing" atoms/molecules
if you can or otherwise monitor carefully what is happening.  A good
test is always to run in both directions and see if the result is within
statistical bounds.
-- 
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 change and volume of the box

[Hannes Loeffler]

Have you had a look at
http://www.alchemistry.org/wiki/GROMACS_4.6_example:_n-phenylglycinonitrile_binding_to_T4_lysozyme

to see if that technique would be applicable to you and look through
the references given?  When you do absolute transformations you will
need to think about standard state (and standard volume).


On Thu, 13 Oct 2016 15:33:47 +0200
Alex  wrote:

> Dear gromacs user,
> 
> As you know the binding free energy is the difference between the free
> energy change inbounded state and unbounded state.
> 
> \Delta (\Delta G)_(binding) = \Delta G_(bounded)  - \Delta
> G_(unbounded).
> 
> Then my question is that if it is necessary to simulate the free
> energy change of bonded and unbounded states in the identical volume
> size?
> 
> In another meaning if the free energy change (either bounded or
> unbounded) is a volume depended quantity?
> 
> I am using FEP method to calculate the binding free energy of amino
> acid to  a solid surface.
> 
> Best 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.

Re: [gmx-users] Gromacs Relative Free Energy Calculation Issue

[Hannes Loeffler]

On Wed, 5 Oct 2016 09:28:16 +
Guanglin Kuang  wrote:

> Dear Hannes,
> 
> Part I:
> 
> I have tried to use FESetup to generate the topology/coordinate files
> for the mutated ligands, but I met some problems.

Also, I do not recognize anything in your archive file that has been
produced by FESetup.
-- 
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 problem with too high dG

[Hannes Loeffler]

On Tue, 4 Oct 2016 15:27:25 +0300
Vlad  wrote:

> I tried to do the decoupling and strangely the result is not the same
> Below is the end of the log
> The total dG is significantly less then previous and the last dG
> value is particularly small I think total dG will get less with more
> precise lambda sampling (by 0.05 instead of 0.1) The thing which
> still is strange is the value of entropy for the last lambda (in the
> previous case it was so high as well for the first lambda)
> 
> So you suggest to use restraints
> Througout the simulation or also gradually switched off?

I suggest that you have a careful look at that web page and related
literature to understand how this type of simulation is supposed to
work. As per my previous explanation I do not see anything particularly
strange at this point.  Making a molecule disappear is just so much
simpler than making it appear.  In the latter case you may get away
with (much) longer equilibration (did you actually do prior
minimisation?).  You may also disrupt your system quite severely esp. in
the case of a protein environment.  Maybe a proper choice of positional
restraints may help.  But I am nor sure if going this route is really
worthwhile except for better understanding the method.  You will
however not get around a careful analysis of your simulations.


Cheers,
Hannes.
-- 
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 problem with too high dG

[Hannes Loeffler]

On Tue, 4 Oct 2016 14:01:52 +0300
Vlad  wrote:
> No reason. I just thought it would give the same result.  

In theory both directions should be symmetrical and give the same
result.


> In the trajectory the ligand diffuse to the end of the cell like a
> ghost through the water molecules. And could the reason lie in the
> fact that in the initial structure during all preparations the ligand
> is "fully existing" and in the first step of simulation suddenly is
> turned off?  


I would suggest to have a look through
http://www.alchemistry.org/wiki/Example:_Absolute_Binding_Affinity and
other material on that web page.  In particular reference 1 on that
page should be of interest to you in how to use restraints within a
general scheme.

In a homogeneous environment like water it doesn't matter that your
"ghost" molecule drifts because on average the interactions with its
environment will be the same regardless of location.  In an
inhomogeneous environment, however, the nearly non-interacting ligand
will have different probabilities of where it is located.  In
principle, you would have to sample _all_ possible locations which in
practice is hard or impossible.  Hence the restraint scheme as
suggested above.

Your ligand should exists as a bonded network throughout the
simulation.  What you are doing however is "decoupling" which means
that the non-bonded interactions between the ligand and its environment
and vice versa a gradually switched off.


Cheers,
Hannes.
-- 
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 Relative Free Energy Calculation Issue

[Hannes Loeffler]

On Mon, 3 Oct 2016 19:57:07 +
Guanglin Kuang  wrote:

> Dear Gromacs users,
> 
> Has any of you managed to use Gromacs to do relative free energy
> calculations?  I have some technical questions that would need your
> suggestions.
> 
> I am trying to reproduce the Amber free energy calculation using
> PMEMD (http://ambermd.org/tutorials/advanced/tutorial9/#overview).

Keep in mind that this is just a tutorial for demonstration purposes
only.  It is not intended to produce absolutely correct results.


> First, I generated the topology and coordinate files of the ligands
> (benzene and phenol) using alchemistry_setup.py, developed by Dr.
> David Mobley (https://github.com/MobleyLab/alchemical-setup), as well
> as antechamber with AM1/BCC charge.

I believe David's tool requires you to provide a mapping (in fact, the
idea is to use Lomap) but at first glance your topology seems fine.
That's why I have suggested to have a look into FESetup which can do
that all for you.  Current downside is that it supports AMBER force
fields only.


> Then I set up the mdp file and topology file for relative free energy
> calculations. I used the single topology together with the one-step
> and three-step transformations, respectively. Please have a look at
> the relevant files  attached for the details.
> 
> Finally, I used g_bar to calculate the free energy
> (alchemistry_analysis.py can give similar results). From the results
> I found that one-step transformation can give good results, only if
> some restraints are applied to the protein-ligand complex, because I
> found that benzene/phenol are not strong binders to lyszome, and
> would drift around in the binding pocket, thus give poor results
> deviating from experimental data.

You should not need to use restraints for this type of relative AFE
simulation.  If you still do I would think you would have to take
the energy cost of the restraint into account as well.  The ligands may
be weak binders but you should see the same effect in other codes as
well.  The tutorial is run only for a very short time.  Making it
converge may have a large effect on the final ddG too.


> However, in the three-step  strategy, the results are much worse,
> this is probably due to the fact that a dual topology is used in
> Amber but a single topology is used in Gromacs (my case), we may not
> be able to do the changes directly (as shown in the attached topology
> files).

The implementation details here should not matter.  With AMBER you can
and you obviously did map atoms which at the end means that you
effectively use a "single topology" (AMBER does energy scaling while
Gromacs may do parameter scaling but I would have to check).  The real
difference is that AMBER does not require you to use explicit "dummy"
atoms and that's how the A9 tutorial has been set up.

I think the real issue with A9 is that for a reason unknown to me the
original author set up the transformation such that a hydrogen (in
benzene) and the -OH group (in phenol) are duplicated (so you _are_
actually using partial dual topology with both codes).  I do not see why
that would be beneficial in this case and would map the hydrogen to the
oxygen to only have a single disappearing atom.  This would also avoid
the awkward creation of partial charges and simplify the protocol.
When there are only disappearing atoms, setup with Gromacs is very
easy: you only need a single topology file and can vary electrostatics
and vdW lambda separately in a single mdp file.  Results should be the
same.


> Even though one-step transformation seems to give a better result in
> this case, it may be due to luck, because it is suggested that
> one-step transformation is usually not as good as three-step
> transformation, and can even give misleading results in some
> occasions (what occasions?).

The one-step protocol may be a problem with AMBER and relative AFEs.
I have tested this with a set of small organic molecules to compute the
free energy of hydration.  It may be that this depends on the number of
dummy atoms or where the dummies are located or both or...  I have not
found out yet why that is.  The other problem is that AMBER has
difficulties reproducing the end-point geometries e.g. in some cases
bond-length involving dummies are too long.  This is particularly
strange given that AMBER does actually use correct end-points.

In a quick test with Gromacs this did not seem a problem i.e. the
one-step protocol appeared to work satisfactorily.  Maybe also that with
AMBER the situation would be different with binding free energies
(maybe any errors would just happen to cancel out in the thermodynamic
cycle).  But that could be easily checked with the A9 tutorial...


Cheers,
Hannes.
-- 
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 

Re: [gmx-users] Gromacs Relative Free Energy Calculation Issue

[Hannes Loeffler]

Hi,

you could have a look into http://www.hecbiosim.ac.uk/fesetup to set up
your system.

It is late evening here so I will answer in more detail tomorrow.

Cheers,
Hannes.


On Mon, 3 Oct 2016 19:57:07 +
Guanglin Kuang  wrote:

> Dear Gromacs users,
> 
> Has any of you managed to use Gromacs to do relative free energy
> calculations?  I have some technical questions that would need your
> suggestions.
> 
> I am trying to reproduce the Amber free energy calculation using
> PMEMD (http://ambermd.org/tutorials/advanced/tutorial9/#overview).
> 
> First, I generated the topology and coordinate files of the ligands
> (benzene and phenol) using alchemistry_setup.py, developed by Dr.
> David Mobley (https://github.com/MobleyLab/alchemical-setup), as well
> as antechamber with AM1/BCC charge.
> 
> Then I set up the mdp file and topology file for relative free energy
> calculations. I used the single topology together with the one-step
> and three-step transformations, respectively. Please have a look at
> the relevant files  attached for the details.
> 
> Finally, I used g_bar to calculate the free energy
> (alchemistry_analysis.py can give similar results). From the results
> I found that one-step transformation can give good results, only if
> some restraints are applied to the protein-ligand complex, because I
> found that benzene/phenol are not strong binders to lyszome, and
> would drift around in the binding pocket, thus give poor results
> deviating from experimental data. However, in the three-step
> strategy, the results are much worse, this is probably due to the
> fact that a dual topology is used in Amber but a single topology is
> used in Gromacs (my case), we may not be able to do the changes
> directly (as shown in the attached topology files).
> 
> Even though one-step transformation seems to give a better result in
> this case, it may be due to luck, because it is suggested that
> one-step transformation is usually not as good as three-step
> transformation, and can even give misleading results in some
> occasions (what occasions?).
> 
> Therefore, could any experienced users of relative free energy
> calculation have a look at my set-up and point out the problems? I
> would appreciate any improvement to the set-up!
> 
> Thank you!
> 
> ---
> Best Regards!
> Guanglin Kuang
> KTH-Royal Institute of Technology
> School of Biotechnology
> Division of Theoretical Chemistry & Biology
> 

-- 
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] typeB in FEP calculations

[Hannes Loeffler]

On Tue, 27 Sep 2016 11:08:45 +0100
Rui Neves  wrote:

> However, what I would like to know is:
>  At the beggining of the topology I call for all the bonded terms
> ('#include ffbonded.itp), and then in the [ atoms ] section, I
> specify the massB, typeB and chargeB for the B-state. Do I have to
> explicitly write the parameters for the B-state of the [ bonds ],
> [ angles ] and [ dihedrals ] section, or are they read from the
> ffbonded.itp file, since I have defined the typeB for the atoms in
> the B-state?

If you wish to have the same parameters (for non-changing bonded terms)
there is not need to fill in the B state. However, if the bonded
parameters change, your topology file will necessarily have to reflect
this and you need to change the B state accordingly.  Grompp doesn't
make any assumption as to what you want esp. in the case of
disappearing/appearing atoms (where the typically approach is to take
the parameters from the other state i.e. the one where the bonded term
exists).  For torsions Gromacs would require a term-by-term matching
(same periodicities).

Cheers,
Hannes.
-- 
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 Binding Question

[Hannes Loeffler]

On Wed, 21 Sep 2016 12:00:29 +
Abdülkadir KOÇAK  wrote:

> In terms of endstates, the state A is the real ligand complexed with
> Protein in water... I did not define dummy atoms for the ligand as
> the state B, which I believe I should have...

I'm not quite sure what you mean here.  I think any recent Gromacs
version (maybe from 4.x?) allows you to use the couple-* parameters in
the mdp file.  In this way you do not need a modified topology.  On the
contrary, you would use a topology just as for standard MD simulation.
-- 
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 Binding Question

[Hannes Loeffler]

On Tue, 20 Sep 2016 13:27:06 +
Abdülkadir KOÇAK <ko...@gtu.edu.tr> wrote:

> What I am doing is exactly as you said; decoupling the ligand once
> from the WT and once from the MUT. And now I see that I am getting
> "absolute" free energies in each case. I did not run MD for the
> proteins in the absence of the ligand. As you mentioned even for the
> absolute binding free energies, the value is too high.

Just to clarify: of course you would need to decouple the ligand in
solution as well if you really wanted to get the absolute binding free
energy.  Otherwise decoupling only from the protein would just mean
transforming the ligand to the vacuum.  But obviously that's the same
procedure for WT and MUT and so would not need to be carried out if you
are really only interested in the _relative_ stability of the two
complexes.

For the absolute BFE you would do the decoupling in solution (only the
ligand in the solvation box) only _once_ and subtract the result from
the protein decoupling step. That's essentially what you mentioned in
your original email.

So the transformation for obtaining the relative BFE in solution would
be

Lig-WT -> WT
Lig-MUT -> MUT

where Lig is decoupled from its protein environment.

For the absolute BFE for either Lig-WT or Lig-MUT also

Lig -> "nothing" (in solution)


The alternative procedure would be

WT -> MUT
Lig-WT -> Lig-MUT

where you do a relative transformation between the affected side-chains.
This is the same as the first cycle but run in vertical direction.  May
be difficult if that is not a point mutation.


> I don't know whether you could help answering for that with the .mdp
> file options used in the production MD step (I attached). I have a
> Zn2+ containing enzyme with charged ligand. The co-crystal structures
> for both WT and MUT with ligand are available.

Attachments don't work with this mailing list.


> As the protocol, I first energy minimized (in two step 1- steep 2-
> cg, 5000 total steps) , NVT (200 ps) and NPT (1 ns) equilibrated the
> systems (For each lamda values of course). Then the equilbrated
> structures were run 10 ns in final MD... For the decoupling
> parameters I used the same lambda values in the .mdp for each
> simulation step.

Ok, but how do you handle the vanishing end states?  See e.g.
http://www.alchemistry.org/wiki/Example:_Absolute_Binding_Affinity


> 
> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of
> Hannes Loeffler <hannes.loeff...@stfc.ac.uk> Sent: Tuesday, September
> 20, 2016 3:45:51 PM To: gromacs.org_gmx-users@maillist.sys.kth.se Cc:
> gmx-us...@gromacs.org Subject: Re: [gmx-users] Free Energy of Binding
> Question
> 
> It usually helps to draw a thermodynamic cycle to inform yourself what
> exactly you are doing.
> 
> What you seem to have been doing is to decouple the ligand, once from
> the WT and once from the MUT.  This should give you the "absolute"
> binding free energy.  Computing the difference between those values
> will give you the relative binding free energy.  It is unclear to me
> what exactly you mean with dG=1134 kJ/mol but if this is the absolute
> binding free energy it seems to be excessively high.  But without
> knowing the details regarding the setup and simulation protocol it
> will not be possible to find out what could have gone wrong.
> 
> Of course, you could also compute the relative free energy by
> transforming WT to MUT, once in presence of the ligand and once
> without the ligand (side-chain mutation). Principally, both
> approaches should give you comparable results.
> 
> 
> On Tue, 20 Sep 2016 11:48:21 +
> Abdülkadir KOÇAK <ko...@gtu.edu.tr> wrote:
> 
> > Dear GMX Community,
> >
> > I am aiming to compare the relative binding energy (BE) of a ligand
> > to wild type (WT) vs mutant (MUT) protein and thus trying to run a
> > Free Energy Calculation for the binding energy of the ligand to both
> > proteins (WT and MUT) using Bennett Acceptance Ratio (BAR).
> >
> > As the first step, I calculated decoupling of the ligand from both
> > proteins in two seperate MD runs by first turning off the coulombic
> > interaction and then the van der waals interaction. The next step
> > would be the solvation free energies of the ligand in order to get
> > the correct BE. However, the ligand is the same for both WT and MUT,
> > so in terms of getting relative BE, should we still calculate the
> > ligand in water or will the ligand solv energies will cancel in the
> > DDG=DG(WT)-DG(MUT)?
> >
> > Besides, from the decoupling of the Protein-ligand complex, I am
> > getting very high DG valu

Re: [gmx-users] Free Energy of Binding Question

[Hannes Loeffler]

It usually helps to draw a thermodynamic cycle to inform yourself what
exactly you are doing.

What you seem to have been doing is to decouple the ligand, once from
the WT and once from the MUT.  This should give you the "absolute"
binding free energy.  Computing the difference between those values
will give you the relative binding free energy.  It is unclear to me
what exactly you mean with dG=1134 kJ/mol but if this is the absolute
binding free energy it seems to be excessively high.  But without
knowing the details regarding the setup and simulation protocol it will
not be possible to find out what could have gone wrong.

Of course, you could also compute the relative free energy by
transforming WT to MUT, once in presence of the ligand and once without
the ligand (side-chain mutation). Principally, both approaches should
give you comparable results.


On Tue, 20 Sep 2016 11:48:21 +
Abdülkadir KOÇAK  wrote:

> Dear GMX Community,
> 
> I am aiming to compare the relative binding energy (BE) of a ligand
> to wild type (WT) vs mutant (MUT) protein and thus trying to run a
> Free Energy Calculation for the binding energy of the ligand to both
> proteins (WT and MUT) using Bennett Acceptance Ratio (BAR).
> 
> As the first step, I calculated decoupling of the ligand from both
> proteins in two seperate MD runs by first turning off the coulombic
> interaction and then the van der waals interaction. The next step
> would be the solvation free energies of the ligand in order to get
> the correct BE. However, the ligand is the same for both WT and MUT,
> so in terms of getting relative BE, should we still calculate the
> ligand in water or will the ligand solv energies will cancel in the
> DDG=DG(WT)-DG(MUT)?
> 
> Besides, from the decoupling of the Protein-ligand complex, I am
> getting very high DG values (1134 kJ/mol). Is this meaningful or not?
> 
> Thx

-- 
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] The role of Jacobian factor in free energy simulations

[Hannes Loeffler]

Hi Stefano,

you might be better off either inspecting the code yourself or
contacting the developers directly (the Gromacs developer list may be
open to this).

Cheers,
Hannes.


On Mon, 12 Sep 2016 18:52:42 +
BOSISIO Stefano  wrote:

> Dear Gromacs staff,
> 
> I am trying to understand some internal issues between my code and
> Gromacs in alchemical free energy calculations.
> 
> Considering a simple alchemical free energy calculation (e.g ethane
> to methanol) with a constraint applied to all the hydrogen bonds
> does Gromacs calculate a Jacobian correction ? (Boresch Stefan, and
> Martin Karplus "The Jacobian factor in free energy simulations." The
> Journal of chemical physics105.12 (1996): 5145-5154  ) Alternatively,
> does Gromacs apply any correction to the computed free energy changes
> if constraints were applied to a solute?
> 
> 
> Thank you
> 
> Best regards,
> 
> Stefano
> 

-- 
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] couple-moltype in FEP

[Hannes Loeffler]

On Tue, 12 Jul 2016 12:44:10 +0200
Alexander Alexander  wrote:

> Yes, Fig.6 comes from the different step in Fig.1. I think the
> easiest way here is to have 6 different "topol-mdp" files
> corresponding to each step. Although unlike the amino acid, defining
> the [ moleculetype ] for ion and a limited number of water molecules
> out of thousand of water molecules in the box would not be that easy.

I think it should be fairly straightforward to define two new
moleculetyps with a new name each, rename and reorder the coordinate
file accordingly and make the few changes in the [molecules] section.


> The harmonic restrants have been applied in step 1 and 3, can you
> please let me know what is the difference between them and also how I
> can apply each of them in md files?

You would need to read the paper.  Since this is a decoupling exercise
my guess would be that the restraints are relative to the Ti surface.
-- 
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] couple-moltype in FEP

[Hannes Loeffler]

On Tue, 12 Jul 2016 10:29:04 +0100
Hannes Loeffler <hannes.loeff...@stfc.ac.uk> wrote:

> On Tue, 12 Jul 2016 10:16:24 +0200
> Alexander Alexander <alexanderwie...@gmail.com> wrote:
> 
> > Hi,
> > Thanks for your response.
> > I want to calculate the binding free energy of a single amino acid
> > to a solid surface in aqueous solution by FEP, alchemical
> > transformation, where perturbations applying on different species in
> > ordered steps. For instance: Step.1 applying harmonic restraints  to
> > the O and H of the amino acid. Step 2. Charges removal from the
> > amino acid and two water molecules and an Ion, happen.
> 
> I think the key figure here is actually Fig. 1.  What that paper
> attempts to do is to replace a molecule on a surface with two water
> molecules through a decoupling scheme.  I guess with Gromacs you would
> have to manipulate the topology files to some extent because the
> decharging and recharging steps involve different species (the sodium
> ion is there only to offset any changes in total charge).
> couple-lamnda0/1 may be helpful for some steps.  There are also both
> insertions and deletions so I think step 3 must be split into two vdW
> (deletion) steps.

Thinking about step 3 you could probably also do it as a relative
transformation i.e. omit couple-moltype and work with explicit A and B
columns.
-- 
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] couple-moltype in FEP

[Hannes Loeffler]

On Tue, 12 Jul 2016 10:16:24 +0200
Alexander Alexander  wrote:

> Hi,
> Thanks for your response.
> I want to calculate the binding free energy of a single amino acid to
> a solid surface in aqueous solution by FEP, alchemical
> transformation, where perturbations applying on different species in
> ordered steps. For instance: Step.1 applying harmonic restraints  to
> the O and H of the amino acid. Step 2. Charges removal from the amino
> acid and two water molecules and an Ion, happen.

I think the key figure here is actually Fig. 1.  What that paper
attempts to do is to replace a molecule on a surface with two water
molecules through a decoupling scheme.  I guess with Gromacs you would
have to manipulate the topology files to some extent because the
decharging and recharging steps involve different species (the sodium
ion is there only to offset any changes in total charge).
couple-lamnda0/1 may be helpful for some steps.  There are also both
insertions and deletions so I think step 3 must be split into two vdW
(deletion) steps.  Probably best to work this out by going through
steps 1 to 6 an pay careful attention to what species have zeroed out
charges.

It's unclear to me why you need to append any data.  What coordinates
you use doesn't really matter because each run is independent from
every other run (unless you see a specific need to take over a
coordinate from an adjacent lambda).
-- 
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] couple-moltype in FEP

[Hannes Loeffler]

On Tue, 12 Jul 2016 07:35:00 +0100
Hannes Loeffler <hannes.loeff...@stfc.ac.uk> wrote:

> On Tue, 12 Jul 2016 01:56:42 +0200
> Alexander Alexander <alexanderwie...@gmail.com> wrote:
> 
> > I was wondering that how I can have for example two different
> > "couple-moltype" in free energy part of my *.mdp file in FEP,
> > alchemical transformation? The reason for having two is to perturb
> > each of them differently in ordered steps.
> 
> You don't quite say what you are trying to do but you can always work
> with more than one topology file.

I got that one wrong.  When you use couple-moltype you obviously want
to make a whole molecule disappear without touching the topology.
Can't you just do two runs with different couple-moltype?

-- 
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] couple-moltype in FEP

[Hannes Loeffler]

On Tue, 12 Jul 2016 01:56:42 +0200
Alexander Alexander  wrote:

> I was wondering that how I can have for example two different
> "couple-moltype" in free energy part of my *.mdp file in FEP,
> alchemical transformation? The reason for having two is to perturb
> each of them differently in ordered steps.

You don't quite say what you are trying to do but you can always work
with more than one topology file.
-- 
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] Software generating gaff for gromacs with am1bcc charges

[Hannes Loeffler]

On Fri, 3 Jun 2016 12:47:34 +0700
Andrian Saputra  wrote:

> Hi all, can anyone suggest me whats softwares can produce gaff
> topology for gromacs with am1bcc charges automatically for 100-200
> atoms ?

All the software I'm aware of wraps around antechamber and so won't do
any magic to fix problems with the input structure.  With that number
of atoms I also wonder what that particular molecule is and whether
you better split that into smaller pieces.


> I tried antechamber and always found error...

That's an entirely useless statement and nobody will be able to help
you with that.  You may want to go through the relevant AMBER tutorial
to learn how to use the tool and ask antechamber specific questions on
the AMBER mailing list.

-- 
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] Clarity on TI free energy terms

[Hannes Loeffler]

On Thu, 2 Jun 2016 08:27:15 +
"Nash, Anthony"  wrote:

> Dear Hannes,
> 
> Thanks for all the help yesterday, it helped. I, hopefully, have just
> the one final question.
> 
> I am still a little confused how Gromacs deals with the interactions
> (vdW & Coul) with the environment when a soft-core potential has been
> used to switch between molecules (typeA and typeB in the topology
> file). I understand how lambda can be used to 1) phase out the charge
> then 2) phase out the vdW interactions for a molecule that is
> disappearing (or reverse for appearing) to capture hydration
> energetics e.g., Justin¹s methane in water FE tutorial. 
> 
> However, in the case of a D2K (aspartic-acid to lysine)
> dual-topology

How is your setup dual-topolgy?  http://dx.doi.org/10.1002/jcc.23804
indicates to me that the pmx approach seems to go for as much single
topology as possible.  Your case seems to fall into this particular
category.


> how are the vdW and charges brought back into play for
> typeB, when the interactions of typeA have been coupled to lambda as
> per your example below? 

In a second setup as explained earlier, assuming you want to separate
vdW from Coulomb.  See attachment for an unrelated example system.
-- 
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] Clarity on TI free energy terms

[Hannes Loeffler]

On Wed, 1 Jun 2016 16:15:31 +
"Nash, Anthony"  wrote:

> In the mean while, do you know of any tutorials (besides the methane
> in water FE tutorial) regarding TI for amino acid substitution? 

I am not aware of one.  You could try
http://www.alchemistry.org/


> And by “q_off” and “vdw_on/off”, I assume you are referring to the
> ‘value’ of the couple-lambda0: parameter?

No.  I'm referring to the respective lambda paths.


Cheers,
Hannes.
-- 
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] Clarity on TI free energy terms

[Hannes Loeffler]

On Wed, 1 Jun 2016 15:00:51 +
"Nash, Anthony"  wrote:

> >  This also assumes that
> >you have vanishing atoms only.  If you have appearing atoms only you
> >would obviously have to revers the order, and when you have both you
> >will have to run with two mdp/tpr setups.
> 
> 
> With aspartic acid transforming into lysine on one polypeptide chain
> and then lysine transforming into aspartic acid polypeptide chain,
> all in the same system (keeps the charges the same and is a complete
> thermodynamic cycle), I will have both appearing and disappearing
> atoms. How do you mean ³to run with two mdp/tpr setups²? Is there an
> example you could give (which I am grateful for) or one in the manual?

Lambda paths are not selective in the sense that you could apply them
to only a subset of the system.  So if you have both disappearing and
appearing atoms you have to:
1) turn off the charges on the disappearing group (or alternatively all
charges to avoid a charged total system), q_off
2) turn off the vdW parameters for the disappearing group, vdw_off
3) turn on the vdW parameters for the appearing group, vdw_on
4) turn on the charges of the appearing group (or all charges), q_on

If you try to do this with Gromacs you will realise that the best you
can do is combine this into 2 mdp steps: 1) q_off and vdw_on/off and 2)
q_on.  Alternatively you can combine the vdw bit with 2) but that
doesn't make any difference.  Of course, if you could assume that a
combined Coulomb/vdW transformation will work, this would be a
non-issue...

Also, I wonder how pmx maps aspartate onto lysine.  I would think that
the gamma-carboxylate should not map onto the gamma-methylene but
rather the residue should be branched off at the C-beta and so
duplicate the rest of the residue.
-- 
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] Clarity on TI free energy terms

[Hannes Loeffler]

On Wed, 1 Jun 2016 12:06:20 +
"Nash, Anthony"  wrote:

> vdw_lambdas  = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
> 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
> mass_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.0 0.0 0.0 0.0
> fep_lambdas =  0.0 0.1 0.1 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.5 0.6 0.6
> 0.7 0.7 0.8 0.8 0.9 0.9 1.0 1.0

This will transform both vdW _and_ Coulomb terms at the same time but
at a different pace.  Maybe you want something like

fep-lambdas = 0.0 0.2 0.4 0.6 0.8 1.0  1.0  1.0  etc.
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0  0.05 0.10 etc.
mass-lambdas as above

This will first transform Coulomb and bonded terms in the first six
lambdas and vdW from the 7th lambda onwards.  This also assumes that
you have vanishing atoms only.  If you have appearing atoms only you
would obviously have to revers the order, and when you have both you
will have to run with two mdp/tpr setups.


> couple-moltype   = protein
> couple-lambda0   = vdw-q
> couple-lambda1   = vdw-q

This will decouple all atoms in the 'protein' selection from the
environment.  This is for absolute transformation and not what you seem
to want to do i.e. a relative transformation of one residue into
another.  So avoid those parameters.
-- 
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] Clarity on TI free energy terms

[Hannes Loeffler]

Set the vector to all-zeroes (or ones).


On Wed, 1 Jun 2016 09:47:59 +
"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:

> Hi Hannes,
> 
> 
> Many thanks for the reply. With regards to your final comment I
> understand conserving mass in theory, but I am a little confused
> regarding, "keep the mass-lambdas at one end-point as they can
> interact badly with constraints". I am testing pmx on a two-molecule
> one-system I.e., G-D2K-G and G-K2D-G in the same system. How ought I
> define the mass-lambdas for this system? (nothing accurate, just an
> example would be great)?
> 
> Thanks
> Anthony
> 
> 
> On 01/06/2016 09:55,
> "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
> Hannes Loeffler" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> on behalf of hannes.loeff...@stfc.ac.uk> wrote:
> 
> >On Wed, 1 Jun 2016 07:54:56 +
> >"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
> >
> >> In the tutorial, charges are off in the topology and the
> >> electrostatic coupling to lambda remains 0 throughout the 20
> >> windows. I assume setting col_lambdas=0 0 0 Š was for that very
> >> reason I.e., the charges were off? Could the charges not have been
> >> left on and col_lambdas defined similar to vdw_lambdas?
> >> (I understand that if charges remain constant, as vdw turns off,
> >> the system will probably blow up as attraction brings molecules
> >> infinitely close).
> >
> >Technically, Gromacs allows you to vary both vdW and Coulomb lambdas
> >simultaneously because Gromacs can apply softcore potentials to both.
> >In practice though it seems that many workers still prefer to
> >separate the two terms from each other.
> >
> > 
> >> If my transition is from a small molecule into a small molecule
> >> e.g., G-D-G to G-K-D, (the PMX paper) should I define all three
> >> lambdas: vdw_lambdas, col_lambdas and bonds_lambdas? Between
> >> states A and B, VdW, charges and bonds are all changing.
> >
> >Lambda paths are only about separating the various force field terms
> >from each other.  If you do not explicitly state any of those lambda
> >vectors they will adopt they same lambdas as specified in
> >fep-lambdas, see manual.  I do not see a reason why you would want
> >to separate out the bonded terms as well.  They are subject to a
> >linear transformation only anyway.
> >
> >What you may want to do is to keep the mass-lambdas at one end-point
> >as they can interact badly with constraints.  In a proper
> >thermodynamic cycle mass contributions must perfectly cancel.
> >
> >
> >Cheers,
> >Hannes.
> >-- 
> >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] Clarity on TI free energy terms

[Hannes Loeffler]

On Wed, 1 Jun 2016 07:54:56 +
"Nash, Anthony"  wrote:

> In the tutorial, charges are off in the topology and the electrostatic
> coupling to lambda remains 0 throughout the 20 windows. I assume
> setting col_lambdas=0 0 0 Š was for that very reason I.e., the
> charges were off? Could the charges not have been left on and
> col_lambdas defined similar to vdw_lambdas? 
> (I understand that if charges remain constant, as vdw turns off, the
> system will probably blow up as attraction brings molecules infinitely
> close).

Technically, Gromacs allows you to vary both vdW and Coulomb lambdas
simultaneously because Gromacs can apply softcore potentials to both.
In practice though it seems that many workers still prefer to separate
the two terms from each other.

 
> If my transition is from a small molecule into a small molecule e.g.,
> G-D-G to G-K-D, (the PMX paper) should I define all three lambdas:
> vdw_lambdas, col_lambdas and bonds_lambdas? Between states A and B,
> VdW, charges and bonds are all changing.

Lambda paths are only about separating the various force field terms
from each other.  If you do not explicitly state any of those lambda
vectors they will adopt they same lambdas as specified in fep-lambdas,
see manual.  I do not see a reason why you would want to separate out
the bonded terms as well.  They are subject to a linear transformation
only anyway.

What you may want to do is to keep the mass-lambdas at one end-point as
they can interact badly with constraints.  In a proper thermodynamic
cycle mass contributions must perfectly cancel.


Cheers,
Hannes.
-- 
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] Release of FESetup 1.2

[Hannes Loeffler]

Hello everyone,

we are delighted to announce the release of FESetup 1.2.

FESetup is a tool to automate the setup of (relative) alchemical free
energy simulations like thermodynamic integration (TI) and free energy
perturbation (FEP) as well as post–processing methods like MM–PBSA and
LIE.  FESetup can also be used for general simulation setup
("equilibration") through an abstract MD engine.  The latest releases
are available from our project web page at
http://www.hecbiosim.ac.uk/fesetup

Release notes can be found at
https://ccpforge.cse.rl.ac.uk/gf/project/ccpbiosim/wiki/?pagename=Releases


Many thanks,
Hannes Loeffler.
-- 
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] Suggested settings for using Amber force field in Gromacs

[Hannes Loeffler]

On Thu, 5 May 2016 13:20:59 +
"Casalini  Tommaso"  wrote:

> Thanks a lot for your quick answer. 
> 
> I am trying to reproduce with Gromacs some results that I obtained
> with Amber software. I put the same thermostat (Langevin) and
> barostat (Berendsen with isotropic scaling) with same time constant,
> the same dt, the same cutoff, and so on. 

You will have to carefully study the manuals for both packages (read
the code actually...), check your settings and see how close you can
get.  But you will need to accept that differences in implementation
for possible lack thereof may exist.  How much your particular system
will depends on such details, I don't know.

BTW, the defaults settings for pmemd/sanders are vdwmeth = 1 (vdW
long-range dispersion correction on) and fswitch = -1 (truncation
cutoff).


Cheers,
Hannes.
-- 
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] Moving this mailing list to Discourse format

[Hannes Loeffler]

Hi Mark,

many thanks.

At the end it is really just a question if the user has flexibility
over how to visualise the communications.  I'm happy to hear that the
old-fashioned email way would still have (some) support :-)

I hope you get sufficient feedback and support from the community.

Cheers,
Hannes.


On Thu, 21 Apr 2016 08:47:56 +
Mark Abraham <mark.j.abra...@gmail.com> wrote:

> Hi Hannes,
> 
> Yes, previous discussion has identified this as a likely key feature,
> and thanks for adding confirmation of that. Discourse can be
> integrated with functionality for replying to threads from email to
> suit people who want that kind of small change to their workflows;
> that will indeed suit simple discussions. Obviously there's
> configurable email activity digests and notifications also, which
> provide links to the forum proper.
> 
> Mark
> 
> On Thu, Apr 21, 2016 at 9:52 AM Hannes Loeffler
> <hannes.loeff...@stfc.ac.uk> wrote:
> 
> > On Thu, 21 Apr 2016 07:11:42 +
> > Mark Abraham <mark.j.abra...@gmail.com> wrote:
> >
> > > Hi people,
> > >
> > > As part of the new BioExcel project (http://bioexcel.eu/) that
> > > supports biomolecular research with codes including GROMACS, we're
> > > exploring migrating the gmx-users mailing list to the Discourse
> > > format. You can see an example at http://try.discourse.org! This
> > > has a number of immediate advantages, including attachments,
> > > single-click sign-up, tagging, voting and contributing to old
> > > discussions. We might also integrate features like replying to
> > > threads via email, or real-time chat. We also know such a change
> > > might be something people don't like, and we would like to get a
> > > feel for that. The community should like its forum, after
> > > all! :-) So please let us know how you feel at
> > > http://goo.gl/forms/k8lhTliffO.
> >
> > I think such a discussion would need a bit more space than just a
> > few pre-canned questions.
> >
> > I am on various mailing lists and happy with my email client as a
> > central "collector" of all those discussions.  I am personally not
> > fond of another web page interface with its very own features,
> > bells and whistles and what not, but typically unlike any other
> > existing one. And this is in essence the major shortcoming unless
> > it would be still possible to receive all messages through a
> > mailing list as it is now.
> >
> > Cheers,
> > Hannes.
> > --
> > 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] Moving this mailing list to Discourse format

[Hannes Loeffler]

On Thu, 21 Apr 2016 07:11:42 +
Mark Abraham  wrote:

> Hi people,
> 
> As part of the new BioExcel project (http://bioexcel.eu/) that
> supports biomolecular research with codes including GROMACS, we're
> exploring migrating the gmx-users mailing list to the Discourse
> format. You can see an example at http://try.discourse.org! This has
> a number of immediate advantages, including attachments, single-click
> sign-up, tagging, voting and contributing to old discussions. We
> might also integrate features like replying to threads via email, or
> real-time chat. We also know such a change might be something people
> don't like, and we would like to get a feel for that. The community
> should like its forum, after all! :-) So please let us know how you
> feel at http://goo.gl/forms/k8lhTliffO.

I think such a discussion would need a bit more space than just a few
pre-canned questions.

I am on various mailing lists and happy with my email client as a
central "collector" of all those discussions.  I am personally not fond
of another web page interface with its very own features, bells and
whistles and what not, but typically unlike any other existing one.
And this is in essence the major shortcoming unless it would be still
possible to receive all messages through a mailing list as it is now.

Cheers,
Hannes.
-- 
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] Thermodynamic integration

[Hannes Loeffler]

TI/FEP are appropriate when you have "small" changes between states.
Whether you delete/create all of a molecule or just part of it doesn't
matter that much.  So a relative simulation between a glycated and
non-glycated sounds quite reasonable to me.  Of course, this will become
increasingly more difficult with increasing chain length of
carbohydrates but you could do in additional steps if necessary.  What
is important is that you draw up a proper thermodynamic cycle to ensure
that you really calculate what you want to.

On Mon, 18 Apr 2016 11:06:40 +
"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:

> Hi Mark,
> 
> I will give that a thorough read. I was wondering if you could
> possibly comment on whether TI is an appropriate tool for calculating
> the free energy difference between two states, A―> non-glycated side
> chain, and b―> glycated side chain? Most examples given focus on the
> inclusion/exclusion of some small molecule in some large system.
> 
> I have tried umbrella sampling and although my results are extremely
> interesting, I’ve had to manipulate the initial systems to take it
> from a crystal structure with defined periodicity in x-y-z
> dimensions, to a slab in the x-y plane, and a water bath in the z
> plane. 
> 
> Thanks
> Anthony 
> 
> 
> 
> 
> On 18/04/2016 11:53,
> "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of Mark
> Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf
> of mark.j.abra...@gmail.com> wrote:
> 
> >Hi,
> >
> >Also you might consider pmx
> >http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4365728/ for such
> >topology generation. There is further work in the pipeline, so do
> >get in touch with Bert if there's something of interest.
> >
> >Mark
> >
> >On Mon, Apr 18, 2016 at 11:56 AM Nash, Anthony <a.n...@ucl.ac.uk>
> >wrote:
> >
> >> From the site, “..or the free energy of a mutation of a side
> >> chain.”
> >>
> >> I think this is what I am after. Many thanks for the link.
> >>
> >> Anthony
> >>
> >>
> >>
> >> On 18/04/2016 10:42,
> >> "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> >>on
> >> behalf of Hannes Loeffler"
> >> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
> >> hannes.loeff...@stfc.ac.uk> wrote:
> >>
> >> >A good starting point is http://www.alchemistry.org/ which has
> >> >quite a lot of detail on relative alchemical free energy
> >> >simulations (not only TI).
> >> >
> >> >On Mon, 18 Apr 2016 09:27:02 +
> >> >"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
> >> >
> >> >> Hi all,
> >> >>
> >> >> I¹m looking for a guide on performing TI between a protein in
> >> >> its crystal periodicity with a particular residue (state A), to
> >> >> the same system but with a different residue (state B).
> >> >>
> >> >> I¹m currently using
> >> >>
> >> >>
> >> 
> >>http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/fr
> >> >>ee
> >> >> _energy/01_theory.html as a guide, however this is more, if I
> >> >> understand having read through it, on the presence-to-absence of
> >> >> methane in a water solvent rather than a replacement with
> >> >> something else.
> >> >>
> >> >> I haven¹t had too much luck googling and I¹m looking piecemeal
> >> >> through the manual with little success.
> >> >>
> >> >> Thanks
> >> >> Anthony
> >> >>
> >> >>
> >> >>
> >> >
> >> >--
> >> >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] Thermodynamic integration

[Hannes Loeffler]

A good starting point is http://www.alchemistry.org/ which has quite a
lot of detail on relative alchemical free energy simulations (not only
TI).

On Mon, 18 Apr 2016 09:27:02 +
"Nash, Anthony"  wrote:

> Hi all,
> 
> I¹m looking for a guide on performing TI between a protein in its
> crystal periodicity with a particular residue (state A), to the same
> system but with a different residue (state B).
> 
> I¹m currently using
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free
> _energy/01_theory.html as a guide, however this is more, if I
> understand having read through it, on the presence-to-absence of
> methane in a water solvent rather than a replacement with something
> else.
> 
> I haven¹t had too much luck googling and I¹m looking piecemeal
> through the manual with little success.
> 
> Thanks
> Anthony  
> 
> 
> 

-- 
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] Relative binding free energy

[Hannes Loeffler]

On Mon, 14 Mar 2016 14:08:32 +
Stefania Evoli  wrote:

> Please, could you give a look at that? I obtained it by using the
> useful tool described in the article 'A Python tool to set up
> relative free energy calculations in GROMACS¹ Pavel V. Klimovich and
> David L. Mobley. Probably I¹m doing something wrong in the MCSS map.

As far as I know Pavel and David's tool reads in a ready-made MCS and
creates a Gromacs topology from it.  You seem to indicate the same.  I
trust that their tool works so you would now need to make sure that
your MCS map is what you think it should be.  I have indicated earlier
what I think it should look like.

I'm writing a tool my self, FESetup
(http://www.hecbiosim.ac.uk/fesetup), which can also create Gromacs
topologies for relative free energy calculations.  I must say though
that at the moment I do not support both appearing and disappearing
atoms in one for Gromacs (and I do no recommend you to that anyway).
Also, the MCSS is based on rather loose rules so may not give the
expected results in your case. But eventually it is based on the
algorithm in RDKit which is fairly easy to use.

If you need a map only once though it shouldn't be too difficult to do
that by hand as indicated earlier.


> You mentioned that I can simultaneously switch off the charges for
> disappearing groups and switch on for the appearing groups with a
> single lambda path. Could you clarify this point and if it is
> possible suggest me someway to do that, please?

No you can't.  That is exactly why you need _two_ topology files e.g. in
the first do q_off followed by vdw on/off.  The second topology file
takes care of q_on for the appearing group (see previous email for
details and the attachment).  You cant switch the electrostatics in the
same topology because you, crucially, have to do it in that specific
order and Gromacs does not allow you to say that a certain lambda path
only applies to a subset of atoms.

But I anyway do not recommend this for your case.  Make your life
easier and just to a maximal map and you will only have disappearing
atoms (probably best to set up your system to mutate from the larger
to the smaller molecule).


Cheers,
Hannes.
-- 
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] Relative binding free energy

[Hannes Loeffler]

I had a quick look at this and I see quite a few problems there.

You don't show the [atomtypes] but I suppose that whatever you call
*_dummy has zero vdW parameters.  BTW, there is no need to invent
different atom type names for every dummy as their non-bonded parameters
are all zero anyway and bonded terms are indexed.

You say this is single topology, well it is, but only almost.  You make
the alpha-H in glycine and the beta-C in serine a dummy which means
that you are duplicating.  A possible problem with this is that any
intermediate state has the alpha-C now five-fold coordinated! You
should always be wary of such setups because this will introduce
additional bonded terms (you don't show them) that are not present in
the end state and can therefore distort the molecule and also lead to
instabilities.

Another problem with this is that you have dummy atoms at both end
states.  With Gromacs you will need _two_ topolgies to ensure that you
change electrostatics and vdW in the correct order.  You will need to
first switch off the charges on the disappearing group followed by
modifying the vdW parameters between the end states (you can do that
simultaneously that is switch off for disappearing groups and switch on
for the appearing groups with a single lambda path).  Eventually, you
switch on the charges for the appearing group.  There are some
varieties in how to do this in practice but that's basically it.

So you have used a MCSS to do the mapping and I guess the algorithm
that you have used actually mapped the alpha-H with beta-C.  I suggest
that you keep it that way to truly have a single topology description
as explained above.  BTW, the typical MCSS algorithm is 2D only and as
such has no idea about stereochemistry.  So make sure that the right
hydrogen maps to the right carbon atom (you say you have D-serine).

Judging from the atom type names it looks to me that those are GAFF
types.  If so I suggest to actually use AMBER peptide parameters (may
be available in the literature or online for single residue zwitterions)
esp. when used together with a protein.  A particular problem here is
that charge derivation will happen in the gas phase and you have
zwitterions! Geometry optimisation of such structures will very likely
lead to rather distorted and unnatural conformations.

HTH,
Hannes.


On Mon, 14 Mar 2016 12:00:28 +
Stefania Evoli  wrote:

> The [atom] section of the single topology of my transformation looks
> like
> 
> [ atoms ]
> ;   nr   type  resnrres   atom   cgnr charge   mass
> typeBchargeB  massB comments
>  1  c  1LIG C1  1   0.932601   12.01000
>   c   0.935601   12.01000 ; MCSS
>  2 c3  1LIGCA1  2  -0.095200   12.01000
>  c3  -0.126500   12.01000 ; MCSS
>  3   c3_dummy  1LIGCB1  30.0   12.01000
>  c3   0.136400   12.01000 ; to be appeared
>  4 hx  1LIGHA1  4   0.0887001.00800
> hx_dummy0.01.00800 ; to be annihilated
>  5 hx  1LIGHA2  5   0.0887001.00800
>  hx   0.0947001.00800 ; MCSS
>  6   h1_dummy  1LIGHB1  60.01.00800
>  h1   0.0682001.00800 ; to be appeared
>  7   h1_dummy  1LIGHB2  70.01.00800
>  h1   0.0682001.00800 ; to be appeared
>  8   ho_dummy  1LIGHG1  80.01.00800
>  ho   0.4480001.00800 ; to be appeared
>  9 n4  1LIG N1  9  -0.835601   14.01000
>  n4  -0.822601   14.01000 ; MCSS
> 10  o  1LIG O1 10  -0.755301   16.0
>   o  -0.759801   16.0 ; MCSS
> 11   oh_dummy  1LIGOG1 110.0   16.0
>  oh  -0.628801   16.0 ; to be appeared
> 12  o  1LIG   OXT1 12  -0.755301   16.0
>   o  -0.759801   16.0 ; MCSS
> 13 hn  1LIG H1 13   0.4438001.00800
>  hn   0.4488001.00800 ; MCSS
> 14 hn  1LIG H2 14   0.4438001.00800
>  hn   0.4488001.00800 ; MCSS
> 15 hn  1LIG H3 15   0.4438001.00800
>  hn   0.4488001.00800 ; MCSS
-- 
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] Relative binding free energy

[Hannes Loeffler]

How does the final structure for lambda=0 minimisation look like?  This
may give you a hint as to what is wrong.

What exactly does "do not work" mean?  That's not very descriptive...

It may help a lot to actually post the relevant A and B columns from
your topologies.  Single topology means in your case that you map both
backbones plus one hydrogen onto each other.  The other hydrogen (the
right one hopefully!) maps onto the C-beta of serine.  The remaining
atoms need to be defined as dummy atoms in the glycine state. Have you
had a look into pmx to have this done for you?


On Mon, 14 Mar 2016 11:05:14 +
Stefania Evoli <stefania.ev...@kaust.edu.sa> wrote:

> The minimization crush after 1182 steps like follow
> 
> Step= 1116, Dmax= 1.2e-02 nm, Epot= -1.76932e+05 Fmax= 9.54828e+03,
> atom= 1 Step= 1117, Dmax= 1.4e-02 nm, Epot= -1.77307e+05 Fmax=
> 2.85691e+04, atom= 8 Step= 1119, Dmax= 8.6e-03 nm, Epot= -1.99060e+05
> Fmax= 3.84768e+07, atom= 8 Step= 1124, Dmax= 6.5e-04 nm, Epot=
> -2.08560e+05 Fmax= 7.94965e+07, atom= 8 Step= 1126, Dmax= 3.9e-04 nm,
> Epot= -2.34096e+05 Fmax= 2.56256e+08, atom= 8 Step= 1128, Dmax=
> 2.3e-04 nm, Epot= -2.39440e+05 Fmax= 3.06063e+08, atom= 8 Step= 1130,
> Dmax= 1.4e-04 nm, Epot= -3.61664e+05 Fmax= 2.64411e+09, atom= 8 Step=
> 1133, Dmax= 4.2e-05 nm, Epot= -1.12275e+06 Fmax= 6.90055e+10, atom= 8
> Step= 1136, Dmax= 1.3e-05 nm, Epot= -1.31351e+06 Fmax= 9.96261e+10,
> atom= 8 Step= 1138, Dmax= 7.5e-06 nm, Epot= -3.71734e+06 Fmax=
> 9.66090e+11, atom= 8 Step= 1141, Dmax= 2.3e-06 nm, Epot= -1.52957e+07
> Fmax= 1.76138e+13, atom= 11020 Step= 1144, Dmax= 6.8e-07 nm, Epot=
> -2.62029e+07 Fmax= 5.21932e+13, atom= 81020
> Step= 1146, Dmax= 4.1e-07 nm, Epot= -4.13002e+07 Fmax= 1.30309e+14,
> atom= 11020
> Step= 1148, Dmax= 2.4e-07 nm, Epot= -7.52155e+07 Fmax= 4.33870e+14,
> atom= 11020
> Step= 1150, Dmax= 1.5e-07 nm, Epot= -1.08171e+08 Fmax= 8.98648e+14,
> atom= 11020
> Step= 1152, Dmax= 8.8e-08 nm, Epot= -2.33175e+08 Fmax= 4.18302e+15,
> atom= 11020
> Step= 1154, Dmax= 5.3e-08 nm, Epot= -2.60654e+08 Fmax= 5.22788e+15,
> atom= 11020
> Step= 1156, Dmax= 3.2e-08 nm, Epot= -9.62197e+08 Fmax= 7.13107e+16,
> atom= 11020
> Step= 1159, Dmax= 9.5e-09 nm, Epot= -2.35812e+09 Fmax= 4.28405e+17,
> atom= 11020
> Step= 1162, Dmax= 2.8e-09 nm, Epot= -6.67786e+10 Fmax= 3.43602e+20,
> atom= 11020
> Step= 1168, Dmax= 1.1e-10 nm, Epot= -6.70934e+11 Fmax= 3.46861e+22,
> atom= 11020
> Step= 1172, Dmax= 1.6e-11 nm, Epot= -1.02118e+12 Fmax= 8.03509e+22,
> atom= 81020
> Step= 1174, Dmax= 9.6e-12 nm, Epot= -1.98975e+12 Fmax= 3.05062e+23,
> atom= 11020
> Step= 1176, Dmax= 5.8e-12 nm, Epot= -2.58715e+12 Fmax= 5.15738e+23,
> atom= 11020
> Step= 1178, Dmax= 3.5e-12 nm, Epot= -6.80566e+12 Fmax= 3.56885e+24,
> atom= 11020
> Step= 1181, Dmax= 1.0e-12 nm, Epot= -7.62947e+13 Fmax= 4.48513e+26,
> atom= 81020
> Step= 1182, Dmax= 1.2e-12 nm, Epot= -5.58957e+12 Fmax= 2.40735e+24,
> atom= 11020
> 
> 
> As you can see the problem seems to rise in correspondence of the
> atom 8 (an HG atom). I tried lambda 0 (works), lambda 1 and 2 (do not
> work) thus increasing lambda value the simulations crush. I¹m using a
> single topology approach to transform a glycine to a D-serine. I
> checked the topology file several times and I visualized the gro file
> too and they seem to e correct.
> 
> Thank you in advance for your help
> 
> ‹Dr. Stefania Evoli
> Post-Doctoral Fellow
> King Abdullah University of Science and Technology
> Catalysis center - Bldg. 3, 4th floor, 4231­WS18
> Thuwal, Kingdom of Saudi Arabia
> stefania.ev...@kaust.edu.sa
> 
> 
> 
> 
> 
> 
> On 3/14/16, 12:21 PM,
> "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
> Hannes Loeffler" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> on behalf of hannes.loeff...@stfc.ac.uk> wrote:
> 
> >At least in the 4.6 series it used to be the case that when
> >couple-moltype was not set (and that's what you want to do) then none
> >of the other couple-* parameters ever had any effect.  So there
> >should be no need for you to set any of those.
> >
> >The more important point though would be to provide much more detail
> >about the crash during minimization.  When does it happen?  At all
> >lambdas, just one?  What was the original structure, what the final
> >one? How have you created the topologies, any errors in there? etc.
> >
> >Cheers,
> >Hannes.
> >
> >
> >On Mon, 14 Mar 2016 08:40:41 +
> >Stefania Evoli <stefania.ev...@kaust.edu.sa> wrote:
> >
> >> Dear Users,
> >>
> >> I¹m performing relative binding free energy by using Gromacs 5.0.5.
> >> As I u

Re: [gmx-users] refining a TI calculation - additional lambda points

[Hannes Loeffler]

On Tue, 26 Jan 2016 15:47:03 +0100
Oliwia Maria Szklarczyk  wrote:

> Dear All,
> 
> I was wondering if you have tips on how to easily create additional
> lambda points in a thermodynamic integration calculation in gromacs
> 5. I have 21 lambda points already simulated, from L_0 to L_20, and
> now I need six more lambda-points between L_16 and L_19. 

As long as you are only interested in TI you can just add points as you
want.  You could, for example, just make a new mdp file with a new
lambda vector only having these points or reuse the old one with adding
the new points in.  That doesn't really matter much.

The practical issue is more how you combine the resulting data such
that your analysis program of choice can use it.  But it shouldn't be
too difficult to quickly script a solution to do that yourself.   You
will want to write the new TI data into a separate xvg file.

If you intend to do a FEP style analysis, however, you need to consider
that you would need to compute energies for at least the nearest
neighbours.  For BAR e.g. you would need to recalculate the windows
bordering L16 and L19 for full MBAR, in principle, all windows...


Cheers,
Hannes.
-- 
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] tutorial for "dual" alchemical transformation in gromacs?

[Hannes Loeffler]

You still have couple-moltype set.

On Tue, 1 Dec 2015 15:39:40 +0100
Sebastian Stolzenberg <s.stolzenb...@fu-berlin.de> wrote:

> OK, thanks a lot, Hannes,
> 
> I combined the "CL-" and "NA+" into one "NC" [moleculetype] and
> switched off the couple-* switches as you suggested (attached),
> but I am still receiving the same warning message:
> "The lambda=0 and lambda=1 states for coupling are identical"
> 
> I guess I am almost there, am I missing something else?
> 
> Thanks again,
> Sebastian
> 
> 
> 
> On 01.12.2015 14:22, Hannes Loeffler wrote:
> > Ok, my attachment seems to have been stripped, so inline only.  Hope
> > it's not too bad.
> > 
> > [ atomtypes ]
> > ;name   bond_type mass charge   ptype   sigma
> > epsilon   Amb
> > DU   DU  0.0  0.0   A   0.0e+00
> > 0.0e+00 ; 0.00  0.
> > OW   OW  0.0  0.0   A   3.15075e-01
> > 6.35968e-01 ; 1.77  0.1520
> > HW   HW  0.0  0.0   A   0.0e+00
> > 0.0e+00 ; 0.00 0.
> > 
> > [ moleculetype ]
> >   ; molname   nrexcl
> >   NA+CL-  1
> > 
> > [ atoms ]
> >   ;   DU->IP
> >   ;   IM->DU
> >   ; = IM->IP
> >   ; id_at type res nr  residu name at name  cg nr  charge
> > mass  typeB chargeB   massB
> > 1   DU  1NA+  NA+   1  0 22.9898
> > IP 0 22.9898
> > 2   IM 1 DU   DU2  0 35.4530
> > DU 0 35.4530
> > 
> > 
> > 
> > On Tue, 1 Dec 2015 13:17:38 +
> > Hannes Loeffler <hannes.loeff...@stfc.ac.uk> wrote:
> > 
> >> Dear Sebastian,
> >>
> >> I have got your attachment.  I do not know all intricacies of
> >> Gromacs but for relative free energies you only use one
> >> moleculetype for the perturbed group, see attachment.  You do not
> >> use the couple-* switches for relative simulations.
> >>
> >> Yes, you are right about typeA and typeB.  If you have larger
> >> molecules you would just combine both molecules in their entirety
> >> into one "pseudo" molecules.  One original molecules has the normal
> >> force field types in the A state and dummy atoms in the final,
> >> while the other molecule has this reversed.
> >>
> >>
> >> Cheers,
> >> Hannes.
> >>
> >>  
> >>
> >> On Tue, 1 Dec 2015 13:34:12 +0100
> >> Sebastian Stolzenberg <s.stolzenb...@fu-berlin.de> wrote:
> >>
> >>> Dear Hannes,
> >>>
> >>>> I have some trouble reading your top file in my web browser
> >>> Apologies, am not sure if this works in this list, but I am now
> >>> also attaching my topol.top and grompp.mdp files.
> >>> I hope this let's you see better what I am doing wrong.
> >>>
> >>>> As your two dummy types appear to be identical
> >>>> you can define just one.
> >>> Thank you, as you can see in topol.top, I now have only one dummy
> >>> atom type called "XX".
> >>>
> >>>> I think your transformations are IPX->Na+, IMX->Cl-.
> >>> my "typeA" columns are "XX" and "IM" (NA+)
> >>> my "typeB" columns are "IM" and "XX" (CL-)
> >>> If I understand correctly (do I?), lambda=0 corresponds to
> >>> "typeA", whereas lambda=1 corresponds to "typeB", so my
> >>> transformation is "XX" -> "Na+", "Cl-" -> "XX"
> >>>
> >>> Thanks a bunch,
> >>> Sebastian
> >>>
> >>>
> >>>
> >>> On 30.11.2015 17:30, hannes.loeff...@stfc.ac.uk wrote:
> >>>> Dear Sebastian,
> >>>>
> >>>> I have some trouble reading your top file in my web browser but I
> >>>> think your transformations are IPX->Na+, IMX->Cl-.  So you would
> >>>> have to reverse the columns for the first atom.  As your two
> >>>> dummy types appear to be identical you can define just one.
> >>>>
> >>>> Cheers,
> >>>> Hannes.
> >>>>
> >>>> 
> >>>> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> >>>> [gromacs.org_gmx-users-boun...@maillist.sys.kth.se

Re: [gmx-users] tutorial for "dual" alchemical transformation in gromacs?

[Hannes Loeffler]

Ok, my attachment seems to have been stripped, so inline only.  Hope
it's not too bad.

[ atomtypes ]
;name   bond_type mass charge   ptype   sigma
epsilon   Amb
DU   DU  0.0  0.0   A   0.0e+00   0.0e+00 ;
0.00  0.
OW   OW  0.0  0.0   A   3.15075e-01   6.35968e-01 ;
1.77  0.1520
HW   HW  0.0  0.0   A   0.0e+00   0.0e+00 ;
0.00 0.

[ moleculetype ]
  ; molname   nrexcl
  NA+CL-  1

[ atoms ]
  ;   DU->IP
  ;   IM->DU
  ; = IM->IP
  ; id_at type res nr  residu name at name  cg nr  charge
mass  typeB chargeB   massB
1   DU  1NA+  NA+   1  0 22.9898
IP 0 22.9898
2   IM 1 DU   DU2  0 35.4530
DU 0 35.4530



On Tue, 1 Dec 2015 13:17:38 +0000
Hannes Loeffler <hannes.loeff...@stfc.ac.uk> wrote:

> Dear Sebastian,
> 
> I have got your attachment.  I do not know all intricacies of Gromacs
> but for relative free energies you only use one moleculetype for the
> perturbed group, see attachment.  You do not use the couple-* switches
> for relative simulations.
> 
> Yes, you are right about typeA and typeB.  If you have larger
> molecules you would just combine both molecules in their entirety
> into one "pseudo" molecules.  One original molecules has the normal
> force field types in the A state and dummy atoms in the final, while
> the other molecule has this reversed.
> 
> 
> Cheers,
> Hannes.
> 
>  
> 
> On Tue, 1 Dec 2015 13:34:12 +0100
> Sebastian Stolzenberg <s.stolzenb...@fu-berlin.de> wrote:
> 
> > Dear Hannes,
> > 
> > > I have some trouble reading your top file in my web browser
> > Apologies, am not sure if this works in this list, but I am now also
> > attaching my topol.top and grompp.mdp files.
> > I hope this let's you see better what I am doing wrong.
> > 
> > > As your two dummy types appear to be identical
> > > you can define just one.
> > Thank you, as you can see in topol.top, I now have only one dummy
> > atom type called "XX".
> > 
> > > I think your transformations are IPX->Na+, IMX->Cl-.
> > my "typeA" columns are "XX" and "IM" (NA+)
> > my "typeB" columns are "IM" and "XX" (CL-)
> > If I understand correctly (do I?), lambda=0 corresponds to "typeA",
> > whereas lambda=1 corresponds to "typeB", so my transformation is
> > "XX" -> "Na+", "Cl-" -> "XX"
> > 
> > Thanks a bunch,
> > Sebastian
> > 
> > 
> > 
> > On 30.11.2015 17:30, hannes.loeff...@stfc.ac.uk wrote:
> > > Dear Sebastian,
> > > 
> > > I have some trouble reading your top file in my web browser but I
> > > think your transformations are IPX->Na+, IMX->Cl-.  So you would
> > > have to reverse the columns for the first atom.  As your two dummy
> > > types appear to be identical you can define just one.
> > > 
> > > Cheers,
> > > Hannes.
> > > 
> > > 
> > > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> > > [gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of
> > > Sebastian Stolzenberg [s.stolzenb...@fu-berlin.de] Sent: 30
> > > November 2015 16:19 To: gromacs.org_gmx-users@maillist.sys.kth.se
> > > Subject: Re: [gmx-users] tutorial for "dual" alchemical
> > > transformation in gromacs?
> > > 
> > > Dear Hannes,
> > > 
> > > Thank you very much for your answer.
> > > 
> > >> Practical question: why would you want to do that?
> > >> Your two molecules will drift apart unless parts of the molecules
> > >> are linked together in a "single" topology fashion or or you
> > >> introduce other restraints/constraints.
> > > I am doing some method development, so my plan is to currently fix
> > > all atoms of A and B.
> > > 
> > >> In practice, you can do that with Gromacs relatively easily.  The
> > >> manual describes how to set up the A and the B state of a
> > >> molecule.
> > > OK, thanks. As a start, I set up a simple "Cl->Na in a water box"
> > > system:
> > > 
> > > gro file:
> > > 1  Na+  NA+1   1.434   1.045   1.049
> > > 2  Cl-  CL-2   1.434   1.045   1.049
> > > 3  WATO3   2.914   1.241   1.997
> > > 3  WAT   H14   2.899   

Re: [gmx-users] tutorial for "dual" alchemical transformation in gromacs?

[Hannes Loeffler]

Dear Sebastian,

I have got your attachment.  I do not know all intricacies of Gromacs
but for relative free energies you only use one moleculetype for the
perturbed group, see attachment.  You do not use the couple-* switches
for relative simulations.

Yes, you are right about typeA and typeB.  If you have larger molecules
you would just combine both molecules in their entirety into one
"pseudo" molecules.  One original molecules has the normal force field
types in the A state and dummy atoms in the final, while the other
molecule has this reversed.


Cheers,
Hannes.

 

On Tue, 1 Dec 2015 13:34:12 +0100
Sebastian Stolzenberg  wrote:

> Dear Hannes,
> 
> > I have some trouble reading your top file in my web browser
> Apologies, am not sure if this works in this list, but I am now also
> attaching my topol.top and grompp.mdp files.
> I hope this let's you see better what I am doing wrong.
> 
> > As your two dummy types appear to be identical
> > you can define just one.
> Thank you, as you can see in topol.top, I now have only one dummy atom
> type called "XX".
> 
> > I think your transformations are IPX->Na+, IMX->Cl-.
> my "typeA" columns are "XX" and "IM" (NA+)
> my "typeB" columns are "IM" and "XX" (CL-)
> If I understand correctly (do I?), lambda=0 corresponds to "typeA",
> whereas lambda=1 corresponds to "typeB", so my transformation is
> "XX" -> "Na+", "Cl-" -> "XX"
> 
> Thanks a bunch,
> Sebastian
> 
> 
> 
> On 30.11.2015 17:30, hannes.loeff...@stfc.ac.uk wrote:
> > Dear Sebastian,
> > 
> > I have some trouble reading your top file in my web browser but I
> > think your transformations are IPX->Na+, IMX->Cl-.  So you would
> > have to reverse the columns for the first atom.  As your two dummy
> > types appear to be identical you can define just one.
> > 
> > Cheers,
> > Hannes.
> > 
> > 
> > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se
> > [gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of
> > Sebastian Stolzenberg [s.stolzenb...@fu-berlin.de] Sent: 30
> > November 2015 16:19 To: gromacs.org_gmx-users@maillist.sys.kth.se
> > Subject: Re: [gmx-users] tutorial for "dual" alchemical
> > transformation in gromacs?
> > 
> > Dear Hannes,
> > 
> > Thank you very much for your answer.
> > 
> >> Practical question: why would you want to do that?
> >> Your two molecules will drift apart unless parts of the molecules
> >> are linked together in a "single" topology fashion or or you
> >> introduce other restraints/constraints.
> > I am doing some method development, so my plan is to currently fix
> > all atoms of A and B.
> > 
> >> In practice, you can do that with Gromacs relatively easily.  The
> >> manual describes how to set up the A and the B state of a molecule.
> > OK, thanks. As a start, I set up a simple "Cl->Na in a water box"
> > system:
> > 
> > gro file:
> > 1  Na+  NA+1   1.434   1.045   1.049
> > 2  Cl-  CL-2   1.434   1.045   1.049
> > 3  WATO3   2.914   1.241   1.997
> > 3  WAT   H14   2.899   1.335   1.984
> > 3  WAT   H25   2.841   1.212   2.052
> > ...
> > 
> > i.e., I transform a "Cl-" into a "Na+" by de-coupling
> > (annihiliating) the non-bonded interactions of "Cl-" and coupling
> > (letting appear) the ones for "Na+". This I attempt by defining
> > "dummy" atom types IMX and IPX for Cl- and Na+, respectively:
> > 
> > top file:
> > ...
> > [ atomtypes ]
> > ;name   bond_type mass charge   ptype   sigma
> > epsilon Amb
> >  IP   IP  0.0  0.0   A 3.32840e-01
> > 1.15897e-02 ; 1.87  0.0028
> >  IM   IM  0.0  0.0   A 4.40104e-01
> > 4.18400e-01 ; 2.47  0.1000
> >  IPX  IP  0.0  0.0   A 0.0e+00
> > 0.0e+00 ; 0.00  0.
> >  IMX  IM  0.0  0.0   A 0.0e+00
> > 0.0e+00 ; 0.00  0.
> > ...
> > [ atoms ]
> >   ; id_at type res nr  residu name at name  cg nr  charge
> > mass typeB chargeB   massB
> > 1   IPX 1  NA+ NA+   1  0
> > 22.9898   IP 0 22.9898
> > ...
> > [ atoms ]
> >   ; id_at type res nr  residu name at name  cg nr  charge
> > mass typeB chargeB   massB
> > 1   IM  1 CL-   CL-  1  0
> > 35.45300  IMX0 35.45300
> > 
> > (Please note that for now, I artificially set the charge of each
> > ion to zero to avoid total charge differences for different lambda
> > values. Also, for "Cl->Na", this strategy is certainly more
> > complicated, but in principle allows for different numbers of atoms
> > between molecule A and B in the "A->B" transformation.)
> > Unfortunately, running grompp with the standard run.mdp file from
> > http://www.gromacs.org/@api/deki/files/196/=fe-tutorial-4.6.pdf
> > (I omitted "couple-lambda0" and "couple-lambda1"
> > and set "couple-moltype   = NA+ CL-")
> > I am puzzled by the following error 

Re: [gmx-users] support mmcif format?

[Hannes Loeffler]

On Mon, 23 Nov 2015 09:57:55 +
Mark Abraham  wrote:

> I think we are unlikely to plan to write mmCIF directly as a
> trajectory format,

It's hard to see for me why someone would even want to use mmCIF as
a trajectory format.  The format serves, primarily, the crystallography
community to store a single set of coordinates together with a huge
amount of meta data, most of which are probably rather irrelevant for
trajectories.   (For archiving purposes this is of course interesting
but no need to write this into every single trajectory, just like the
topology is kept separate).

I think for the simulation community the mmCIF format should be what the
PDB format should always have been: a read-only file format to be
converted into a format developers think is best for their particular
piece of software.  Any future setup program should be able to read
mmCIF but several libraries are already available to do exactly that.

Cheers,
Hannes.
-- 
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] solvation free energy, electrostatic transformation

[Hannes Loeffler]

What happens when you do both transformations simultaneously, i.e.
sc-coul=yes?


On Mon, 27 Jul 2015 15:55:04 +0200
Daniele Veclani danielevecl...@gmail.com wrote:

 I'm trying to calculate the  solvation free energy of a molecule (M).
 
 I have done a VdW.  transformation.
 
 I have done also a  electrostatic transformation. But I get an
 incorrect value (large overestimation of DG solvation).
 
 I use gromacs 5.0.4 with gromos 54a7 force fields and my .mdp file
 (for lambda 00) is:
 integrator   = sd
 tinit= 0
 dt   = 0.002
 nsteps   = 250
 comm_mode= angular
 nstcomm  = 100
 nstxout  = 500
 nstvout  = 500
 nstfout  = 0
 nstlog   = 500
 nstenergy= 500
 nstxout-compressed   = 0
 cutoff-scheme= group
 nstlist  = 0
 ns_type  = simple
 pbc  = no
 rlist= 0
 coulombtype  = cutoff
 rcoulomb = 0
 epsilon_r= 1
 vdwtype  = cutoff
 rvdw = 0
 DispCorr  = no
 fourierspacing   = 0.12
 pme_order= 6
 ewald_rtol   = 1e-06
 epsilon_surface  = 0
 tcoupl   = berendsen
 tc_grps  = system
 tau_t= 0.1
 ref_t= 300
 Pcoupl   = no
 tau_p= 1.0
 compressibility  = 4.5e-05
 ref_p= 1.0
 free_energy  = yes
 init_lambda_state= 0
 delta_lambda = 0
 calc_lambda_neighbors= 1
 vdw_lambdas  = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 1.00 1.00 1.00
 coul_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70
 0.80 0.90 1.00
 bonded_lambdas   = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 restraint_lambdas= 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 temperature_lambdas  = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 sc-alpha = 0.5
 sc-coul  = no
 sc-power = 1.0
 sc-sigma = 0.3
 couple-moltype   = QUI
 couple-lambda0   = vdw-q
 
 couple-lambda1   = vdw
 couple-intramol  = no
 nstdhdl  = 10
 gen_vel  = no
 constraints  = all-bonds  ; we only have C-H bonds here
 constraint-algorithm = lincs
 continuation = no
 lincs-order  = 12
 
 
 What's wrong?
 
 Best regards
 D.V.

-- 
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] solvation free energy

[Hannes Loeffler]

Any particular reason why you want to do it the hard way?  I suspect
the answer to this is somewhere in section 5.3 where [pairs_nb] is
discussed.  You would have to explicitly set A and B states and not
use couple-*. You may find some old-style input files on
http://www.alchemistry.org/wiki/Main_Page .

On Fri, 17 Jul 2015 12:34:42 +0200
Daniele Veclani danielevecl...@gmail.com wrote:

 I read the manual and also some examples of  free energy calculation.
 
 when I do the simulation:
 M(vacuo) -- nothing(vacuo)
 
 I leave only intramolecular interactions.
 
 How can I just leave these interactions? I don't find the answer in
 the manual.
 
 
 Best regards
 D.V.
 
 2015-07-17 10:10 GMT+02:00 hannes.loeff...@stfc.ac.uk:
 
  I do not know how your top file looks like and what you have done
  to it.
 
  The setup procedure to create the topolog file would be as for a
  standard MD simulation with M.  You would make sure that the QUI
  label referes to the right indexes in your index file and run TI
  with you preferred protocol.  This simplified setup procedure
  ensures that you do not have to modify the topology file in any way
  esp. not to have to set up an explicit B state.
 
  I suggest to read all parts relating to couple-* in the manual.
  This is scattered over several places in the text.
 
  Cheers,
  Hannes.
 
  
  From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se [
  gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of
  Daniele Veclani [danielevecl...@gmail.com]
  Sent: 17 July 2015 08:57
  To: gmx-us...@gromacs.org
  Subject: Re: [gmx-users] solvation free energy
 
  Thank you for your answer.
 
  I read the manual. As I understand it, I should be deleted,
  from .top file, the
  section of non-bonded parameters. It's right?
 
  Of course, after the VdW I'll change the electrostatic
  transformation.
 
  Best regards
  D.V.
 
  2015-07-16 16:55 GMT+02:00 Hannes Loeffler
  hannes.loeff...@stfc.ac.uk:
 
   The couple-* parameters take already care of including the
   non-bonded terms internal to your molecule to correctly describe
   the transfer of M to vacuum.  That's the point of those
   parameters so that you would not have to run an additional
   correction in vacuo.  See the discussion in the manual (it's
   section 5.3 for Gromacs 4.6.x).  You would also need to compute
   the electrostatic transformation.
  
   On Thu, 16 Jul 2015 16:26:44 +0200
   Daniele Veclani danielevecl...@gmail.com wrote:
  
I'm trying to calculate the  solvation free energy of a
molecule (M). I have done:
M+water --- dum+water
   
Now I have to do:
M(vacuo) -- dum(vacuo)
   
In this case I have a problem, in fact I find a DG = 0.0 and
within the .xvg file there are only zeros.
   
Where is the problem?
   
I use gromacs 5.0.4 and my .mdp file (for lambda 00) is:
integrator   = sd
tinit= 0
dt   = 0.002
nsteps   = 250
comm_mode= angular
nstcomm  = 100
nstxout  = 500
nstvout  = 500
nstfout  = 0
nstlog   = 500
nstenergy= 500
nstxout-compressed   = 0
cutoff-scheme= group
nstlist  = 0
ns_type  = simple
pbc  = no
rlist= 0
coulombtype  = cutoff
rcoulomb = 0
epsilon_r= 1
vdwtype  = cutoff
rvdw = 0
DispCorr  = no
fourierspacing   = 0.12
pme_order= 6
ewald_rtol   = 1e-06
epsilon_surface  = 0
tcoupl   = berendsen
tc_grps  = system
tau_t= 0.1
ref_t= 300
Pcoupl   = no
tau_p= 1.0
compressibility  = 4.5e-05
ref_p= 1.0
free_energy  = yes
init_lambda_state= 0
delta_lambda = 0
calc_lambda_neighbors= 1
vdw_lambdas  = 0.00 0.10 0.20 0.30 0.40 0.50 0.60
0.70 0.80 0.90 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
bonded_lambdas   = 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
restraint_lambdas= 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
temperature_lambdas  = 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
sc-alpha = 0.5
sc-coul  = no
sc-power = 1.0
sc-sigma = 0.3
couple-moltype   = QUI

Re: [gmx-users] solvation free energy

[Hannes Loeffler]

The couple-* parameters take already care of including the non-bonded
terms internal to your molecule to correctly describe the transfer of M
to vacuum.  That's the point of those parameters so that you would not
have to run an additional correction in vacuo.  See the discussion in
the manual (it's section 5.3 for Gromacs 4.6.x).  You would also need to
compute the electrostatic transformation.

On Thu, 16 Jul 2015 16:26:44 +0200
Daniele Veclani danielevecl...@gmail.com wrote:

 I'm trying to calculate the  solvation free energy of a molecule (M).
 I have done:
 M+water --- dum+water
 
 Now I have to do:
 M(vacuo) -- dum(vacuo)
 
 In this case I have a problem, in fact I find a DG = 0.0 and
 within the .xvg file there are only zeros.
 
 Where is the problem?
 
 I use gromacs 5.0.4 and my .mdp file (for lambda 00) is:
 integrator   = sd
 tinit= 0
 dt   = 0.002
 nsteps   = 250
 comm_mode= angular
 nstcomm  = 100
 nstxout  = 500
 nstvout  = 500
 nstfout  = 0
 nstlog   = 500
 nstenergy= 500
 nstxout-compressed   = 0
 cutoff-scheme= group
 nstlist  = 0
 ns_type  = simple
 pbc  = no
 rlist= 0
 coulombtype  = cutoff
 rcoulomb = 0
 epsilon_r= 1
 vdwtype  = cutoff
 rvdw = 0
 DispCorr  = no
 fourierspacing   = 0.12
 pme_order= 6
 ewald_rtol   = 1e-06
 epsilon_surface  = 0
 tcoupl   = berendsen
 tc_grps  = system
 tau_t= 0.1
 ref_t= 300
 Pcoupl   = no
 tau_p= 1.0
 compressibility  = 4.5e-05
 ref_p= 1.0
 free_energy  = yes
 init_lambda_state= 0
 delta_lambda = 0
 calc_lambda_neighbors= 1
 vdw_lambdas  = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70
 0.80 0.90 1.00
 coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 bonded_lambdas   = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 restraint_lambdas= 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 temperature_lambdas  = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00
 sc-alpha = 0.5
 sc-coul  = no
 sc-power = 1.0
 sc-sigma = 0.3
 couple-moltype   = QUI
 couple-lambda0   = vdw
 
 couple-lambda1   = none
 couple-intramol  = no
 nstdhdl  = 10
 gen_vel  = no
 constraints  = all-bonds  ; we only have C-H bonds here
 constraint-algorithm = lincs
 continuation = no
 lincs-order  = 12
 
 Best regards
 Daniele Veclani.

-- 
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

[Hannes Loeffler]

On Wed, 15 Jul 2015 16:20:25 +0200
Julian Zachmann frankjulian.zachm...@uab.cat wrote:

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

To add to what has been said already:  increasing the output within a
fixed set of time won't give you any more information.  Your data is
probably highly correlated.  You may want to look into a tool like
https://github.com/mobleylab/alchemical-analysis if you haven't done
already and have a look how much data survives the correlation analysis
implemented in there.

Also, check for convergence to see if you have still strongly varying
gradients/energies.  Something as simple as a block analysis can be
very effective.


Cheers,
Hannes.
-- 
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] [pairs] func type 2

[Hannes Loeffler]

Hi,

do I interpret table 5.5 (Gromacs 4.6.7) correctly that function type 2
in section [pairs] can be used to realize mixed 1-4 scaling as it is
necessary with mixing e.g. an AMBER protein force field with the GLYCAM
force field?  GLYCAM doesn't scale 1-4 interactions.

Thanks,
Hannes.
-- 
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

[Hannes Loeffler]

Hi,

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.

One minor point is that your duplication of atomtypes (with i and m
prefixes) seems pretty redundant to me.

Cheers,
Hannes.
-- 
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

[Hannes Loeffler]

On Tue, 19 May 2015 17:48:17 +0200
Julian Zachmann frankjulian.zachm...@uab.cat wrote:

 Concerning the changes of the atom types: the changes in the charges
 are really small and the results will probably be the same if I would
 not convert the whole ligand, just the 4 atoms which are changing.
 However, as i had both ligands parametrised, I thought it would be a
 good idea to really convert the whole ligand from one to the other,
 and change all charges on all atoms.

That's not what I meant.  I was wondering why you define a, say, ic3
type and an mc3 type.  They are identical!  The charges are explicitly
listed for both states so I think the atom types only affect the vdW
terms.

Cheers,
Hannes.
-- 
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

[Hannes Loeffler]

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.

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

[Hannes Loeffler]

On Tue, 19 May 2015 12:03:41 -0400
Michael Shirts mrshi...@gmail.com wrote:

 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.

This reminds me: there was recently a discussion around couple-moltype
on this list.  My understanding is that this can be used for absolute
free energies.  But what I really want to make sure is this: as far as
I can see in the code and I think the manual says that too is that when
couple-moltype is not set in the .mdp none of the other couple-
parameters have any effect should I set any of the to something.  Do I
get this right?

Many thanks for the answers,
Hannes.



 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] difference of calculation time between the numbers of fep-lambdas values

[Hannes Loeffler]

On Mon, 27 Apr 2015 07:43:51 -0400
Justin Lemkul jalem...@vt.edu wrote:

 Note that without couple-moltype, you're going to be decoupling the
 whole system, which is (1) not what you want for calculating
 solvation free energy and (2) extremely slow.  There is an inherent
 slowdown when running the free energy code, but it should not be so
 large.

Really? I just looked into the code yesterday and I think that if you do
not set couple-moltype, none of the coupling code is ever triggered.  I
have only played a bit with relative free energies (Gromacs 4.6) where I
think you don't need any of the couple- parameters at all.

Cheers,
Hannes.
-- 
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] difference of calculation time between the numbers of fep-lambdas values

[Hannes Loeffler]

On Mon, 27 Apr 2015 09:05:05 -0400
Justin Lemkul jalem...@vt.edu wrote:

 
 
 On 4/27/15 9:02 AM, Hannes Loeffler wrote:
  On Mon, 27 Apr 2015 07:43:51 -0400
  Justin Lemkul jalem...@vt.edu wrote:
 
  Note that without couple-moltype, you're going to be decoupling the
  whole system, which is (1) not what you want for calculating
  solvation free energy and (2) extremely slow.  There is an inherent
  slowdown when running the free energy code, but it should not be so
  large.
 
  Really? I just looked into the code yesterday and I think that if
  you do not set couple-moltype, none of the coupling code is ever
  triggered.  I have only played a bit with relative free energies
  (Gromacs 4.6) where I think you don't need any of the couple-
  parameters at all.
 
 
 Maybe I got that backwards.  For relative free energies, this is true
 (we're doing these now); topological differences are all that are
 needed.  For an absolute solvation free energy, one needs
 couple-moltype unless the B-state (dummy) is explicitly defined in
 the topology, AFAIK.  If there is no couple-moltype, and no B-state
 defined in the topology, I'm not sure what the code would even be
 doing, actually.

I think it runs a normal MD simulation for each lambda.  A quick test
with a A-state only topology and no couple- parameters gives only
gradients with zero value.  In fact, the output file says

There are 0 atoms and 0 charges for free energy perturbation


Cheers,
Hannes.

-- 
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] Announcement: release of FESetup 1.1

[Hannes Loeffler]

Dear Gromacs community,

We are pleased to announce the release of FESetup 1.1.

FESetup is a tool to automate setup of alchemical relative free energy
simulations like thermodynamic integration (TI) and can also be used
for general simulation setup (equilibration). The tool will
automatically parameterise ligands (AM1/BCC) and map all atoms of each
mutational pair (for codes supporting the single topology paradigm).
Supported molecular simulation packages implementing free energy
simulation are currently GROMACS, AMBER and Sire (all these codes are
hybrid single/dual topology).

General simulation setup through an abstract MD engine is available for
AMBER, GROMACS, NAMD and DL_POLY. Supported force fields are all modern
AMBER force fields including GAFF. Future plans include extending the
code to support other popular biomolecular simulation software,
additional force fields and parameterisation schemes. We particularly
aim at automatisation where it makes sense and is possible, ease of use
and robustness of the code.

Please find the software at
http://www.hecbiosim.ac.uk/repo/download/2-software/3-fesetup

Kind regards,
Hannes Loeffler.
-- 
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.