Optimization in MOPAC
Dear Netters:
Many thanks for the suggestions on handling PAHs in MOPAC. For the benefit
of anyone else who didn't get the direct replies to me, here is a list of
the replies.
Jack Houser
Department of Chemistry
The University of Akron
Akron, Ohio 44325-3601
-----------------------------
Try keyword XYZ. It works in most cases.
Ya-Jun Zheng
=========================================================================
I've encountered this problem with MOPAC ( although I run MOPAC 6 ),
occassionally and I've found that using an optimizer other than the
default (BFGS) usually allows the computation to proceed. I usually
use the eigen following method which is unfortunately not available
in Mopac 5, but you might try NLLSQ or one of the others.
Jim Schmidt
University of Rochester
Pharmacology Dept.
jims %-% at %-%
duce.medicine.rochester.edu
=========================================================================
This is a typical problem when you create ring structures via a Z matrix.
If the atoms in the ring are defined such that each is positioned at some
angle to the previous one, then the collective changes of them add, and
can couse large motions of the last one defines. Try using a dummy atom
at the center of the rings and defining the atoms likes spokes on a wheel
relative to some fixed point of reference (like another dummy atom).
This usually alievaiates the problem and speeds up geometry
optimizations, as the variable are less coupled.
_______________________________________________________________________
Steven K. Pollack |
University of Cincinnati |
Department of Materials Science | I notice that you are still
& Engineering | using polymers!
498 Rhodes Hall ML#12 | - Lt. Commander Montgomery Scott
Cincinnati, OH 45221-0012 | U.S.S. Enterprise
spollack %-% at %-% ester.mse.uc.edu |
________________________________________________________________________
=========================================================================
=========================================================================
I found the same problems you describe with MOPAC when
I worked up my buckminsterfullerene calculation. I solved them in
the following manner.
First, make certain you specify the carbons as aromatic carbons
in PCModel rather than alternating double and single bonds. The latter
actually gives alternating bond lengths even if you request the
pi calculation.
Secondly, give the key word XYZ. This specifies the calculation be
done in cartessian coordinates rather than with the Z-matrix.
In case it is not already obvious, the problem is the significant
difference between the PCModel input geometry and the final MOPAC
geometry. I think one could also solve this problem by asking MOPAC
to take smaller steps in the iteration but I did not figure out how
to do that. Good luck.
Jim Gano
=========================================================================
I used to have this problem alot, especially in systems with a lot of rings.
I've found that it was helpfull to choose a z-matrix in which there were a
minimum of dependencies. By that i mean that the def. of atom 'n' was
dependent on as few other atoms as reasonably possible. This was usually
done by choosing the first atom defined to be placed at a strategically
centralized location.
This didn't allways solve my problem, and i found that if i took the z-matrix
of the system immediatly prior to the crash, and used that for a new
".DAT"
file, things would get better, (but usually only after quite a few iterations
of this rather labor intensive procedure).
A few times none of these methods solved the problem. But in those
cases i discovered that the problem was with a bad starting position. I
now create all my initial conformations from a minimized structure using some
empirical force field approach.
Hope this is helpfull.
I've also noticed that using MOPAC6 with the EF option finds better minimums
faster. Unfortunaly i not sure i trust the new parameters in the newer AM1
...
i've recently heard some other people say they don't trust them, but they
gave no reason. (I don't think i believe some of the results i get in
the 3rd derivative of energy vs. distance/angle). If you happen to know
of any documented problems with MOPAC6/AM1 could you please drop me a line?
phil
Philip E. Klunzinger
------------------------------------------------------------------
email : philk %-% at %-% eby.polymer.uakron.edu phone: (216) 972-5810
usmail: Institute of Polymer Science FAX : (216) 972-5190
Univsrity of Akron,
Akron, OH 44325-3909
------------------------------------------------------------------
=========================================================================
The fixed point exception sounds like a real bug in your version, or
else the symmetry input is wrong. An easier way to solve the problem
is to do optimiztion in cartesian coordinates. Add the keyword XYZ to the
command line. For internal coordinates, it really is possible that
"SMALL CHANGES IN INTERNAL COORDINATES ARE CAUSING A LARGE CHANGE IN THE
DISTANCE". Cartesian optimizations are generally more stable in this
respect.
Another option is to get MOPAC 6 and use the EF geometry optimization
option. This is quite an efficient optimizer.
George Fitzgerald
Cray Research, Inc.
655E Lone Oak Dr.
Eagan, MN 55121
612-683-3676
=========================================================================
Have you tried optimizing in cartesians using the keyword XYZ.
I have often seen gradients "explode" when they are transformed from
cartesians to internals, which usually indicates an unfortunate choice
of internals. After the gradient has settled down in cartesian space,
the same z-matrix will often work again. Just an idea.
Regards, Jan Jensen
North Dakota State, Dept. of Chem.
Fargo, ND
=========================================================================
When I received your first error message, I have had some success
rerunning the job using the keyword XYZ, which runs the calculation
using cartesian coordinates instead of internal. The job usually runs
correctly, but it takes 2-4 times more cpu time. If you have input good
coordinates (bond lengths, angles) for the rigid aromatic portions of
your molecules, perhaps you could not optimize the lengths, angles,
and torsions for these atoms, and only optimize the flexible regions;
then you might be able to run in internal coords. and much faster
than optimizing the entire structure. Hope this helps,
Mark Bures
Abbott Laboratories
Abbott Park, IL 60064-3500
======================================================================== 42
I know your problem all too well, but it is easy to solve. I
solved it with the help of J. Stewart's group, and am really just passing
information along...
I was looking at perylene derivatives as a warm-up to doing C60
calculations and ran into the same errors. They derive from the way
the z matrix is set up in that small changes at one end of the molecule
rattle around to cause big problems at the other end of the molecule. Some
of this is dependent on the order in which you input the atoms.
There are two fixes. One is to reduce DMAX to prevent large
changes in the geometry from one SCF to the next. The method I settled
on is just to use the XYZ keyword. Use the z matrix you have, but then
all optimizations are done in cartesian coordinates (no symmetry!!!), and
the whole problem just goes away. This is how we've done all our
polycyclic hydrocarbons and C60 and derivatives.
AM1 is probably very good for your application. Avoid MNDO because
we found that perylene is twisted by this method! I use PM3 alot. Also,
be sure RHF calculations are adequate, evidence we have on C60 suggests
that UHF calculations are much more appropriate.
Paul Cahill
Organization 1811
Sandia National Laboratories
Albuquerque, NM 87185
voice 505-844-5754
fax 505-844-9624
pcahill %-% at %-% sandia.gov
=========================================================================
You might like to try using the XYZ keyword to run the
optimization in Cartesian space...?
--Dave Ricketts
=========================================================================
Jack, I am a regular user of MOPAC and I have written my
own molecular geometry creation programs to help solve
problems such as yours. I have never used PC-Model so what
follows is simply a suggestion based upon past experience with
large ring systems (porphorins, enzymes, etc.)
The default way most users of internal coordinates (Z-matrix) think
of the reference atoms (often referred to as the connectivity) is:
NA represents a bonded neighbor to the current atom
NB represents an atom bonded to NA, and
NC is an atom bonded to NB.
This works ok for non-annular systems and for small rings. It fails
for large rings because as MOPAC (or most any other geometry optimizing
code) varies the geometric parameters it actually alters bonded interactions
in the far part of the ring. For example, when MOPAC increases the
bond angle between atoms 1-2-3 of, eg., benzene, it moves the remainder
of the ring (atoms 4,5, and 6) away from atom 1. This breaking of a bond
looks like a very large gradient change and forces the optimization code to
quit. While the distances involved in benzene are small and the
optimizer handles them, for a large ring the distances can be enormous
(think of playing "crack-the-whip" on ice skates.)
There are two solutions: (1) optimize the geometry in cartesian space
using the XYZ keyword, or (2) use different NB and NC reference atoms.
I prefer option (2) because it allows ME to choose and I can use symmetry
or keep some functions constant. I will often use atoms 2 and 1 as NB
and NC respectively. (You will also have to watch out for reference atoms
which are colinear with atoms 1 and 2.)
There is a larger discussion of this problem in the manual for my
graphic program DRAW which is available from QCPE.
I hope this helps.
Donn M. Storch, USAF
Wright Laboratory
WPAFB, OH 45433-6563
=========================================================================
we have found the same type of problems here (I am a post-doc at EPA
NC), we use Gaussian with the AM1 HAMILTONIAN, and we are also
running some PAHs (B[a]p, anthracene, etc..). After redefining
the Z-matrix we got rid of the problem, so I would suggest you to
redefine your Z-matrix. Hope that help.
Ruben
sa %-% at %-% linus.rtpnc.epa.gov
=========================================================================
Jack Houser asked anout some error messages that were generated when
running an AM1 calculation. While many of these types of errors
are relics of older codes, this one sounds pretty common, so I'll
address it to the net in general.
Large molecules sometimes have gradient problems, and this sounds a
great deal like on of those. The FIRST thing that should be done
is to check the actaul values of the individual gradients and
the gnorm itself. (This is done in AMPAC by using the keyword
GRAD.) If one of these gradients is enormous, there is the problem.
The molecule can be redefined in order that the offending parameter
is removed. In the case of such a large system, I am leery of
using an automatic z-matrix generation routine. (Topical question:
Why use it at all if you don't use it for big guys?) Defining
symmetry to a large gradient parameter only worsens the situation
in that the bad parameter is even more important now than when it
caused generation of the first error!
There are several ways to approach one of these. I'll list but a few
below in no particular order. (Others please participate as well.
There needs to be more discussion of "real" chemistry (i.e. quantum
(haha)) on here anyhow.)
1. Redefine the molecule to get rid of the bad geomtric variable.
This is not always possible and may lead to other problems.
2. Optimize a biggy in pieces. Divide the system into easily
handled chunks and only set optimization flags on those
pieces. After this is accomplished, put the optimized
chunks back in and let the whole thing go again. This is
a bit tedious, but I've had good success with it.
3. Use dummy atoms to redefine the parameters with big gradients.
4. The wavefunction for the molecule may not be stable. It may
be that the system is a biradical and you are trying to define
it via RHF. Polyenes also suffer from nonconvergence along
this line as well. UHF or limited (or extensive) CI may be
required if this happens. [ Side note: ALWAYS use DERINU
keyword within AMPAC with CI. ]
5. Check the geometry for any wacky things, like atoms too close.
AMPAC will warn you about this within limits.
6. Draw the final geometry that is dumped and see what has happened
to the geometry. This may give you a clue about what is causing
the optimization problem.
A few more ideas. Users of AMPAC (yes, Virginia there is a MOPAC as
well) MUST check final geometries by computing frequencies at the
level of theory used to determine the geometry. This requires either
use of FORCE (force constants in Cartesian coordinates) or LTRD
(AMPAC only, force constants in internal coordinates). There
should be no negative force constants for a ground state, one for
a transition state and so forth (see page 13 of the AMPAC 2.1
manual). Before this can be done, however, the gnorm must be
reduced. You must check this at the end of your optimization.
I always use the PRECISE keyword when doing any AMPAC calculation.
This requires better refinement of the wave function and the
gnorms. It is almost required when doing FORCE analysis. For
large systems, the geometry may need to be preminimized prior to
application of PRECISE.
Andy
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
DR. ANDREW HOLDER
Assistant Professor of Computational/Organic Chemistry
Department of Chemistry || BITNET Addr: AHOLDER %-% at %-%
UMKCVAX1
University of Missouri - Kansas City || Internet Addr: aholder %-% at %-%
vax1.umkc.edu
Spencer Chemistry, Room 502 || Phone Number: (816) 235-2293
Kansas City, Missouri 64110 || FAX Number: (816) 235-1717
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
=========================================================================
It is easy to fix this problem...
1. Use "EF" in mopac...we've never had a problem with this in
MOPAC.
2. In FLEPO you can fix the problem by reducing the variable DEL
by a factor of 10 andrestore the previous geometry. This
has ALWAYS worked for us.
John McKelvey
Eastman Kodak
=========================================================================
c...This contains patches to fix the recently discussed problem
c...that can occur with FLEPO optimisation in AMPAC/MOPAC
c...between the "cNEW........" lines
C
C *
C RESTART SECTION
C *
C
cNEW..............................
89 continue
cNEW..............................
RESET=.TRUE.
DO 90 I=1,NVAR
C
C MAKE THE FIRST STEP A WEAK FUNCTION OF THE GRADIENT
C
STEP=ABS(GRAD(I))*0.0002D0
STEP=MAX(0.01D0,MIN(0.04D0,STEP))
C# XD(I)=XPARAM(I)-SIGN(STEP,GRAD(I))
XD(I)=XPARAM(I)-SIGN(DEL,GRAD(I))
90 CONTINUE
C# WRITE(6,'(10F8.3)')(XD(I)-XPARAM(I),I=1,NVAR)
C
C THIS CALL OF COMPFG IS USED TO CALCULATE THE SECOND-ORDER MATRIX IN H
C IF THE NEW POINT HAPPENS TO IMPROVE THE RESULT, THEN IT IS KEPT.
C OTHERWISE IT IS SCRAPPED, BUT STILL THE SECOND-ORDER MATRIX IS O.K.
C
C# WRITE(6,*)' RESET HESSIAN'
CALL COMPFG (XD,.TRUE.,FUNCT2,.TRUE.,GD,.TRUE.)
IF(.NOT. GEOOK .AND. SQRT(DOT(GD,GD,NVAR))/GNORM.GT.10.
1 AND.GNORM.GT.20.AND.JCYC.GT.2)THEN
C
C THE GEOMETRY IS BADLY SPECIFIED IN THAT MINOR CHANGES IN INTERNAL
C COORDINATES LEAD TO LARGE CHANGES IN CARTESIAN COORDINATES, AND THESE
C LARGE CHANGES ARE BETWEEN PAIRS OF ATOMS THAT ARE CHEMICALLY BONDED
C TOGETHER.
WRITE(IPRT,'('' GRADIENTS OF OLD GEOMETRY, GNORM='',F13.6)')
1 GNORM
WRITE(IPRT,'(6F12.6)')(GRAD(I),I=1,NVAR)
GDNORM=SQRT(DOT(GD,GD,NVAR))
WRITE(IPRT,'('' GRADIENTS OF NEW GEOMETRY, GNORM='',F13.6)')
1 GDNORM
WRITE(IPRT,'(6F12.6)')(GD(I),I=1,NVAR)
cNEW..............................
del=del/10.0
c...You may wish to play with thresh.
thresh=0.00005d0
if(del.ge.0.thresh)then
go to 89
else
cNEW..............................
WRITE(IPRT,'(///20X,''CALCULATION ABANDONED AT THIS POINT!''
1)')
WRITE(IPRT,'(//10X,'' SMALL CHANGES IN INTERNAL COORDINATES
1ARE '',/10X,'' CAUSING A LARGE CHANGE IN THE DISTANCE BETWEEN'',
2/ 10X,'' CHEMICALLY-BOUND ATOMS. THE GEOMETRY OPTIMIZATION'',/
3 10X,'' PROCEDURE WOULD LIKELY PRODUCE INCORRECT RESULTS'')')
CALL GEOUT(1)
STOP
cNEW..............................
endif
cNEW..............................
=========================================================================
In MOPAC another way around the problem is to use the keyword XYZ.
This instructs the optimization to be carried out in cartesian coordinates.
brian yates.
=========================================================================
MOPAC GEOM OPT A QUICK SOLUTION
Whenever I ran into a problem with large floppy cyclic molecules
using Mopac i did the same calc using the AM1 or MNDO option on
the latest set of Gaussian prog. The math tecniques are entirely differe
nt and probably more precise since the tolerences for convergence are
much greater.(ANDY H. correct me if I am wrong). Anyway whatever the
reason it worked. In one or two cases especially when the semiemp.
predicted a geometry which I thought was unreasonable the G. prog
gave me a more reasonable result. The difference in computer time was ma
maybe a factor of two. (Since 2*0=0 I decided this was more efficient
with respect to my time than fiddling with the Z matrix. However, I
believe there is a recent paper on how to obtain a
Z matrix for cyclic molecules probably in J comp. Chem 1991
(John E. Bloor, pa13808 at UTKVM1)
=========================================================================
One problem which must always be avoided in the definition of a Z matrix
is defining your rings simply as a ring. In other words, if you have a
six membered ring with atoms numbered 1-6 (clockwise) so that the ring
closure is defined by 6-1, you should not define your Z matrix such that
the atomic distance column includes 1-2, 2-3, 3-4, 4-5 and 5-6, nor should
the bond or dihedral angles be defined using such a scheme (i.e. 3-2-1,
4-3-2, or 4-3-2-1, etc.). This can lead to small changes at the beginning
of the sequence causing large distance or angle changes at the end of the
sequence. (Anybody remember playing crack the whip as a child?) I am not
sure how PCModel defines rings, but if it does what I described above, it
might lead to one of the errors originally posted here.
We now use the AVS Chemistry Viewer to set up many of our jobs, and do not
get poorly defined Z matrices.
Doug Smith
Assistant Professor of Chemistry
The University of Toledo
Toledo, OH 43606-3390
voice 419-537-2116
fax 419-537-4033
email fax0236 %-% at %-% uoft02.utoledo.edu
=========================================================================
In reply to John Bloor,
The criteria are not more stringent in GAUSSIAN than in AMPAC.
What is different is the optimization algorithm. I believe
that the default optimizer in GAUSSIAN is what is known as EF
(eigenvector following) on the semiempirical side. (Someone
please correct me if I am wrong. This is only an impression
not a fact.) If this is the case, John McKelvey pointed this
out before as a solution. I am looking at his other fix and
trying to figure out just what overall effect it will have. I
think that Henry Rzepa's implementation of EF will find its way
into the next version of AMPAC.
Andy Holder
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
DR. ANDREW HOLDER
Assistant Professor of Computational/Organic Chemistry
Department of Chemistry || BITNET Addr: AHOLDER %-% at %-%
UMKCVAX1
University of Missouri - Kansas City || Internet Addr: aholder %-% at %-%
vax1.umkc.edu
Spencer Chemistry, Room 502 || Phone Number: (816) 235-2293
Kansas City, Missouri 64110 || FAX Number: (816) 235-1717
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
=========================================================================
There seems to be some confusion about what Gaussian does during
semi-empirical optimizations. Here are some details:
1. Gaussian defaults to accuracy comparable to that produced by the
PRECISE keyword in MOPAC.
2. The step-size used in numerically differentiating the integrals as
part of the gradient calculation is small in Gaussian, because a
smaller step can be taken reliably if it is known in advance that
PRECISE is always turned on. This provides somewhat greater
accuracy in the derivatives. I think something like this (but
with different values for the step size) is done in MOPAC version
6 but not in earlier versions.
3. Only the energy and cartesian gradient are computed using code from
MOPAC; the rest of the optimization uses the standard Gaussian
routines in exactly the same manner as an ab initio optimization
would.
4. The cartesian forces are converted to internal coordinates using
different algorithms. I don't know the details of how MOPAC does
the conversion, but this appears to make a difference. Gaussian
uses analytic expressions and doesn't suffer from the sensitivity
to internal coordinate definition that MOPAC apparently does,
since MOPAC gives the warning message about small changes in
internal coordinates causing large displacements for cases where
Gaussian does not have difficulties.
5. Gaussian's default optimizer is not the same as either EF in
Gaussian or EF in MOPAC. It does use an RFO (Rational Function
Optimization) step for the quadratic portion of the step as in the
EF method, but it also does gradient-based linear corrections at
every step as suggested by Schlegel, and it uses an improved
version of Schlegel's update scheme for the Hessian. EF (Baker's
algorithm) is available as an option, but is usually inferior to
the default, since the default is a mixture of the best parts of
Baker's and Schlegel's methods.
Michael Frisch
Gaussian, Inc.
-------
=========================================================================
Thanks to Mike Frisch for clearing up the confusion on G90 optimi-
izers. As I noted in my earlier post, I would reccomend all final
optimizations employ the keyword PRECISE for refinement of the
gradients and accuracy in the wavefunction.
As regards some of the ideas in Doug Smith's notes on optimization
of undefined bonds in rings, I have never experienced this diff-
iculty with AMPAC optimizers. Because the time required for an
energy calculation is so much shorter than ab initio methods, it
is my impression that semiempiricl optim"G~zers are a bit more
"laid back" for want of a better term. This type of event was
a classic case with early versions of BERNY I believe. I think
that GAUSSIAN, Inc. Has repaired this though.
Andy Holder
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
DR. ANDREW HOLDER
Assistant Professor of Computational/Organic Chemistry
Department of Chemistry || BITNET Addr: AHOLDER %-% at %-%
UMKCVAX1
University of Missouri - Kansas City || Internet Addr: aholder %-% at %-%
vax1.umkc.edu
Spencer Chemistry, Room 502 || Phone Number: (816) 235-2293
Kansas City, Missouri 64110 || FAX Number: (816) 235-1717
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
=========================================================================
By now you've probably had a lot of replies (I've been away a week). I did
some MOPAC work with pyrene at one point. I found that I had to define the
z-matrix in such a way as to minimize the lengths of chains of atoms which
end up referring back to each other. If that's not clear, let me know &
I'll
surface mail you a drawing of what I mean...
Irene Newhouse
Neely %-% at %-% ducvax.auburn.edu
=========================================================================
Try a redefinition of the internal coordinates. Perhaps you could place a
dummy atom in the middle of the rings and define the atoms in reference to that
dummy atom. An alternate approach is to run the calculation in cartesian
coordinates.
Eve Zoebisch
zoebisch %-% at %-% crvax.sri.com
=========================================================================
Date: Thu, 23 Apr 92 14:43:03 EDT
From: "Dr./CPT Christopher J. Cramer" <cjcramer %-% at %-%
crdec6.apgea.army.mil>
To: chemistry %-% at %-% ccl.net
Subject: Huckel MO Theory software
Status: RO
For use in a course on applied theory at all levels, I am interested in
finding some friendly software which does Huckel theory. Macintosh, UNIX
mainframe, or IBM PC clone all acceptable as platforms in decreasing
order of desirability.
I am aware of an HMO program distributed by Trinity at moderate cost, and
I'm told of some shareware by Farrell and Haddon (although short of
writing them I've no idea yet how to get it). I'd appreciate input from
authors/users of such programs. Happy to summarize for net.
Thanks,
Chris Cramer
(cjcramer %-% at %-% crdec.apgea.army.mil or mf12131 %-% at %-% msc.edu [milnet
preferred])