From dchin@polygenh.prospect.com Thu Apr 23 02:58:31 1992 Date: Thu, 23 Apr 92 02:47:11 EDT From: dchin@polygenh.prospect.com (Donovan Chin) Status: RO Dear Dr. Heller, I would like to be considered as a reviewer for software in question. I am a post-doctoral scientist in chemistry at harvard university, and have been a computational chemist throughout my graduate and present post-doc terms. I have worked part time on and off for many years at polygen (now molecular simulations inc). Our lab at harvard contains several quadras, indigo workstations and high end pc's (386 and 486). I am also speaking for several of my collegues (post docs and grad students) who also have expressed interest in reviewingyour software. Thank you for your consideration. Dr. Donovan Chin dept of chemistry harvard university 12 oxford st cambridge, MA 02138 From chemistry-request@ccl.net Thu Apr 23 15:31:30 1992 Date: Thu, 23 Apr 92 09:53:54 EST From: Jack Houser Subject: Optimization in MOPAC To: Supercomputer bulletin board Status: RO 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@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@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@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@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@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@UMKCVAX1 University of Missouri - Kansas City || Internet Addr: aholder@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@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@UMKCVAX1 University of Missouri - Kansas City || Internet Addr: aholder@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@UMKCVAX1 University of Missouri - Kansas City || Internet Addr: aholder@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@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@crvax.sri.com ========================================================================= Date: Thu, 23 Apr 92 14:43:03 EDT From: "Dr./CPT Christopher J. Cramer" To: chemistry@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@crdec.apgea.army.mil or mf12131@msc.edu [milnet preferred])