Dear Jason,
Thank you for your answer...
>> If one loads the same number of _improper_ dihedrals (FF parameters)
>> in LEaP using either the generic (i.e. improper defined with 'X'
>> letters) or the specific 'way' (no 'X' letter) LEaP does not
>> rigorously apply the same impropers...
>>
>> I take an example:
>> - I first executed leaprc.ff99SB (so there are generic improper
>> dihedrals in parm99.dat) and saved the prmtop/prmcrd files for
>> ACE-TYR-NME: 10 impropers are generated.
>>
>> - Then I loaded my own frcmod file with only specific impropers; see below:
>> CA-CA-C -OH 1.10000000e+00 180.0 2. taken from parm99.dat
>> CA-CA-CA-CT 1.10000000e+00 180.0 2. taken from parm99.dat
>> C -CT-N -H 1.10000000e+00 180.0 2. taken from parm99.dat
>> CT-N -C -O 1.05000000e+01 180.0 2. adapted from
>> parm99.dat
>> CA-CA-CA-HA 1.10000000e+00 180.0 2. adapted from
>> parm99.dat
>> C -CA-CA-HA 1.10000000e+00 180.0 2. adapted from
>> parm99.dat
>> and saved the prmtop/prmcrd files for the same ACE-TYR-NME dipeptides
>>
>> I compare the generated prmtop files:
>>
>> diff generic.top specific.top
>> 230c230
>> < 39 48 -42 -45 18 54 42 -48 -51 18
>> ---
>> > 39 48 -42 -45 18 42 54 -48 -51 18
>> -> One observes differences in the %FLAG DIHEDRALS_INC_HYDROGEN
>> section of the prmtop file
>>
>> Then, I used 'rdparm' to study these 10 impropers; saved the 'rdparm'
>> output in two text files, that I compared:
>>
>> diff generic.txt specific.txt
>> 63c63
>> < B 62: 1.100 3.14 2.0 :1.CZ :1.CD1 :1.CE1 :1.HE1 (19,15,17,18)
>> ---
>> > B 62: 1.100 3.14 2.0 :1.CD1 :1.CZ :1.CE1 :1.HE1 (15,19,17,18)
>>
>> A single CA-CA-CA-HA improper is different; however, there are 4
>> CA-CA-CA-HA impropers in the Tyrosine residue; this means 3 are
>> applied using an homogeneous way, and for one the two first atoms are
>> inverted.
>>
>> Obviously I understand that the impact on the energy value is
>> negligible, but this is quite painful when one tries to generate
>> identical prmtop files:
>>
>> Any idea where/how to correct that in the LEaP source code?
>
> I had noticed that myself not too long ago. In 'real' impropers, the only
> atom whose position 'matters' is the third (central) atom -- the order of
> the others is arbitrary. However, since we treat impropers the same way as
> we treat proper dihedrals, the overall order will impact the energy a
> little bit. Therefore, the convention is applied to order the atoms in
> alphabetical order by type.
Yes, for all; however, in case of ambiguity (the case I reported) the
alphabetic order by atom name has to be added besides/after_? the
alphabetic order by atom type.
I think this was correctly implemented in the (old) PLEP program.
> However, there is no rigorously correct way of sorting 2 identical types
> (CA in your case). As a result, there is no well-defined ordering in this
> particular case (note that one of the CAs will always be the third atom, it
> is only the other 2 that can swap).
I also found other cases I did not report... This means well
designed/agreed empirical rules have to be added to select these 3
improper atoms (1st, 2nd & 4th ones).
> A more rigorous way of fixing the ordering is to order based on alphabetic
> order of types (to maintain 'classical' behavior of all force fields), then
> order degenerate types based on their atom names. While it is fairly
> arbitrary, it will at least yield reliably consistent topology files. And
> as you mentioned, improper ordering will have a minor effect on energies,
> but an even smaller effect on dynamics, I'd imagine.
I 100% agree ;-)
> A good place to start may be the zParmSetOrderImproperAtoms function
> (parmSet.c).
ok
regards, Francois
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jun 06 2013 - 00:00:03 PDT
