[gmx-users] BAR and non-bonded tables revisited
rustadjr at corning.com
Wed May 23 17:24:11 CEST 2012
Thanks much for your previous help. I now have the calculation running.
However I get suspicious results that make me think that BAR may not be
integrated with tables correctly.
In this problem I am trying to change an Mg ion into an Fe ion (really its a
ferrous iron to a ferric iron, but I’m calling the ferrous iron “Mg” for
Here are the relevant fragments of my current 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
Mg 55.30500 1.20000 A 0.0 0.0
Fe 55.30500 1.80000 A 0.0 0.0
;(The tables “table_O_Fe.xvg” and “table_O_Mg.xvg” exist and are read in
[ nonbond_params ] ; these now give the multipliers for the tables columns
; 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 Mg 1 7.542 1.92800E-10
O Fe 1 40.425524240 1.92800E-10
. . .
Now I have defined a molecule type for the changing ion labeled “Cr” in the
.gro file, which transforms from atom type Mg to atom type Fe
[ moleculetype ]
; name nrexcl
[ atoms ]
; nr type resnr residue atom cgnr charge mass TypeB ChargeB
1 Mg 1 Mg Mg 2 1.2 55.3050 Fe 1.8
. . .
And the relevant portion of my .mdp file for lambda=0
coulombtype = PME-Switch
rcoulomb = 1.0
vdw_type = User
. . .
free-energy = yes
couple-moltype = Cr
init-lambda = 0.0
foreign-lambda = 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
couple-lambda0 = vdw-q
couple-lambda1 = vdw-q
It runs fine, but when I look at the RDF for the changing ion (Cr), it does
not go smoothly between Mg and Fe as lambda goes from 0 to 1, but, at
lambda=1, the peak position is pushed way further out than either Mg (~0.2
nm) and Fe (~0.19 nm) (it’s not an aqueous system b.t.w.) and sits at 0.23
nm or more.
At lambda=0, everything looks OK (i.e. the "Cr" ion gives the same g(r) as
the "Mg" ion)
So it makes me wonder if the short-ranged interactions in the tables are
simply being added together without regard for what the lambda value is?
(i.e. at lambda =1, we still get the entire short range interaction for the
"A" state added to that of the "B" state).
Is it possible that there is a bug in the integration of non-bonded tables
into the free energy capability?
Or do I misunderstand how the lambda is applied in this "mutating" type free
Thanks for all the help already, and for any further insight anyone can
provide. I can show a picture if the RDFs if you want to see them.
James R. Rustad, Ph.D.
Research Associate, Corning, Inc.
Professor Emeritus, UC Davis
One Science Center Drive
SP TD 01-1
Corning, NY 14831
View this message in context: http://gromacs.5086.n6.nabble.com/BAR-and-non-bonded-tables-revisited-tp4997686.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
More information about the gromacs.org_gmx-users