[gmx-users] alchemical free energy changes with OPLS
Alexandre Suman de Araujo
asaraujo at if.sc.usp.br
Tue Jul 12 12:38:58 CEST 2005
Are you using single or double precision?
In single precision the noise is general high generating this kind of
problem.
--
Alexandre Suman de Araujo
asaraujo at if.sc.usp.br
UIN: 6194055
IFSC - USP - São Carlos - Brasil
Brian Stephenson wrote:
> Hi,
>
> I've recently been computing free energies of hydration by converting
> all-atoms in a molecule to non-interacting dummy atoms in water and
> vacuo. DeltaG hydration then equals the difference between this
> free-energy
> change in vacuo and this free-energy change in water. I have been
> getting good
> results (in comparison with experiment) for my hydration free energy
> calculations using dummy atoms which I define in my .top file.
>
> However, I've been trying to set up some alchemical simulations lately
> and these
> aren't working as well. As you can see in the attached topology file, I'm
> morphing ethane into ethane, which should (obviously) have a net free
> energy
> change of zero. I'm trying to mimic a study published by David
> Pearlman ("A
> Comparison of Alternative Approaches to Free Energy Calculations" J. Phys.
> Chem. 1994, 1487-1493). I start with a propane molecule and use pdb2gmx
> to get a .top file. Then I change some atoms to be dummies such that
> at lambda
> = 0 and 1, only ethane is present additional dummy atoms attached.
> The dummy
> groups switch from one side of the molecule to the other over the
> course of the
> simulation.
>
> I've been using windows of lambda to estimate dG/dlambda for lambda
> between 0
> and 1, but unfortunately, my results for dG/dlambda are only roughly,
> and not
> exactly, symmetric. Also the results are very noisy, even for
> simulation times
> MUCH longer than reported by Pearlman in his paper (a factor of 1000
> longer).
> Instead of getting a free energy change of zero I get results between
> -30 and
> 30 kJ/mol (depending on implementation details and simulation time).
> I get similar
> problems in vacuo and water, although transformations in both
> environments should
> have a deltaG of zero.
>
> Is the way I am implementing dummy atoms with OPLS correct? Because I was
> getting OK results for my hydration free energy calculations I thought
> it would
> be, but maybe there is something I'm doing incorrectly for alchemical
> changes?
>
> Is there something else I am missing?
>
> Thanks very much for your help,
> Brian
>
>
>;
>; File
>'ethane_al.top' was generated
>; By user:
>bstephenson (503)
>; At date: Sun
>Jul 10 16:46:37 2005
>;
>; This is your
>topology file
>; "Shake
>Yourself" (YES)
>;
>; Include forcefield parameters
>#include "ffoplsaa.itp"
>
>[ atomtypes ]
>;type
>mass
>charge ptype
>c6
>c12
>DUMH 1.008
>0.0000 A
>0.0 0.0
>
>[ bondtypes ]
>HC
>DUMH 1 0.070 284512.0
>; H to dummy H using CT to H bond distance and potential
>
>[ angletypes ]
> DUMH HC
>DUMH 1
>107.800 276.144 ; CHARMM 22 parameter file
>for HC CT HC used for dumh h dumh interaction
> CT
>HC
>DUMH 1
>110.700 313.800 ; CHARMM 22 parameter file
>for CT CT HC used for c h dumh interaction
>
>[ dihedraltypes ]
> DUMH
>HC CT HC 3
>0.00000 0.00000 0.00000
>0.000000 0.00000 0.00000 ; hydrocarbon *new*
>11/99
>
>CT CT HC DUMH
>3 0.00000 0.00000
>0.00000 0.000000 0.00000 0.00000 ;
>hydrocarbon all-atom
>
>
>[ moleculetype ]
>; Name
>nrexcl
>ETH
>3
>
>[ atoms ]
>; nr type resnr
>residue atom cgnr
>charge mass
>typeB chargeB massB
> 1
>opls_135 1
>ETH CH1
>1 -0.18
>12.011
>opls_140 0.06 1.008;
>qtot -0.18
> 2
>opls_140 1
>ETH OH11
>1 0.06
>1.008
>DUMH 0.0 1.008;
>qtot -0.12
> 3
>opls_140 1
>ETH OH12
>1 0.06
>1.008
>DUMH 0.0 1.008;
>qtot -0.06
> 4
>opls_140 1
>ETH OH13
>1 0.06
>1.008
>DUMH 0.0 1.008;
>qtot 0
> 5
>opls_135 1
>ETH CAA
>1 -0.18
>12.011 ; qtot -0.18
> 6
>opls_140 1
>ETH HAA1
>1 0.06
>1.008 ; qtot -0.12
> 7
>opls_140 1
>ETH HAA2
>1 0.06
>1.008 ; qtot -0.06
> 8
>opls_140 1
>ETH CH2
>1 0.06
>1.008 opls_135 -0.18 12.011; qtot 0
> 9
>DUMH
>1 ETH OH21
>1
>0 1.008
>opls_140 0.06 1.008;
>qtot 0
> 10
>DUMH
>1 ETH OH22
>1
>0 1.008
>opls_140 0.06 1.008;
>qtot 0
> 11
>DUMH
>1 ETH OH23
>1
>0 1.008
>opls_140 0.06 1.008;
>qtot 0
>
>[ bonds ]
>; ai aj
>funct
>c0
>c1
>c2 c3
> 1 2 1
> 1 3 1
> 1 4 1
> 1 5 1
> 5 6 1
> 5 7 1
> 5 8 1
> 8 9 1
> 8 10 1
> 8 11 1
>
>[ pairs ]
>; ai aj
>funct
>c0
>c1
>c2 c3
> 1 9 1
> 1 10 1
> 1 11 1
> 2 6 1
> 2 7 1
> 2 8 1
> 3 6 1
> 3 7 1
> 3 8 1
> 4 6 1
> 4 7 1
> 4 8 1
> 6 9 1
> 6 10 1
> 6 11 1
> 7 9 1
> 7 10 1
> 7 11 1
>
>[ angles ]
>; ai aj ak
>funct
>c0
>c1
>c2 c3
> 2 1
>3 1
> 2 1
>4 1
> 2 1
>5 1
> 3 1
>4 1
> 3 1
>5 1
> 4 1
>5 1
> 1 5
>6 1
> 1 5
>7 1
> 1 5
>8 1
> 6 5
>7 1
> 6 5
>8 1
> 7 5
>8 1
> 5 8
>9 1
> 5 8
>10 1
> 5 8
>11 1
> 9 8
>10 1
> 9 8
>11 1
> 10 8
>11 1
>
>[ dihedrals ]
>; ai aj ak al
>funct
>c0
>c1
>c2
>c3
>c4 c5
> 2 1
>5 6 3
> 2 1
>5 7 3
> 2 1
>5 8 3
> 3 1
>5 6 3
> 3 1
>5 7 3
> 3 1
>5 8 3
> 4 1
>5 6 3
> 4 1
>5 7 3
> 4 1
>5 8 3
> 1 5
>8 9 3
> 1 5
>8 10 3
> 1 5
>8 11 3
> 6 5
>8 9 3
> 6 5
>8 10 3
> 6 5
>8 11 3
> 7 5
>8 9 3
> 7 5
>8 10 3
> 7 5
>8 11 3
>
>; Include Position restraint file
>#ifdef POSRES
>#include "posre.itp"
>#endif
>
>; Include water topology
>#include "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 generic topology for ions
>#include "ions.itp"
>
>[ system ]
>; Name
>ETH
>
>[ molecules ]
>; Compound #mols
>ETH
>1
>
>------------------------------------------------------------------------
>
>_______________________________________________
>gmx-users mailing list
>gmx-users at gromacs.org
>http://www.gromacs.org/mailman/listinfo/gmx-users
>Please don't post (un)subscribe requests to the list. Use the
>www interface or send it to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list