[gmx-users] tutorial for "dual" alchemical transformation in gromacs?
Sebastian Stolzenberg
s.stolzenberg at fu-berlin.de
Tue Dec 1 15:39:57 CET 2015
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.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+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 1 NA+ NA+ 1 0 22.9898
> IP 0 22.9898
> 2 IM 1 DU DU 2 0 35.4530
> DU 0 35.4530
>
>
>
> On Tue, 1 Dec 2015 13:17:38 +0000
> Hannes Loeffler <Hannes.Loeffler at 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.stolzenberg at 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.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 --------------
; 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
; and set the free energy parameters
free-energy = yes
couple-moltype = NC
; 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
-------------- 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
DU DU 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
NC 1
[ atoms ]
; id_ at type res nr residu name at name cg nr charge mass typeB chargeB massB
1 DU 1 NA+ NA+ 1 0 22.9898 IP 0 22.9898
2 IM 1 CL- CL- 2 0 35.45300 DU 0 35.45300
[ exclusions ]
1 2
[ 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
NC 1
WAT 414
More information about the gromacs.org_gmx-users
mailing list