[gmx-users] Free energy calculation of multiple solutes

[Jeff Chen]

Thank you very much for the quick reply.

Without using tpbconv and mdrun -rerun, I'm not sure how to get the
solute-solute energies from gromacs (please let me know it there's a way to
do this) -- since there are many more solvent molecules, the total
potential energy is dominated by solvent-solvent and solvent-solute
interactions. However, the total potential energy from a simulation (not a
free energy simulation) with the solute charges set to zero is different
than that obtained from a free energy simulation with lambda=1 -- the
differences are fairly small, but big enough to be worrisome.

I'll run a system with fewer solvent molecules to get a better handle on
what's going on -- thanks again for your suggestions.

It seems most people do free energy calculations with a single solute
molecule -- in your opinion, is it very unlikely that there is any problem
in the way solute-solute interactions are handled in the free energy
simulations?

Thanks again,

Jeff


On Wed, Feb 19, 2014 at 7:34 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 2/19/14, 1:05 PM, Jeff Chen wrote:
>
>> Dear All,
>>
>> I am trying to calculate the change in free energy that occurs when
>> charging multiple type B molecules in an A-B mixture. For example, I have
>> a
>> mixture of 900 molecules of type A and 100 of type B. When coulomb
>> interactions are on, B is identical to A (something like TIP4P water), but
>> when coulomb interactions are off, each B is simply a Lennard-Jones
>> particle (in the topology file, I've labelled B differently and called the
>> molecule SOLB, but all interaction information is identical to A). To do
>> the free energy calculations, I perform different simulations at each
>> lambda point. For example, for a simulation with lambda=0, I have the
>> following in my .mdp file:
>>
>> free_energy = yes
>> init_lambda = 0.0
>> delta_lambda = 0
>> foreign_lambda= 0.05
>> couple-lambda0=vdw-q
>> couple-lambda1=vdw
>> couple-moltype=SOLB
>>
>> and when lamba=1,
>>
>> free_energy = yes
>> init_lambda = 1.0
>> delta_lambda = 0
>> foreign_lambda= 0.95
>> couple-lambda0=vdw-q
>> couple-lambda1=vdw
>> couple-moltype=SOLB
>>
>> I expected that when lambda=0, B molecules would behave identically to A
>> molecules with the coulombic energy fully turned on
>> (couple-lambda0=vdw-q),
>> and when lambda=1, B molecules would behave like LJ particles, with no
>> coulombic interactions, either with each other, or with A molecules.
>>
>> However, when I look at the coulombic energy for B-B interactions (I used
>> tpbconv -o to make a new *.tpr for only B molecules, then used mdrun
>> -rerun
>> with a trajectory containing only information about B molecules), I find
>> that both the short-range and reciprocal coulombic energies are NOT zero
>> as
>> I had expected for lambda = 1, and are in fact exactly equal to the
>> negative of the respective coulombic energies obtained using the same
>> trajectory data when lambda=0.
>>
>
> What are the energies without using tpbconv and mdrun -rerun?  In other
> words, it's necessary to figure out whether this is a problem in the
> original simulation with the free energy code (unlikely, given how simple
> the system is) or a problem using mdrun -rerun in conjunction with the free
> energy code (far more likely).
>
> -Justin
>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==================================================
> --
> 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.
>