Custom Search
---
> 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<<
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Sep 25 2010 - 06:00:04 PDT