[gmx-users] BAR calculations with tabulated non-bonded potentials

[Justin A. Lemkul]

On 5/18/12 4:21 PM, Rustad, James R wrote:
> Dear GMX users,
>
> I cannot figure out from the manual whether it is possible to use tabulated
> non-bonded potentials in a BAR calculation.
>
> I want to estimate the free energy for changing Fe2+ to Fe3+. These have
> different charges, of course, and also different non-bonded interactions. These
> non-bonded interactions are neither LJ nor Buckingham, so I am using a lookup
> table (table_O_Fe2.xvg and table_O_Fe3.xvg) for the “vdw” part (vdw_type=User).
>
> I think I should be able to define Fe2 and Fe3 atoms types, and associate the
> tables (table_O_Fe2.xvg and table_O_Fe3.xvg) with those types in the mdp file:
>
> energygrps = O Si Na Ca Al Fe2 Fe3
>
> energygrp_table = O Al O Ca O Fe2 O Fe3 O Na O O O Si Si Al Si Ca Si Fe2 Si Fe3
> Si Na Si Si Na Al Na Ca Na Fe2 Na Fe3 Na Na Ca Al Ca Ca Ca Fe2 Ca Fe3 Al Al Al
> Fe2 Al Fe3 Fe2 Fe2 Fe2 Fe3 Fe3 Fe3
>
> and in the .top file:
>
> [ atomtypes ]
>
> ;name mass charge ptype A B
>
> O 15.99940 -1.20000 A 0.0 0.0
>
> Na 22.99000 0.60000 A 0.0 0.0
>
> Si 28.08000 2.40000 A 0.0 0.0
>
> Ca 40.07800 1.20000 A 0.0 0.0
>
> Al 26.98200 1.80000 A 0.0 0.0
>
> Fe2 55.30500 1.20000 A 0.0 0.0
>
> Fe3 55.30500 1.80000 A 0.0 0.0
>
> [ nonbond_params ]
>
> ; i j func C6 C12
>
> O O 1 4.0904959800 2.12268E-09
>
> O Si 1 32.858468400 9.64853E-11
>
> O Na 1 2.2541869930 4.82427E-10
>
> O Ca 1 2.9149186000 4.82427E-10
>
> O Al 1 34.887265640 8.68368E-11
>
> O Fe3 1 40.425524240 1.92800E-10
>
> O Fe2 1 7.542 1.92800E-10
>
> ...
>
> Also in the .top file, I define a moleculetype Fe, which is Fe2 in state A and
> Fe3 in state B:
>
> [ moleculetype ]
>
> ; name nrexcl
>
> Fe 0
>
> [ atoms ]
>
> ; nr type resnr residue atom cgnr charge mass typeB chargeB massB
>
> 1 Fe2 1 Fe Fe 2 1.2 55.3050 Fe3 1.8 55.305
>
> I expect that I could specify this atom in the gro file like this:
>
> ...
>
> 1764Na Na 1764 3.448 3.335 0.801 0.9202 0.1009 -0.6229
>
> 1765Fe Fe 1765 1.123 1.300 0.562 0.3254 -1.0409 -0.6704
>
> 1766Ca Ca 1766 2.536 3.320 1.379 -0.5737 -0.4492 -1.8149
>
> ...
>
> The problem is that grompp fails, telling me
>
> “
>
> Fatal error:
>
> Group Fe2 not found in indexfile
>
> “
>
> So, while I want to associate the tables with atom types, it seems GMX is
> expecting that I am associating the tables (through the energy groups) with a
> particular atom, which, in my case, moves between two atom types in the BAR
> calculation.
>
> Can anyone tell me whether I am just suffering from a misconception about how
> the tables (or energy groups) are used, or whether it is in fact just not
> possible to do BAR (or other free energy) calculations in this way?
>

Group names must be valid [moleculetype] names or whatever is listed in a 
suitable index file.  You've used the name Fe2 as a group, when in fact it is an 
atom type.  If you want to set this as an energygrp, then the appropriate name 
is Fe, which is the name of the [moleculetype] shown in the topology above.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================