[gmx-users] BAR and non-bonded tables revisited

[jrustad]

All

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
convenience

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
without problems)


[ 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
Cr     0
[ atoms ]

; nr   type   resnr  residue   atom   cgnr   charge  mass     TypeB  ChargeB 
MassB
1       Mg     1     Mg        Mg      2       1.2   55.3050   Fe    1.8      
55.305

 . . .

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
energy problem?

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

Corning Incorporated
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.