First of all, and most important, before any modifications are done a back-up copy of the standard MOPAC should be made. This will prove invaluable in pinpointing deviations from the standard working. This point cannot be over-emphasized--make a back-up before modifying MOPAC!
Clearly, a bug can occur almost anywhere, and a logical search sequence is necessary in order to minimize the time taken to locate it.
If possible, perform the debugging with a small molecule, in order to save time (debugging is, of necessity, time consuming) and to minimize output.
The two sets of subroutines in MOPAC, those involved with the electronics and those involved in the geometrics, are kept strictly separate, so the first question to be answered is which set contains the bug. If the heats of formation, derivatives, I.P.s, and charges, etc., are correct, the bug lies in the geometrics; if faulty, in the electronics.
Use formaldehyde for this test. Use keywords 1SCF, DEBUG, and any others necessary.
The main steps are:
If an error is found, use MOLDAT to verify the orbital character, etc.
If faulty the error lies in READMO, GETGEO or MOLDAT.
Otherwise the error lies in HCORE, H1ELEC or ROTATE.
If the starting matrices are correct, go on to step (2).
If an error is found, check the first Fock matrix. If faulty, the bug lies in ITER, probably in the Fock subroutines FOCK1 or FOCK2, or in the (guessed) density matrix (MOLDAT). An exception is in the UHF closed-shell calculation, where a small asymmetry is introduced to initiate the separation of the α and β UHF wavefunctions.
If no error is found, check the second Fock matrix. If faulty, the error lies in the density matrix DENSIT, or the diagonalization RSP.
If the Fock matrix is acceptable, check all the Fock matrices. If the error starts in iterations 2 to 4, the error probably lies in CNVG, if after that, in PULAY, if used.
If SCF is achieved, and the heat of formation is faulty, check HELECT. If C.I. was used check MECI.
If the derivatives are faulty, use DCART to verify the Cartesian derivatives. If these are faulty, check DCART and DHC. If they are correct, or not calculated, check the DERIV finite difference calculation. If the wavefunction is non-variationally optimized, check DERNVO.
If the geometric calculation is faulty, use FLEPO or PRNT=5 to monitor the optimization, DERIV may also be useful here.
For the FORCE calculation, DCART or DERIV are useful for variationally optimized functions, COMPFG for non-variationally optimized functions.
For reaction paths, verify that FLEPO is working correctly; if so, then PATHS is faulty.
For saddle-point calculations, verify that FLEPO is working correctly; if so, then REACT1 is faulty.
Keep in mind the fact that MOPAC is a large calculation, and, while intended to be versatile, many combinations of options have not been tested. If a bug is found in the original code, please communicate details to Dr. James J. P. Stewart, Stewart Computational Chemistry, 15210 Paddington Circle, Colorado Springs, CO 80921-2512; E-mail: openmopac@gmail.com.