[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 +0000
Stefania Evoli <stefania.evoli at kaust.edu.sa> wrote:

> The [atom] section of the single topology of my transformation looks
> like
> 
> [ atoms ]
> ;   nr       type  resnr    res   atom   cgnr     charge       mass
> typeB    chargeB      massB comments
>      1          c      1    LIG     C1      1   0.932601   12.01000
>   c   0.935601   12.01000 ; MCSS
>      2         c3      1    LIG    CA1      2  -0.095200   12.01000
>  c3  -0.126500   12.01000 ; MCSS
>      3   c3_dummy      1    LIG    CB1      3    0.00000   12.01000
>  c3   0.136400   12.01000 ; to be appeared
>      4         hx      1    LIG    HA1      4   0.088700    1.00800
> hx_dummy    0.00000    1.00800 ; to be annihilated
>      5         hx      1    LIG    HA2      5   0.088700    1.00800
>  hx   0.094700    1.00800 ; MCSS
>      6   h1_dummy      1    LIG    HB1      6    0.00000    1.00800
>  h1   0.068200    1.00800 ; to be appeared
>      7   h1_dummy      1    LIG    HB2      7    0.00000    1.00800
>  h1   0.068200    1.00800 ; to be appeared
>      8   ho_dummy      1    LIG    HG1      8    0.00000    1.00800
>  ho   0.448000    1.00800 ; to be appeared
>      9         n4      1    LIG     N1      9  -0.835601   14.01000
>  n4  -0.822601   14.01000 ; MCSS
>     10          o      1    LIG     O1     10  -0.755301   16.00000
>   o  -0.759801   16.00000 ; MCSS
>     11   oh_dummy      1    LIG    OG1     11    0.00000   16.00000
>  oh  -0.628801   16.00000 ; to be appeared
>     12          o      1    LIG   OXT1     12  -0.755301   16.00000
>   o  -0.759801   16.00000 ; MCSS
>     13         hn      1    LIG     H1     13   0.443800    1.00800
>  hn   0.448800    1.00800 ; MCSS
>     14         hn      1    LIG     H2     14   0.443800    1.00800
>  hn   0.448800    1.00800 ; MCSS
>     15         hn      1    LIG     H3     15   0.443800    1.00800
>  hn   0.448800    1.00800 ; MCSS