- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

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.

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:

This is not an Amber thing, it is a general concept in the simulation of

periodic systems. Just google "minimum image convention".

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/amberReceived on Wed Oct 01 2014 - 08:30:02 PDT

Custom Search