QMMM keyblock options¶

Introduction¶

Chapter 2 explains setting up an ADF QM/MM simulation. This section describes the available options that you can define in the QMMM key block. This section is essentially a reference source. The main components of the QMMM key block, the MM_CONNECTION_TABLE, the FORCE_FIELD_FILE, and the LINK_BONDS key blocks that are described in detail in the previous chapter are repeated here (with some additional notes). Please note that the global optimization options are not well tested and are prone to crashing the run. If one is interested in using these options, please be aware of this fact. We appreciate reports of any failures.

Example Input¶

In this section we simply provide a few examples of the QM/MM key block. In some examples, the MM_CONNECTION_TABLE and LINK_BONDS subkey blocks are not filled.

Example 3.1 This example depicts a global optimization of the MM region with the simulated annealing-like optimizer available in the QM/MM program. The global search involves 100 ps of MD at 1000 K with 100 structures sampled in regular intervals during the simulation. Each of the 100 structures is then partially optimized, and then the 10 best are fully optimized. At the end of this, the lowest energy structure is used for the QM/MM run. Note that in this MD search, the QM atoms, including the link atoms are frozen.

QMMM
  FORCE_FIELD_FILE sybyl.ff
  OUTPUT_LEVEL 1
  WARNING_LEVEL 1
  ELSTAT_COUPLING_MODEL 1
  OPTIMIZE
    GLOBAL
      METHOD MD_SEARCH
      FREQUENCY ONCE
    SUBEND
  SUBEND
  MD_SEARCH
    TIME{PS} 100.0
    N_STRUCTURES 100
    TEMPERATURE 1000.0
  SUBEND
  MM_CONNECTION_TABLE
  SUBEND

  LINK_BONDS
  SUBEND
END

Example 3.2 In this example, custom charges are assigned to some of the atoms. Charges for atoms that were not given specific charges in the QMMM key block are assigned on a per atom-type basis from the force field file. Also note that this example has no LINK_BONDS subkey block. This is only allowed if there are not link bonds, as in the example in Figure 1-2a.

QMMM
 FORCE_FIELD_FILE sybyl.ff
 OUTPUT_LEVEL 1
 WARNING_LEVEL 1
 ELSTAT_COUPLING_MODEL 1
 MM_CONNECTION_TABLE
 SUBEND

 CHARGES
   0.4
   0.3
   3 -0.1
 SUBEND
END

Description of Options¶

FORCE_FIELD_FILE

Keyword (required, default = amber95.ff) This keyword simply defines the full path of the force field file to be used for the molecular mechanics potential. The location of the force field file is given after the keyword. The full path can be give, or just the file name. In the latter case, the program checks the current directory that ADF is executing in. Examples:

FORCE_FIELD_FILE /home/username/sybyl.ff

FORCE_FIELD_FILE sybyl.ff

NEWQMMM

Keyword (Default is not to use the subkey NEWQMMM) Key to be used for more efficient QM/MM calculations, work in progress. It also allows more QM/MM atoms than in a default calculation. This key should be used ONLY with the amber force field. This key also offers the possibility to use the new QM/MM input format, which can, for example, be made with the utility pdb2adf. The old input format remains working if one includes this NEWQMMM subkey.

OUTPUT_LEVEL

Keyword (Default = 1) The integer following this keyword specifies the amount of output to be printed to the ADF output file. 0: minimal output 1: normal output 2: troubleshooting output OUTPUT_LEVEL 2 is recommending for initially setting up a job. However, once the job is set up properly this output level is probably too verbose.

WARNING_LEVEL

Keyword (Default = 1) The performs some checking of the input, ranging from examining all interatomic distances to examining the input order of the QM, LI and MM atoms. The integer following this keyword specifies how many warnings to report and when to stop the run due to the warning. -1: Report only the most severe warnings, and never stop the run. Useful when user is knowingly violating the ‘rules’. 0: Report severe warnings and only stop at ‘fatal’ errors. 1: Report all warnings, stop at severe and fatal errors. This is the default. 2: Report all warnings and stop at any of them. Useful when initially setting up a job.

MDC_LEVEL

Keyword (Default = 2) The integer following this keyword specifies the level of the Multipole Derived Charge analysis [7] used in conjunction with ELSTAT_COUPLING_MODEL=4. Note that the default value has been changed compared to ADF2007.01 1: MDC-m charges are used to update charges of QM system. 2: MDC-d charges used. 3: MDC-q charges used.

MM_CONNECTION_TABLE

Subkey block (required) This key block defines the connection table, the force field atom types and the partitioning of the full system into QM and MM regions. It is critical that the atoms specified in this key block are in the same order as in the ATOMS key block. This is important, because the program may not detect this type of input error and you would get ridiculous results.

MM_CONNECTION_TABLE
 n FF_LABEL MM_TYPE connection numbers
 ...
END

The labels are defined in the following table

Input column
1	N	atom number
2	FF_LABEL	Force field atom type.
		These labels correspond to the atom types defined in the force field file.
		They can be up to four characters long. Xx defines dummy atoms.
3	MM_TYPE	QM, MM or LI
4-	connection	These define to which atoms the current atom has a covalent bond.
	numbers	These connections are used to generate the molecular mechanics potential.
		Currently, a maximum of 6 connections is allowed per atom.

The connection table should be a fully redundant one. In other words, if atom #1 is bonded to atom #5, they each should have the other atom listed in their connections.

Example:

1 C_2 QM 2 3 4 5
2 O_2 QM 1
3 H QM 1
4 C_3 QM 1 5 6 7
5 Cu QM 4 1
6 H QM 4

A fully non-redundant connection table is also supported. In such a connection table, once a bond is mentioned, it is not mentioned again. In other words, the connection list for any atom cannot contain an atom, which precedes it in the atom numbering.

Example:

1 C_2 QM 2 3 4 5
2 O_2 QM
3 H QM
4 C_3 QM 5 6 7
5 Cu QM
6 H QM

These two connection tables are equivalent. Connection tables that are semi-redundant might cause problems. We recommend the fully redundant connection table.

LINK_BONDS

Required for systems with LINK bonds This key block required for systems with covalent bonds that cross the QM/MM boundary. These bonds are referred to in this document as the link-bonds. Each link bond has a constant parameter ‘alpha‘ associated with it, which is defined as the ratio of the bond length in the real system and of the capping bond in the model QM system. See Section 1 for more details. To determine the alpha parameters for each link bond, one can take the capping atom bond distance in a ‘pure QM’ calculation of the QM model system and ratio it to the corresponding bond distance in the real system. These ratios are typically around 1.30 to 1.50 when hydrogen as a capping atom. The LINK_BONDS subkey block has the following format:

LINK_BONDS
 atom_a - atom_b alpha replacement_fragment addremove_force_field_type
SUBEND

Example:

LINK_BONDS
  15 - 3 1.42 H H1
  8 - 1 1.40 Cl.dzp Cl
SUBEND

The integers atom_a and atom_b refer to the numbering of the two atoms involved in the link bond. One of the atoms will be a LI type atom whereas the other will be a QM type atom. Atom_a and atom_b must be separated by ” - ” with at least one space between the integer and the hyphen. In other words ‘3 - 4’ is correct, but not ‘3- 4’ ‘3 -4’. Atoms need not be in any particular order, and the order of the link bonds is also not important. Following this is the alpha parameter for that specific bond. The replacement_fragment is the ADF atom used for the capping atom. Often the capping atom is a hydrogen atom, however, it need not be. The replacement_fragment must be present in the FRAGMENT key block in the ADF input file. The addremove_force_field_type need only be present for the AddRemove model [3], and indicates the force field type of the capping atom (similar to FF_LABEL in the MM_CONNECTION_TABLE block). Important note: It is very important to realize that the Hamiltonian depends on the a parameters used. When comparing relative energies for example, one needs to take care that the a’s corresponding to the same bonds are identical.

CHARGES

Key block (optional)** This key block defines the initial charges on each atom based on their atom **number. Atom numbers must be carefully specified because the program does not assume any order. Charges can also be assigned based on their atom type from the force field file in the CHARGE PARAMETERS key block. If this key block does not specify charges, the program looks for charge assignments from the force field file. If charges are not assigned in either this key block or the force field file, then a charge of 0.0 is assigned. When polarizable electrostatic coupling is invoked, the charges for QM and LI atoms are not read, because the MM point charges interact with the QM charge distribution. With ELSTAT_COUPLING_MODEL=4, these charges (for the QM system) are used only in the first cycle of the geometry optimization: after each cycle, the QM charges are replaced with the Multipole Derived Charges. Example:

CHARGES
 3 -0.10000
 1 0.05000
 2 0.05000
SUBEND

ELSTAT_COUPLING_MODEL

Keyword (optional, default=1) This keyword controls the type of electrostatic model that is used, including whether true electrostatic coupling between the QM and MM regions is evoked.

Coupling Model	Description
0	Electrostatics OFF
1	Simple electrostatic coupling, where there is no polarization of the QM wave function, i.e. pure MM coupling.
4	Simple electrostatic coupling, like option 1. But the point charges of the QM region are updated throughout the geometry optimization using the Multipole Derived Charge analysis [7](MDC-x level depending on MDC_LEVEL).

OPTIMIZE

Key block (optional) This key block allows the user to modify the geometry optimization settings. An example of a key block with many available options is shown below:

OPTIMIZE
 MAX_STEPS 1000
 MAX_GRAD 0.001
 PRINT_CYCLES 20
 METHOD BFGS

 GLOBAL
 METHOD GRID
 FREQUENCY ONCE
 SUBEND

 GRID
 INCREMENT 20.0
 BOND 2 - 4
 BOND 2 - 3
 SUBEND
SUBEND

Sub-options to this key are described next.

OPTIMIZE: MAX_STEPS

Keyword (optional, default = 1000) This keyword defines the maximum number of optimization steps allowed before the optimization is discontinued.

OPTIMIZE: MAX_GRADIENT

Keyword (optional, default = (0.01 kcal/mol)/Angstrom) This keyword allows the user to change the convergence criteria. For now, the optimization is considered converged when the maximum gradient on any MM atom is less than MAX_GRADIENT. The default value will provide gradients that are very small, especially when compared to the convergence criteria specified most electronic structure codes. NOTE: The gradients on the QM atoms due to the MM potentials are not accounted for in the convergence criteria. Lrge MM forces can exist on the QM atoms after the optimization.

OPTIMIZE: ENERCVG

Keyword (optional, default = 0.001 kcal/mol) This keyword allows the user to change the convergence criteria for the energy between successive cycles. Can only be used in case of NEWQMMM.

OPTIMIZE: METHOD

Keyword (optional, default = BFGS, available: BFGS, STEEPEST_DESCENT, CONJGRAD) For the most part, the default quasi-Newton optimizer with BFGS Hessian update scheme is very stable, and converges well. Other optimizers available are the steepest descent method (STEEPEST_DESCENT) and conjugate gradient (CONJGRAD). STEEPEST_DESCENT and CONJGRAD are almost always less efficient than the BFGS optimizer (particularly close to the minimum). It is notable that the Hessian based BFGS method requires more memory than the STEEPEST_DESCENT method and so for very large systems may be problematic to use. In that case, it is best to use the CONJGRAD method.

OPTIMIZE: MM_NOTCONVERGED

Keyword (optional, default = 1) This keyword defines what should happen if the MM geometry is not fully optimized after MAX_STEPS steps; set this to zero for large (biochemical) systems where it may be problematic to get the optimization to converge fully in a limited number of steps (1000): the QMMM run will continue as if the MM-optimization had converged.

OPTIMIZE: FIX_MM_GEOMETRY

Keyword (optional, default = .false.) If this keyword is specified in the OPTIMIZE subblock, the MM system will be frozen, i.e. no geometry optimization will be done on any of its atoms. If NEWQMMM is included use ‘IRUNTYPE_QMMM 0’ as separate keyword (outside the OPTMIZE subblock).

OPTIMIZE: PRINT_CYCLES

Keyword (optional, default = 100) PRINT_CYCLES represents the number of optimization cycles between which the optimization status is printed and the MM restart file is written.

OPTIMIZE: GLOBAL

Sub key block (optional) - CURRENTLY IN A BETA STATE This subkey block controls the global optimization options in the program. Currently the global optimization option has not been thoroughly tested and should be considered to be in beta form. The normal optimizers are designed only to locate the “nearest” local minimum and therefore you are not guaranteed to find the best overall structure, which is termed the global minimum structure. Currently, only two global optimization algorithms have been implemented: - Molecular dynamics based optimizer related to a simulated annealing algorithm - A grid search, which generates conformations by rotations about bonds specified by the user. Both optimizers generate a number of structures (100s to 1000s), which are all partially optimized. The partially optimized structures are then sorted based on their energies. The best 10 of these structures are then fully optimized. The best of these fully optimized structures is kept and assumed to be the global minimum. The global optimization is only applied to the MM region with the QM atoms frozen. Therefore, the structure can only be considered the global minimum structure on the constrained surface where the QM atoms and QM charge density are frozen.

OPTIMIZE
  GLOBAL
    METHOD MD_SEARCH
    FREQUENCY ONCE
  SUBEND
  MD_SEARCH
    TIME{PS} 100.0
    N_STRUCTURES 100
    TEMPERATURE 1000.0
  SUBEND
SUBEND

Global optimization is not the default. Therefore, to invoke a global optimization, the GLOBAL subkey block must exist. It is important to note that the subkey blocks that control the global optimization schemes are subkey blocks of the OPTIMIZE key block and not sub-sub key blocks within the GLOBAL subkey block. The above example demonstrates this.

OPTIMIZE: GLOBAL: METHOD

Keyword (optional, default = MD_SEARCH) This key block specifies the global optimization method to be used. To date there are only two methods, MD_SEARCH which is the default and GRID. More detail one how these methods work is given in the description of the MD_SEARCH and GRID subkey blocks.

OPTIMIZE: GLOBAL: FREQUENCY

Keyword (optional, default = ONCE) This key block specifies how often the global optimization algorithm, if it is specified, is called. Since the global optimization is very time consuming it is not recommended that it be used every QM iteration. The default is that it is done only on the first iteration. The options available are tabulated below.

ONCE	Only at the first iteration
EVERY_TIME	At every iteration
N_CYCLES X	At each X-th iteration including the first.
	Here X is the integer following the “N_CYLES” keyword.
	e.g. N_CYCLES 4

OPTIMIZE: MD_SEARCH

Subkey block (optional, default settings specified below) The MD_SEARCH method involves performing molecular dynamics on the MM subsystem at a high temperature. The high temperature dynamics allows the MM subsystem to “get out of” the local minimum of the initial structure and explore other regions phase space, potentially leading to lower energy structures. During the molecular dynamics, structures are sampled at specified intervals and stored. When the dynamics is complete, the stored structures are optimized and sorted in terms of their energy. This procedure is similar to simulated annealing, except that the temperature of the dynamics is not ramped up and down in a cyclic fashion. At the beginning, the dynamics is immediately pulsed up to the specified temperature with a random excitation on each of the free MM degrees of freedom. An example of the key block, with good settings is given below.

MD_SEARCH
 TIME{PS} 100.0
 N_STRUCTURES 100
 TEMPERATURE 700.0
SUBEND

In the above example, the MM subsystems are heated up to a temperature of 700 Kelvin. Dynamics is run for a total of 100 picoseconds, with a total of 101 structures sampled (100 plus the initial structure). Each structure is sampled every 1.0 picoseconds. The default timestep is 0.5 femtoseconds, and therefore in the above example 200,000 timesteps will be performed. This global search technique is the most general and robust of the two methods implemented. It is therefore the default global optimization method. This subkey block is optional, since the default settings should work reasonably with most systems.

OPTIMIZE: GRID

Subkey block (optional, required if method selected) The GRID method provides a systematic search for global minimum by rotating about specified covalent bonds in the MM subsystem. This method is only efficient for small systems or systems where the conformational variability is confined to torsions involving a few bonds. The user must specify the bonds that are to be rotated in the search, up to a maximum of 10, and the increment (in degrees) by which the bonds are to be rotated between subsequent structures. The program does not allow bonds that are completely within the QM subsystem (link bonds are allowed, however, or part of a ring system. Finally, since QM atoms cannot be rotated, at least one of the two fragments resulting from splitting the specified bond must contain no QM atoms. An example of the key block is shown below where three bonds are rotated, in 60° increments. 216 structures (6x6x6) will be generated corresponding to a full 360° rotation about the three bonds in 60° increments and all combinations thereof. :: GRID INCREMENT 60.0 BOND 7 - 6 BOND 8 - 7 BOND 9 - 8 SUBEND NOTE: It is important to realize that the program uses the connection table specified in the input to determine which atoms to rotate. MASSES Key block (optional) This is used to assign custom masses to individual atoms. If no custom masses are specified, then the default masses defined in the force field file are used. Below is example input.

MASSES
 15 32.066
 8 2.0
SUBEND

The first column is an integer specifying the atom number and the second column is a real specifying the custom mass of that atom in atomic mass units. The atoms need not be in any particular order and it is not necessary to specify custom masses for all atoms. It should be noted that only masses of the MM atoms and the link atoms could be customized. Masses of the QM atoms and the capping atoms are taken from the QM borderleft.

PARTITION

Keyword (Default = 5) Parameter for partitioning QM and MM systems (capping method). Relevant if LINK_BONDS is present. Possible values: 5 (corresponding to the AddRemove scheme), 3 (corresponding to the IMOMM scheme)

Electrostatic Embedding¶

The ELSTAT_COUPLING_MODEL=2 option can be used for both the SCF (bonding and total) energies and gradients and for the TDDFT energies and gradients. The MM point charges are considered during the Becke grid setup and are affected by the QPNEAR keyword of the BECKEGRID block. Symmetry should be explicitly set to NOSYM, because the MM atoms are not considered when determining the symmetry and actually might break it. In other ELSTAT_COUPLING_MODEL modes, the QMMM driver is usually called twice: once during geometry setup and once before the SCF. In the electrostatic embedding mode, the QMMM driver may also be called after the SCF finishes and after the TDDFT gradient was computed because these steps provide a new electron density necessary to compute the complete QMMM energies/forces.

The excited state geometry optimization (with EXCITEDGO) is possible, although a few restrictions apply due to the fact that QM and MM atoms are optimized separately by different parts of the program. If only the QM atoms should be optimized, the MM atoms can be kept frozen with the MM optimizer method SKIP:

QMMM
   ...
   OPTIMIZE
      METHOD SKIP
      ...
   SUBEND
END

This option is also useful when only the energy gradient needs to be computed without actually moving the atoms (e.g., external optimizer/dynamics driver). If only MM atoms should be optimized, do not use the GEOMETRY/EXCITEDGO keywords but use only the QMMM/OPTIMIZE. If the QM and MM atoms should be optimized simultaneously, the internal workflow is as follows:

1. compute the energy, energy gradient, and the density matrix for the given geometry;
2. optimize the MM atoms using the density matrix from step 1;
3. compute a quasi-Newton step for the QM atoms using the energy gradient from step 1;
4. repeat steps 1-3 until convergence.

This means that the QM electron density is constant during the MM optimization, so if the MM atoms have moved far away from their original positions the QM gradient may be not completely consistent anymore. Therefore it is advisable to not take too many steps in the MM optimization part:

QMMM
   ...
   OPTIMIZE
      MAX_STEPS 5
      MM_NOTCONVERGED 0
      ...
   SUBEND
END

The MM optimization might still take relatively many iterations to converge because it is Hessian-free. Hence, it might be worthwhile to preoptimize the MM part without electrostatic embedding first.