On Wed, Oct 15, 2014 at 10:54 PM, Arun Srikanth <askforarun.gmail.com>
wrote:
> That was an awesome explanation for the impropers (I should have thought
> about bonds angles taking precedence) thanks. This fact has not been
> adressed properly anywhere atleast as far I have searched. I couldn't find
> this in cornell's paper or GAFF (junmei) or may be it is obvious for people
> who have experience in MD.
>
> I rely on parmed.py extensively for my simulations, often to check the
> dihedral connectivity information and also for exclusions.
>
> Your statment :
>
> "However, since the
> atom types defined in that torsion *may* be used elsewhere in your
> molecule, you need to make sure that the nonbonded scaling factors for that
> parameter is "correct".
>
I've thought of a clearer way of saying this. The scaling factors are
defined in the parameter files (parm10.dat and frcmod files), with the
default values of 1.2 and 2.0 (EEL and vdW) assigned for any torsion whose
scaling factors are not defined.
The only way you can know for *certain* that the end-group interactions
for a particular torsion term will *always* be excluded from the 1-4
nonbonded terms is if it is one of the first N-1 terms of an N-term
torsion. Only then can you assign scaling factors of 0 in the parameter
file. [1] The torsions you mentioned previously that was excluded (by
virtue of being in a ring) was a single, normal torsion term defined
between 4 atom types, so it needs to assign the "normal" scaling factors
(normal being whatever is normal for that force field). So even if a
torsion is only used in ring systems where those interactions are excluded,
you don't know that for certain when you write the force field parameter
database, so you can't set the scaling factors to 0.
To take into consideration the above statement, in my simulation I have the
> number of dihedral types equal to number of dihedrals that is each
> dihedral has its own type (as defined by parmed.py).
Note that tleap does not do this by default (nor do I think there is an
option). There is a DIHEDRAL_PHASE, DIHEDRAL_PERIODICITY, and
DIHEDRAL_FORCE_CONSTANT array (and SCXX_SCALE_FACTOR
when the scale factors are assigned). Each dihedral listed in
DIHEDRALS_INC_HYDROGEN and DIHEDRALS_WITHOUT_HYDROGEN contain 5 numbers --
the atom indexes and a pointer into the 3 (or 5) arrays I previously
listed. In my experience, the combined lengths of the 2 dihedral arrays
has never equalled the lengths of the 3 (or 5) parameter arrays (although
that's not impossible).
Although I am afraid
> that the end group interactions (1-4 non bonded) of a multiterm dihedral
> (atoms with the same type ) may be counted twice ( I do not know yet how
> lammps calculates the 1-4 interaction for a multi term dihedral, if the
> multiterm dihedral is defined twice as two new type of dihedrals)
>
>
>
> Your opinion on my above imeplementation would be helpful
>
I
would hope that LAMMPS does it correctly. I would bet that they do,
since the other common biomolecular force fields work the same way (like
CHARMM, for instance). Note that using the dihedral list to build the 1-4
nonbonded pair list (that have to be treated *specially*) is an
implementation detail. Amber does it because it's convenient, and they
"decorate" the dihedral lists with negatives in order to make the code
easier to write. But you certainly don't *have* to do it this way. You
can just look at all of the bonds, and build the 1-4 nonbonded arrays from
unique collections of atoms that are separated by at least 3 bonds (note in
5-member rings, torsions are separated by 3 bonds in 1 direction and 2
bonds in another, so they don't get added to this list).
That is an equally valid, and definitely more general, way to define the
1-4 nonbonded pair list. Both will give you the same answer for the Amber
force field (again, unless you used ParmEd to mess around with the dihedral
or bond list).
What I hope you've gathered from my tl;dr answers so far is that many of
the file formats in Amber are intricately linked to implementation details
in the code or relics from a code that once was (for instance, why are all
atom indexes represented as multiples of 3 indexed from 0? why are there
negatives in torsion terms? ...etc.) The 1-4 pair list is derived from
the torsion list, which is why some things in the prmtop are the way they
are. This makes some things easy (like coding the Amber force field), but
can make other things more difficult (like supporting coarse-grained force
fields and even the CHARMM force field). Programs that were designed to
support multiple force fields and file types (like NAMD, OpenMM, maybe
LAMMPS ??) may make very different implementation choices than Amber,
because they were guided by different design principles and different
goals. [2]
HTH,
Jason
[1] This is actually done in tleap, but based on the same principle.
[2] I've found that reading Amber code helps decode certain file formats --
see
http://ambermd.org/doc/OFF_file_format.txt for example -- and vice
versa (e.g., I knew what to look for in the code based on how a file format
was laid out, like with the 1-4 scaling terms defined in the prmtop). In
many cases these file formats and some old code relics help tell the story
of how the Amber codes have developed.
--
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 15 2014 - 21:30:02 PDT
