Hi All,
I am currently trying to update some of the tutorials for amber 10 /
ambertools. One of these is the Antechamber tutorial for sustiva
http://amber.scripps.edu/tutorials/basic/tutorial4/index.htm
There appear to be some serious issues with the use of mopac in the charge
derivation. My installation is good and all the test cases pass including
the sustiva test case. However if we take a close look at the mopac output
generated in this test case (the same system that is actually used in the
tutorial) the mopac output doesn't seem to make any sense.
amber10/test/antechamber/sustiva/mopac.out
CYCLE: 1 TIME: 0.18 TIME LEFT: 3599.7 GRAD.: 355.525 HEAT:-107.1790
CYCLE: 2 TIME: 0.18 TIME LEFT: 3599.5 GRAD.: 145.809 HEAT:-112.8673
CYCLE: 3 TIME: 0.17 TIME LEFT: 3599.3 GRAD.: 152.735 HEAT:-114.4029
CYCLE: 4 TIME: 0.23 TIME LEFT: 3599.1 GRAD.: 109.340 HEAT:-116.7827
CYCLE: 5 TIME: 0.16 TIME LEFT: 3598.9 GRAD.: 97.314 HEAT:-117.0437
CYCLE: 6 TIME: 0.16 TIME LEFT: 3598.8 GRAD.: 68.803 HEAT:-117.2610
CYCLE: 7 TIME: 0.18 TIME LEFT: 3598.6 GRAD.: 58.979 HEAT:-117.3911
CYCLE: 8 TIME: 0.15 TIME LEFT: 3598.4 GRAD.: 42.097 HEAT:-117.5702
CYCLE: 9 TIME: 0.18 TIME LEFT: 3598.3 GRAD.: 38.033 HEAT:-117.6394
HEAT OF FORMATION TEST SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
CYCLE: 10 TIME: 0.10 TIME LEFT: 3598.2 GRAD.: 54.640 HEAT:-117.6395
CYCLE: 11 TIME: 0.17 TIME LEFT: 3598.0 GRAD.: 56.588 HEAT:-117.8282
CYCLE: 12 TIME: 0.25 TIME LEFT: 3597.7 GRAD.: 101.377 HEAT:-117.9484
CYCLE: 13 TIME: 0.24 TIME LEFT: 3597.5 GRAD.: 59.564 HEAT:-118.2088
CYCLE: 14 TIME: 0.26 TIME LEFT: 3597.2 GRAD.: 77.089 HEAT:-118.2228
CYCLE: 15 TIME: 0.16 TIME LEFT: 3597.1 GRAD.: 54.505 HEAT:-118.3628
CYCLE: 16 TIME: 0.16 TIME LEFT: 3596.9 GRAD.: 51.136 HEAT:-118.4075
CYCLE: 17 TIME: 0.16 TIME LEFT: 3596.8 GRAD.: 42.986 HEAT:-118.4901
CYCLE: 18 TIME: 0.09 TIME LEFT: 3596.7 GRAD.: 56.844 HEAT:-118.7318
CYCLE: 19 TIME: 0.24 TIME LEFT: 3596.4 GRAD.: 63.560 HEAT:-118.7464
TEST ON X SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
CYCLE: 20 TIME: 0.24 TIME LEFT: 3596.2 GRAD.: 70.835 HEAT:-118.7517
CYCLE: 21 TIME: 0.25 TIME LEFT: 3596.0 GRAD.: 48.988 HEAT:-118.8432
CYCLE: 22 TIME: 0.10 TIME LEFT: 3595.8 GRAD.: 47.115 HEAT:-118.8828
CYCLE: 23 TIME: 0.08 TIME LEFT: 3595.8 GRAD.: 21.254 HEAT:-118.9406
CYCLE: 24 TIME: 0.24 TIME LEFT: 3595.5 GRAD.: 21.108 HEAT:-118.9553
CYCLE: 25 TIME: 0.08 TIME LEFT: 3595.4 GRAD.: 29.882 HEAT:-118.2161
CYCLE: 26 TIME: 0.08 TIME LEFT: 3595.4 GRAD.: 27.310 HEAT:-114.3819
CYCLE: 27 TIME: 0.08 TIME LEFT: 3595.3 GRAD.: 36.062 HEAT:-115.9819
CYCLE: 28 TIME: 0.08 TIME LEFT: 3595.2 GRAD.: 37.861 HEAT:-108.4102
CYCLE: 29 TIME: 0.08 TIME LEFT: 3595.1 GRAD.: 53.848 HEAT:-105.8221
CYCLE: 30 TIME: 0.19 TIME LEFT: 3594.9 GRAD.: 29.006 HEAT:-48.31796
CYCLE: 31 TIME: 0.17 TIME LEFT: 3594.8 GRAD.: 38.843 HEAT: 25.54461
CYCLE: 32 TIME: 0.17 TIME LEFT: 3594.6 GRAD.: 45.101 HEAT:-86.86664
CYCLE: 33 TIME: 0.16 TIME LEFT: 3594.4 GRAD.: 30.893 HEAT:-71.26195
CYCLE: 34 TIME: 0.19 TIME LEFT: 3594.2 GRAD.: 44.271 HEAT:-73.27885
CYCLE: 35 TIME: 0.16 TIME LEFT: 3594.1 GRAD.: 53.310 HEAT:-81.22255
HEAT OF FORMATION IS ESSENTIALLY STATIONARY
Well this last line is a blatant lie - how is the heat of formation here
essentially stationary? It changed by almost 8 Kcal/mol on the last step
along. This of course makes my statement in the tutorial that one should
check that the QM calculation worked successfully somewhat confusing...
Also taking a look at the final heat of formation that is reported:
FINAL HEAT OF FORMATION = -41.13867 KCAL
This is very different to the divcon energy from Amber 9:
HEAT OF FORMATION = -119.65798436 KCAL/MOL
In fact the mopac.out.save file looks pretty good:
CYCLE: 37 TIME: 0.11 TIME LEFT: 3595.5 GRAD.: 11.335 HEAT:-119.3189
TEST ON X SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
CYCLE: 38 TIME: 0.19 TIME LEFT: 3595.3 GRAD.: 15.148 HEAT:-119.3208
TEST ON X SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
CYCLE: 39 TIME: 0.14 TIME LEFT: 3595.2 GRAD.: 16.906 HEAT:-119.3210
HEAT OF FORMATION TEST SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
CYCLE: 40 TIME: 0.17 TIME LEFT: 3595.0 GRAD.: 20.318 HEAT:-119.3218
HEAT OF FORMATION TEST SATISFIED
HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN 0.20
With a heat of formation of -119.3218 which is very close to the values one
used to get from divcon in amber 9. However, the sustiva test does not
actually check if the mopac.out file matches the mopac.out.save file.
So it looks initially like something is very wrong in mopac here. The
structure also comes out quite a bit different as well:
Mopac.out:
1 C 0.0000 0.0000 0.0000
2 H 1.1033 0.0000 0.0000
3 C -0.6913 1.2084 0.0000
4 C 0.0207 2.5361 0.0263
Mopac.out.save
1 C 0.0000 0.0000 0.0000
2 H 1.1040 0.0000 0.0000
3 C -0.6541 1.2344 0.0000
4 C 0.1141 2.5253 0.0110
So the structures are obviously quite different. And the resulting AM1/BCC
charges produced by Mopac are understandably different:
Mopac.out
ATOM NO. TYPE CHARGE ATOM ELECTRON DENSITY
1 C -0.0465 4.0465
2 H 0.1670 0.8330
3 C -0.1667 4.1667
4 C 0.2767 3.7233
5 C -0.2160 4.2160
6 C -0.0109 4.0109
7 C -0.0971 4.0971
8 H 0.1463 0.8537
Mopac.out.save
ATOM NO. TYPE CHARGE ATOM ELECTRON DENSITY
1 C -0.0460 4.0460
2 H 0.1658 0.8342
3 C -0.1734 4.1734
4 C 0.2691 3.7309
5 C -0.2163 4.2163
6 C -0.0115 4.0115
7 C -0.0945 4.0945
8 H 0.1443 0.8557
However, the resulting sustiva.mol2 file is identical to the
sustiva.mol2.save file. How can this be?
Upon close inspection it appears that the charges in the mol2.save file
match those from the mopac.out file I get and not the mopac.out.save file.
Thus it seems that the mol2.save file was produced by just blindly running
antechamber with this copy of mopac and not checking the actual mopac.out
file. This is evident since the mopac.out.save file was not updated when the
sustiva.mol2.save file was updated.
This is evident by looking at the cvs logs for mopac.out.save:
revision 10.0
date: 2008/04/15 23:31:14; author: case; state: Exp; lines: +0 -0
bring tags to 10.0
----------------------------
revision 9.1
date: 2007/12/15 22:42:14; author: case; state: Exp; lines: +400 -400
updates from Junmei
----------------------------
And for sustiva.mol2.save:
revision 10.0
date: 2008/04/15 23:31:14; author: case; state: Exp; lines: +0 -0
bring tags to 10.0
----------------------------
revision 9.4
date: 2008/02/07 19:44:45; author: junmei_wang; state: Exp; lines: +29
-29
update
Junmei
----------------------------
revision 9.3
date: 2007/12/15 22:42:14; author: case; state: Exp; lines: +4 -4
updates from Junmei
----------------------------
revision 9.2
date: 2007/06/26 18:59:03; author: case; state: Exp; lines: +33 -33
change to mopac charges
----------------------------
This is very worrying since it means ALL the charges being generated by
antechamber with AmberTools are incorrect and incorrect to the point that
the test cases are bogus since they contain the incorrect values.
For example if we look at the charges for this test case for amber8, amber9
and ambertools1.0 we get
Amber 8 (mopac) Amber 9 (divcon) AmberTools 1.0 (mopac)
1 C1 -0.04968 -0.0496 -0.046500
2 H1 0.16937 0.1694 0.167000
3 C2 -0.17826 -0.1785 -0.174000
4 C3 0.31502 0.3150 0.328900
5 C4 -0.19720 -0.1967 -0.189100
6 C5 0.01368 0.0135 0.016000
7 C6 -0.07873 -0.0787 -0.084700
8 H2 0.10860 0.1086 0.107000
Thus the ambertools 1.0 charges are clearly incorrect.
One end user has definitely already come across this problem:
http://archive.ambermd.org/200806/0419.html
Only he was fobbed off being told that the difference was likely not
important and to be expected between different AM1 implementations. However
looking at the message closely you can see that the difference is way higher
than what one would expect from using say a working mopac vs divcon for
example.
So comments people? This looks pretty serious... Plus I have to give a
workshop in London in a week and I'm not sure what to do here - as it stands
right now I can't see an option except to have them have both amber 9 and
amber 10 / ambertools installed and for the antechamber tutorial tell them
they have to use amber 9 since otherwise I'll end up updating the tutorial
with incorrect values and get loads of questions relating to the mopac
calculation not actually having converged... Unless we can rush out a bugfix
in time.
Comments?
/\
\/
|\oss Walker
| Assistant Research Professor |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross.rosswalker.co.uk |
|
http://www.rosswalker.co.uk | PGP Key available on request |
Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.
- application/octet-stream attachment: mopac.out
Received on Sun Jul 13 2008 - 06:07:40 PDT
