Re: [AMBER] LIE FUNCTION

From: Jason Swails <jason.swails.gmail.com>
Date: Wed, 01 Oct 2014 11:17:30 -0400

Please keep responses on the mailing list.

On Wed, 2014-10-01 at 15:53 +0200, Loïc Etheve wrote:
> Hi,
>
> Could you explain what is rhe minimum image convention because i'm not
> so experienced in AMBER.

This is not an Amber thing, it is a general concept in the simulation of
periodic systems. Just google "minimum image convention".

> When I did calculations between 2 atoms closer (less than 5Å) results
> are not equal

You can check the code in AmberTools/src/cpptraj/src/Action_LIE.cpp
(specifically lines 141 to 183), which does a double-loop over all atoms
in the ligand mask and all atoms in the surroundings mask (which is
everything _except_ the ligand mask by default), computes the square of
their minimum image distance, and then accumulates the LJ interactions
as A*r6*r6-B*r6 (where r6 is 1/r^6).

The calculation you are doing by hand is not the same as the calculation
that cpptraj is doing -- that is the only explanation I have for why you
can't get the energies to match. Make sure you know _exactly_ what
atoms are selected in each mask in the LIE action and that the distance
you are using really is the minimum image distance.

Also make sure you are extracting the correct A and B coefficients from
the prmtop file. You need to use the atom type pointers for each,
compute a pointer into the NONBONDED_PARM_INDEX section, and use that
pointer as the index for the LENNARD_JONES_ACOEF and BCOEF sections
(indexes in the prmtop file begin from 1). Or save yourself some
headache and extract the LJ parameters using ParmEd and re-combine them
to get the A and B coefficients.

HTH,
Jason

P.S.

I just ran a quick calculation using the following input file:

parm sol.top
trajin test.mdcrd
distance :54.CH2 :143.CD out dist.dat
lie :54.CH2 :143.CD out lie.dat cutvdw 1000 cutelec 1000

in the directory $AMBERHOME/AmberTools/test/mmpbsa_py/EstRal_Files/.
This will specifically compute the LJ interaction between only the two
atoms selected. I extracted the LJ parameters from ParmEd using the
script:

parm sol.top
printDetails :54.CH2|:143.CD

and got the output:

The mask :54.CH2|:143.CD matches 2 atoms:

   ATOM RES RESNAME NAME TYPE LJ Radius LJ Depth Mass Charge GB Radius GB Screen
    827 54 TRP CH2 CA 1.9080 0.0860 12.0100 -0.1134 1.7000 0.7200
   2270 143 LYS CD CT 1.9080 0.1094 12.0100 -0.0479 1.7000 0.7200

Using the Lorentz-Berthelot combining rules, Rmin=2*1.9080 and
epsilon=sqrt(0.0860*0.1094). So the A coefficient is epsilon*Rmin^12 and
the B coefficient is 2*epsilon*Rmin^6. cpptraj computes the distance
for the first frame as 3.8995 A (the same value I get with VMD).
Plugging that into the formula A*r^12-B*r^6, I get a vdW energy of
-0.095558088 kcal/mol. The value in the lie.dat file (-0.0956) matches.

-- 
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Oct 01 2014 - 08:30:02 PDT
Custom Search