Later, in 1977 MNDO [1] appeared. Initially, parameters were available for C, H, N, and O only [1]. Results were much better [98] than for MINDO/3, and MNDO was soon extended to Be[150], B [99], F [100],  Al [99], Si [101], P [102], S [103], Cl [104], K [135], Zn [105], Ge [106], Br [107], Sn [108], I [109], Hg [110], and Pb [111]. In 2004, the remaining main-group elements Na, Mg, K, Ca, Ga, As Se, Rb, Sr, In, Sb, Te, Cs, Ba, Tl and Bi were published [138].  Earlier versions of MOPAC used parameters for
Na [134], and K [135], but these have been replaced by newer parameters that are intended to be more accurate.
Many of the faults in MNDO were corrected in AM1 [3], which also was initially available only for H, C, N, and O [3]. AM1 has since been extended to F [113], Al [114], Si [115], P [116], S [117], Cl [113], Zn [118], Ge [119], Br [113], Mo [133] , I [113], and Hg [120]. In 2004, the remaining main-group elements Li, Be, Na, Mg, K, Ca, Ga, As, Se, Rb, Sr, In, Sn, Sb, Te, Cs, Ba, Pb, and Bi were published [138]. Sparkle models for La, Ce, Pr[142], Nd[141], Pm, Sm, Eu[137],Gd[137], Tb[137], Dy[139], Ho[140], Er, Tm[143], Yb[144], and Lu and Th[145] have been developed.
In 1985, PM3 [4] was developed. Initially, 12 elements were available: H, C, N, O, F, Al, Si, P, S, Cl, Br, and I [121]. This set was expanded in 1991 to include Be, Mg, Zn, Ga, Ge, As, Se, Cd, In, Sn, Sb, Te, Hg, Tl, Pb, and Bi [122]. In 1993 lithium was parameterized [123], and in 2004, the remaining main-group elements B, Na, K, Ca, Rb, Sr, Cs, and Ba were published [138].


        A parameter set was developed to give improved results for systems containing the elements H, C, N, O, P, S, F, Cl, Br, and I only   For details, see Rocha, G.B., R.O. Freire, A.M. Simas, and J.J.P. Stewart., RM1: A Reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I. J. Comp. Chem., 27(10): 1101-1111 (2006).       

When d orbitals are added to the basis set [6], the accuracy of the method rises. At present (1998), parameters are available for Al, Si, P, S, Cl, Br, and I [6].
PM6 parameters[148] are available for all elements from H to Bi, except the lanthanides from Ce to Yb. All the lanthanides[149] can be used, as their 3+ sparkles.  To use them, add keyword SPARKLE.  The sparkle versions of Ln give good geometries.

All methods 

Eight point-charges, called "sparkles" are provided, these have charges of +3 +2, +1, +1/2, -1/2, -1, -2, and -3.  They were developed for use as counterions.  No journal reference exists, the only description of the sparkles is in this manual.


By default, geometries are optimized using Baker's EF routine [7]. If this is not desired, then the Broyden [19]-Fletcher [20]-Goldfarb [21]-Shanno [22] method can be used. This replaces the older Davidon [124]-Fletcher [125]-Powell [126] method. The Powell line search in both the BFGS and DFP methods has been upgraded with Thiel's FSTMIN technique [44].

The energy minimum for intersystem crossing can be calculated using the method of Anglada and Bofill [28].

Reaction paths can be followed, but sometimes (unavoidable) numerical instability causes difficulty [127].

Various methods can be used to locate transition states. In the region of a transition state, TS [7] is the best method to use. Other options are NLLSQ [8] and SIGMA [9,10].

However, getting to this region can often be difficult. One effective strategy is to use the reaction path option. A more costly method is the SADDLE technique [18]. Once the transition state is located, the geometry must be refined

All systems can be characterized by determining the number of imaginary vibrations using FORCE [73]. Ground state systems should have none, and transition states should have exactly one imaginary frequency. The normal modes can then be used to calculate heat capacities and entropies and other thermodynamic quantities [50].

Molecular dynamics can be followed via DRC [29].


Semiempirical methods use approximations to the Roothaan [52]-Hall [53] equations. The four approximate methods are MNDO [1], AM1 [3], PM3 [4] , RM1*, and PM6 [4] For RHF open shell systems, a Multi-Electron-Configuration-Interaction procedure [39] is available. SCF convergence can be assisted by an inexpensive SHIFT technique [49], by a slightly more complicated Direct Inversion of the Iterative Subspace method (DIIS) [48], or by an expensive but sophisticated interpolation procedure [128].
For large systems, an alternative to conventional methods is to use localized molecular orbitals [129]. This is usually more efficient, particularly for geometry optimization [130].
For variationally-optimized wavefunctions, two derivative methods are available. The default is to use finite difference, but under user request (ANALYT), analytical derivatives can be used [16]. When analytical derivatives are used, STO-6G [17] Gaussian functions are used instead of Slater orbitals.
When non-variationally optimized wavefunctions are used, the derivatives are much more complicated [43].

*: Rocha, G.B., R.O. Freire, A.M. Simas, and J.J.P. Stewart., RM1: A Reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I. J. Comp. Chem., 27(10): 1101-1111 (2006).


Fundamental constants are taken from the CODATA report [71]. A good introduction to MOPAC can be found in Tim Clark's book [131].

The SCF M.O.s, which diagonalize the Fock matrix, can be localized [38] to give M.O.s which can be identified with the conventional picture of two-electron bonds and lone pairs. The localization scheme is faster at the semiempirical level than the Edmiston-Ruedenberg [62] or Boys [61] methods. Associated with each conventional M.O. is a bond-index [23], which represents the contribution to the bond-order matrix due to each M.O. Bond orders and valencies can be displayed by use of BONDS [15]. Other phenomena relating to bonding can also be calculated [24,25,26]. An alternative to the normal Coulson density matrix is the Mulliken [132],  the center of mass is used. Higher terms, e.g., polarizability, and first and second hyperpolarizability, can be calculated [47] by POLAR.
Ionization potentials [59] can be corrected using Green's Functions [34,35,36,33,37,32]

Solvent and Electrostatics
Solvent phenomena can be studied. The COSMO technique [30], unlike the self-consistent reaction fields [74], allows geometries to be optimized. Although the Miertus-Scrocco-Tomasi model [75,76] cannot optimize geometries, is more sophisticated in that it allows cavitation effects. This model has been modified [77,78,79,80,81] to allow NDDO methods to be used. In this, optimized VdW radii [83,84] are used to construct [82] a cavity.
The free energy of hydration is computed as the addition of three contributions:
1. The electrostatic term, which is computed from the linear free energy response theory [75,76,77,78,79,80,81].
2. The cavitation contribution, which is computed from Pierotti's scaled particle theory [85].
3. The van der Waals terms, which is computed using a linear relation with the solute accessible surface, and optimized "hardness" parameters [83,84].

In addition to the free energy of hydration a "solvent-adapted" wavefunction is obtained. Such a wavefunction can be used to determine changes in solute properties due to the solvent [86,87,88,89].

Electrostatic potentials can be used with the MST method both by deorthogonalizing the wavefunction [90,91,40] and by keeping the wavefunction orthogonal [92,93]

Other ESP methods available are the Merz-Bessler technique [31] and the Ford-Wang procedure [45,46]. The Ford-Wang is much faster and more accurate than the Merz-Bessler method, but is limited to AM1 calculations on systems containing H, C, N, O, F, and Cl, only.