6.15. ORCA MM Module¶
Since version 4.2 ORCA features its own independent MM implementation.
The minimum input necessary for a MM treatment looks as follows.
!MM
%mm
ORCAFFFilename "UBQ.ORCAFF.prms"
end
In this section we discuss the basic keywords and options, i.e.
the basic structure of the ORCA Forcefield File,
how to generate the ORCA Forcefield File,
how to manipulate the ORCA Forcefield File,
how to speed up MM calculations,
further MM options and keywords.
Further options important for QM/MM calculations will be discussed in section ORCA Multiscale Implementation.
6.15.1. ORCA Forcefield File¶
For the MM part of the QM/MM calculation force-field parameters are necessary. ORCA has its own parameter file format (ORCA forcefield file - ORCAFF.prms), which includes the atom specific parameters for nonbonded interactions:
partial charge
LJ coefficients
and parameters for bonded interactions:
bonds
angles
Urey-Bradley terms
dihedrals
impropers
CMAP terms (cross-terms for backbone, currently not used)
Individual parameters, like e.g. atomic charge, equilibrium bond length and force constant, …, can be conveniently modified directly within the ORCA Forcefield File.
6.15.1.1. How to generate the ORCA Forcefield File¶
The easiest way to generate a ORCAFF.prms file is currently to convert from psf (protein structure file) files. Psf files are specific to the CHARMM forcefield and its application via NAMD. Psf files for a specific protein system can easily be generated by the popular molecular visualization program VMD and its extension QwikMD, but also with other extensions in the VMD program (e.g. psfgen or fftk). The psf file contains information on the atom types and on the bonded interactions of all atoms. It does, however, not contain the parameters that belong to these interactions. These parameteres are stored in specific files, often ending with prm, but also par or str. The CHARMM parameter files come with VMD installation, can be directly downloaded, or can be generated with the VMD extension fftk (forcefield toolkit).
Once a ORCAFF.prms file is available, it can be manipulated, i.e. split up into several parts for individual molecules, new ORCAFF.prms files can be generated for non-standard molecules, and individual ORCAFF.prms files can be merged, as described in the following:
6.15.1.1.1. Conversion from psf or prmtop files to ORCAFF.prms: convff¶
The orca_mm module can convert psf and prm files (CHARMM), prmtop files (AMBER) or xml files (open force field from the openff toolkit, compatible to AMBER) to the ORCAFF file with the -convff flag. Input options are:
orca_mm -convff <optional:-verbose> <FFInput> <PSFFILE> <PRMFILE(S)>
Keywords:
<FFInput> = -CHARMM
<FFInput> = -AMBER
<FFInput> = -OPENMM
For CHARMM topologies, when a psf file is available for a system with standard residues, prepared by e.g. QwikMD, psfgen or other vmd tools, the conversion needs the psf plus the prm files as input:
CHARMM example:
orca_mm -convff -CHARMM 1C1E.psf par_all36_prot.prm toppar_water_ions_namd.str
ORCA can also convert Amber topologies to the ORCAFF file. Here, only the prmtop file is required:
AMBER example:
orca_mm -convff -AMBER complex.prmtop
ORCA can also convert xml files from the openff toolkit (AMBER compatible) to the ORCAFF file. Here, only the xml file is required:
OPENFF example:
orca_mm -convff -OPENMM complex.xml
6.15.1.1.2. Divide a forcefield file: splitff¶
If a ORCAFF.prms file should be subdivided into several files, e.g. if the psf file stems from QWikMD with non-standard molecules included, e.g. a ligand. In that case first the parameters of the ligand are split from the remaining system, next the ligand needs to be protonated, then a simple ORCAFF.prms file is generated via orca_mm’s makeff option, and finally the ligand’s new ORCAFF.prms file is added to the main systems file via the above described mergeff option. Note that the file can only be split into files for nonbonded fragments.
Input options:
orca_mm -splitff <optional:-verbose> <ORCAFFFILE> <A1> <optional:A2> ...
Keywords:
<ORCAFFFILE> = ORCA forcefield file.
<A1> = Atom number of first atom of fragment that should belong
to a new ORCA forcefield file
<A2> = Atom number of first atom of fragment that should belong
to a new ORCA forcefield file
... = More split atoms possible
Note that atoms start counting at 1.
Example:
orca_mm -splitff 1C1E_substrate_noH.ORCAFF.prms 7208
6.15.1.1.3. Merge forcefield files: mergeff¶
If several ORCAFF.prms files are available and should be merged for an ORCA calculation, e.g. for a standard plus a non-standard molecule.
Input options:
orca_mm -mergeff <optional:-verbose> <ORCAFFFILE1> <ORCAFFFILE2> ...
Keywords:
<ORCAFFFILE1> = First ORCA forcefield file
<ORCAFFFILE2> = Second ORCA forcefield file
... = More ORCA forcefield files possible
Example:
orca_mm -mergeff 1C1E.ORCAFF.prms substrate_withH.ORCAFF.prms
6.15.1.1.4. Repeat forcefield files: repeatff¶
In case the same ORCAFF.prms file needs to be repeated multiple times, the repeatff functionality is available.
Input options:
orca_mm -repeatff <optional:-verbose> <ORCAFFFILE> <A>
<ORCAFFFILE> = ORCA forcefield file.
<A> = Factor (integer) defining how often this forcefield file should be repeated.
Example:
orca_mm -mergeff methanol.ORCAFF.prms 580
This feature can be useful e.g. in the case of solvating a molecule, i.e. adding multiple copies of a solvent to a solute. First the solvent can be repeated N times, and subsequently the solute’s prms file can be merged together with the solvent prms file.
6.15.1.1.5. Divide a forcefield file: splitpdb¶
When splitting a ORCAFF.prms file, also splitting of the pdb file is required. The file can be split into an arbitrary number of individual files.
This can be useful together with the splitff command.
Input options:
orca_mm -splitpdb <optional:-verbose> <PDBFILE> <A1> <optional:A2> ...
Keywords:
<PDBFILE> = PDB file.
<A1> = Atom number of first atom of fragment that should belong
to a new ORCA forcefield file
<A2> = Atom number of first atom of fragment that should belong
to a new ORCA forcefield file
... = More split atoms possible
Note that atoms start counting at 1.
Example:
orca_mm -splitpdb 1C1E_substrate_noH.pdb 7208
6.15.1.1.6. Merge PDB files: mergepdb¶
If several PDB files are available and need to be merged for an ORCA calculation, e.g. a protein and a ligand or multiple ligands, or a ligand that was first removed from a complex, then modified, and finally should get back into the complex PDB file.
This can be useful together with the mergeff command.
Input options:
orca_mm -mergepdb <optional:-verbose> <1PDBFILE> <2PDBFILE> ...
Keywords:
<1PDBFILE> = First PDB file
<2PDBFILE> = Second PDB file
... = More PDB files possible
Example:
orca_mm -mergepdb 1C1E.pdb substrate_withH.pdb
6.15.1.1.7. Create simple force field: makeff¶
The orca_mm tool can generate an approximate forcefield for a molecule, storing it in ORCAFF.prms format. Here, the LJ coefficients are based on UFF parameters and the partial charges are based on a simple PBE or XTB calculation. The resulting topology is certainly not as accurate as an original CHARMM topology, but can still be used for an approximate handling of the molecule. Herewith, the molecule can be part of the QM region (having at least the necessary LJ coefficients), or part of the MM region as a non-active spectator - being not too close to the region of interest. In the latter case it is important that the molecule is not active, since bonded parameters are not available. However, it can still be optimized as a rigid body, see sections Geometry Optimizations using the L-BFGS optimizer and the usage in QM/MM calculations in section Optimization with the Cartesian L-BFGS Minimizer, on MM level, optimizing its position and orientation with respect to the specific environment.
Input options:
orca_mm -makeff <optional:-verbose> <XYZ/PDBFILE> <optional:-C CHARGE>
<optional:-M MULT> <optional:-nproc N> <optional:-CHARGE_OPTIONS>
Keywords:
<CHARGE> = charge of system
<MULT> = multiplicity of system
<-nproc N> = number of processors (Default 1)
<-CHARGE_OPTIONS> = Structure Charge calculation
<-PBE> input PBE
<-PBEOpt> PBE opt. PBE
<-PBEOptH> PBE H-opt. PBE
<-XTB> input XTB
<-XTBOpt> XTB opt. XTB
<-XTBOptH> XTB H-opt. XTB
<-XTBOptPBE> XTB opt. PBE
<-noChargeCalc> distribute net charge evenly
PBE Opt and SP level: RI-PBE/def2-SVP CPCM(Water), CHELPG charges
XTB Opt and SP level: GFN2-XTB GBSA(Water), Mulliken partial charges
Example:
orca_mm -makeff substrate_withH.xyz -M 2 -XTBOptPBE
Note that ORCA generates bonds based on simple distance rules, which enables ORCA to detect where to add link atoms between QM and MM atoms, see also section QM-MM, QM-QM2 and QM2-MM Boundary. But the user is advised to treat a molecule, for which the ORCAFF.prms file was generated with the makeff option, either fully in the QM, or fully in the MM region, unless the charge distribution has been properly taken care of (due to the need of integer charges in QM and MM system).
6.15.1.1.8. Get standard hydrogen bond lengths: getHDist¶
For the definition of the link atoms standard bond lengths between C, N and O and hydrogen are directly set by ORCA but can be modified by the user, see section QM-MM, QM-QM2 and QM2-MM Boundary. If other atom types are on the QM side of the QM-MM boundary, their distance to the link atom has to be defined. In this case a file can be provided to ORCA which defines the standard bond length to hydrogen for all possible atoms. Such a file can be generated via the following command:
Input options:
orca_mm -getHDist <optional:-verbose> <XYZ/PDBFILE>
Example:
orca_mm -getHDist 1C1E.xyz
This file can then be modified, the required values can be added, and the resulting file can be defined as input for the QMMM calculation.
6.15.1.1.9. Create ORCAFF.prms file for IONIC-CRYSTAL-QMMM¶
For IONIC-CRYSTAL-QMMM calculations, section IONIC-CRYSTAL-QMMM, an ORCAFF.prms file with initial charges and connectivities is required. If you are not using the orca_crystalprep tool for setting up such calculations, see section orca_crystalprep, you can directly prepare the ORCAFF.prms file with the command:
orca_mm -makeff <XYZ/PDBFILE> -CEL <ELEMENT1> <OXIDATION_STATE1>
-CEL <ELEMENT2> <OXIDATION_STATE2>
Keywords:
<ELEMENT1> = element name
<OXIDATION_STATE1> = formal oxidation state of element 1
... = More elements possible
Example:
orca_mm -makeff na4cl4.xyz -CEL Na 1.0 -CEL Cl -1.0
Here, na4cl4.xyz is the supercell structure file (it can contain tens of thousands of atoms).
Note
For supercells with more complex oxidation states’, e.g. Co\(_3\)O\(_4\), the ORCAFF.prms file can be generated conveniently via the orca_crystalprep tool, orca_crystalprep.
6.15.2. Speeding Up Nonbonded Interaction Calculation¶
For MM calculations of very large systems with hundreds of thousands of atoms, and for QMMM calculations with fast QM methods (e.g. XTB, AM1) and / or very small QM systems, the computation of the nonbonded interactions can become a bottleneck. Different schemes for speeding up the calculation of nonbonded interactions are available within the ORCA MM implementation. Two schemes truncate long-range interaction, another scheme can be used for calculations with active regions, i.e. calculations where only a part of the system is active or optimized. For more information on active regions see section Active and Non-Active Atoms - Optimization, Frequency Calculation, Molecular Dynamics and Rigid MM Water.
6.15.2.1. Force Switching for LJ Interaction¶
With force switching for the LJ interaction (described in reference [818]) a smooth switching function is used to truncate the LJ potential energy smoothly at the outer cutoff distance LJCutOffOuter. If switching is set to false, the LJ interaction is not truncated at LJCutOffOuter. The parameter LJCutOffInner specifies the distance at which the switching function starts taking an effect to bring the van der Waals potential to 0 smoothly at the LJCutOffOuter distance, ensuring that the force goes down to zero at LJCutOffOuter, without introducing discontinuities. Note that LJCutOffInner must always be smaller than LJCutOffOuter.
%mm
SwitchForceLJ true # Use the switch force scheme for the LJ interaction.
# Default: true.
LJCutOffInner 10. # Distance (in Ang). Default: 10 Ang.
LJCutOffOuter 12. # Distance (in Ang). Default: 12 Ang.
end
6.15.2.2. Force Shifting for Electrostatic Interaction¶
With force shifting for the electrostatic interaction (described in reference [818]) the electrostatic potential is shifted to zero at the cutoff distance CoulombCutOff. If shifting is set to false, the electrostatic interaction is not truncated at CoulombCutOff.
%mm
ShiftForceCoulomb true # Use the shift force scheme for the Coulomb interaction.
# Default: true.
CoulombCutOff 12. # Distance (in Ang). Default: 12 Ang.
end
6.15.2.3. Neglecting Nonbonded Interactions Within Non-Active Region¶
When using active regions (see section Active and Non-Active Atoms - Optimization, Frequency Calculation, Molecular Dynamics and Rigid MM Water) for optimizations and MD runs, the nonbonded interactions at the MM level can be neglected for those atom pairs, which are both non-active, without loss of accuracy for the results. Even relative energies between two structures are correct, if the atom positions of the non-active atoms are identical. For all other cases, i.e. if the positions of atoms in the non-active region differ, the full nonbonded interaction should be computed in the final single-point energy calculation. By default this option is switched off.
%mm
Do_NB_For_Fixed_Fixed true # Compute MM-MM nonbonded interaction also for
# non-active atom pairs. Default true.
end
6.15.3. Rigid Water¶
As TIP3P water might have to be treated as rigid bodies due to its parametrization - please check out the specifics of the applied force field parametrization - we offer a keyword for the automated rigid treatment of all active MM waters. The following keyword applies bond and angle constraints to active MM waters in optimization as well as MD runs:
%mm
Rigid_MM_Water false # Default: false.
end
6.15.4. Available Keywords for the MM module¶
Here we list all keywords that are accessible from within the mm block and that are relevant to MM, but also QM/MM calculations. Some of the MM keywords can also be accessed via the qmmm block, see section Additional Keywords.
!MM # or QMMM, as the MM keywords will also affect the MM part of the QMMM calculation
%mm
# Schemes for the truncation of long-range
# Coulomb and LJ interaction:
# The Shift Force scheme for the Coulomb interaction shifts the Coulomb potential
# such that it becomes zero at the cutoff distance.
ShiftForceCoulomb true # Use the shift force scheme for the Coulomb interaction.
# Default: true.
CoulombCutOff 12. # Distance (in Ang). Default: 12 Ang.
# With the Switch Force scheme for the LJ interaction is unchanged up to
# LJCutOffInner. Between LJCutOffInner and LJCutOffOuter a smooth switching function
# is applied onto the LJ potential so that the force goes down to zero at
# LJCutOffOuter, without introducing discontinuities.
SwitchForceLJ true # Use the switch force scheme for the LJ interaction.
# Default: true.
LJCutOffInner 10. # Distance (in Ang). Default: 10 Ang.
LJCutOffOuter 12. # Distance (in Ang). Default: 12 Ang.
DielecConst 1. # dielectric constant used in calc. of electrostatic
# interaction. Default: 1.
Coulomb14Scaling 1. # Scaling factor for electrostatic interaction between
# 1,4-bonded atoms. Default: 1.
PrintLevel 1 # Printing options: Can be 0 to 4, 0=nothing, 1=normal, ...
# Keywords that can be accessed from the mm as well as the qmmm block.
# For a description see qmmm block.
# Information about topology and force field
ORCAFFFilename "UBQ.ORCAFF.prms"# If available, e.g. from a previous run, or after
# modification, the ORCA Forcefield can be provided.
# Computing MM nonbonded interactions within non-active region.
Do_NB_For_Fixed_Fixed true # Compute MM-MM nonbonded interaction also for
# non-active atom pairs. Default true.
# Optimization and MD of active MM waters
RIGID_MM_WATER false # Default: false.
# Extended active region
ExtendActiveRegion distance # rule to choose the atoms belonging to activeRegionExt.
# no - do not use activeRegionExt atoms
# cov_bonds - add only atoms bonded covalently to
# active atoms
# distance (default) - use a distance criterion (VDW
# distance plus Dist_AtomsAroundOpt)
Dist_AtomsAroundOpt 1. # in Angstrom. Default 1 Ang.
OptRegion_FixedAtoms { 2 9} end # Default: empty list.
# The following keywords will affect the behavior of MM (without QMMM) calculations,
# but have to be provided via the qmmm block
PrintOptRegion true # Additionally print xyz and trj for opt region
PrintOptRegionExt true # Additionally print xyz and trj for extended opt region
PrintQMRegion true # Additionally print xyz and trj for QM region
PrintPDB true # Additionally print pdb file for entire system, is
# updated every iteration for optimization
end
*pdbfile 0 1 ubq.pdb # structure input via pdb file, but also possible via xyz file
6.16. ORCA Multiscale Implementation¶
With ORCA 5.0 ORCA ‘s multiscale functionality has been extensively expanded. ORCA 5 features five different multiscale methods for
proteins, DNA, large molecules, explicit solvation:
additive QMMM (Additive QMMM)
subtractive QM1/QM2 methods (2-layered ONIOM) (Subtractive QM/QM2 Method)
QM1/QM2/MM methods (3-layered ONIOM) (QM/QM2/MM Method)
CRYSTAL-QMMM for crystals:
MOL-CRYSTAL-QMMM for molecular crystals (MOL-CRYSTAL-QMMM)
IONIC-CRYSTAL-QMMM for semiconductors and insulators (IONIC-CRYSTAL-QMMM)
The multiscale features are optimally connected to all other modules and tools available in ORCA allowing the user to handle multiscale calculations from a QM-centric perspective in a simple and efficient way, with a focus on simplifying the process to prepare, set up and run multiscale calculations.
From the input side all methods share a common set of concepts and keywords, which will be outlined in the first part of this chapter. In the subsequent parts of this chapter, the different methods are described and further input options are discussed.
6.16.1. General Settings and Input Structure¶
Some of the keywords in this section are common to all five multiscale features, and some are not. If keywords are not available for one of the multiscale features, this will be mentioned.
6.16.1.1. Overview on Combining Multiscale Features with other ORCA Features¶
The multiscale features can be used together with all other possible ORCA methods:
- Single Point Calculations
Use all kinds of available electronic structure methods as QM method.
- Optimization
Use all kinds of geometry optimizations with all kinds of constraints, TS optimization, relaxed surface scans, and the ScanTS feature. Use the L-Opt and L-OptH features including the combination of all kinds of fragment optimizations (fix fragments, relax fragments, relax only specific elements in fragments, treat a fragment as a rigid body).
- Transition States and Minimum Energy Paths
Use all kinds of Nudged-Elastic Band calculations (Fast-NEB-TS, NEB, NEB-CI, NEB-TS, including their ZOOM variants) and Intrinsic Reaction Coordinate calculations. (not implemented for MOL-CRYSTAL-QMMM and IONIC-CRYSTAL-QMMM)
- Frequency Calculations
Use regular frequency calculations. If required, ORCA automatically switches on the Partial Hessian Vibrational Analysis (PHVA) calculation. (not tested for IONIC-CRYSTAL-QMMM)
- Molecular Dynamics
Use the Molecular Dynamics (MD) module for Born-Oppenheimer MD (BOMD) with QM/MM in combination with all kinds of electronic structure methods. (not implemented for MOL-CRYSTAL-QMMM and IONIC-CRYSTAL-QMMM)
- Property Calculation
All kinds of regular property calculations are available. For electrostatic embedding the electron density is automatically perturbed by the surrounding point charges.
- Excited State Calculations
Use all kinds of excited state calculations (TD-DFT, EOM, single point calculations, optimizations, frequencies). (For the ONIOM calculations the low-level calculations are carried out in the ground state)
6.16.1.2. Overview on Basic Aspects of the Multiscale Feature¶
In the following, the basic concepts are introduced.
- QM atoms
The user can define the QM region either directly, or via flags in a pdb file. See QM Atoms.
- QM2 atoms
Only applicable for QM/QM2/MM. For the QM/QM2/MM method the low level QM region (QM2) is defined via the input or via flags in a pdb file. See QM2 Atoms. For QM/QM2 the low level region consists of all atoms but the QM atoms.
- Active atoms
The user can choose an active region, e.g. for geometry optimizations the atoms that are optimized, for a frequency calculation the atoms that are allowed to vibrate for the PHVA, or for an MD run the atoms that are propagated. See Active and Non-Active Atoms - Optimization, Frequency Calculation, Molecular Dynamics and Rigid MM Water and Fig. 6.57.
- Forcefield
ORCA has its own forcefield file format (stored in files called basename.ORCAFF.prms). For a convenient setup the orca_mm module offers the option to convert from other forcefield formats. Currently, the following formats can be converted to the ORCA forcefield file format:
- CHARMM psf files
protein structure file from the CHARMM forcefield. The psf files can be easily prepared with the popular molecular visualizer VMD, together with its extensions (psfgen, QwikMD, fftk, which works together with ORCA ).
- AMBER prmtop files
topology files from the AMBER force field. Tutorials on how to generate AMBER prmtop files (for standard and non-standard molecules) can be found here.
- Open Force Field
xml files from the openforcefield initiative. With the openff-toolkit xml topology files (compatible with the AMBER force field) can be easily generated for small and medium-sized non-standard molecules. For a tutorial see here.
- Simple forcefield for small to medium-sized molecules
Alternatively, the orca_mm module can generate a simple approximate ORCAFF.prms file. For more options, see ORCA Forcefield File.
This concept has the following advantages:
- Modification of forcefield parameters
Atom and bond specific parameters can be easily modified within the ORCA forcefield file, allowing the user maximum flexibility in modifying the forcefield, which might be particularly useful for chemical systems like transition metal complexes. See ORCA Forcefield File.
- Standard and Non-Standard Ligands
Ligands can be easily and flexibly exchanged or added to the system, see ORCA Forcefield File.
- Boundary Treatment
ORCA automatically detects QM-MM boundaries, i.e. bonds that have to be cut between QM and MM region. ORCA automatically generates the link atoms and keeps them at their relative position throughout the run, even allowing to optimize the bond along the boundary. See QM-MM, QM-QM2 and QM2-MM Boundary. Not applicable for CRYSTAL-QMMM.
- Treatment of overpolarization
ORCA automatically adapts the charges at the QM-MM boundary. See QM-MM, QM-QM2 and QM2-MM Boundary. Not applicable for both CRYSTAL-QMMM.
- Embedding types
The electrostatic and mechanical embedding schemes are available. See Embedding Types.
Detailed information on the different available input and runtime options and additionally available keywords (see Additional Keywords) are given below.
6.16.1.3. QM Atoms¶
QM atoms can be defined either directly
!QMMM
%qmmm
QMAtoms {0 1 2 27 28} end
end
or via the occupancy column of a pdb file.
%qmmm # use either
QMAtoms {0:4} end # 1. list of atoms (counting starts from 0) or
Use_QM_InfoFromPDB true # 2. get the definition from the pdb file. Default false.
end # If (2) is set to true, (1) is ignored
*pdbfile 0 1 ubq.pdb
If Use_QM_InfoFromPDB is set to true, a pdb file should be used for
the structural input. QM atoms are defined via 1 in the occupancy
column, MM atoms via 0. QM2 atoms (for QM/QM2/MM calculations, see
QM2 Atoms) can be defined via 2 in the occupancy
column. The IONIC-CRYSTAL-QMMM method can have even further entries in
the PDB file, see Different QM and MM regions Stored in the PDB file. Note that the Use_QM_InfoFromPDB
keyword needs to be written before the coordinate section.
ubq.pdb:
...
ATOM 327 N ASP A 21 29.599 18.599 9.828 0.00 0.00 P1 N
ATOM 328 HN ASP A 21 29.168 19.310 9.279 0.00 0.00 P1 H
ATOM 329 CA ASP A 21 30.796 19.083 10.566 0.00 0.00 P1 C
ATOM 330 HA ASP A 21 31.577 18.340 10.448 0.00 0.00 P1 H
ATOM 331 CB ASP A 21 31.155 20.515 10.048 2.00 0.00 P1 C
ATOM 332 HB1 ASP A 21 30.220 21.082 9.865 2.00 0.00 P1 H
ATOM 333 HB2 ASP A 21 31.754 21.064 10.801 2.00 0.00 P1 H
ATOM 334 CG ASP A 21 31.923 20.436 8.755 1.00 0.00 P1 C
ATOM 335 OD1 ASP A 21 32.493 19.374 8.456 1.00 0.00 P1 O
ATOM 336 OD2 ASP A 21 31.838 21.402 7.968 1.00 0.00 P1 O
ATOM 337 C ASP A 21 30.491 19.162 12.040 0.00 0.00 P1 C
ATOM 338 O ASP A 21 29.367 19.523 12.441 0.00 0.00 P1 O
...
Note that contrary to the hybrid36 standard of PDB files, ORCA writes non-standard pdb files as:
atoms 1-99,999 in decimal numbers
atoms 100,000 and beyond in hexadecimal numbers, with atom 100,000 corresponding to index 186a0.
This ensures a unique mapping of indices. If you want to select an atom with an idex in hexadecimal space (index larger than 100,000), convert the hexadecimal number into decimals when choosing this index in the ORCA input file. Note also, that in the pdb file, counting starts from 1, while in ORCA counting starts from zero.
6.16.1.4. Active and Non-Active Atoms - Optimization, Frequency Calculation, Molecular Dynamics and Rigid MM Water¶
The systems of multiscale calculations can become quite large with tens and hundreds of thousands of atoms. In multiscale calculations the region of interest is often only a particular part of the system, and it is common practice to restrict the optimization to a small part of the system, which we call the active part of the system. Usually this active part consists of hundreds of atoms, and is defined as the QM region plus a layer around the QM region. The same definition holds for frequency calculations, in particular since after optimization non-active atoms are not at stationary points, and a frequency calculation would lead to artifacts in such a scenario. MD calculations on systems with hundreds of thousands of atoms are not problematic, but there are applications where a separation in active and non-active parts can be important (e.g. a solute in a solvent droplet, with the outer shell of the solvent frozen).
Note
If no active atoms are defined, the entire system is treated as active.
The active region definitions also apply to MM calculations, but have to be provided via the qmmm block.
6.16.1.4.1. Input Format¶
Active atoms can be defined either directly or via the B-factor column of a pdb file.
%qmmm # use either
ActiveAtoms {0:5 16 21:30} end # 1. list of atoms (counting starts from 0) or
Use_Active_InfoFromPDB true # 2. get the definition from the pdb file.
# Default false.
end # If (2) is set to true, (1) is ignored
*pdbfile 0 1 ubq.pdb
If Use_Active_InfoFromPDB is set to true, a pdb file should be used
for the structural input. Active atoms are defined via 1 in the B-factor
column, non-active atoms via 0. Note that the Use_Active_InfoFromPDB
keyword needs to be written before the coordinate section.
ubq.pdb:
...
ATOM 327 N ASP A 21 29.599 18.599 9.828 0.00 0.00 P1 N
ATOM 328 HN ASP A 21 29.168 19.310 9.279 0.00 0.00 P1 H
ATOM 329 CA ASP A 21 30.796 19.083 10.566 0.00 1.00 P1 C
ATOM 330 HA ASP A 21 31.577 18.340 10.448 0.00 1.00 P1 H
ATOM 331 CB ASP A 21 31.155 20.515 10.048 1.00 1.00 P1 C
ATOM 332 HB1 ASP A 21 30.220 21.082 9.865 1.00 1.00 P1 H
ATOM 333 HB2 ASP A 21 31.754 21.064 10.801 1.00 1.00 P1 H
ATOM 334 CG ASP A 21 31.923 20.436 8.755 1.00 1.00 P1 C
ATOM 335 OD1 ASP A 21 32.493 19.374 8.456 1.00 1.00 P1 O
ATOM 336 OD2 ASP A 21 31.838 21.402 7.968 1.00 1.00 P1 O
ATOM 337 C ASP A 21 30.491 19.162 12.040 0.00 0.00 P1 C
ATOM 338 O ASP A 21 29.367 19.523 12.441 0.00 0.00 P1 O
...
Note that in the above example also the QM atoms are defined along with the active atoms.
6.16.1.4.2. Optimization in redundant internal coordinates¶
In ORCA’s QM/MM geometry optimization only the positions of the active atoms are optimized. The forces on these active atoms are nevertheless influenced by the interactions with the non-active surrounding atoms. In order to get a smooth optimization convergence for quasi-Newton optimization algorithms in internal coordinates, it is necessary that the Hessian values between the active atoms and the directly surrounding non-active atoms are available. For that reason the active atoms are extended by a shell of surrounding non-active atoms which are also included in the geometry optimization, but whose positions are constrained, see Fig. 6.57. This shell of atoms can be automatically chosen by ORCA. There are three options available:
- Distance
(Default) The parameter
Dist_AtomsAroundOptcontrols which non-active atoms are included in the extension shell, i.e. non-active atoms that have a distance of less than the sum of their VDW radii plusDist_AtomsAroundOptare included.- Covalent bonds
All (non-active) atoms that are covalently bonded to active atoms are included.
- No
No non-active atoms are included.
The user can also provide the atoms for the extension shell manually. This will be discussed in section Frequency Calculation.
Fig. 6.57 Active and non-active atoms. Additionally shown is the extension shell, which consists of non-active atoms close in distance to the active atoms. The extension shell is used for optimization in internal coordinates and PHVA.¶
The set of active atoms is called the ‘activeRegion’, and the set of active atoms plus the surrounding non-active atoms is called ‘activeRegionExt’. During geometry optimization the following trajectories are stored (which can be switched off):
- basename_trj.xyz
Entire QM/MM system
- basename.QMonly_trj.xyz
Only QM region
- basename.activeRegion_trj.xyz
Only active atoms
- basename.activeRegionExt_trj.xyz
Active atoms plus extension
The following files are stored containing the optimized structures - if requested:
- basename.pdb
Optimized QM/MM system in pdb file format
- basename.xyz
Optimized QM/MM system
- basename.QMonly.xyz
Only QM region
- basename.activeRegion.xyz
Only active atoms
- basename.activeRegionExt.xyz
Active atoms plus extension
6.16.1.4.3. Optimization with the Cartesian L-BFGS Minimizer¶
For very large active regions the quasi-Newton optimization in internal coordinates can become costly and it can be advantageous to use the L-Opt or L-OptH feature, see section Geometry Optimizations using the L-BFGS optimizer. For the L-Opt(H) feature there exist two ways to define the active region:
via the ActiveAtoms keyword (or the Use_Active_InfoFromPDB flag) or
via fragment definition and the different keywords for fragment optimization. Available options are:
- FixFrags
Freeze the coordinates of all atoms of the specified fragments.
- RelaxHFrags
Relax the hydrogen atoms of the specified fragments. Default for all atoms if !L-OptH is defined.
- RelaxFrags
Relax all atoms of the specified fragments. Default for all atoms if !L-Opt is defined.
- RigidFrags
Treat each specified fragment as a rigid body, but relax the position and orientation of these rigid bodies.
Note
The L-Opt(H) option together with the fragment optimization can be used in order to quickly preoptimize your system at MM level. E.g. you can optimize the hydrogen positions of the protein and water molecules, and at the same time relax non-standard molecules, for which no exact bonded parameters are available, as rigid bodies.
!MM L-OptH
%mm
ORCAFFFilename "DNA_hexamer.ORCAFF.prms"
end
*pdbfile 0 1 protein_ligand.pdb
%geom
Frags # all other atoms belong to fragment 1 by default
2 {22307:22396} end # cofactor
3 {22397:22423} end # ligand
end
RigidFrags {2 3} end # treat cofactor and ligand as individual rigid bodies
end
6.16.1.4.4. Frequency Calculation¶
If all atoms are active, a regular frequency calculation is carried out when requesting !NumFreq. If there are also non-active atoms in the QM/MM system, the Partial Hessian Vibrational Analysis (PHVA, see section Partial Hessian Vibrational Analysis) is automatically selected. Here, the PHVA is carried out for the above defined activeRegionExt, where the extension shell atoms are automatically used as ‘frozen’ atoms. Note that the analytic Hessian is not available due to the missing analytic second derivatives for the MM calculation. Note that in a new calculation after an optimization it might happen that the new automatically generated extended active region is different compared to the previous region before optimization. This means that when using a previously computed Hessian (e.g. at the end of an optimization or a NEB-TS run) the Hessian does not fit to the new extended active region. ORCA tries to solve this problem by fetching the information on the extended region from the hess file. If that does not work (e.g. if you distort the geometry after the Hessian calculation) you should manually provide the list of atoms of the extended active region. This is done via the following keyword:
%qmmm
OptRegion_FixedAtoms {27 1288:1290 4400} end # manually define the
# activeRegionExt atoms.
end
Note that ORCA did print the necessary information in the output of the calculation in which the Hessian was computed:
Fixed atoms used in optimizer ... 27 1288 1289 1290 4400
6.16.1.4.5. Nudged Elastic Band Calculations¶
NEB calculations (section Nudged Elastic Band Method) can be carried out in combination with the multiscale features in order to e.g. study enzyme catalysis. When automatically building the extension shell at the start of a Multiscale-NEB calculation, not only the coordinates of the main input structure (‘reactant’), but also the atomic coordinates of the ‘product’ and, if available, of the ‘transition state guess’ are used to determine the union of the extension shell of the active region. For large systems it is advised to use the active region feature for the NEB calculation. Note that the atomic positions of the non-active atoms of reactant and product and, if available, transition state guess, should be identical.
6.16.1.4.6. Molecular Dynamics¶
If there are active and non-active atoms in the multiscale system, only the active atoms are allowed to propagate in the MD run. If all atoms are active, all atoms are propagated. For more information on the output and trajectory options, see section Regions.
6.16.1.4.7. Rigid MM Water¶
As TIP3P water might have to be treated as rigid bodies due to its parametrization - please check out the specifics of the applied force field parametrization - we offer a keyword for the automated rigid treatment of all active MM water molecules. The following keyword applies bond and angle constraints to active MM water molecules in optimizations as well as MD runs:
Rigid_MM_Water false # Default: false.
6.16.1.5. Forcefield Input¶
For the MM part of the QM/MM calculation forcefield parameters are necessary. Internally, ORCA uses the ORCA forcefield. For a description on the format, how to obtain and manipulate the forcefield parameters, see section ORCA Forcefield File.
Note
ORCAFF.prms files only need to be provided for QM/MM, QM/QM2/MM and IONIC-CRYSTAL-QMMM calculations.
For QM/QM2 and MOL-CRYSTAL-QMMM calculations there is no need to provide a ORCAFF.prms file.
The ORCAFF.prms file for the IONIC-CRYSTAL-QMMM calculation can be conveniently set up with the orca_crystalprep tool, see section orca_crystalprep.
For IONIC-CRYSTAL-QMMM and MOL-CRYSTAL-QMMM calculations the self-consistently optimized MM point charges of the entire supercell are stored in an ORCAFF.prms file, see section Charge Convergence between QM and MM region. This ORCAFF.prms file can then be used in subsequent calculations with larger QM regions, different methods and basis sets, excited state calculations, etc.
The force field filename is provided via the keyword ORCAFFFilename:
%qmmm
ORCAFFFilename "UBQ.ORCAFF.prms"
end
6.16.1.6. QM-MM, QM-QM2 and QM2-MM Boundary¶
This section is important for the QM/MM, QM/QM2 and QM/QM2/MM methods. For the latter method two boundary regions are present (between QM and QM2 as well as between QM2 and MM region), and both can go through covalent bonds. In the following we will only discuss the concept for the boundary between QM and MM, but the same holds true for the other two boundaries.
6.16.1.6.1. Link Atoms¶
ORCA automatically generates link atoms based on the information of the QM region and on the topology of the system (based on the ORCAFF.prms file). ORCA places link atoms on the bond between QM and MM atoms.
Important
When defining the QM, QM2 and MM regions, make sure that you only cut through single bonds (not aromatic, double, triple bonds, etc.).
6.16.1.6.2. Bond Length Scaling factor¶
The distance between QM atom and link atom is determined via a scaling factor (in relation to the QM-MM bond length) that is computed as the ratio of the equilibrium bond length between QM and hydrogen atom (d0_QM-H) and the equilibrium bond length between QM and MM atom (d0_QM-MM).
6.16.1.6.3. Standard QM-H Bond Length¶
For the equilibrium bond lengths to hydrogen, ORCA uses tabulated standard values for the most common atoms involved in boundary regions (C, N, O), which can be modified via keywords as defined further below. ORCA stores these values on-the-fly in a simple file (basename.H_DIST.prms). If necessary, the user can modify these values atom-specific or add others to the file and provide this file as input to ORCA (see also paragraph Get standard hydrogen bond lengths: getHDist). For QM/QM2 and QM/QM2/MM methods the equilibrium bond lengths to hydrogen are explicitly calculated.
%qmmm
# standard equilibrium bond lengths with hydrogen can be modified
Dist_C_HLA 1.09 # d0_C-H
Dist_O_HLA 0.98 # d0_O-H
Dist_N_HLA 0.99 # d0_N-H
# file can be provided which provides the used d0_X-H values specific to all atoms
H_Dist_FileName "QM_H_dist.txt"
end
6.16.1.6.4. Bonded Interactions at the QM-MM Boundary¶
The MM bonded interactions within the QM region are neglected in the additive coupling scheme. However, at the boundary, it is advisable to use some specific bonded interactions which include QM atoms. By default ORCA neglects only those bonded interactions at the boundary, where only one MM atom is involved, i.e. all bonds of the type QM1-MM1, bends of the type QM2-QM1-MM1, and torsions of the type QM3-QM2-QM1-MM1. Other QM-MM mixed bonded interactions (with more than two MM atoms involved) are included. The user is allowed to include the described interactions, which is controlled via the following keywords:
%qmmm
DeleteLADoubleCounting true # Neglect bends (QM2-QM1-MM1) and torsions
# (QM3-QM2-QM1-MM1). Default true.
DeleteLABondDoubleCounting true # Neglect bonds (QM1-MM1)
end
6.16.1.6.5. Charge Alteration¶
If QM and MM atoms are connected via a bond (defined in the topology file), the charge of the close-by MM atom (and its neighboring atoms) has to be modified in order to prevent overpolarization of the electron density of LA and QM region. This charge alteration is automatically taken care of by ORCA. ORCA provides the most popular schemes that can be used to prevent overpolarization:
- CS
Charge Shift - Shift the charge of the MM atom away to the MM2 atoms so that the overall charge is conserved
- RCD
Redistributed Charge and Dipole - Shift the charge of the MM atom so that the overall charge and dipole is conserved
- Z0
Keep charges as they are. This MM atom will probably lead to overpolarization
- Z1
Delete the charge on the MM1 atom (no overall charge conservation)
- Z2
Delete the charges on the MM1 atom and on its first (MM2) neighbors (no overall charge conservation)
- Z3
Delete the charges on the MM1 atom and on its first (MM2) and second (MM3) neighbors (no overall charge conservation)
6.16.1.7. Embedding Types¶
The following embeding schemes are available:
- Electrostatic
The electrostatic interaction between QM and MM system is computed at the QM level. Thus, the charge distribution of the MM atoms can polarize the electron density of the QM region. The LJ interaction between QM and MM system is computed at the MM level.
%qmmm
Embedding Electrostatic # Electrostatic (Default)
# Mechanical
end
In the scheme of electrostatic embedding, the evaluation of the electrostatic potential generated by the MM part can be accelerated by using the FMM algorithm (described in reference [317]). This will speed up the building of the Fock Matrix. The default recommended setup can be called using the FMM-QMMM keyword directly in the keywords line. However, more details about the algorithm parameters and all input options can be found in (see also paragraph Fast Multipole Method)
! FMM-QMMM
It is recommended to use that option whenever the MM part is composed of more than 10,000 atoms (or point charges for ECM methods).
6.16.2. Additive QMMM¶
The minimum input necessary for an additive QM/MM calculation looks as follows.
!QMMM
%qmmm
QMAtoms {0 1 2 27 28} end
ORCAFFFilename "UBQ.ORCAFF.prms"
end
6.16.3. ONIOM Methods¶
For the simulation of large systems with up to 10000 atoms, or for large QM regions in biomolecules, ORCA provides the QM/QM2 and QM/QM2/MM methods. The specifics of the two different methods are discussed further below. Here we are presenting the common concepts and keywords of both methods. For subtractive methods, we use a high level (QM) and a low level (QM2) of accuracy for different parts of the system. The advantages of this - in contrary to QM-MM methods - are as follows:
QM2 methods are polarizable, the interaction with the high level region is more accurate.
No MM parameters are needed for the atoms that are described at the QM2 level.
6.16.3.1. Available QM2 Methods¶
ORCA has several built in QM2 methods:
Semiempirical methods (AM1, PM3)
Tight-binding DFT (XTB0, XTB1, XTB (or XTB2))
Composite methods (HF-3c, PBEh-3c, r2SCAN-3c)
User-defined QM2 methods
The individual keywords for these methods are:
!QM/XTB or !QM/XTB/MM
!QM/XTB1 or !QM/XTB1/MM
!QM/XTB0 or !QM/XTB0/MM
!QM/HF-3C or !QM/HF-3C/MM
!QM/PBEH-3C or !QM/PBEH-3C/MM
!QM/R2SCAN-3C or !QM/R2SCAN-3C/MM
!QM/PM3 or !QM/PM3/MM
!QM/AM1 or !QM/AM1/MM
Users can define their own low-level methods in the following way
!QM/QM2 or !QM/QM2/MM
%qmmm
QM2CUSTOMMETHOD "B3LYP"
QM2CUSTOMBASIS "def2-SVP def2/J"
end
Alternatively, a custom QM2 method / basis set file can be provided:
!QM/QM2 or !QM/QM2/MM
%qmmm
QM2CustomFile "myQM2Method.txt" # File should be available in working directory.
end
The custom QM2 method file can contain any desired input, as e.g. the file myQM2Method.txt:
!cc-pVDZ HF TightSCF NOSOSCF KDIIS
%basis
NewAuxJKGTO Mg "AutoAux" end
end
Note
Only add methods (including convergence settings) and basis sets for the QM2Custom options. Everything else (parallelization, memory, solvation, etc.) is taken care of by ORCA itself.
6.16.3.2. Electrostatic Interaction between high and low level¶
By default we are using electrostatic embedding, i.e. the high level system sees the atomic point charges of the low level (QM2) system. These point charges are derived from the full system low level (QM2) calculation. The following methods for determining these charges are available:
Charge_Method Hirshfeld # Hirshfeld (default)
# MBIS
# CHELPG
# Mulliken
# Loewdin, default for QM2 = AM1 or PM3
The QM2 point charges can be scaled with the following keyword.
%qmmm
Scale_QM2Charges_MMAtom 1. # default is 1.
end
6.16.3.3. Boundary Region¶
The boundary between high and low level part of the system can contain covalent bonds. For the detection and realistic treatment of these covalent bonds, a topology of the large QM2 system is generated using the following keyword.
AutoFF_QM2_Method XTB # XTB (default)
# XTB1
# XTB0
# GFNFF
# HF3C
# PBEH3C
# R2SCAN3C
# PM3
# AM1
Note
By default ORCA uses the XTB method for the preparation of the QM2 topology. In order to use the default you need to make sure to have the otool_xtb binary in your ORCA PATH, see Semi-empirical tight-binding methods: Grimme’s GFN0-xTB, GFN-xTB and GFN2-xTB.
6.16.3.4. Subtractive QM/QM2 Method¶
The QM/QM2 method is a very convenient way of running multiscale calculations without the need to prepare any parameters. This method is a subtractive QM-QM method, in which we treat a part of interest on a higher level of accuracy, and the remainder of the system on lower level of accuracy. The implementation follows similar works as e.g. described in reference [571].
The method can be used in a similar way as a regular QM calculation. Let us have a look at the proton transfer in propionic acid, which can be modeled as follows:
!QM/XTB BP86 def2-TZVP def2/J
!Fast-NEB-TS NumFreq
!pal8
%qmmm
QMAtoms {0:3} end
end
%neb
product "propionicAcid_prod.xyz"
preopt true
end
*xyz 0 1
H -0.738352472 0.000000000 -5.836214279
O -0.738352472 -0.587240971 -5.061536853
O -0.738352472 1.434717404 -4.069730302
C -0.738352472 0.227304724 -3.975502162
C -0.738352472 -0.566448428 -2.687358498
H 0.133951528 -1.231202352 -2.710760176
H -1.610656472 -1.231202352 -2.710760176
C -0.738352472 0.318369069 -1.443687014
H -0.738352472 -0.294739868 -0.538164669
H 0.142397528 0.965221387 -1.423275731
H -1.619102472 0.965221387 -1.423275731
*
with the product structure file (propionicAcid_prod.xyz):
11
C3H6O2
H -0.738352472 1.628728096 -5.020130139
O -0.738352472 -0.587240971 -5.061536853
O -0.738352472 1.434717404 -4.069730302
C -0.738352472 0.227304724 -3.975502162
C -0.738352472 -0.566448428 -2.687358498
H 0.133951528 -1.231202352 -2.710760176
H -1.610656472 -1.231202352 -2.710760176
C -0.738352472 0.318369069 -1.443687014
H -0.738352472 -0.294739868 -0.538164669
H 0.142397528 0.965221387 -1.423275731
H -1.619102472 0.965221387 -1.423275731
As can be seen from the input, the only difference to a regular calculation is the necessity to define the high level region via the QMAtoms keyword.
6.16.3.4.1. System charges and multiplicities¶
The two subsystems can have different (integer) charges and multiplicities. Defining the correct charges and multiplicities is important. The charge and multiplicity defined via the coordinate section defines the charge and multiplicity of the high level region (QMAtoms). The user still needs to define the charge and multiplicity of the total system (corresponding to the sum of the charge of the high level and low level parts, and corresponding to the overall multiplicity).
%qmmm
QMAtoms {0:3} end # high level region
Charge_Total 0 # charge of the full system. Default 0.
Mult_Total 1 # multiplicity of the full system. Default 1.
end
*xyz 0 1 # charge and mult. of the high level region, i.e. atoms 0 to 3
6.16.3.4.2. Available low level methods¶
The following QM2 (low level) methods are available:
!QM/XTB
!QM/XTB1
!QM/XTB0
!QM/HF-3C
!QM/PBEH-3C
!QM/R2SCAN-3C
!QM/PM3
!QM/AM1
!QM/QM2
For information on how to specify the custom QM/QM2 method please see Available QM2 Methods.
6.16.3.4.3. Solvation¶
Implicit Solvation effects can be included in QM/QM2 calculations. On
the one hand, for QM/XTB calculations, one can adopt the
analytical linearized Poisson-Boltzmann (ALPB) solvation model, the
domain decomposition COSMO (ddCOSMO), or the extended conductor-like
polarizable continuum model (CPCM-X), and on the other hand, if no XTB
is requested, ORCA uses the C-PCM. The user just needs to add the
following tags in the ORCA input file,
XTB for the QM2 region:
!QM/XTB ALPB(Water)
or
!QM/XTB DDCOSMO(Water)
or
!QM/XTB CPCMX(Water)
No XTB for the QM2 region:
!QM/HF-3c CPCM(Water)
If the ddCOSMO (XTB) or the C-PCM (non-XTB) are requested, there are two possible ONIOM/implicit-solvation methods:[874]
C-PCM/B: The effect of the solvent is, in the first place, included in the calculation for the large QM2 system. Once this calculation finishes, the solvation charges located on the surface of the cavity for the large system are used as point charges for the subsequent low-level and high-level calculations for the small system.
C-PCM/C: The effect of the solvent is only included in the calculation for the large QM2 system.
The user can choose one scheme or the other via the tag “solv_scheme”
in the “qmmm” block:
%qmmm
solv_scheme CPCM_B # CPCM_B (default)
# CPCM_C
end
If the ALPB model or the CPCM-X are requested (within QM/XTB methods), the solvation effect is just included in the calculation for the large QM2 system (as one does for the C-PCM/C scheme).
6.16.3.5. QM/QM2/MM Method¶
The QM/QM2/MM method uses a combination of the subtractive scheme for the QM-QM2 part, and the additive scheme for the (QM-QM2) - (MM) interaction. It can be used if very large QM regions are required for biomolecules and explicitly solvated systems. The system is divided into a high level (QM), low level (QM2), and MM region (MM).
6.16.3.5.1. QM2 Atoms¶
QM2 atoms need to be defined for QM/QM2/MM calculations. They can be defined either directly
%qmmm
QMAtoms {0:4} end # list of QM atoms (counting starts from 0) or
QM2Atoms {5:22} end # list of QM2 atoms
end # an atom should not occur in both lists
*pdbfile 0 1 ubq.pdb
or via the occupancy column of a pdb file (see QM Atoms).
6.16.3.5.2. System charges and multiplicities¶
The high and low level subsystems can have different (integer) charges and multiplicities. Defining the correct charges and multiplicities is important. The charge and multiplicity defined via the coordinate section defines the charge and multiplicity of the high level region (QMAtoms). The user still needs to define the charge and multiplicity of the medium system (corresponding to the sum of the charge of the high level and low level regions, and corresponding to the overall multiplicity of the combined high and low level region). The charge of the MM region is determined based on the MM parameters provided by the forcefield.
%qmmm
QMAtoms {0:3} end # high level region
Charge_Medium 0 # charge of the medium system. Default 0.
Mult_Medium 1 # multiplicity of the medium system. Default 0.
end
*xyz 0 1 # charge and mult. of the high level region, i.e. atoms 0 to 3
6.16.3.5.3. Available low level methods¶
The following QM2 (low level) methods are available:
!QM/XTB/MM
!QM/XTB1/MM
!QM/XTB0/MM
!QM/HF-3C/MM
!QM/PBEH-3C/MM
!QM/R2SCAN-3C/MM
!QM/PM3/MM
!QM/AM1/MM
!QM/QM2/MM
For information on how to specify the custom QM/QM2/MM method please see Available QM2 Methods.
6.16.3.5.4. Example Input¶
An example for a QM/QM2/MM calculation looks as follows:
!QM/HF-3c/MM Opt B3LYP def2-TZVP def2/J NumFreq
%qmmm
ORCAFFFilename "peptideChain.ORCAFF.prms"
QMAtoms {16:33 68:82} end
QM2Atoms {0:12 83:104} end
ActiveAtoms { 0:38 65:120} end
Charge_Medium 0
end
*pdbfile -1 1 peptideChain.pdb
6.16.4. CRYSTAL-QMMM¶
For the simulation of extended systems ORCA provides the CRYSTAL-QMMM features:
- MOL-CRYSTAL-QMMM
for molecular crystals.
- IONIC-CRYSTAL-QMMM
for semiconductors and insulators.
In this section we present the concepts and keywords that are common to both methods. ORCA is a molecular code and cannot carry out periodic calculations, but ORCA has been used for modeling selected properties for ionic solids using embedded cluster models already in the past (see e.g reference [213]). With ORCA 5 we provide the CRYSTAL-QMMM method, which simplifies setting up and running these types of embedded cluster calculations. The user needs to provide the structure with a (large enough) supercell representative for the crystal, and provide a list of QM atoms. The QM region should be located in the central part of the supercell. The QM atoms are embedded in the remainder of the supercell, the surrounding MM atoms, which are represented by atom-centered point charges and influence the QM calculations. How these charges are obtained, is outlined in the next paragraph.
6.16.4.1. Charge Convergence between QM and MM region¶
The core idea of the CRYSTAL-QMMM method is to reach charge convergence between the QM and the MM atoms. For this purpose, the MM charges are updated self-consistently with the atomic (CHELPG) charges of the QM atoms. This idea follows the work of reference [213] for the IONIC-CRYSTAL-QMMM and of reference [109] for the MOL-CRYSTAL-QMMM.
!IONIC-CRYSTAL-QMMM or MOL-CRYSTAL-QMMM
%qmmm
Conv_Charges true # default true for both CRYSTAL-QMMM variants
Conv_Charges_MaxNCycles 10 # default 30 for MOL-CRYSTAL-QMMM
# default 10 for IONIC-CRYSTAL-QMMM
Conv_Charges_ConvThresh 0.01 # threshold (default 0.01) for maximum charge change
#for atom type between two subsequent charge convergence
# cycles
Scale_FormalCharge_MMAtom 1. # default is 1. MM atomic charges used in QM part of
# the CRYSTAL-QMMM calculation are scaled by this factor
end
During optimizations, the charge convergence scheme can be used (i) only at the very beginning, (ii) every N geometry steps, or (iii) for each geometry step, using the following keyword:
%geom
ReConvCharge 1 # default is 1. Redo charge convergence scheme every N steps.
end
After charge convergence, the converged charges are stored in the file basename.convCharges.ORCAFF.prms. It can be used for a later calculation (with the same or different electronic structure settings) with the following keyword combination:
!IONIC-CRYSTAL-QMMM or !MOL-CRYSTAL-QMMM
%qmmm
Conv_Charges false
ORCAFFFilename "firstjob.convCharges.ORCAFF.prms"
end
6.16.4.2. MM-MM Interaction¶
Unlike with the other multiscale methods (QMMM, QM/QM2, QM/QM2/MM) the MM region is only treated as a perturbation to the QM region. The (bonded and nonbonded) MM-MM interaction is not computed - since it would not affect the QM-MM interaction - and thus the active region (optimizations, frequencies, …) in CRYSTAL-QMMM calculations should be restricted to atoms of the QM region.
6.16.4.3. Preparing CRYSTAL-QMMM Calculations¶
The input structures and input files for the CRYSTAL-QMMM calculations can be prepared with the orca_crystalprep module (see section orca_crystalprep).
6.16.4.4. MOL-CRYSTAL-QMMM¶
Conceptually we use an additive QM/MM calculation, in which the QM region interacts with the MM region only via nonbonded interactions. Lennard-Jones interaction between QM and MM region is computed using van der Waals parameters from the UFF force field. The charge convergence treatment between QM and MM region starts with zero charges on the MM atoms, and is iterated until convergence of the QM atomic (CHELPG) and MM point charges is achieved.
The MOL-CRYSTAL-QMMM method expects as structural input a supercell that consists of repeating, identical molecular subunits, i.e. the atom ordering of the molecules should always be the same, and one molecular subunit should follow the next one. The minimum input necessary for a MOL-CRYSTAL-QMMM run looks as follows.
%qmmm
NUnitCellAtoms 16 # provide the number of atoms per molecular subunit
QMAtoms {160:175} end # QM atoms, should be in central part of the supercell.
# Must be at least one complete molecular subunit.
end
Note
For molecular crystals it is assumed that the entire supercell consists of repeating units which each have zero charge. QM regions should consist of one or multiple of such charge-neutral units.
6.16.4.4.1. Extending the QM Region¶
ORCA can extend the QM region (which we call QM core region, defined by
the original QMAtoms input) by one or multiple layers of molecular
subunits using the HFLayers keyword. If the HFLayers keyword is
used, the molecular subunits of the defined number of layers around the
QM core region is added to the QM region. The layers of molecular
subunits around the QM core region are detected via distance criteria.
%qmmm
HFLayers 0 # default 0
HFLayerGTO "LANL2DZ" # default. Use this basis set for the HFLayer atoms
HFLayerECP "HayWadt" # default. Use these ECPs for the HFLayer atoms.
Conv_Charge_UseQMCoreOnly true # only use the QM core region for the charge
# convergence scheme
end
The HFLayer can be seen as a buffer region between the molecular subunit of interest (original QM atoms) and the surrounding subunits. The different possible reasons for HFLayers are as follows:
For DLPNO calculations with HFLayers, the DLPNO multilevel feature is automatically switched on. The molecules of the HFlayer form an own fragment which is treated on HF level only, i.e. these molecules are not included in the Post-HF treatment.
The HFLayers can also be used for non-DLPNO calculations. In such cases, the HFLayer molecules are in the same way included in the QM calculation as the other QM molecules. But their definition is a convenient way to choose a different basis set (and ECPs) for those molecules.
It can be assumed that the QM charges computed for the QM core region are more realistic than the charges of the HFLayer atoms near the MM atomic charges. Thus, using only the QM atomic charges of the QM core region represent more realistic charges for the charge convergence scheme.
Note
The term HFLayers might be misleading. Strictly speaking, only for MOL-CRYSTAL-QMMM calculations with DLPNO methods (e.g. DLPNO-CCSD(T), DLPNO-MP2, DLPNO-B2PLYP) the HF layer atoms are treated on HF level. For other methods (e.g. DFT) the HF layer atoms are treated with the same electronic structure method as the QM core region atoms.
If necessary, the atoms of the HFLayer can be defined by the user. Make sure that complete molecular subunits are used here.
%qmmm
HFLayerAtoms {0:15} end
end
6.16.4.4.2. Example Input¶
An example for a MOL-CRYSTAL-QMMM calculation looks as follows:
! MOL-CRYSTAL-QMMM
! PBE def2-SVP Opt NumFreq
%qmmm
NUnitCellAtoms 13
QMAtoms {221:233} end
ActiveAtoms {221:233} end
end
*xyzfile 0 1 Ala_SC.xyz
6.16.4.5. IONIC-CRYSTAL-QMMM¶
Conceptually we use an embedded cluster calculation, in which the QM region interacts only with the MM atomic point charges of the surrounding. Unlike with MOL-CRYSTAL-QMMM, the Lennard-Jones interaction between QM and MM region is not computed. The charge convergence treatment between QM and MM region starts with the initial MM charges (the charges initially read from the ORCAFF.prms file) and is iterated until convergence of the QM atomic (CHELPG) and MM point charges is achieved.
6.16.4.5.1. Boundary Region¶
For ionic crystals a boundary region needs to be introduced between the QM region and the atomic point charges of the MM atoms in order to avoid spurious electron leakage from the clusters and to treat dangling bonds on the surface of the QM region. This boundary region consists of capped effective core potentials (cECPs). These cECPs are placed on the position of the MM atoms that are directly adjacent to the QM region. ORCA analyzes the connectivity of the atoms of the supercell and can thus detect the layers of atoms around the QM region, where the first layer consists of the atoms directly bonded to the QM atoms, the second layer consists of the atoms directly bonded to the atoms of the first layer and so on. The charges of these cECPs are determined in the same way as the MM atomic charges during the charge convergence scheme.
%qmmm
ECPLayerECP "SDD" # cECPs used for the boundary region
ECPLayers 3 # Number of cECP layers around the QM region.
# Default is 3.
Scale_FormalCharge_ECPAtom 1. # default is 1. Charges on cECPs are scaled by
# this factor
end
6.16.4.5.2. Extending the QM Region¶
ORCA can extend the QM region (which we call QM core region, defined by
the original QMAtoms input) by one or multiple layers of atoms using the
HFLayers keyword. If the HFLayers keyword is used, the atoms of the
defined number of layers around the QM core region is added to the QM
region. The layers of atoms around the QM core region are detected in
the same way as defined above for the ECPLayers.
%qmmm
HFLayers 0 # default 0
HFLayerGTO "LANL2DZ" # default. Use this basis set for the HFLayer atoms
HFLayerECP "HayWadt" # default. Use these ECPs for the HFLayer atoms.
Conv_Charge_UseQMCoreOnly true # only use the QM core region for the charge
# convergence scheme
end
The HFLayer can be seen as a buffer region between the actual atoms of interest (original QM atoms) and the surrounding cECP and MM point charge environment.The different possible reasons for HFLayers are as follows:
For DLPNO calculations with HFLayers, the DLPNO multilevel feature is automatically switched on. The atoms of the HFLayer form an own fragment which is treated at HF level only, i.e. these atoms are not included in the Post-HF treatment.
It can be assumed that the QM charges computed for the QM core region are more realistic than the charges of the HFLayer atoms near the cECPs and MM atomic charges, in particular for highly charged systems. Thus, using only the QM atomic charges of the QM core region represent more realistic charges for the charge convergence scheme.
Note
The term HFLayers might be misleading. Strictly speaking, only for IONIC-CRYSTAL-QMMM calculations with DLPNO methods (e.g. DLPNO-CCSD(T), DLPNO-MP2, DLPNO-B2PLYP) the HF layer atoms are treated on HF level. For other methods (e.g. DFT) the HF layer atoms are treated with the same electronic structure method as the QM core region atoms.
If necessary, the atoms of the HFLayer can be defined by the user:
%qmmm
HFLayerAtoms {29:35} end
end
The charge of the HFLayer is automatically computed based on the formal charges of its atoms. If necessary, the HFLayer charge can be provided by the user.
%qmmm
Charge_HFLayer 10 # by default it is computed automatically based on the formal
# charges from the ORCAFF.prms file
end
6.16.4.5.3. Net Charge of the Entire Supercell¶
For ionic crystals, the QM region (as well as the entire supercell) can consist of an unequal number of atoms of each atom type. As a result of that, the QM region may have non-zero net charge. Consequently, when assigning the atomic point charges to the MM atoms during the charge convergence scheme, the sum of the charge of the QM region, of the ECP layer, and of the atomic charges of the MM atoms, can deviate from the ideal charge of the supercell. In order to enforce the ideal charge of the supercell, the total charge of the entire system is enforced by equally shifting the charges of all MM (including boundary) atoms.
%qmmm
Charge_Total 0 # default 0. Total charge of the supercell
EnforceTotalCharge true # enforce the total charge by shifting MM charges
end
The charge shifting procedure can be restricted to the MM atoms on the
outer parts of the supercell by defining the number of OuterPCLayers. If
OuterPCLayers is set to 1, only the atoms of the most outer layer of
the supercell (recognized based on the connectivity of the atoms) are
included in the charge shifting procedure. If OuterPCLayers is larger
than 2, the atoms connected to the most outer layer are additionally
included in the charge shifting procedure, etc.
%qmmm
OuterPCLayers 0 # default 0, i.e. charge shifting for all MM atoms
end
Note
The charge of the QM core region can be chosen to be assigned automatically. This overwrites the charge provided in the xyz or pdb section.
%qmmm
AutoDetectQCCharge false # default is false
end
6.16.4.5.4. Example Input¶
A minimal example for an IONIC-CRYSTAL-QMMM calculation looks as follows:
! IONIC-CRYSTAL-QMMM
! NMR
! PBE pcSseg-2 def2/JK def2-svp/C
%qmmm
QMAtoms {0:6} end
ORCAFFFilename "nacl6.ORCAFF.prms"
CHARGE_TOTAL 0
end
*xyzfile -5 1 nacl6.xyz
6.16.4.5.5. Different QM and MM regions Stored in the PDB file¶
The stored PDB file contains the following entries in its occupancy column.
- 0
MM atoms mimicked by surrounding point charges.
- 1
QM core region atoms
- 2
HFLayer atoms
- 3
cECPs
- 4
OuterPCLayer atoms (subset of MM atoms) if defined, are the only atoms that are used in the charge shift procedure for enforcing the total charge.
6.16.5. Additional Keywords¶
Here we list additional keywords and options that are accessible from
within the %qmmm block and that are relevant to QM/MM calculations.
Some of these keywords can also be accessed via the %mm block, see
section Available Keywords for the MM module.
!QMMM
%qmmm
# Charge alteration scheme preventing overpolarization.
ChargeAlteration CS # CS (Default)
# RCD
# Z0
# Z1
# Z2
# Z3
# Parameters for placing the shifted and redistributed charges for RCD and CS schemes.
Scale_RCD 0.5 # Defines where on the bond between MM1 and MM2 the
# shifted charge is positioned. Default: midpoint.
Scale_CS 0.06 # Defines where on the bond between MM2 and MM1/MM3 the
# shifted charge is positioned. Default: 0.06 x bond.
# The extended active region, activeRegionExt, contains the atoms surrounding the
# active atoms.
ExtendActiveRegion distance # rule to choose the atoms belonging to activeRegionExt.
# no - do not use activeRegionExt atoms
# cov_bonds - add only atoms bonded covalently to
# active atoms
# distance (default) - use a distance criterion (VDW
# distance plus Dist_AtomsAroundOpt)
Dist_AtomsAroundOpt 0.0 # in Angstrom (Default 0). Meaning only freeze atoms
# which touch the active atoms by their VDW spheres.
# only needed for ExtendActiveRegion distance
OptRegion_FixedAtoms {27 1288:1290 4400} end # manually define the activeRegionExt
# atoms. Default: empty list.
# Printing options. All are true by default.
PrintLevel 1 # Can be 0 to 4, 0=nothing, 1=normal, ...
PrintOptRegion true # Additionally print xyz and trj for opt region
PrintOptRegionExt false # Additionally print xyz and trj for extended opt region
PrintQMRegion true # Additionally print xyz and trj for QM region
PrintFullTRJ true # Print trajectory of full system. Default true.
PrintPDB true # Additionally print pdb file for entire system, is
# updated every iteration for optimization
# Computing MM nonbonded interactions within non-active region.
Do_NB_For_Fixed_Fixed true # Compute MM-MM nonbonded interaction also for
# non-active atom pairs. Default true.
# Treats all active water molecules that are treated on MM level as rigid bodies
# in optimization and MD simulation, see section "Rigid Water".
Rigid_MM_Water false # Default false.
end
*pdbfile 0 1 ubq.pdb # structure input via pdb file, but also possible via xyz file
6.17. QM/MM via Interfaces to ORCA¶
ORCA is easy to interface as a QM engine in pretty much any QM/MM environment, as it will accept a set of point charges from an external file (see section Inclusion of Point Charges) and it will return, in a transparent format, all the required information for computing energies and gradients to the driving program. In our research group we have experience with four different QM/MM environments: ChemShell, Gromacs, pDynamo and NAMD. In the following each of the interfaces are described. Is beyond the scope of the manual to be exhaustive, and only minimal working examples are going to be presented. These will cover mainly the technical aspects with respect to the QM part of the QM/MM calculation. For the initial preparation of the system the user is referred to the documentation of the driving program.
6.17.1. ORCA and Gromacs¶
In cooperation with the developers of Gromacs software package we developed an interface to ORCA. The interface is available starting with Gromacs version 4.0.4 up to version 4.5.5.
As mentioned above, the initial setup has to be done with the Gromacs.
In the QM/MM calculation Gromacs will call ORCA to get the energy and
the gradients of the QM atoms. The interface can be used to perform all
job types allowed by the Gromacs software package, e.g. optimizations,
molecular dynamics, free energy. In addition, for geometry optimizations
we have implemented a microiterative scheme that can be requested by
setting the keyword bOpt = yes in the Gromacs .mdp file. This will
cause ORCA to perform a QM geometry optimization every time it is called
by Gromacs. During this optimization ORCA will also compute the
Lennard-Jones interaction between the QM and MM atoms, and freeze the
boundary atoms. This microiterative scheme can also be used to perform a
transition state optimization. If you are looking for a transition state
and have a good initial guess for the structure, you can carry out an
optimization of the MM system and at the same time perform a transition
state optimization of the QM system with ORCA via the microiterative
scheme. Since it is expensive to calculate the Hessian for each
microiterative microiterative step, the user can tell ORCA to use the
(updated) Hessian matrix of the previous step via read_temp_Hess in
the ORCA input.
In order to allow full flexibility to the user, the information for the
QM run are provided in a separate plain text file, called
GromacsBasename.ORCAINFO. When Gromacs writes the input for the
ORCAcalculation, it will merge this file with the information on the
coordinates, point charges, Lennard-Jones coefficients and type of
calculation (EnGrad, Opt, TSOpt).
Below is a typical example of an input file created by Gromacs, where for each Gromacs optimization step, a full optimization of the QM-part will be performed by ORCA, instead of just doing the energy and gradient calculation.
# Optimization step in the Lennard-Jones and point charges field of the MM atoms
! QMMMOpt
# file containing the Lennard Jones coefficients for the Lennard-Jones interaction
%LJCoefficients "temp.LJ"
# file containing the point charges for the electrostatic interaction
%pointcharges "temp.pc"
%geom
# calculate the exact Hessian before the first optimization step
Calc_Hess true
# in case of a TS optimization the updated Hessian of the previous
# TS optimization run is read instead of calculating a new one
read_temp_Hess true
end
Note
If you are using
bOptorbTSyou have to take care of additional terms over the boundary. Either you can zero out the dihedral terms of the Q2-Q1-M1-M2 configurations, or you can fix the respective Q2 atoms.If you want to use the ORCA constraints, you have to also put them in the Gromacs part of the calculation.
If there are no bonds between the QM and MM systems, the
bOptoptimization of the QM system might have convergence problems. This is the case if the step size is too large and thus QM atoms come too close to MM atoms. The convergence problems can be circumvented by adding extra coordinates and Hessian diagonal values for Cartesian coordinates of selected QM atoms to the set of internal coordinates. This should be done for only few atoms that are in the boundary region.
%geom
# add the Cartesian position of atoms 2 and 4 to the
# set of internal coordinates with a diagonal Hessian value of 0.1
Hess_Internal
{C 2 D 0.1}
{C 4 D 0.1}
end
end
6.17.2. ORCA and pDynamo¶
The interface with the pDynamo library is briefly documented here. It uses the same plain text files to exchange information between the pDynamo library and ORCA. In order to have simiar functionality as presented above, we have extended the interface to generate more complex ORCA input files by accepting verbatim blocks of text. We have also implemented in pDynamo the capability of writing files containing Lennard-Jones parameters.
A simple example which calculates the QM/MM energy for a small designed peptide in which the side chain of one amino acid is treated QM is ilustrated below. For this example, a complete geometry optimization is going to be performed during the ORCA call.
import os
import pCore
import pBabel
import pMolecule
import pMoleculeScripts
# Read the initial coordinates from the .pdb file.
system = pMoleculeScripts.PDBFile_ToSystem(
"1UAO.pdb", modelNumber=1, useComponentLibrary=True)
# Instantiate the required models.
mmmodel = pMolecule.MMModelOPLS("protein")
nbmodel = pMolecule.NBModelORCA()
qcmodel = pMolecule.QCModelORCA(
command=os.getenv("ORCA_COMMAND"),
deleteJobFiles=False, header="bp86 def2-svp qmmmopt/pdynamo",
job="chignolin", run=True)
# Assign the models to the system.
system.DefineMMModel(mmmodel)
system.DefineQCModel(
qcmodel, qcSelection=pCore.Selection([35, 36, 37, 34, 40, 41]))
system.DefineNBModel(nbmodel)
system.electronicState = pMolecule.ElectronicState(
charge=-1, multiplicity=1)
# Print a summary and calculate the energy.
system.Summary()
system.Energy()
After the execution of the above Python program, a series of files are
going to be created chignolin.inp, chignolin.pc, chignolin.lj and
ORCA is going to be called. The generated ORCAinput file is listed
below.
! bp86 def2-svp qmmmopt/pdynamo
% geom
constraints
{C 0 C}
{C 1 C}
end
end
% pointcharges "chignolin.pc"
% ljcoefficients "chignolin.lj"
* xyz -1 1
H -1.0637532468 1.1350324675 2.4244220779
C -0.5230000000 0.6870000000 3.2490000000
C 0.4180000000 1.7240000000 3.8660000000
O -0.0690000000 2.7620000000 4.2830000000
O 1.6090000000 1.4630000000 3.9110000000
H -1.2240000000 0.3460000000 3.9970000000
H 0.0550000000 -0.1510000000 2.8890000000
*
There are few points that have to be raised here. Because the keyword
qmmm/pdynamo was specified in the header variable, the pDynamo
library will automatically add the constraint block in the ORCA input,
which will freeze the link atoms and the QM atoms to which they are
bound. It will also generate the chignolin.lj file containing all the
Lennard-Jones parameters. The important parts of this file, which is
somewhat different than the one generated by Gromacs, are listed next.
# number of atoms combination rule
138 0
# x y z sigma epsilon id
-6.778000 -1.424000 4.200000 3.250000 0.711280 -1
-6.878000 -0.708000 2.896000 3.500000 0.276144 -1
-5.557000 -0.840000 2.138000 3.750000 0.439320 -1
...
0.433000 0.826000 0.502000 2.960000 0.878640 -1
-0.523000 0.687000 3.249000 3.500000 0.276144 1
0.418000 1.724000 3.866000 3.750000 0.439320 2
-0.069000 2.762000 4.283000 2.960000 0.878640 3
1.609000 1.463000 3.911000 2.960000 0.878640 4
-2.259000 -0.588000 1.846000 0.000000 0.000000 -1
-1.795000 2.207000 2.427000 2.500000 0.125520 -1
-1.224000 0.346000 3.997000 2.500000 0.125520 5
0.055000 -0.151000 2.889000 2.500000 0.125520 6
-0.311000 2.922000 0.557000 3.250000 0.711280 -1
...
-1.387000 -2.946000 5.106000 2.500000 0.125520 -1
# number of special pairs
22
# atom1 atom2 factor
34 32 0.000000
35 39 0.500000
40 31 0.000000
41 30 0.500000
41 32 0.500000
36 31 0.500000
40 32 0.500000
40 39 0.500000
34 31 0.000000
35 30 0.500000
34 11 0.500000
34 38 0.500000
41 31 0.000000
37 31 0.500000
34 33 0.500000
34 39 0.000000
40 30 0.500000
41 39 0.500000
34 30 0.000000
35 31 0.000000
34 42 0.500000
35 32 0.500000
The second number on the first line refers to the type of combination
rule used to calculate the Lennard-Jones interaction. It is 0 if a
geometric average is used (OPLS force field), or 1 for the
Lorentz-Berthelot rules (AMBER force field). The id on the last column
is -1 for MM atoms and is equal to the atom number for the QM atoms. In
this case the hydrogen link atom is atom 0. The last block of the file
is composed of atom pairs and a special factor by which their
Lennard-Jones interaction is scaled. In general this factor is equal to
1, but for atoms one or two bonds apart is zero, while for atoms three
bonds apart depends on the type of force field, and in this case is 0.5.
After successful completion of the ORCA optimization run, the information will be relayed back the pDynamo library, which will report the total QM/MM energy of the system. At this point the type QM/MM of calculation is limited only by the capabilities of the pDynamo library, which are quite extensive.
6.17.3. ORCA and NAMD¶
Since version 2.12, NAMD is able to perform hybrid QM/MM calculations. A more detailed explanation of all available key words, setting up the calculation and information on tutorials and on the upcoming graphic interface to VMD are available on the NAMD website.
Similar to other calculations with NAMD, the QM/MM is using a pdb file
to control the active regions. An example is shown below, where the
sidechain of a histidine protonated at N\(\epsilon\) is chosen to be the
QM region. Either the occupancy column or the b-factor column of the
file are used to indicate which atom are included in a QM area and which
are treated by the forcefield. In the other column, atoms which are
connecting the QM area and the MM part are indicated similarly. To
clarify which column is used for which purpose, the keywords qmColumn
and qmBondColumn have to be defined in the NAMD input.
...
ATOM 1737 CA HSE P 117 14.762 47.946 31.597 1.00 0.00 PROT C
ATOM 1738 HA HSE P 117 14.751 47.579 32.616 0.00 0.00 PROT H
ATOM 1739 CB HSE P 117 14.129 49.300 31.501 1.00 1.00 PROT C
ATOM 1740 HB1 HSE P 117 14.407 49.738 30.518 0.00 1.00 PROT H
ATOM 1741 HB2 HSE P 117 13.024 49.194 31.509 0.00 1.00 PROT H
ATOM 1742 ND1 HSE P 117 13.899 51.381 32.779 0.00 1.00 PROT N
ATOM 1743 CG HSE P 117 14.572 50.261 32.582 0.00 1.00 PROT C
ATOM 1744 CE1 HSE P 117 14.615 52.043 33.669 0.00 1.00 PROT C
ATOM 1745 HE1 HSE P 117 14.356 53.029 34.064 0.00 1.00 PROT H
ATOM 1746 NE2 HSE P 117 15.678 51.318 33.982 0.00 1.00 PROT N
ATOM 1747 HE2 HSE P 117 16.369 51.641 34.627 0.00 1.00 PROT H
ATOM 1748 CD2 HSE P 117 15.706 50.183 33.335 0.00 1.00 PROT C
ATOM 1749 HD2 HSE P 117 16.451 49.401 33.388 0.00 1.00 PROT H
ATOM 1750 C HSE P 117 13.916 47.000 30.775 0.00 0.00 PROT C
ATOM 1751 O HSE P 117 12.965 46.452 31.334 0.00 0.00 PROT O
...
NOTES:
If one wants to include more than one QM region, integers bigger than 1 can be used to define the different regions.
Charge groups cannot be split when selecting QM and MM region. The reason is that non-integer partial charges may occur if a charge group is split. Since the QM partial charges are updated in every QM iteration, this may lead to a change in the total charge of the system over the course of the MD simulation.
The occupancy and b-factor columns are used for several declarations in NAMD. If two of these come together in one simulation, the keyword
qmParamPDBis used to define which pdb file contains the information about QM atoms and bonds.To simplify the selection of QM atoms and writing the pdb file a set of scripts is planned to be included in future releases of NAMD.
To run the calculation, the keyword qmForces on must be set. To select
ORCA qmSoftware "orca" must be specified and the path to the
executables must be given to qmExecPath, as well as a directory where
the calculation is carried out (qmBaseDir). To pass the method and
specifications from NAMD to ORCA qmConfigLine is used. These lines
will be copied to the beginning of the input file and can contain both
simple input as well as block input. To ensure the calculation of the
gradient, the engrad keyword should be used.
The geometry of the QM region including the selected links as well as
the MM point charges are copied to the ORCA inputfile automatically.
Multiplicity and charge can be defined using qmMult and qmCharge,
although the latter can be determined automatically by NAMD using the MM
parameters. It should be noted at this point that NAMD is capable to
handle more than one QM region per QM/MM calculation. Therefore for each
region, charge and multiplicity are expected. In the case of only one QM
region, the input looks like the following:
qmMult "1 1"
qmCharge "1 0"
Currently, two charge modes are available: Mulliken and CHELPG. They
have to be specified in the NAMD input using QMChargeMode and in the
qmConfigLine, respectively. Different embedding schemes, point charge
schemes and switching functions are available, which will be not further
discussed here. Another useful tool worth mentioning is the possibility
to call secondary executables before the first or after each QM software
execution using QMPrepProc or QMSecProc, respectively. Both are
called with the complete path and name to the QM input file, allowing
e.g. storage of values during an QM/MM-MD.
It is strongly enphasized that at this points both programs are constantly developed further. For the latest information, either the ORCA forum or the NAMD website should be consulted.