[gmx-users] tutorial for "dual" alchemical transformation in gromacs?
Sebastian Stolzenberg
s.stolzenberg at fu-berlin.de
Tue Dec 1 13:34:16 CET 2015
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.loeffler at 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-bounces at maillist.sys.kth.se [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Sebastian Stolzenberg [s.stolzenberg at fu-berlin.de]
> Sent: 30 November 2015 16:19
> To: gromacs.org_gmx-users at 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 WAT O 3 2.914 1.241 1.997
> 3 WAT H1 4 2.899 1.335 1.984
> 3 WAT H2 5 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.00000 0.00000 A 3.32840e-01 1.15897e-02
> ; 1.87 0.0028
> IM IM 0.00000 0.00000 A 4.40104e-01 4.18400e-01
> ; 2.47 0.1000
> IPX IP 0.00000 0.00000 A 0.00000e+00 0.00000e+00
> ; 0.00 0.0000
> IMX IM 0.00000 0.00000 A 0.00000e+00 0.00000e+00
> ; 0.00 0.0000
> ...
> [ 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 IMX 0 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 message:
> "The lambda=0 and lambda=1 states for coupling are identical"
>
> Is it possible from this error message to see what I am doing wrong?
>
> Many Thanks,
> Sebastian
>
>
>
>
> On 29.11.2015 14:36, hannes.loeffler at stfc.ac.uk wrote:
>> 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.
>>
>> 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.
>>
>> ________________________________________
>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Sebastian Stolzenberg [s.stolzenberg at fu-berlin.de]
>> Sent: 29 November 2015 13:11
>> To: gromacs.org_gmx-users at maillist.sys.kth.se
>> Subject: [gmx-users] tutorial for "dual" alchemical transformation in gromacs?
>>
>> Dear Gromacs experts,
>>
>> I am looking for a Gromacs tutorial to perform a "dual" alchemical
>> transformation in a water box, i.e. a mutation or any other alchemical
>> transformation, in which:
>> the nonbonded interactions between the environment and a molecule A
>> disappear as lambda goes from 0 to 1, and
>> the nonbonded interactions between the environment and a molecule B
>> appears simultaneously.
>> Interactions between atoms A and B are to be excluded.
>>
>> So far, I have been trying to use the Gromacs user guide to extend the
>> ethanol tutorial
>> http://www.gromacs.org/@api/deki/files/196/=fe-tutorial-4.6.pdf
>> into such a "dual" transformation, but also because I am currently faced
>> with manually manipulating my topology file, this option is highly prone
>> to my own errors.
>>
>> Would you be able to help?
>>
>> With Thanks and Best Regards,
>> Sebastian
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
>>
>
>
> --
> ----------------------------------------------
> Sebastian Stolzenberg, Ph.D.
> PostDoc, Computational Molecular Biology group
> Freie Universität Berlin
>
> Phone: (+49) (0)30 838 75153
> Mail: Arnimallee 6, 14195 Berlin, Germany
> ----------------------------------------------
>
>
--
----------------------------------------------
Sebastian Stolzenberg, Ph.D.
PostDoc, Computational Molecular Biology group
Freie Universität Berlin
Phone: (+49) (0)30 838 75153
Mail: Arnimallee 6, 14195 Berlin, Germany
----------------------------------------------
-------------- next part --------------
; in.amber_GMX.top created by acpype (Rev: 8810) on Mon Nov 30 15:20:20 2015
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.8333
[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon Amb
IP IP 0.00000 0.00000 A 3.32840e-01 1.15897e-02 ; 1.87 0.0028
IM IM 0.00000 0.00000 A 4.40104e-01 4.18400e-01 ; 2.47 0.1000
XX IP 0.00000 0.00000 A 0.00000e+00 0.00000e+00 ; 0.00 0.0000
OW OW 0.00000 0.00000 A 3.15075e-01 6.35968e-01 ; 1.77 0.1520
HW HW 0.00000 0.00000 A 0.00000e+00 0.00000e+00 ; 0.00 0.0000
[ moleculetype ]
; molname nrexcl
NA+ 1
[ atoms ]
; id_ at type res nr residu name at name cg nr charge mass typeB chargeB massB
1 XX 1 NA+ NA+ 1 0 22.9898 IP 0 22.9898
[ moleculetype ]
; molname nrexcl
CL- 1
[ 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 XX 0 35.45300
[ moleculetype ]
; molname nrexcl ; TIP3P model
WAT 2
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 OW 1 WAT O 1 -0.834 16.00000
2 HW 1 WAT H1 1 0.417 1.00800
3 HW 1 WAT H2 1 0.417 1.00800
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 462750.4 0.09572 462750.4
1 3 1 0.09572 462750.4 0.09572 462750.4
[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.520 836.800 104.520 836.800
#else
[ settles ]
; i j funct length
1 1 0.09572 0.15139
[ exclusions ]
1 2 3
2 1 3
3 1 2
#endif
[ system ]
in.amber
[ molecules ]
; Compound nmols
NA+ 1
CL- 1
WAT 414
-------------- next part --------------
; we'll use the sd integrator with 100000 time steps (200ps)
integrator = sd
nsteps = 100000
dt = 0.002
nstenergy = 1000
nstlog = 5000
; turn off trajectory writing
nstxout = 0
nstvout = 0
; We use the old group scheme because as of writing, the Verlet scheme
; does not support free energy calculations with coupled molecules.
cutoff-scheme = group
; cut-offs at 1.0nm
rlist = 1.0
dispcorr = EnerPres
vdw-type = cut-off
rvdw = 1.0
; Coulomb interactions
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.12
; Constraints
constraints = all-bonds
; set temperature to 300K
tcoupl = v-rescale
tc-grps = system
tau-t = 0.2
ref-t = 300
; set pressure to 1 bar with a thermostat that gives a correct
; thermodynamic ensemble
pcoupl = parrinello-rahman
ref-p = 1
compressibility = 4.5e-5
tau-p = 5
; energygrp-excl = Cys Ser
; freezegrps = Cys Ser
; freezedim = Y Y Y Y Y Y
; and set the free energy parameters
free-energy = yes
couple-moltype = NA+ CL-
; these 'soft-core' parameters make sure we never get overlapping
; charges as lambda goes to 0
sc-power = 1
sc-sigma = 0.3
sc-alpha = 1.0
; we still want the molecule to interact with itself at lambda=0
couple-intramol = yes
; couple-lambda1 = vdwq
; couple-lambda0 = vdwq
init-lambda-state = 0
; These are the lambda states at which we simulate
; for separate LJ and Coulomb decoupling, use
fep-lambdas = 0.0 0.2 0.4 0.6 0.8 0.9 1.0
More information about the gromacs.org_gmx-users
mailing list