The Fragment Molecular Orbital method

coded by D.G. Fedorov, M.Chiba, T. Nagata and K. Kitaura at
    Research Institute for Computational Sciences (RICS)
     National Institute of Advanced Industrial Science
                   and Technology (AIST)
           AIST Tsukuba Central 2, Umezono 1-1-1,
                 Tsukuba, 305-8568, Japan.
with code contributions by:
       N. Asada (Kyoto U.), C. H. Choi (Kyungpook U.),
            C. Steinmann (U. Copenhagen).

The method was proposed by Professor Kitaura and coworkers 
in 1999, based on the Energy Decomposition Analysis (EDA, 
sometimes called the Morokuma-Kitaura energy 
decomposition). The FMO method is completely independent of 
and bears no relation to:
      1. Frontier molecular orbitals (FMO),
      2. Fragment molecular orbitals (FMO).
The latter name is often used for the process of 
construction of full molecular orbitals by combining MO 
diagrams for parts of a molecule, ala Roald Hoffmann.
The effective fragment molecular orbital method (EFMO) is 
closely related to but also bears significant difference to 
FMO, and discussed below.

The FMO program was interfaced with GAMESS and follows 
general GAMESS guidelines for code distribution and usage.  
The users of the FMO program are requested to cite the 
FMO3-RHF paper as the basic FMO reference,
        D.G. Fedorov, K. Kitaura,
        J. Chem. Phys. 120, 6832-6840(2004)
and other papers as appropriate (see below).
 
The basic idea of the method is to acknowledge the fact the 
exchange and self-consistency are local in most molecules 
(and clusters and molecular crystals), which permits 
treating remote parts with Coulomb operators only, ignoring 
the exchange.  This idea further evolves into doing 
molecular calculations, piecewise, with Coulomb fields due 
to the remaining parts.  In practice one divides the 
molecule into fragments and performs n-mer calculations of 
these in the Coulomb field of other fragments (n=1,2,3).  
There are no empirical parameters, and the only departure 
from ab initio rigor is the subjective fragmentation.  It 
has been observed that if performed physically reasonably, 
the fragmentation scheme alters the results very little.  
What changes the accuracy the most is the fragment size, 
which also determines the computational efficiency of the 
method.
 
The first question is how to get started.  The easiest way 
to prepare an FMO input file for GAMESS is to use the free 
GUI software Facio, developed by M. Suenaga at Kyushu 
University. It can do molecular modeling, automatic 
fragmentation of peptides, nucleotides and saccharides and 
create GAMESS/FMO input files:
    http://www1.bbiq.jp/zzzfelis/Facio.html
A web bsed interface to FMO is maintained by Y. Alexeev 
(Argonne National Lab):
    http://www.fmo-portal.info (shut down at present)

Alternatively, if you prefer a command line interface, and 
your molecule is a protein found in the PDB
    http://www.rcsb.org/pdb
you can simply use the fragmentation program "fmoutil" that 
is provided with GAMESS in tools/fmo, or the FMO home page
    http://staff.aist.go.jp/d.g.fedorov/fmo/main.html
If you have a cluster of identical molecules, you can 
perform fragmentation with just one keyword ($FMO NACUT=). 

Computationally, it is always better to partition in a 
geometrical way (close parts together), so that the 
distance-based approximations are more efficient. The 
accuracy depends mainly upon the locality of the density 
distribution, and the appropriateness of partitioning it 
into fragments. There is no simple connexion between the 
geometrical proximity of fragmentation and accuracy.

Supposing you know how to fragment, you should choose a 
basis set and fragment size.  We recommend 2 amino acid 
residues or 2-4 water molecules per fragment for final 
energetics (or, even better, three-body with 1 molecule or 
residue per fragment).  For geometry optimizations one may 
be able to use 1 res/mol per fragment, especially if 
gradient convergence to about 0.001 is desired.  Note that 
although it was claimed that FMO gradient is analytic 
(Chem. Phys. Lett., 336 (2001), 163.) it is not so. Neither 
theory nor program for fully analytic gradient has been 
developed, to the best of our knowledge up to this day 
(December 21, 2006).  The gradient implementation is nearly 
analytic, meaning three small terms are missing, one which 
can now be included using MODGRD=8+2.  The magnitude of 
these small terms depends upon the fragment size (larger 
fragments have smaller errors).  It has been our experience 
that in proteins with 1 residue per fragment one gets 1e-
3...1e-4 error in the gradient, and with 2 residues per 
fragment it is about 1e-4...1e-5. If you experience energy 
rising during geometry optimizations, you can consider two 
countermeasures:
1. increase approximation thresholds, e.g. RESPPC from
   2.0->2.5, RESDIM from 2.0 -> 2.5.
2. increase fragment size (e.g. by merging very small
   fragments with their neighbors). 
Finally a word of caution: optimizing systems with charged 
fragments in the absence of solvent is frequently not a 
good idea: oppositely charged fragments will most likely 
seek each other, unless there is some conformational 
barrier.

For basis sets you should use general guidelines and your 
experience developed for ab initio methods. There is a file 
provided (HMOs.txt) that contains hybrid molecular orbitals 
(HMO) used to divide the MO space along fragmentation 
points at covalent bonds. If your basis set is not there 
you need to construct your own set of HMOs. See the example 
file makeLMO.inp for this purpose.

Next you choose a wave function type. At present one can 
use RHF, DFT, MP2, CC, and MCSCF (all except MCSCF support 
the 3-body expansion).  Geometry optimization can be 
performed with all of these methods, except CC.

Note that presence of $FMO turns FMO on.

Surfaces and solids

Until 2008, for treating covalently connected fragments, 
FMO had fully relaxed electron density of the detached 
bonds. This method is now known as FMO/HOP (HOP=hybrid 
orbital projection operator). It allows for a full 
polarization of the system and is thus well suited to very 
polar systems, such as proteins with charged residues. In 
2008, an alternative fragmentation was suggested, based on 
adaptive frozen orbitals (AFO), FMO/AFO. In it, the 
electron density for each detached bond is first computed 
in the automatically generated small model system (with the 
bond intact), and in the FMO fragment calculations this 
electron density is frozen. It was found that FMO/AFO works 
quite well for surfaces and solids, where there is a dense 
network of bonds to be detached in order to define 
fragments (and the detached bonds interact quite strongly). 
In addition, by restricting the polarization, FMO/AFO was 
found to give a more balanced properties for large basis 
sets (triple-zeta with polarization or larger), or in 
comparing different isomers. However, for proteins with 
charged residues the original FMO/HOP scheme has a better 
accuracy (except large basis sets). At this point, FMO/AFO 
was applied to zeolites only, and some more experience is 
needed to give more practical advice to applications. 
FMO/AFO is turned on by a nonzero rafo(1) parameter (rafo 
array provides the thresholds to build model systems).


FMO variants

In 2007, Dahlke et al. introduced the Electrostatically 
Embedded Many-Body Expansion method (see E. E. Dahlke and 
D. G. Truhlar, J. Chem. Theory Comput. 4, 1-6 (2008) for 
more recent work). This method is essentially FMO with the 
RESPPC approximation (point charges for the electrostatic 
field) applied to all fragments, with the further provision 
that these charges may be defined at will (whereas RESPPC 
uses Mulliken charges), and they are kept frozen (not 
optimized, as in FMO). Next, Kamiya et al. suggested a fast 
electron correlation method (M. Kamiya, S. Hirata, M. 
Valiev, J. Chem. Phys. 128, 074103 (2008)), where again FMO 
with the RESPPC approximation to all fragments is applied 
with the further provision that the charges are derived 
from the electrostatic potential (so called ESP charges), 
and BSSE correction is added. The Dahlke's method was 
generalized in GAMESS with the introduction of an arbitrary 
hybrid approach, in which some fragments may have fixed and 
some variationally optimized charges. This implementation 
was employed in FMO-TDDFT calculations of solid state 
quinacridone (see Ref. 16 below) by using DFT/PBC frozen 
charges.  The present energy only implementation is mostly 
intended for such cases as that (i.e., TDDFT), and some 
more work is needed to finish it for general calculations. 
To turn this on, set RESPPC=-1 and define NOPFRG for frozen 
charge fragments to 64, set frozen charges in ATCHRG. 
Another FMO-like method is EFMO, see its own subsection 
below. EFMO itself is related to several methods (PMISP: P. 
Soederhjelm, U. Ryde, J. Phys. Chem. A 2009, 113, 617?627; 
another is G. J. O. Beran, J. Chem. Phys. 2009, 130, 
164115).

Effective fragment molecular orbital method (EFMO)

EFMO has been formulated by combining the physical models 
in EFP and FMO, namely, in EFMO, fragments are computed 
without the ESP (of FMO), and the polarization is estimated 
using EFP models of fragment polarizabilities, which are 
computed on the fly, so this can be thought of as 
automatically generated potentials in EFP. Consequently, 
close dimers are computed quantum-mechanically (without 
ESP) and far dimers are computed using the electrostatic 
multipole models of EFP. At present, only vacuum closed-
shell RHF and DFT are supported, for energy and gradient; 
and only molecular clusters can be computed (no systems 
with detached bonds). From the user point of view, EFMO 
functionality is very intensively borrowed from FMO, and 
the calculation setup is almost identical. Most additional 
physical models such as PCM are not supported in EFMO. EFMO 
should not be confused with FMO/EFP. The latter uses FMO 
for some fragments and EFP for others. EFMO uses the same 
model (EFMO), which is neither FMO nor EFP. For 
approximations, EFMO at present has only RESDIM.

EFMO references are:
1. Effective Fragment Molecular Orbital Method: A Merger of
   the Effective Fragment Potential and Fragment Molecular
   Orbital Methods.
       C. Steinmann, D. G. Fedorov, J. H. Jensen
       J. Phys. Chem. A 114, 8705-8712 (2010).
2. The Effective Fragment Molecular Orbital Method for
   Fragments Connected by Covalent Bonds.
       C. Steinmann, D. G. Fedorov, J. H. Jensen
       PLoS One, 7, e41117(2012).
3. Mapping enzymatic catalysis using the effective fragment 
molecular orbital method: towards all ab initio 
biochemistry.
       C. Steinmann, D. G. Fedorov, J. H. Jensen
       PLoS One 8, e60602 (2013).

Guidelines for approximations with FMO3

Three sets are suggested, for various accuracies:
  low:    resppc=2.5 resdim=2.5  ritrim(1)=0.001,-1,1.25
  medium: resppc=2.5 resdim=3.25 ritrim(1)=1.25,-1,2.0
  high:   resppc=2.5 resdim=4.0  ritrim(1)=2,2,2
For correlated runs, add one more value to ritrim, equal to 
the third element (i.e., 1.25 or 2.0).  Note that gradient 
runs do not support nonzero RESDIM and thus use RESDIM=0 if 
gradient is to be computed.  The "low" level of accuracy 
for FMO3 has an error versus full ab initio similar to 
FMO2, except for extended basis sets (6-311G** etc) where 
it is substantially better than FMO2. Thus the low level is 
only recommended for those large basis sets, and if a 
better level cannot be afforded.  The medium level is 
recommended for production FMO3 runs; the high level is 
mostly for accuracy evaluation in FMO development.  The 
cost is roughly: 3(low), 6(medium), 12(high). This means 
that FMO3 with the medium level takes roughly six times 
longer than FMO2.

Some of the default tolerances were changed as of January 
2009, when FMO 3.2 was included in GAMESS.  In general, 
stricter parameters are now enforced when using FMO3, which 
of course is intended to produce more accurate results.  If 
you wish to reproduce earlier results with the new code, 
use the input to revert to the earlier values:
        former -> FMO2 or FMO3 (as of 1/2009)
  RESPPC: 2.0      2.0    2.50
  RESDIM: 2.0      2.0    3.25
  RCORSD: 2.0      2.0    3.25
  RITRIM: 2.0,2.0,2.0,2.0 -> 1.25,-1.0,2.0,2.0 (FMO3 only)
  MODESP: 1         0       1
  MODGRD: 0        10       0
and two other settings which are not strictly speaking FMO 
keywords may change FMO results:
  MTHALL: 2 -> 4  (FMO/PCM only, see $TESCAV)
  DFT grid: spherical -> Lebedev (FMO-DFT only, see $DFT)
Note that FMO2 energies printed during a FMO3 run will 
differ from those in a FMO2 run, due to the different 
tolerances used.

How to perform FMO-MCSCF calculations

Assuming that you are reasonably acquainted with ab initio 
MCSCF, only FMO-specific points are highlighted. The active 
space (the number of orbitals/electrons) is specified for 
the MCSCF fragment. The number of core/virtual orbitals for 
MCSCF dimers will be automatically computed.  The most 
important issue is the initial orbitals for the MCSCF 
monomer.  Just as for ab initio MCSCF, you should exercise 
chemical knowledge and provide appropriate orbitals. There 
are two basic ways to input MCSCF initial orbitals:
   A) through the FMO monomer density binary file
   B) by providing a text $VEC group.
The former way is briefly described in INPUT.DOC (see 
orbital conversion). The latter way is really identical to 
ab initio MCSCF, except the orbitals should be prepared for 
the fragment (so in many cases you would have to get them 
from an FMO calculation). Once you have the orbitals, put 
them into $VEC1, and use the IJVEC option in $FMOPRP (e.g., 
if your MCSCF fragment is number 5, you would use $VEC1 and 
ijvec(1)=5,0).  For two-layer MCSCF the following 
conditions apply.  Usually one cannot simply use F40 
restart, because its contents will be overwritten with RHF 
orbitals and this will mess up your carefully chosen MCSCF 
orbitals. Therefore, two ways exist. One is to modify A) 
above by reordering the orbitals with something like
 $guess guess=skip norder=1 iorder(28)=29,30,31,32,28 $end
Then the lower RHF layer will converge RHF orbitals that 
you reorder with iorder in the same run (add 512 to nguess 
in $FMO). This requires you know how to reorder before 
running the job so it is not always convenient.  Probably 
the best way to run two-layer MCSCF is verbatim B) above, 
so just provide MCSCF monomer orbitals in $VEC1. Finally, 
it may happen that some MCSCF dimer will not converge. 
Beside the usual MCSCF tricks to gain convergence as the 
last resort you may be able to prepare good initial dimer 
orbitals, put them into $VEC2 ($VEC3 etc) and read them 
with ijvec.  SOSCF is the preferred converger in FMO, and 
the other one (FULLNR) has not been modified to eradicate 
the artefacts of convergence (due to detached bonds).  In 
the bad cases you can try running one or two monomer SCF 
iterations with FULLNR, stop the job and use its orbitals 
in F40 to do a restart with SOSCF.  We also found useful to 
set CASDII=0.005 and nofo=10 in some cases running FOCAS 
longer to get better orbitals for SOSCF.
 
How to perform multilayer runs

For some fragments you may like to specify a different 
level of electron correlation and/or basis set. In a 
typical case, you would use high level for the reaction 
center and a lower level for the remaining part of the 
system.  The set up for multilayer runs is very similar to 
the unilayer case.  You only have to specify to what layer 
each fragment belongs and for each layer define DFTTYP, 
MPLEVL, SCFTYP as well as a basis set.  If detached bonds 
are present, appropriate HMOs should be defined.  See the 
paragraph above for multilayer MCSCF.  Currently geometry 
optimizations of multilayer runs require adding 128 to 
NGUESS, if basis sets in layers differ from each other.

How to mix basis sets in FMO

You can mix basis sets in both uni and multilayer cases. 
The difference between a 2-layer run with one basis set per 
layer and a 1-layer run with 2-basis sets is significant: 
in the former case the lower level densities are converged 
with all fragments computed at the lower level. In the 
latter case, the fragments are converged simultaneously, 
each with its own basis set. In addition, dimer corrections 
between layers will be computed differently: with the lower 
basis set in the former case and with mixed basis set in 
the latter.  The latter approach may result in unphysical 
polarization, so mixing basis sets is mainly intended to 
add diffuse functions to anionic (e.g., carboxyl) groups, 
not as a substitute for two-layer runs.

How to perform FMO/PCM calculations

Solvent effects can be taken into account with PCM.  PCM in 
FMO is very similar to regular PCM.  There is one basic 
difference: in FMO/PCM the total electron density that 
determines the electrostatic interaction is computed using 
the FMO density expansion up to n-body terms.  The cavity 
is constructed surrounding the whole molecule, and the 
whole cavity is used in each individual m-mer calculation. 
There are several levels of accuracy (determined by the "n" 
above), and the recommended level is FMO/PCM[1(2)], 
specified by:

 $pcm ief=-10 icomp=2 icav=1 idisp=1 ifmo=2 $end
 $fmoprp npcmit=2 $end
 $tescav ntsall=240 $end
 $pcmcav radii=suahf $end

Many PCM options can be used as in the regular PCM. The 
following restrictions apply:
   IEF may be only -3 or -10, IDP must be 0.
Multilayer FMO runs are supported.  Restarts are limited to 
IREST=2, and in this case PCM charges (the ASCs) are not 
recycled. However, the initial guess for the charges is 
fairly reasonable, so IREST=2 may be useful although 
reading the ASCs may be implemented in future.

Note for advanced users. IFMO < NBODY runs are permitted. 
They are denoted by FMOm/PCM[n], where m=NBODY and n=IFMO. 
In FMOm/PCM[n], the ASCs are computed with n-body level. 
The difference between FMO2/PCM[1] and FMO2/PCM[1(2)] is 
that in the former the ASCs are computed at the 1-body 
level, whereas for the former at the 2-body level, but 
without self-consistency (which would be FMO2/PCM[2]). 
Probably, FMO3/PCM[2] should be regarded as the most 
accurate and still affordable (with a few thousand nodes) 
method.  However, FMO3/PCM[1(2)] (specified with NBODY=3, 
IFMO=2 and NPCMIT=2) is much cheaper and slightly less 
accurate than FMO3/PCM[2].  FMO3/PCM[3] is the most 
accurate and expensive level of all.

How to perform FMO/EFP calculations

Solvent effects can also be taken into account with the 
Effective Fragment Potential model.  The presence of both 
$FMO and $EFRAG groups selects FMO/EFP calculations.  See 
the $EFRAG group and the $FMO group for details.

In the FMO/EFP method, the Quantum Mechanical part of the 
calculation in the usual EFP method is replaced by the FMO 
method, which may save time for large molecules such as 
proteins.

In the present version, only FMOn/EFP1 (water solvent only) 
is available for RHF, DFT and MP2.  One can use the MC 
global optimization technique for FMO/EFP by RUNTYP=GLOBOP. 
Of course, the group DDI (GDDI) parallelization technique 
for the FMO method can be used.

Geometry optimization or saddle point search for FMO

The standard optimizers in GAMESS are now well 
parallelized, and thus recommended to be used with FMO up 
to the limit hardwired in GAMESS (2000 atoms). In practice, 
if more than about 1000 atoms are present, numeric Hessian 
updates often result in the improper curvature and 
optimization stops. One can either do a restart, or use 
RUNTYP=OPTFMO (which does not diagonalize the Hessian).

RUNTYP=OPTIMIZE applies to Cartesian coordinates or DLC.  
RUNTYP=OPTFMO works only with Cartesian coordinates.  If 
your system has more than 2000 atoms you can consider 
RUNTYP=OPTFMO, which can now use Hessian updates and 
provides reasonable way to optimize although it is not as 
good as the standard means in RUNTYP=OPTIMIZE.

A transition state search for FMO can be performed with 
RUNTYP=SADPOINT using either Cartesian coordinates or DLC.

IRC calculations can be performed.



FMO hessian calculations

Analytic FMO Hessian with RUNTYP=HESSIAN may be computed 
for RHF, ROHF, UHF, RDFT, and UDFT in the gas phase (no 
PCM, EFP etc), provided that RESPCC is set to 0.

Molecular dynamics with FMO

MD can be run for any FMO method, which has the gradient 
implemented. However, in many cases the approximations in 
the gradient for a particular method may lead to large 
discrepancies in MD.  The following methods have a fully 
analytic gradient (which has to be turned on with $FMO 
keyword MODGRD=42):  FMO-RHF, FMO-MP2, FMO-RHF/EFP; the 
following condition should be satisfied: no ESP 
approximations, RESPPC=0.

Pair interaction energy decomposition analysis (PIEDA)

PIEDA can be performed for the PL0 and PL states.  The PL0 
state is the electronic state in which fragments are 
polarised by the environment in its free (noninteracting) 
state.  The simplest example is that in a water cluster, 
each molecule is computed in the electrostatic field 
exerted by the electron densities of free water molecules.  
The PL state is the FMO converged monomer state, that is, 
the state in which fragments are polarised by the self-
consistently converged environment. Namely, following the 
FMO prescription, fragments are recomputed in the external 
field, until the latter converges.  Using the PL0 state 
requires a series of separate runs; and it also relies on a 
"free state" which can be defined in many ways for 
molecules with detached covalent bonds.

What should be done to do the PL0 state analysis?
1. run FMO0.
This computes the free state for each fragment, and those 
electron densities are stored on file 30 (to be renamed 
file 40 and reused in step 3).
2. compute BDA energies (if detached bonds are present), 
using sample files in tools/fmo/pieda.  This corrects for 
artifacts of bond detaching, and involves running a model 
system like H3C-CH3, to amend for C-C bond detaching.
3. Using results of (1) and (2), one can do the PL0 
analysis.  In addition to pasting the data from the two 
punch files in steps 1,2 and the density file in step 1 
should be provided.

What should be done to do the PL state analysis?  The PL 
state itself does not need either the free state or PL0 
results. However, if the PL0 results are available, 
coupling terms can be computed, and in this case IPIEDA is 
set to 2; otherwise to 1.

So the easiest and frequently sufficient way to run PIEDA 
is to set IPIEDA=1 and do not provide any data from 
preceding calculations.  The result of a PIEDA calculation 
is a set of pair interaction energies (interfragment 
interaction energies), decomposed into electrostatic, 
exchange-repulsion, charge transfer and dispersion 
contributions.

Finally, PIEDA (especially for the PL state) can be thought 
of as FMO-EDA, EDA being the Kitaura-Morokuma decomposition 
(RUNTYP=MOROKUMA).  In fact, PIEDA (for the PL state) in 
the case of just two fragments of standalone molecules is 
entirely equivalent to EDA, which can be easily verified, 
by running the full PIEDA analysis (ipieda=2).  Note that 
PIEDA can be run as direct SCF, whereas EDA cannot be, and 
for large fragments PIEDA code can be used to perform EDA.  
Also, EDA in GAMESS has no explicit dispersion.

In 2012, PIEDA/PCM was developed describing the solvent 
screening.  RO-PIEDA based on RO-(HF, MP2 or CC) may be 
used for radicals.  Grimme's dispersion models may be used 
in PIEDA.

Excited states

At present, one can use CI, MCSCF, or TDDFT to compute 
excited states in FMO.  MCSCF is discussed separately 
above, so here only TDDFT and CI are explained.  They are 
enabled by setting the IEXCIT option (EXCIT(1) defines the 
excited state's fragment ID).

Two levels are implemented for TDDFT (FMO1-TDDFT and FMO2-
TDDFT).  In the former, only monomer TDDFT calculations are 
performed, whereas the latter adds pair corrections from 
TDDFT dimers.  PCM may be used for solvent effects with 
TDDFT (PCM[1] is usually sufficient).

CI can only be done for CIS at the monomer level (nbody=1), 
FMO1-CIS. The set-up for CI is similar to that for TDDFT.

Selective and sussystem FMO

Sometimes, one is interested only in some pair 
interactions, for example, between ligand and protein, or 
the opposite, only pair interactions within ligand. This 
saves a lot of CPU time by omitting all other pair 
calculations, but does not give the total properties. To 
use this feature, define MOLFRG and MODMOL. RUNTYP=ENERGY 
only is implemented.

In the subsystem analysis, one can divide fragments into                        
subsystems and obtain various properties of subsystems.

Frozen domain

To accelerate geometry optimisations, one can specify that 
the electronic state of the first layer in a 2-layer FMO 
can be computed at the initial geometry and consequently be 
frozen. One can define the polarizable buffer (equal to 
layer 2) and frozen domain (layer 1).  Fragments in the 
polarizable buffer which contain the atoms active in 
geometry optimisation form the active domain.  The 
fragments in the active domain should have a nonzero 
separation from the frozen domain. In FMO/FD all dimers in 
the polarizable buffer are computed; in FMO/FDD only those 
dimers which have at least one monomer in the active domain 
are computed. FMO/FD and FNI/FDD are only implemented for 
RUNTYP=OPTIMIZE.  MODFD and IACTAT in $FMO specify 
FMO/FD(D) atop of the usual multilayer FMO setup with some 
atoms frozen in geometry optimization by the standard means 
(i.e., IACTAT in $STATPT). Note that in FMO/FD(D) the 
Hessian as used in RUNTYP=OPTIMIZE is formed only for the 
atoms in the second layer, so this upper layer should not 
have more than the GAMESS limit (currently, 2000 atoms).

IMOMM with FMO

IMOMM (namely, SIMOMM) calculations can be performed with
the "MO" in IMOMM treated using FMO, i.e., this is like 
QM/MM but without electronic embedding of QM by MM.
This calculation uses Tinker, a plug-in source code, 
available from the GAMESS web site.  You should compile and 
link in the Tinker plug-in by changing a single line in 
comp/compall/lked,
  set TINKER=false
into
  set TINKER=true
In addition, you should change MAXATM=10 (the maximum 
number of atoms in the whole system, as used by Tinker) in 
several GAMESS source files into MAXATM=12000 (this number 
is used inside Tinker).  If you need a larger number, 
change it within Tinker as well.  After changing this, 
recompile and link GAMESS.

The input file style is in general like that of SIMOMM 
(q.v.).  Different from regular FMO, the atomic coordinates 
are given in $TINXYZ, not in $FMOXYZ. The fragmentation in 
FMO applies to QM atoms only, selected by IQMATM, and 
numbered consequently in FMO, so that INDAT in $FMO applies 
to the atoms renumbered from 1 (defined in IQMATM). Other 
than $FMOXYZ being superceded by $TINXYZ, the rest of FMO 
options is like in normal FMO.  IMOMM based on FMO is 
usually referred to as FMO/MM for short, although "FMO-
based SIMOMM" is probably easier to understand.  The 
somewhat tautological FMO-IMOMM has also been used by some.  
Covalent boundaries between FMO and MM are supported (via 
link atoms).  FMO/MM can be used to run geometry 
optimizations, whichis really what it is designed for.

Analyzing and visualizing the results

Annotated outputs provided in tools/fmo have matching 
mathematical formulae added onto the outputs, for easier 
reading.

Facio (http://www1.bbiq.jp/zzzfelis/Facio.html) can plot 
various FMO properties such as interaction energies, using 
interactive GUI viewers. 
To plot orbitals for an n-mer, set NPUNCH=2 in $SCF and 
PLTORB=.T.  There are several ways to produce cube files 
with electron densities.  They are described in detail in 
tools/fmo/fmocube/README.  To plot pair interaction maps, 
use tools/fmo/fmograbres to generate CSV files from GAMESS 
output, which can be easily read into Gnuplot or Excel.

FMO portal offers tools for visialising FMO results:
http://www.fmo-portal.info/ (shut down at present).

Parallelization of FMO runs with GDDI

The FMO method has been developed within a 2-level 
hierarchical parallelization scheme, group DDI (GDDI), 
allowing massively parallel calculations.  Different groups 
of processors handle the various monomer, dimer, and maybe 
trimer computations.  The processor groups should be sized 
so that GAMESS' innate scaling is fairly good, and the 
fragments should be mapped onto the processor groups in an 
intelligent fashion.

This is a very important and seemingly difficult issue. It 
is very common to be able to speed up parallel runs at 
least several times just by using GDDI better.  First of 
all, do not use plain DDI and always employ GDDI when 
running FMO calculations. Next, learn that you can and 
should divide nodes into groups to achieve better 
performance. The very basic rule of thumb is to try to have 
several times fewer groups than jobs. Since the number of 
monomers and dimers is different, group division should 
reflect that fact. Ergo, find a small parallel computer 
with 8-32 nodes and experiment changing just two numbers: 
ngrfmo(1)=N1,N2 and see how performance changes for your 
particular system.

Limitations of the FMO method in GAMESS

1. Dimensions: in general none, except that the standard 
GAMESS engines RUNTYP=OPTIMIZE and IRC are limited to 2000 
atoms (for FD(D), domain B may not exceed this limit). The 
limit can be increased by changing the source and 
recompiling GAMESS (see elsewhere).
2. CHARMM may not be combined with FMO, and some other 
extensions may not work.  Not every illegal combination is 
trapped, caveat emptor!
3. RUNTYP is limited to ENERGY, GRADIENT, OPTIMIZE, OPTFMO, 
IRC, FMO0, MD, GLOBOP, SADPOINT, FMOHESS and RAMAN only:
        Do not even try other ones!
4. Three-body FMO-MCSCF and FMO-TDDFT are not implemented.
5. No MOPAC semiempirical methods may be used, but DFTB was 
interfaced with FMO..

What will work the same way as ab initio:
The various SCF convergers, all DFT functionals, in-core 
integrals, direct SCF.

Restarts with the FMO method

RUNTYP=ENERGY can be restarted from anywhere before 
trimers. To restart monomer SCF, copy file F40 with monomer 
densities to the grandmaster node.  To restart dimers, 
provide file F40 and monomer energies ($FMOENM).  
Optionally, some dimer energies can be supplied ($FMOEND) 
to skip computation of corresponding dimers.

RUNTYP=GRADIENT can be easily restarted from monomer SCF 
(which really means it is a restart of RUNTYP=ENERGY, since 
gradient is computed at the end of this step).  Provide 
file F40. There is another restart option (1024 in $FMOPRP 
irest=), supporting full gradient restart, requiring, 
however, that you set this option in the original run 
(whose results you use to restart). To use this option, you 
would also need to keep (or save and restore) F38 files on 
each node (they are different).

RUNTYP=OPTIMIZE can be restarted from anywhere within the 
first RUNTYP=GRADIENT run (q.v.).  In addition, by 
replacing FMOXYZ group, one can restart at a different 
geometry.

RUNTYP=OPTFMO can be restarted by providing a new set of 
coordinates in $FMOXYZ and, optionally, by transferring 
$OPTRST from the punch into the input file.
 
Note on accuracy

The FMO method is aimed at computation of large molecules. 
This means that the total energy is large, for example, a 
6646 atom molecule has a total energy of -165,676 Hartrees. 
If one uses the standard accuracy of roughly 1e-9 (that 
should be taken relatively), one gets an error as much as 
0.001 hartree, in a single calculation.  FMO involves many 
ab initio single point calculations of fragments and their 
n-mers, thus it can be expected that numeric accuracy is 1-
2 orders lower than that given by 1e-9. Therefore, it is 
compulsory that accuracy should be raised, which is done by 
default.

The following default parameters are reset by FMO:
  ICUT/$CONTRL (9->12), ITOL/$CONTRL(20->24), 
  CONV/$SCF(1e-5 -> 1e-7), 
  CUTOFF/$MP2 (1e-9->1e-12), CUTTRF/$TRANS(1e-9->1e-10).
  CVGTOL/$DET,$GUGDIA (1e-5 -> 1e-6)
This to some extent slows down the calculation (perhaps on 
the order of 10-15%). It is suggested that you maintain 
this accuracy for all final energetics.  However, you may 
be able to drop the accuracy a bit for the initial part of 
geometry optimization if you are willing to do manual work 
of adjusting accuracy in the input.  It is recommended to 
keep high accuracy at the flat surfaces (the final part of 
optimizations) though.  For DFT the numeric grid's accuracy 
may be increased in accordance with the molecule size, e.g. 
extending the default grid of 96*12*24 to 96*20*40.  
However, some tests indicate that energy differences are 
quite insensitive to this increase.


FMO References

I. Basic FMO papers


A book chapter contains an introduction to FMO basics:
   Theoretical development of the fragment molecular
   orbital (FMO) method, D. G. Fedorov, K. Kitaura,
   in "Modern methods for theoretical physical chemistry of
   biopolymers", E. B. Starikov, J. P. Lewis, S. Tanaka,
   Eds., pp 3-38, Elsevier, Amsterdam, 2006.

There is now a full FMO book (11 chapters), which contains 
an introduction to FMO aimed at general application 
chemists, and a wealth of practical advice on doing FMO 
calculations:
   The Fragment Molecular Orbital Method: Practical
   Applications to Large Molecular System,
   D. G. Fedorov, K. Kitaura, Eds.,
   CRC Press, Boca Raton, FL, 2009.

FMO reviews:
   D. G. Fedorov, K. Kitaura  (Feature Article)
     J. Phys. Chem. A 111, 6904-6914 (2007).
   D. G. Fedorov, T. Nagata, K. Kitaura (Perspective)
     Phys. Chem. Chem. Phys., 14, 7562-7577 (2012)

A review of FMO in the context of other fragment-based 
methods is
    M. S. Gordon, D. G. Fedorov, S. R. Pruitt,
    L. V. Slipchenko Chem. Rev. 112, 632-672 (2012).

A very concise and detailed mathematical formulation of FMO 
including various extensions and property calculations is 
published as
   Mathematical formulation of the fragment molecular
   orbital method.
   T. Nagata, D. G. Fedorov, K. Kitaura.
   In "Linear-Scaling Techniques in Computational Chemistry
   and Physics". R. Zalesny, M. G. Papadopoulos,
   P. G. Mezey, J. Leszczynski, Eds., pp. 17-64,
   Springer, New York, 2011.


1. Fragment molecular orbital method: an approximate 
computational method for large molecules"
   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi
   Chem. Phys. Lett., 313, 701(1999).
2. Fragment molecular orbital method: application to 
polypeptides
   T. Nakano, T. Kaminuma, T. Sato, Y. Akiyama,
   M. Uebayasi, K. Kitaura Chem.Phys.Lett. 318, 614(2000).
3. Fragment molecular orbital method: analytical energy 
gradients
   K. Kitaura, S.-I. Sugiki, T. Nakano, Y. Komeiji, 
   M. Uebayasi, Chem. Phys. Lett., 336, 163(2001).
4. Fragment molecular orbital method: use of approximate 
electrostatic potential
   T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa,
   Y. Akiyama, M. Uebayasi, K. Kitaura
   Chem. Phys. Lett., 351, 475(2002).
5. The extension of the fragment molecular orbital method 
with the many-particle Green's function,
   K. Yasuda, D. Yamaki, J. Chem. Phys. 125, 154101(2006).
6. The role of the exchange in the embedding electrostatic 
potential for the fragment molecular orbital method.
   D. G. Fedorov, K. Kitaura
   J. Chem. Phys.  131, 171106(2009).
7. Analytic second derivatives of the energy in the 
fragment molecular orbital method.
   H. Nakata, T. Nagata, D. G. Fedorov, S. Yokojima, K.                         
   Kitaura, S. Nakamura, J. Chem. Phys. 138 (2013) 164103.

II. FMO in GAMESS

1. A new hierarchical parallelization scheme: generalized 
distributed data interface (GDDI), and an application to 
the fragment molecular orbital method (FMO).
   D. G. Fedorov, R. M. Olson, K. Kitaura, M. S. Gordon,
   S. Koseki  J. Comput. Chem.  25, 872-880(2004).
2. The importance of three-body terms in the fragment 
molecular orbital method.
   D. G. Fedorov and K. Kitaura
   J. Chem. Phys.  120, 6832-6840(2004).
3. On the accuracy of the 3-body fragment molecular orbital
method (FMO) applied to density functional theory
   D. G. Fedorov and K. Kitaura
   Chem. Phys. Lett.  389, 129-134(2004).
4. Second order Moeller-Plesset perturbation theory based 
upon the fragment molecular orbital method.
   D. G. Fedorov and K. Kitaura
   J. Chem. Phys. 121, 2483-2490(2004).
5. Multiconfiguration self-consistent-field theory based 
upon the fragment molecular orbital method.
   D. G. Fedorov and K. Kitaura
   J. Chem. Phys. 122, 054108/1-10(2005).
6. Multilayer Formulation of the Fragment Molecular Orbital 
Method (FMO).
    D. G. Fedorov, T. Ishida, K. Kitaura
    J. Phys. Chem. A. 109, 2638-2646(2005).
7. Coupled-cluster theory based upon the Fragment Molecular 
Orbital method.
    D. G. Fedorov, K. Kitaura
    J. Chem. Phys. 123, 134103/1-11 (2005)
8. The polarizable continuum model (PCM) interfaced with 
the fragment molecular orbital method (FMO).
    D. G. Fedorov, K. Kitaura, H. Li, J. H. Jensen,
    M. S. Gordon, J. Comput. Chem., 27, 976-985(2006)
9. The three-body fragment molecular orbital method for 
accurate calculations of large systems,
    D. G. Fedorov, K. Kitaura
    Chem. Phys. Lett. 433, 182-187(2006).
10. Pair interaction energy decomposition analysis,
    D. G. Fedorov, K. Kitaura
    J. Comp. Chem. 28, 222-237(2007).
11. On the accuracy of the three-body fragment molecular 
orbital method (FMO) applied to Moeller-Plesset 
perturbation theory,
    D. G. Fedorov, K. Ishimura, T. Ishida, K. Kitaura,
    P. Pulay, S. Nagase
    J. Comput. Chem., 28, 1476-1484 (2007).
12. The Fragment Molecular Orbital method for geometry 
optimizations of polypeptides and proteins,
    D.G.Fedorov, T. Ishida, M. Uebayasi, K. Kitaura
    J.Phys.Chem.A, 111, 2722-2732(2007).
13. Time-dependent density functional theory with the 
multilayer fragment molecular orbital method
    M. Chiba, D. G. Fedorov, K. Kitaura
    Chem. Phys. Lett. 444, 346-350 (2007).
14. Time-dependent density functional theory based upon the 
fragment molecular orbital method
     M. Chiba, D. G. Fedorov, K. Kitaura
     J. Chem. Phys. 127, 104108(2007).
15. Polarizable continuum model with the fragment molecular 
orbital-based time-dependent density functional theory.
     M. Chiba, D. G. Fedorov, K. Kitaura
     J. Comput. Chem. 29, 2667-2676 (2008).
16. Theoretical Analysis of the Intermolecular Interaction 
Effects on the Excitation Energy of Organic Pigments: Solid 
State Quinacridone. 
     H. Fukunaga, D.G.Fedorov, M. Chiba, K. Nii, K. Kitaura
     J. Phys. Chem. A 112, 10887-10894 (2008).
17. Covalent Bond Fragmentation Suitable To Describe Solids 
in the Fragment Molecular Orbital Method.
     D. G. Fedorov, J. H. Jensen, R. C. Deka, K. Kitaura
     J. Phys. Chem. A 112, 11808-11816 (2008).
18. Excited state geometry optimizations by time-dependent 
density functional theory based on the fragment molecular 
orbital method.
     M. Chiba, D. G. Fedorov, T. Nagata, K. Kitaura
     Chem. Phys. Lett. 474, 227-232 (2009).
19. Derivatives of the approximated electrostatic 
potentials in the fragment molecular orbital method.
     T. Nagata, D. G. Fedorov, K. Kitaura,
     Chem. Phys. Lett. 475, 124-131 (2009).
20. A combined effective fragment potential - fragment 
molecular orbital method. I. The energy expression and 
initial applications.
     T. Nagata, D. G. Fedorov, K. Kitaura, M. S. Gordon,
     J. Chem. Phys. 131, 024101 (2009).
21. Analytic gradient for the adaptive frozen orbital bond 
detachment in the fragment molecular orbital method.
     D. G. Fedorov, P. V. Avramov, J.H. Jensen, K. Kitaura,
     Chem. Phys. Lett. 477, 169-175 (2009).
22. Fragment molecular orbital study of the electronic 
excitations in the photosynthetic reaction center of 
Blastochloris viridis. 
    T. Ikegami, T. Ishida, D. G. Fedorov, K. Kitaura,
    Y. Inadomi, H. Umeda, M. Yokokawa, S. Sekiguchi,
    J. Comp. Chem. 31, 447-454 (2010).
23. Open-Shell Formulation of the Fragment Molecular 
Orbital Method.
     S. R. Pruitt, D. G. Fedorov, K. Kitaura, M. S. Gordon
     J. Chem. Theor. Comp.  6, 1-5 (2010)
24. Energy gradients in combined fragment molecular orbital 
and polarizable continuum model (FMO/PCM) calculation.
    H. Li, D. G. Fedorov, T. Nagata, K. Kitaura,
    J. H. Jensen, M. S. Gordon
    J. Comput. Chem. 31, 778-790 (2010).
25. Nuclear-Electronic Orbital Method within the Fragment 
Molecular Orbital Approach.
    B. Auer, M. V. Pak, S. Hammes-Schiffer,
    J. Phys. Chem. C 114, 5582-5588 (2010).
26. Importance of the hybrid orbital operator derivative 
term for the energy gradient in the fragment molecular 
orbital method. 
    T. Nagata, D. G. Fedorov, K. Kitaura, 
    Chem. Phys. Lett. 492, 302-308 (2010).
27. Systematic Study of the Embedding Potential Description 
in the Fragment Molecular Orbital Method. 
    D. G. Fedorov, L. V. Slipchenko, K. Kitaura,
    J. Phys. Chem. A 114, 8742-8753 (2010).
28. A combined effective fragment potential - fragment 
molecular orbital method. II. Analytic gradient and 
application to the geometry optimization of solvated 
tetraglycine and chignolin.
    T. Nagata, D. G. Fedorov, T. Sawada, K. Kitaura,
    M. S. Gordon, J. Chem. Phys. 134, 034110 (2011).
29. Geometry optimization of the active site of a large 
system with the fragment molecular orbital method.
    D. G. Fedorov, Y. Alexeev, K. Kitaura,
    J. Phys. Chem. Lett. 2, 282-288 (2011).
30. Fully analytic energy gradient in the fragment 
molecular orbital method.
    T. Nagata, K. Brorsen, D. G. Fedorov, K. Kitaura,
    M. S. Gordon, J. Chem. Phys. 134, 124115(2011).
31. Analytic energy gradient for second-order Moeller-
Plesset perturbation theory based on the fragment molecular 
orbital method.
    T. Nagata, D. G. Fedorov, K. Ishimura, K. Kitaura,
    J. Chem. Phys. 135, 044110 (2011).
32. Large-Scale MP2 Calculations on the Blue Gene 
Architecture Using the Fragment Molecular Orbital Method.
    G. D. Fletcher, D. G. Fedorov, S.R.Pruitt, T.L.Windus,
    M. S. Gordon, J. Chem. Theory Comput. 8, 75-79(2012).
33. Energy decomposition analysis in solution based on the 
fragment molecular orbital method.
    D.G.Fedorov, K.Kitaura
    J.Phys.Chem. A 116, 704-719(2012).
34. Analytic gradient and molecular dynamics simulations 
using the fragment molecular orbital method combined with 
effective potentials.
    T. Nagata, D. G. Fedorov, K. Kitaura
    Theor. Chem. Acc. 131, 1136 (2012).
35. Geometry Optimizations of Open-Shell Systems with the 
Fragment Molecular Orbital Method.
    S. R. Pruitt, D. G. Fedorov, M. S. Gordon,
    J. Phys. Chem. A, 116, 4965-4974 (2012).
36. Analytic gradient for second order Moeller-Plesset 
perturbation theory with the polarizable continuum model
based on the fragment molecular orbital method.
    T. Nagata, D. G. Fedorov, H. Li, K. Kitaura,
    J. Chem. Phys., 136, 204112 (2012).
37. Reducing scaling of the fragment molecular orbital 
method using the multipole method.
    C. H. Choi, D. G. Fedorov
    Chem. Phys. Lett. 543, 159-165(2012).
38. Unrestricted Hartree-Fock based on the fragment 
molecular orbital method: energy and its analytic gradient.
    H. Nakata, D. G. Fedorov, T. Nagata, S. Yokojima, K.
    Ogata, K. Kitaura, S. Nakamura
    J. Chem. Phys. 137, 044110 (2012).
39. Analytic gradient for the embedding potential with 
approximations in the fragment molecular orbital method.
    T. Nagata, D. G. Fedorov, K. Kitaura
    Chem. Phys. Lett. 544, 87-93 (2012).
40. Analysis of solute-solvent interactions in the fragment 
molecular orbital method interfaced with the effective 
fragment potentials: theory and application to solvated 
griffithsin-carbohydrate complex.
    T. Nagata, D. G. Fedorov, T. Sowada, K. Kitaura
    J. Phys. Chem. A, 116, 9088-9099 (2012).
41. Open-shell pair interaction energy decomposition 
analysis (PIEDA):  Formulation and application to the 
hydrogen abstraction in tripeptides.
    M.C.Green, D.G.Fedorov, K.Kitaura, J.S.Francisco,
    L.V.Slipchenko, J.Chem.Phys. 138 (2013) 074111.

Other FMO references including applications can be found 
at:
    http://staff.aist.go.jp/d.g.fedorov/fmo/main.html

EFMO references are given in its own subsection.




created on 7/7/2017