Hello Alan,
This was indeed a very detailed description. One thing I have to comment on
-- the fact that you're using antechamber I think makes all the difference.
The problem is, for all intents and purposes, the CX and CT atom types are
completely identical chemically. They are sp3 carbons with 4 bonds, so
there's really no way for antechamber to distinguish between them.
Historically, there was never any reason to make a difference between these
2 atom types. However, it was found that the backbone torsion fixes were
affecting the sidechains in an undesired way, so the specific sp3 carbon at
the C-alpha position was re-typed to CX so that those backbone dihedrals
could be uniquely specified. It would be difficult for antechamber to
automatically insert CX for some sp3 carbons and CT for the others, I would
imagine.
The main point is this study changed too many things. You're not comparing
ff99SB vs. ff10, you're comparing vanilla ff99SB vs. antechamber-ed ff10.
To really compare the 2, you should run the same tleap scripts for ff99SB
and ff10 and omit the antechamber step, and see if those 2 results are
identical. They should be, even up to the non-bonded terms for 2 steps.
The problem is that you're assuming that antechamber will assign the same
atom types to the various atoms in the residues as the creators of ff10 did,
which I think is invalid. I'm also not sure it's the goal of antechamber to
do so, since it's designed primarily to find charge derivations and force
field parameters for small ligands that have not yet been parameterized.
Hope this helps,
Jason
On Sat, Sep 25, 2010 at 8:30 AM, Alan <alanwilter.gmail.com> wrote:
> First, I am not using ff10 for simulations. My interest is Antechamber (as
> it is the core of ACPYPE). I am testing ff10 to validate according my
> criteria (still working on a paper to be published) and doing so I hope to
> help improve AmberTools.
>
> One point of this exercise is to know how good is Antechamber in
> reproducing
> results for which we know the results.
>
> Ok, I will try to put in examples and then more questions will come along.
> However be prepared, it's long and although I did to the best of my
> attention, as careful and detailed as possible, I kindly ask you to pay
> very
> careful attention as well.
>
> Case 1: Using ff99SB/ff99bsc0 (they gave the same result)
>
> #------------------------------------------------------------------------
> # Generate input PDB using tleap and sander with ff99SB and parm99.dat
> cat << EOF >| mdin
> Minimization
> &cntrl
> imin=1, maxcyc=200, ntmin=2, ntb=0, igb=0,cut=999,
> //
> EOF
>
> cat << EOF >| leap.in
> source leaprc.ff99SB
> mol = sequence {NHIE HIE CHIE}
> savepdb mol HHH.pdb
> saveamberparm mol prmtop inpcrd
> quit
> EOF
>
> tleap -f leap.in
> sander -O
> ambpdb < restrt >| HHH.pdb
>
> # Single Point Energy for Reference
> cat << EOF >| leap.in
> source leaprc.ff99SB
> hhh = loadpdb HHH.pdb
> saveamberparm hhh prmtop inpcrd
> savepdb hhh hhh.pdb
> quit
> EOF
> tleap -f leap.in
>
> cat << EOF >| mdin
> Single point
> &cntrl
> imin=0, maxcyc=0,
> ntmin=2,
> ntb=0,
> igb=0,
> cut=999,
> //
> EOF
>
> sander -O
> mv mdinfo mdinfo_ref
>
> # Now using Antechamber
> babel -ipdb HHH.pdb -omol2 HHH.mol2
> antechamber -i HHH.mol2 -fi mol2 -o HHH_gas_amber.mol2 -fo mol2 -c gas -nc
> 0
> -m 1 -s 2 -df 2 -at amber -pf n
> parmchk -i HHH_gas_amber.mol2 -f mol2 -o HHH_AC.frcmod -p
> $AMBERHOME/dat/leap/parm/parm99.dat
>
> cat << EOF >| tleap.in
> source leaprc.ff99SB
> mods = loadamberparams HHH_AC.frcmod
> HIE = loadmol2 HHH_gas_amber.mol2
> check HIE
> saveamberparm HIE HHH_AC.prmtop HHH_AC.inpcrd
> quit
> EOF
>
> tleap -f tleap.in
>
> sander -O -i mdin -o mdout -p HHH_AC.prmtop -c HHH_AC.inpcrd
> diff mdinfo mdinfo_ref
>
> < Etot = 27.3750 EKtot = 0.0000 EPtot =
> 27.3750
> < BOND = 1.2774 ANGLE = 62.2606 DIHED =
> 17.9565
> < 1-4 NB = 7.6326 1-4 EEL = 43.3055 VDWAALS =
> -6.6478
> < EELEC = -98.4098 EHBOND = 0.0000 RESTRAINT =
> 0.0000
> ---
> > Etot = -17.0210 EKtot = 0.0000 EPtot =
> -17.0210
> > BOND = 1.2774 ANGLE = 62.2606 DIHED =
> 17.9163
> > 1-4 NB = 7.6326 1-4 EEL = 235.4812 VDWAALS =
> -6.6478
> > EELEC = -334.9414 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> Question 1: ignore non-bonded energies; why DIHED differs? Running Case 1
> with ALA and I have no differences for bonded energies (see Case 2).
> Using 'rdparm' to shed a light here:
>
> cat << EOF >| rdparm.in
> atoms
> bonds
> angles
> dihedrals
> pertbonds
> pertangles
> pertdihedrals
> printExcluded
> printLennardJones
> printTypes
> info
> quit
> EOF
>
> rdparm prmtop < rdparm.in > rdparm_ref.out
> rdparm HHH_AC.prmtop < rdparm.in > rdparm.out
>
> Using an text editor and some operations I found out that only these
> dihedrals differ:
>
> from rdparm.out:
> (10,17,16,14)
> (13,14,12,11)
> (27,34,33,31)
> (30,31,29,28)
> (44,51,50,48)
> (47,48,46,45)
>
> from rdparm_ref.out:
> (10,14,16,17)
> (11,14,12,13)
> (27,31,33,34)
> (28,31,29,30)
> (44,48,50,51)
> (45,48,46,47)
>
> But, AFAIK, those dihedral with the 2 central atoms equal are equivalent,
> i.e., (11,14,12,13) = (13,14,12,11).
> However (10,14,16,17), (27,31,33,34) and (44,48,50,51) seems to be wrong.
>
> To validate what I found above, since I have no idea of how fixing this
> supposed wrong dihe parameters in a prmtop file, I used then gromacs with
> ffamber (in this case using ff99sb) and acpype. Once calculating the
> dihedral energies, correcting the respective class of dihedrals from
> CG-HD2-CD2-NE2 to CG-NE2-CD2-HD2 (for the 3 residues in the gromacs top
> file) and the dihedral energy matched the one done for reference with
> tleap/sander.
>
> So I believe antechamber is generating this particular class of dihedrals
> wrongly.
>
> End of Case 1
>
> #------------------------------------------------------------------------
> Now Case 2: Running example above for ALA (AAA) and comparing ff99SB x ff10
> NOTE: every time I run an example I delete all the content of the folder,
> just in case.
>
> # AAA
> # Generate input PDB using tleap and sander with ff99SB and parm99.dat
> cat << EOF >| mdin
> Minimization
> &cntrl
> imin=1, maxcyc=200, ntmin=2, ntb=0, igb=0,cut=999,
> //
> EOF
>
> cat << EOF >| leap.in
> source leaprc.ff99SB
> mol = sequence {NALA ALA CALA}
> savepdb mol AAA.pdb
> saveamberparm mol prmtop inpcrd
> quit
> EOF
>
> tleap -f leap.in
> sander -O
> ambpdb < restrt >| AAA.pdb
>
> # Single Point Energy for Reference
> cat << EOF >| leap.in
> source leaprc.ff99SB
> aaa = loadpdb AAA.pdb
> saveamberparm aaa prmtop inpcrd
> savepdb aaa aaa.pdb
> quit
> EOF
> tleap -f leap.in
>
> cat << EOF >| mdin
> Single point
> &cntrl
> imin=0, maxcyc=0,
> ntmin=2,
> ntb=0,
> igb=0,
> cut=999,
> //
> EOF
>
> sander -O
> mv mdinfo mdinfo_ref
> cat mdinfo_ref
> # Etot = -18.5791 EKtot = 0.0000 EPtot =
> -18.5791
> # BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 13.4340
> # 1-4 NB = 7.9882 1-4 EEL = 232.7906 VDWAALS =
> -3.1222
> # EELEC = -273.3628 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> # Now using Antechamber
> babel -ipdb AAA.pdb -omol2 AAA.mol2
> antechamber -i AAA.mol2 -fi mol2 -o AAA_gas_amber.mol2 -fo mol2 -c gas -nc
> 0
> -m 1 -s 2 -df 2 -at amber -pf n
> parmchk -i AAA_gas_amber.mol2 -f mol2 -o AAA_AC.frcmod -p
> $AMBERHOME/dat/leap/parm/parm99.dat
>
> cat << EOF >| tleap.in
> source leaprc.ff99SB
> mods = loadamberparams AAA_AC.frcmod
> ALA = loadmol2 AAA_gas_amber.mol2
> check ALA
> saveamberparm ALA AAA_AC.prmtop AAA_AC.inpcrd
> quit
> EOF
>
> tleap -f tleap.in
>
> sander -O -i mdin -o mdout -p AAA_AC.prmtop -c AAA_AC.inpcrd
> diff mdinfo mdinfo_ref
> cat mdinfo
> # Etot = -27.0552 EKtot = 0.0000 EPtot =
> -27.0552
> # BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 13.4340
> # 1-4 NB = 7.9882 1-4 EEL = 54.2986 VDWAALS =
> -3.1222
> # EELEC = -103.3469 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> # only non-bonded energies diffs
> # Anyway, let's have rdparm results for further comparison against Case 2
> AAA with ff10
>
> cat << EOF >| rdparm.in
> atoms
> bonds
> angles
> dihedrals
> pertbonds
> pertangles
> pertdihedrals
> printExcluded
> printLennardJones
> printTypes
> info
> quit
> EOF
>
> rdparm prmtop < rdparm.in >| rdparm_ff99SB_ref.out
> rdparm AAA_AC.prmtop < rdparm.in >| rdparm_ff99SB_AC.out
>
> # Identical rdparm out (except for non-bonded)
>
> Now the same above but using ff10
>
> [replace ff99SB by ff10 and skipping to the more important results]
>
> sander -O
> mv mdinfo mdinfo_ff10_ref
> cat mdinfo_ff10_ref
> # Etot = -18.5791 EKtot = 0.0000 EPtot =
> -18.5791
> # BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 13.4340
> # 1-4 NB = 7.9882 1-4 EEL = 232.7906 VDWAALS =
> -3.1222
> # EELEC = -273.3628 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> # This result is absolutely identical to the one from AAA with ff99SB ref,
> great, since we expected that.
>
> # Now using Antechamber
> [snip]
> sander -O -i mdin -o mdout -p AAA_AC.prmtop -c AAA_AC.inpcrd
> diff mdinfo mdinfo_ff10_ref
>
> #< Etot = -30.3944 EKtot = 0.0000 EPtot =
> -30.3944
> #< BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 10.0948
> #< 1-4 NB = 7.9882 1-4 EEL = 54.2986 VDWAALS =
> -3.1222
> #< EELEC = -103.3469 EHBOND = 0.0000 RESTRAINT =
> 0.0000
> #---
> #> Etot = -18.5791 EKtot = 0.0000 EPtot =
> -18.5791
> #> BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 13.4340
> #> 1-4 NB = 7.9882 1-4 EEL = 232.7906 VDWAALS =
> -3.1222
> #> EELEC = -273.3628 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> # Ops! Ignoring non-bonded, only DIHED from AAA ff10 Antechamber differs
> from AAA ff99SB/ff10 ref
> # Could it be a similar problem as in Case 1, I doubt since we don't have
> dihedrals CG-NE2-CD2-HD2.
>
> Like before, using rdparm:
>
> rdparm prmtop < rdparm.in >| rdparm_ff10_ref.out
> rdparm AAA_AC.prmtop < rdparm.in >| rdparm_ff10_AC.out
>
> Comparing rdparm_ff10_ref.out x rdparm_ff10_AC.out, and I found in ref we
> have atom type CX instead of CT. Fine. But then, having that in mind I will
> find that there's 9 missing dihedrals in rdparm_ff10_AC.out and they are:
>
> (11,13,15,21) * 3 : Type C - N - CT - C
> (13,15,21,23) * 3 : Type N - CT - C - N
> (21,23,25,31) * 3 : Type C - N - CT - C
>
> So, AAA ff10_ref works fine because when creating AAA, it gets CX atom type
> instead of CT.
> However, when using Antechamber, atom type identified is CT and dihedrals
> involving C - N - CT - C and N - CT - C - N are not present in parm10.dat.
> If I amend parm10.dat (by including the missing C - N - CT - C and N - CT -
> C - N) and rerun Case 2 again:
>
> mdinfo_ff10_ref doesn't change: Good
>
> diff mdinfo mdinfo_ff10_ref
> # only non-bonded energies diffs: Great!
>
> cat mdinfo
> # Etot = -27.0552 EKtot = 0.0000 EPtot =
> -27.0552
> # BOND = 0.6320 ANGLE = 3.0611 DIHED =
> 13.4340
> # 1-4 NB = 7.9882 1-4 EEL = 54.2986 VDWAALS =
> -3.1222
> # EELEC = -103.3469 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> So, concluding, amending parm10.dat and then Antechamber generate a prmtop
> that really reproduce the same DIHED = 13.4340 seen in mdinfo_ff10_ref.
>
> Of course, another solution would be to make Antechamber smarter enough to
> recognise when to use CX instead of CT.
>
> End of Case 2
> #------------------------------------------------------------------------
>
> Bottom line:
> Case 1: Antechamber (with ff99SB/bsc0, not tested with ff10) is generating
> wrong dihedral parameters?
> Case 2: It seems to be missing dihed parameters C - N - CT - C and N - CT -
> C - N in parm10.dat
>
> I hope all this divagation may help in anyway, if you had the patience to
> read up here.
>
> Best regards,
>
> Alan
>
> On 24 September 2010 17:05, case <case.biomaps.rutgers.edu> wrote:
>
> > On Fri, Sep 24, 2010, Alan wrote:
> > >
> > > However, when running with ff10, I noticed that C -N -CT-C was getting
> > zero
> > > everywhere.
> >
> > Can you give an example? I would think that the C-N-CT-C dihedral would
> > not
> > be needed. This sounds like an error in one of the prep files. (And
> > another
> > reason not to have wild-cards, but that is a separate issue.)
> >
> > For the rest of the list:
> >
> > (a) Don't use ff10: there is a reason it is not documented.
> >
> > (b) If you do feel you absolutely have to use ff10, please report any
> ways
> > in
> > which it is different than ff99SB in the energies it produces. (Don't
> > report
> > diffs for RNA: we already know that ff10 for RNA differs from that in
> > ff99SB/ff99bsc0.)
> >
> > ...thanks...dac
> >
> >
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
>
>
>
> --
> Alan Wilter S. da Silva, D.Sc. - CCPN Research Associate
> Department of Biochemistry, University of Cambridge.
> 80 Tennis Court Road, Cambridge CB2 1GA, UK.
> >>http://www.bio.cam.ac.uk/~awd28 <http://www.bio.cam.ac.uk/%7Eawd28><<
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Graduate Student
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Sep 25 2010 - 10:30:05 PDT
