MOLCAS manual:

Next: 6.15 expbas Up: 6. Programs Previous: 6.13 embq


6.14 espf (+ QM/MM interface)

6.14.1 Description

The ElectroStatic Potential Fitted (ESPF) method adds contributions to the one-electron Hamiltonian for computing the interaction between the charge distribution in MOLCAS and any external electrostatic potential, field, and field derivatives. The approximate interaction energy writes:

\Delta E^{\mathrm{ESPF}} = \left [ \sum_a \left < \Psi\left \vert Q^a\right\vert \Psi \right > + Z_a \right ] V^a
\end{displaymath} (6.1)

with Qa a multipole-like operator which matrix elements are fitted to the electron potential integrals (determined on a grid surrounding the QM atoms) and Va the external electrostatic potential (and derivatives) at nucleus a. Both energy and gradient computations are available. A call to ESPF right after SEWARD is required to carry on such calculations. NOTE: always run SEWARD + ESPF. If not, very strange results may happen, due to interactions counted twice or more! NOTE: symmetry is ignored since the external potential usually breaks the one given in GATEWAY.

If no external potential is given, the ESPF module can be used to compute atomic point charges fitted to the electrostatic potential produced by the nuclei and electrons.

6.14.2 ESPF and QM/MM

Whereas the ESPF method can be used standalone, it has been developed for hybrid quantum mechanics/molecular mechanics (QM/MM) computations, in which an extended molecular system is divided into two subsystems: the 'active' center described by any QM method and its surroundings qualitatively treated with an empirical atomic forcefield. The current implementation uses a modified version of the Tinker program as MM code. In order to obtain the modified Tinker code, you must run the "molcas get_tinker" command.

The current patched version of Tinker6.1 is 6.3.2.

IMPORTANT: the environment variable TINKER must point to the directory in which the Tinker executable binaries are stored (usually in $MOLCAS/tinker/bin).

The most convenient way to define (i) the QM and MM subsystem and (ii) which atoms are to be known by MOLCAS (all the QM ones and some MM ones, see below) requires to simply add the keyword Tinker in GATEWAY. This way, GATEWAY will ask Tinker to pass him all informations needed.

Alternatively, an old input style can be used in GATEWAY, using the following syntax:

Basis  set
X.....  /  MM
  name  x  y  z  Angstrom
End  of  basis
where q gives its MM point charge value, name is its name and x, y, z are its coordinates given in au or in Å if the Angstrom keyword is given. This way, a MM atom is simply a QM atom without basis function and with a non-integer atomic charge. Usually only MM atoms needed to define some constrained link atoms (see below) positions are included, however this is not mandatory. The QM/MM method.

The Hamiltonian of the full QM/MM system is divided into three terms:


the first one describes the QM part as it would be in vacuo, the second one describes the surroundings using a classical MM forcefield and the last one deals with the interactions between the QM and the MM subsystems. In its usual formulation, the last term is (for q point charges interacting with N nuclei and n electrons):
\end{displaymath} (6.3)

The first two terms deal with the electrostatic interactions between the QM charge distribution and the MM electrostatic potential. In MOLCAS the ESPF method is used for this purpose. A short-range van der Waals term is added (van der Waals parameters are assigned to all the atoms - both QM and MM). If the frontier between the two subsystems involves a bond, some empirical bonded terms may also be used. For the sake of simplicity, the standard MM parameters are kept unchanged for the MM atoms but should be modified (or calculated) for the QM atoms (e.g. it may be necessary to fit the QM van der Waals parameters).

The usual forcefields use the "1-4 condition" to separate the bonded interactions (stretching, bending, torsion) from the non-bonded ones (electrostatic and vdw). This means than the non-bonded potentials are applied only if atoms are separated by 3 bonds or more. As for the QM/MM interactions, this procedure is kept with the exception that all the QM atoms experience the electrostatic potential generated by all the MM point charges (the QM/MM frontier case is considered later).

NOTE that, starting with MOLCAS 8, all the MM point charges interact with the QM charge distribution using the ESPF method (at variance with previous MOLCAS versions in which the few MM atoms defined in GATEWAY were interacting directly with the QM electrons and nuclei). Link atoms.

When no bond are involved between the QM and the MM parts, the QM/MM frontier definition is obvious and only the electrostatic and vdw interactions are taken into account. However, if one or several chemical bonds exist, the definition of a smooth but realistic frontier is needed. Several schemes, more or less sophisticated, have been proposed. In the current implementation, only the most basic one, the Link atom (LA) approach is included. In the LA approach, each QM/MM bond that should be cut is saturated with a monovalent atom - most often a hydrogen atom - on the QM side. The position of a link atom is often restrained: frozen distance from the corresponding QM frontier atom and always on the segment defined by the two frontier atoms (Morokuma's method select by the LAMorokuma keyword).

From the macromolecular point of view, link atoms do not exist, ie. they should not interact with the MM part. However, this leads to severe overpolarization of the frontier, due to unbalanced interactions. Hence interactions between the link atoms and the MM potential is kept. To remove problems that may arise from too strong interactions between a link atom and the closest MM point charges, these point charges may be spread in the MM neighborhood. For instance, in a protein, this procedure is mainly justified if the MM frontier atom is an $\alpha$ carbon (Amber or Charmm-typed forcefields usually set these point charges close to zero). Geometry optimizations - microiterations.

In a QM/MM geometry optimization job, a MOLCAS step costs as hundreds of TINKER steps. Thus it is very convenient to use the microiterations technique, that is converging the MM subsystem geometry every MOLCAS step. This can be requested directly within the TINKER keywords file. However, in order to improve the optimization convergence, an improved QM/MM hessian can be built in SLAPAF using the rHidden keyword (note that adding the keyword Cartesian may help too).

6.14.3 Dependencies

The ESPF program depends on SEWARD for modifying the core Hamiltonian matrix and on ALASKA for computing the extra contributions to the gradient.

6.14.4 Files

ESPF will use the following input files: RYSRW, ABDATA,RUNFILE, ONEINT (for more information see [*]). In addition, ESPF uses ESPFINP (the ESPF input file) and SEWARINP (the Seward input file).

Please note that the external potential can be given within a file, separated from the ESPF input file. Intermediate files

All the intermediate files are related to the use of ESPF together with a MM code (i.e. Tinker) and allow for communication between the two programs. MOLCAS uses one file to pass the QM atoms coordinates and ESPF-derived point charges to Tinker. Tinker uses the same file to pass the external potential, the MM-only energy and gradient components to MOLCAS.

TINKER.LOGThe log file of the Tinker run
$Project.xyzThe coordinates file for Tinker.
$Project.keyThe keywords file for Tinker.
$Project.qmmmThe communication file between MOLCAS and Tinker. Output files

ONEINTOne-electron integral file generated by the SEWARD program.
RUNFILECommunication file for subsequent programs.
ESPF.DATAAscii file containing some specific informations needed for subsequent calls to the ESPF module.

6.14.5 Input

Below follows a description of the input to ESPF.

In addition to the keywords and the comment lines the input may contain blank lines. The input for each module is preceded by its name like:


Compulsory keywords
ExternalSpecify how the external potential is given. This keyword is compulsory in the first run of ESPF. On the next line, one integer or a text string must be given:
  • One integer n is given. If n is 0, the next lines give the numbering, the values for the external potential, the field and field gradients for each atom. If n is greater than 0, the n next lines specify the sources of the external potential, each line gives three cartesian coordinates, one point charge, and (optionally) three dipole components. If Å is used as the length unit, the Angstrom keyword must be given right after n.
  • The NONE word means that no external potential is given. Accordingly, the ESPF module will compute the atomic point charges (and optionally dipoles) deriving from the electrostatic potential due to all electrons and nuclei.
  • The word is TINKER, it means that the current job is a QM/MM job using the MOLCAS/Tinker interface. Accordingly the external potential will be computed directly by Tinker. Note Tinker requires at least two input files, ending with .xyz (coordinates) and .key (keywords). These files must share the name of the current MOLCAS project. Optionally, you can add the MULLIKEN or LOPROP keyword after TINKER: it indicates what kind of charges are passed to Tinker. These charges may be used during the MM microiterations. If no keyword is given, the ESPF multipoles are selected.
  • Any other word. The following characters up to the next space are taken as a file name and the rest of the line is ignored. Instead, the full input (including the first line) is read from the specified file and must follow the syntax specified above.

Optional keywords
TitleTitle of the job.
MultipoleOrderGive the multipolar order of the ESPF operators. Only 0 (charge-like) or 1 (charge- and dipole-like) are allowed and should be written on the next input line. Default value is 0.
GridModify the grid specifications. The grid is made of points belonging to molecular surfaces defined according to the van der Waals radii of each quantum atom. Two schemes are available. The first one is the GEPOL procedure, as implemented into the PCM SCRF method. The other one is called PNT and is the default. On the next line, first select the method with the GEPOL or PNT option. On the same line, one integer number and one real number are given if PNT is selected. The first one gives the maximum number of shells around the van der Waals surface of the quantum atoms. The second one gives the distance between the shells. Note that all points within the van der Waals envelope are discarded to avoid the penetration effects. Default values are 4 shells separated by 1 Å. Alternatively, if GEPOL is selected, the same line must contain 1 integer indicating the number of surfaces to be computed (must be < 6).
ShowRequires the printing of the ESPF.DATA file.
LAMorokumaActivate the Morokuma's scheme for scaling the link atom positions (QM/MM calculation). Note the scaling factor is currently hard-coded and is actually determined from the radii of the atoms involved in the QM/MM frontier bond.

6.14.6 Examples ESPF example

This is a typical input for the calculation of the energy and the gradient of glycine molecule feeling the external potential of 209 TIP3P water molecules.

Basis  set
  C1  1.11820  0.72542  -2.75821  Angstrom
  C2  1.20948  0.66728  -1.25125  Angstrom
End  of  basis
Basis  set
  O1  2.19794  1.10343  -0.67629  Angstrom
End  of  basis
Basis  set
  H1  2.02325  1.18861  -3.14886  Angstrom
  H2  0.25129  1.31794  -3.04374  Angstrom
  H3  1.02458  -0.28460  -3.15222  Angstrom
End  of  basis
Basis  set
  N1  0.17609  0.12714  -0.61129  Angstrom
End  of  basis
Basis  set
  C3  0.09389  -0.01123  0.84259  Angstrom
  C4  -1.21244  -0.67109  1.28727  Angstrom
End  of  basis
Basis  set
  O2  -2.06502  -1.02710  0.48964  Angstrom
End  of  basis
Basis  set
  H4  -0.61006  -0.21446  -1.14521  Angstrom
  H5  0.92981  -0.61562  1.19497  Angstrom
  H6  0.16338  0.97444  1.30285  Angstrom
End  of  basis
Basis  set
  N2  -1.41884  -0.85884  2.57374  Angstrom
End  of  basis
Basis  set
  H7  -0.73630  -0.57661  3.25250  Angstrom
  H8  -2.28943  -1.29548  2.82140  Angstrom
End  of  basis


MultipoleOrder  =  0
External  =  0
1  -0.048  -0.002  -0.006  -0.001  0.007  -0.009  0.002  -0.001  0.001  -0.001
2  -0.047  -0.002  0.001  -0.002  0.003  0.000  -0.004  0.000  -0.001  0.000
3  -0.053  0.004  0.000  -0.011  0.002  0.002  -0.004  0.002  0.003  -0.007
4  -0.046  0.011  -0.009  -0.001  0.006  -0.005  -0.001  0.003  0.003  -0.004
5  -0.042  -0.016  -0.011  -0.006  0.005  -0.007  0.003  -0.004  -0.001  -0.005
6  -0.050  0.000  0.008  0.001  0.006  -0.006  0.000  -0.002  0.000  -0.001
7  -0.039  -0.008  0.001  0.000  0.001  -0.002  0.001  -0.001  -0.001  -0.001
8  -0.032  -0.007  -0.002  0.004  0.002  -0.003  0.001  -0.002  0.002  -0.001
9  -0.011  -0.009  0.004  0.001  0.002  0.000  -0.002  -0.001  0.001  0.001
10  0.000  -0.011  0.003  0.004  0.001  0.002  -0.003  0.001  -0.001  0.001
11  -0.028  -0.008  0.004  -0.001  -0.001  -0.002  0.002  -0.001  0.001  -0.002
12  -0.026  0.003  -0.008  0.014  0.002  -0.001  -0.001  -0.008  0.006  -0.009
13  -0.037  -0.008  -0.003  0.004  -0.007  0.007  0.000  0.001  0.007  -0.001
14  -0.016  -0.007  0.007  -0.008  0.003  0.003  -0.006  0.000  0.002  0.002
15  -0.025  0.003  0.012  -0.007  0.003  -0.001  -0.002  -0.006  0.005  0.009
16  -0.010  -0.011  0.000  -0.014  0.001  0.007  -0.008  0.001  0.000  -0.001

Charge  =  0

  &alaska QM/MM example

A typical start for a QM/MM calculation is given in the following input. It is quite general since all the informations related to the QM and MM subsystems definitions are already included into the Tinker key file.

>  EXPORT  TINKER=$MOLCAS/tinker/bin_qmmm
>  COPY  $PATH_TO/$  $WorkDir/$
>  COPY  $PATH_TO/$Project.key  $WorkDir/$Project.key

Basis  =  STO-3G
Group  =  Nosym


External  =  Tinker

This can be used, eg. with the following Tinker files. In this example, the asparate anion is cut into two pieces, the QM subsystem contains the end of the side-chain until the $\beta$ carbon atom. There is a link atom between the QM $\beta$ and MM $\alpha$ carbon atoms.

16  ASP
 1 N3    -0.040452    0.189961    0.173219   448     2     6    14    15
 2 CT    -0.011045   -0.060807    1.622395   449     1     3     7    11
 3 C      1.446535   -0.110535    2.028518   450     2     4     5
 4 O      1.902105    0.960982    2.409042   452     3
 5 O      2.137861   -0.898168    1.387158   452     3
 6 H      0.559257   -0.496270   -0.262338   451     1
 7 CT    -0.789906   -1.336520    1.982558   216     2     8    12    13
 8 C     -2.256402   -1.184505    1.571038   218     7     9    10
 9 O2    -2.460769   -0.949098    0.356151   219     8
10 O2    -3.120135   -1.188969    2.465678   219     8
11 H1    -0.478878    0.773493    2.145163   453     2
12 HC    -0.356094   -2.194944    1.466324   217     7
13 HC    -0.720511   -1.505463    3.058628   217     7
14 H     -0.996208    0.061130   -0.151911   451     1
15 H      0.304306    1.116522   -0.018698   451     1
16 HLA   -0.283317   -0.506767    1.748300  2999     2     7


parameters $PATH_TO_TINKER/params/amber99.prm
QM -8 10 7 12 13
MM 2
LA 16
* Add the atom type for the LA
atom   2999    99    HLA     "Hydrogen Link Atom"        1      1.008     0
charge -2  0.0
charge -11 0.0

next up previous contents
Next: 6.15 expbas Up: 6. Programs Previous: 6.13 embq