amber-developers: Serious problems with Antechamber runs due to faulty Mopac in AmberTools?

From: Ross Walker <ross.rosswalker.co.uk>
Date: Thu, 10 Jul 2008 20:40:14 -0700

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.





Received on Sun Jul 13 2008 - 06:07:40 PDT
Custom Search