QM/MM Theory
QM/MM is a type of multilevel method (see Hybrid Theory) where a QM theorylevel and an MMtheory level are used for different parts of the system to give a combined description of the system. The purpose is usually to approximate the QMtheory description of the whole system (usually out of reach) by replacing most of the system by a classical MMtheory while retaining a quantum description of the important part of the system.
The method can be described by the equations:
where \(E_{QM}\) is the QMenergy of the QMregion, \(E_{MM}\) is the MMenergy of the MMregion and \(E_{coupling}\) is the interaction between the QM and MM region. The coupling terms account for electrostatic, vdW and covalent (bonded) interactions between the QM and MM regions.
QM/MM in ASH is highly flexible as one can combine any QMtheory in ASH with an MMTheory object of class NonbondedTheory (see MM Interfaces) or OpenMMTheory (see OpenMM interface). One simply needs to combine a defined QMtheory object (QM Interfaces) and an MMtheory object (usually OpenMMTheory, see OpenMM interface) within a QMMMTheory object. Additionally one specifies which atoms are QM and which are MM and the type of QMMM coupling (either mechanical or electrostating embedding). Note that in contrast to a QMtheory object or an MMtheory object, we also pass a fragment object to a QMMMTheory object as the atominformation is needed for defining the division of QMregion and MMregions.
QMMMTheory class
class QMMMTheory:
def __init__(self, qm_theory=None, qmatoms=None, fragment=None, mm_theory=None, charges=None,
embedding="elstat", printlevel=2, numcores=1, actatoms=None, frozenatoms=None, excludeboundaryatomlist=None,
unusualboundary=False, openmm_externalforce=False, TruncatedPC=False, TruncPCRadius=55, TruncatedPC_recalc_iter=50,
qm_charge=None, qm_mult=None, chargeboundary_method="shift",
dipole_correction=True, linkatom_method='simple', linkatom_simple_distance=None,
linkatom_forceproj_method="adv", linkatom_ratio=0.723):
QMMMTheory options:
Keyword 
Type 
Default value 
Details 


ASHTheory 
None 
Required: The theory level for the QMregion. Must be a valid ASH Theory
that supports electrostatic embedding.


list 
None 
Required: List of QMatom indices that defined the QMregion. Atoms not
in list are treated as MM atoms.


ASHTheory 
None 
Required: The theory level for the MMregion. Must be an object of class
OpenMMTheory or NonbondedTheory.


ASH Fragment 
None 
Required: ASH fragment, needed for setting up QMregion and MMregion. 

integer 
None 
Optional: Specify the charge of the QMregion. This takes precedence
over other charge specifications.


integer 
None 
Optional: Specify the spin multiplicity of the QMregion. This takes
precedence over other mult specifications.


list 
None 
Optional: list of atom charges. If not defined then charges will be read
from mm_theory.


integer 
2 
Optional: The printlevel setting. If printlevel >= 3 then more printing
and gradient files are written to disk.


integer 
1 
Optional: Number of CPU cores to use for qm_theory. If defined, takes
precedence over QMTheory setting.


list 
None 
Optional: List of atoms that are excluded from adding linkatoms to. 

Boolean 
False 
Optional: Boundaryoption: overrides ASH from quitting if an unusual QMMM boundary is found. 

Boolean 
False 
Optional: Option for passing QM/MM force as an external force to OpenMMTheory. 

Boolean 
False 
Optional: Truncated Pointcharge Option on or off. 

float 
55 
Optional: Truncated PC option; Radius (Å) for the truncated PC region. 

integer 
50 
Optional: Truncated PC option; frequency for recalculating with full PC field. 

list 
None 
Optional: List of active atoms in QM/MM. NOTE: Only compatible if
mm_theory is of NonBondedTheory class.


list 
None 
Optional: List of frozen atoms in QM/MM, alternative to actatoms.
NOTE: Only compatible if mm_theory is of NonBondedTheory class.


string 
shift 
What chargeboundary method to use for covalent QMMM boundary.
Default option: shift' . Other option: 'rcd'


Boolean 
True 
For chargeboundary='shift', whether to add additional charges to preserve dipole


string 
'simple' 
What linkatom method to use. Options: 'simple', 'ratio'


float 
None 
For linkatom_method='simple', what QM1L linkatom distance to use. Default setting is 1.09 Å.


float 
0.723 
For linkatom_method='ratio', what ratio to use. Default is 0.723.


string 
'adv' 
What linkatom force projection method to use. Options: 'adv', 'lever'

Dummy example showing how to combine a QMTheory and MMTheory object into a QMMMTheory object:
#QM theory: xTB
qm = xTBTheory(xtbmethod='GFN1')
#Creating new OpenMM object from OpenMM XML files (builtin CHARMM36 and a userdefined one)
omm = OpenMMTheory(xmlfiles=["charmm36.xml", "charmm36/water.xml", "./specialresidue.xml"], pdbfile="topology.pdb",
periodic=True, platform='CPU', numcores=numcores, autoconstraints=None, rigidwater=False)
#QM/MM theory object. QMregion defined as atom indices 500,501,502 and 503
qmmm = QMMMTheory(qm_theory=qm, mm_theory=omm, fragment=fragment, embedding="elstat",
qmatoms=[500,501,502,503], printlevel=2, qm_charge=1, qm_mult=6)
Defining the charge and spin multiplicity of the QMregion
To define the charge and spin multiplicity of the QMregion in QM/MM calculations you can choose between 3 options:
 Define qm_charge and qm_mult attributes when defining the QMMMTheory object (recommended):
qmmm = QMMMTheory(qm_theory=qm, mm_theory=omm, fragment=frag, qm_charge=1, qm_mult=6)
 Define as input to the jobfunction (e.g. Singlepoint):
Singlepoint(theory=qmmm, fragment=frag, charge=1, mult=6)
 Provide the information in the fragment definition:
frag=Fragment(xyzfile="system.xyz", charge=1, mult=6)
This information will be passed onto the QMprogram when called. The qm_charge/qm_mult option takes precedence over the other options, followed by the jobtype keyword.
Note that the specified charge and multiplicity of the QMregion needs to be consistent with what chemical groups are present in the QMregion.
Defining QMregion and active region
The QMregion needs to be defined in the QMMMTheory object by specifying which atom indices (of the full system) should be QMatoms (everything else is MM).
qmmm = QMMMTheory(qm_theory=qm, mm_theory=omm, fragment=fragment,
qmatoms=[500,501,502,503], qm_charge=1, qm_mult=6)
Similarly the activeregion (when performing a geometry optimization) needs to be defined by specifying which atoms are allowed to move. This information should be provided to the Optimizer instead (not the QMMMTheory object).
Optimizer(fragment=fragment, theory=QMMMobject, ActiveRegion=True, actatoms=[400,450,500,501,502,503,550,600,700])
Note that for MD simulations one should use the frozenatoms option instead of actatoms.
Both the QMregion and Active regions are thus defined as simple Python lists of integers (corresponding to atom indices). This approach allows you considerable flexibility in defining the QM/MM job. The QMregions and active regions can be the same or different (quite common).
Definition of the QMregion when part of a larger molecule (e.g. a protein) requires a bit of insight into the system and knowledge of how the QM/MM boundary works (see next section). It is usually best to define the QMregion by manually creating the list of atoms. One can doublecheck whether the region is correct by using the fragedit.py script (see Coordinates and fragment tools) or check the QMregion coordinates printed in the ASH output.
The active region is typically much larger than the QMregion (for a protein, an active region of approx. 1000 atoms is common) and it is usually inconvenient to define it manually. ASH provides a convenient function actregiondefine (see Coordinates and fragment tools) to define such a large list of atom indices. This function can also be used to define the QMregion.
As these lists can be large it is convenient to read them from a file. ASH provides a function read_intlist_from_file (see Coordinates and fragment tools) to read a list of integers from a file and return a Python list. The file should contain integers separated by spaces or newlines.
qmatoms = read_intlist_from_file("qmatoms")
actatoms = read_intlist_from_file("active_atoms")
qmmm = QMMMTheory(qm_theory=qm, mm_theory=omm, fragment=fragment,
qmatoms=qmatoms, qm_charge=1, qm_mult=6)
Optimizer(fragment=fragment, theory=QMMMobject, ActiveRegion=True, actatoms=actatoms)
Note
Note that if one wants to use an active region in MD simulations at the QM/MM level one would have to define frozenatoms inside the OpenMMTheory object.
QM/MM coupling: mechanical vs. electrostatic embedding
QM/MM typically comes in 2 flavours: mechanical embedding and electrostatic embedding. The approaches differ in how the the QM/MM energy expression is actually constructed:
Mechanical embedding is the simplest QM/MM coupling scheme where the \(E_{elstat}\) term is calculated at the MMlevel as a classical Coulomb term of pointcharge interactions between the QM and MM regions. Choose embedding = 'mechanical' when defining the QMMMTheory object to use mechanical embedding. Mechanical embedding requires pointcharges to be defined for each atom inside the QMregion which can introduce problem if the QMregion contains exotic entities such as metal complexes or clusters, and the QMregion atom charge definitions will requires some care. The main drawback of mechanical embedding is that the QMenergy of the QMregion (\(E_{QM}\)) is calculated entirely without any environment present. For systems with strong polarization effects between regions this can be a major drawback.
Electrostatic embedding is a more sophisticated QM/MM coupling scheme where the \(E_{elstat}\) term is calculated at the QMlevel, by calculating it at the same time as the (\(E_{QM}\)) term via the QMprogram. Choose embedding = 'elstat' when defining the QMMMTheory object to use electrostatic embedding (it is the default). By including all MM pointcharges as additional nucleilike terms in the 1electron Hamiltonian of the QMenergy expression, the QMenergy is calculated in a field of the MM pointcharges, i.e. the QM electron density is polarized by the environment. The \(E_{vdW}\) is in contrast calculated at the MMlevel as a classical LennardJones term between the QM and MM regions and is calculated at the same time as the MMenergy of the MMregion. The covalent bonded term (\(E_{covalent}\)) also gets incorporated in the MMenergy calculation (though the linkatom part is handled by the QMpart). This means that in electrostatic embedding the QM/MM energy expression is actually calculated like this:
where the \(E_{QM}^{pol} = (E_{QM} + E_{elstat})\) term is calculated simultaneously as one term by the QMprogram while the \(E_{MM}^{mod} = (E_{MM} + E_{vdW} + E_{covalent})\) term is calculated as one term by the MMprogram. The presence of the MM pointcharges during the QMcalculation has the effect of the QMcalculation sensing the electrostatic part of the environment, the QMdensity will be (mostly) correctly polarized and hence QM properties will also be more realistic. Electrostatic embedding is considered the standard QM/MM coupling scheme and is the default in ASH. It is more sophisticated than mechanical embedding and is usually the preferred choice for QM/MM calculations. The drawbacks of electrostatic embedding are :
It can only be used if the QMprogram supports pointcharge embedding (including gradients on pointcharges). ASH currently supports pointcharge embedding for programs: ORCA, CFour, MRCC, xTB, pySCF, NWChem, QUICK, CP2K, MNDO, TeraChem.
The presence of a large number of MM pointcharges in the QMcalculation can slow down the QMcalculation considerably. Especially the QMpointcharge gradient can be slow to calculate. See the TruncatedPC option below for a way to deal with this issue.
It requires some care in the handling of the covalent QM/MM boundary (see next sections on linkatoms, chargeshifting etc.)
More sophisticated polarized embedding approaches are not yet available in ASH.
QM/MM boundary treatment: linkatoms
If the QMregionMMregion boundary is between two bonded atoms, then a boundary correction needs to be applied as the QMregion will otherwise have a dangling bond, which would result in artifacts. In ASH this is treated by the popular linkatom method where a hydrogenlinkatom is added to cap the QMsubsystem. The hydrogen linkatoms are only visible to the QM theory, not the MM theory. The linkatoms are only used temporarily (automatically created and deleted) during the calculation of the QMpart and are never part of the system.
The need for a linkatom is automatically detected by ASH by noticing that 2 boundary atoms (QM1 and MM1) are bonded to each other according to connectivity information (determined by distances). ASH next places a hydrogen atom (H) along the bond axis between QM1 and MM1. The linkatom distance is determined according to which linkatom_method has been chosen. The standard linkatom_method = 'simple' option uses a fixed linkatom distance which is by default 1.09 Å (corresponds to a CH bond length). The default distance can be changed by setting the linkatom_simple_distance keyword in the QMMMTheory object. This simple fixedlinkatom distance method is simplistic but works well in most cases as long as the QMMM boundary is chosen well, i.e. the QMMM boundary is not through a polar bond but rather a nonpolar CC bond. See Tutorial: QM/MM boundary for more information on how to define a good QM/MM boundary for proteins.
An alternative linkatom_method option is linkatom_method = 'ratio' which calculates the linkatom position by scaling the difference between the QM1 and MM1 positions:
\(r_{L} = r_{QM1} + ratio*(r_{MM1}r_{QM1})\)
The linkatom_ratio is by default 0.723 but can be changed (linkatom_ratio keyword).
Finally, during a QM/MM gradient calculation there will be a gradient/force calculated on the (fictious) linkatom. This force is projected onto the QM1 and MM1 atoms to give the correct gradient for the QM/MM system. Two different forceprojections are available in ASH, controlled by the linkatom_forceproj_method keyword. The default is linkatom_forceproj_method = 'adv' which is an advanced projection of the linkatom force onto the QM and MM atoms while the alternative is linkatom_forceproj_method = 'lever' utilizes the simple lever rule to determine how the force should be projected onto QM1 and MM1. Both approaches give similar results.
Overall, the recommended way of using link atoms is to define the QMMM boundary for two carbon atoms that are as nonpolar as possible. In the CHARMM forcefield one should additionally make sure that one does not make a QMMM boundary through a chargegroup (check topology file). By default ASH will exit if you try to define a QMMM covalent boundary between two atoms that are not carbon atoms (since this is almost never desired). To override this behaviour add "unusualboundary=True" as keyword argument when creating QMMMTheory object.
In rare cases you may want to prevent ASH from adding a linkatom for a specific QMatom, e.g. if you are making unusual QMMM boundaries. This can be accomplished like below. Note, however, that the QMMM bonded terms will still be included.
#Excluding QMatom 5785 from linkatomcreation.
qmmmobject = QMMMTheory(qm_theory=orcaobject, mm_theory=openmmobject, fragment=frag, embedding="Elstat",
qmatoms=qmatoms, excludeboundaryatomlist=[5785])
General recommendations for biomolecular systems:
Special care should be taken when defining a QMregion for a biomolecular system
Always cut a CC bond that is as nonpolar as possible.
Focus on including nearby sidechains of residues that are charged (e.g. Arg, LYS, ASP, GLU) or are involved in important hydrogen bonding.
Amino acid sidechains are straighforward but make sure to not cut through CHARMM charge groups.
Including protein backbone is more involved and needs careful inspection. The only good option is typically to cut the CC bond between the C=O and the Calpha.
See Tutorial: QM/MM boundary for more information on how to define a good QM/MM boundary for proteins.
QM/MM boundary treatment: mechanical vs. electrostatic embedding
The chosen coupling scheme (mechanical vs. electrostatic) influences the treatment of the QM/MM boundary, including the linkatom handling. For mechanical embedding there is nothing besides the linkatomtreatment (see above) that needs to be done: the linkatoms are present during the QMcalculation but invisible to the MMpart and the linkatom force is projected onto the QM1 and MM1 atoms.
However, in electrostatic embedding, the presence of the linkatom, as well as a bonded MM atom being so close, created problems, that if not treated this would lead to some artifical overpolarization. To prevent this overpolarization, the atom charge of the MMatom is traditionally shifted towards its bonded neighbours (MM2 atoms) with a possible dipole correction also applied.
ASH includes 2 different chargeboundarymethods for preventing overpolarization at the QMMM boundary which are controlled by the chargeboundary_method keyword in the QMMMTheory object:
Chargeshift method
The chargeboundary_method = 'shift' option employs the popular chargeshifting strategy by Paul Sherwood and coworkers. See de Vries et al. J. Phys. Chem. B 1999, 103, 61336141. This method is the default in ASH when electrostatic embedding is chosen.
The charge of the MM1 atom is set to 0.0 and is shifted towards the MM2 atoms. Effectively, the original chargevalue of the MM1 is divided by the number of MM2 atoms bonded to the MM1 atom and each MM2 atom receives a fraction of the original MM1 charge. This chargeshifting has the effect of avoiding the overpolarization that would have occured in the QM1L and MM1 region while maintaining the overall charge of the system. The drawback, however, is that the MM1MM2 dipole is no longer correct which is why a dipole correction is also applied. The dipole correction adds extra pointcharges around the MM2 atom to compensate. In ASH the dipole correction is applied automatically by default but can be turned off ( dipole_correction=False).
RCD: Redistributed charge and Dipole scheme
The chargeboundary_method = 'rcd' option employs the RCD method by Donald Truhlar and coworkers. See Lin et al. J. Phys. Chem. A 2005, 109, 39914004.
The RCD method is similar to the chargeshift method above but has some additional flexibility and can sometimes give better results. It also involves setting the charge of the MM1 pointcharge to 0.0 and redistributing the charge away. However, instead of placing a fraction of the MM1 charge on the MM2 atoms the charges are instead placed along the MM1MM2 bond midpoints. This defines the RC (redistributed charge) method. The RCD method involves in addition, changing the values of the charges placed on the MM1MM2 bond midpoints to be twice as large as the divided MM1 chargefraction. Additionally the pointcharge on each MM2 atom is reduced by the same amount as the original MM1 chargefraction. This redistribution in the RCD method has the effect of preserving the MM1MM2 bond dipoles.
How QM/MM works behind the scenes in ASH
During a QM/MM energy+gradient calculation in ASH the following steps take place:
QM/MM program reads in the full model of the system, containing all atoms (no fake atoms,dummy atoms or linkatoms).
ASH determines connectivity of the system, i.e. finds what atoms are bonded to each other. E.g. atom no. 17 and atom no. 18 may be close enough that ASH thinks they are bonded.
The program parses the qmatoms list. The qmatoms list only contains real atoms (not linkatoms because they don‘t exist yet). The qmatoms list may e.g. contain atoms 1,2,3,4,14,15,16,17) and the program next checks if there is a covalent QMMM boundary. Since atom no. 17 was (according to step 2) bound to atom no. 18 (which is an MMatom) then that means we have a covalent QMMM boundary.
ASH will next automatically calculate the need for a linkatom for all QMMM boundaries. Any required linkatom (H) coordinates will be calculated and MM charges will be modified to account for the QMMM boundary.
The Cartesian coordinates of the QMatoms are taken and passed to the QMprogram. A hydrogen link atom is automatically added to this list of QMcoordinates (so that the QMsystem will not have a dangling boundary). Additionally the MM pointcharges (also passed to the QMprogram) are modified so that overpolarization will not occur (electron density at atom no. 17 and linkatom would be overpolarized by the closeness to MMatom no. 18 ). Additional MM charges are also added so that the dipole is more realistic.
The QMprogram calculates the energy and gradient of the QMatoms (with linkatoms) with the electrostatic effect of the MMatoms included (enters the 1electron Hamiltonian). The QM energy and gradient (of the QMatoms and also the PC gradient of the MM atoms) is passed back to ASH.
An MM calculation is performed for the whole system. The pointcharges of the atoms that have been labelled QMatoms have been set to zero to avoid calculating the electrostatic energy (would doublecount otherwise). Bonded MM terms for the same QMatoms are removed (if present, to avoid doublecounting). Bonded terms at the QMMM boundary are dealt with in a special way. The MM program never sees any linkatoms,chargeshifted MM charges or dipole charges.
The QM/MM energy is calculated by combining the contribution calculated by the QMprogram and the MMprogram. This will include the coupling energy of the QM and MM subsystems. Correction for artificial linkatom energy could be done here (not done in practice in ASH). The QM/MM gradient of the full system is assembled from the QMgradient, MMgradient and PCgradient. The gradient calculated on the dummy linkatom during the QMatoms (which does not exist in the real system) is taken and it is projected onto the relevant MM and QMatoms instead.
The complete QM/MM gradient of the whole system is used to make a step in the relevant job.
if part of a geometry optimization then the step is taken so as to minimize the QM/MM gradient
or: if part of a dynamics simulation is taken according to Newton‘s equations (the QM/MM gradient or force is used to calculate an acceleration which results in a change in velocity and positions of all the (real) atoms).
Note
Neither the geometry optimization or dynamics algorithms see or experience any linkatoms, only real atoms of the system.
QM/MM Truncated PC approximation
For large systems (e.g. > 50 000 atoms) the evaluation of the QMpointcharge interaction (calculated by the QMcode) will start to dominate the cost of the calculation in each QM/MM calculation step. The QMpointcharge gradient calculation is the main culprit and it depends on the QMcode how efficiently this step is carried out for a large number of pointcharges. ASH features a convenient workaround for this problem in QM/MM geometry optimizations. Instead of reducing the system size, ASH can temporarily reduce the size of the PC field (MM calculation size remains the same) during the geometry optimization which can speed up the calculation a lot. The size of the truncated PC field is controlled by the TruncPCRadius variable (radius in Å) which results in a truncated spherical PC field.
The algorith works like this:
Opt cycle 1:
Calculate truncated and full pointcharge field. Calculate gradient and energy correction.
Opt cycle n:
if Opt cycle n is a multiple of TruncatedPC_recalc_iter then:
Recalculate correction using both full pointcharge field and truncated.
else:
Use truncated PC field (defined by TruncPCRadius) in each QM run. Combine with energy and gradient corrections.
Final Opt cycle:
Recalculate final geometry using full pointcharge field.
In a typical truncatedPC QM/MM optimization, the full pointcharge field (e.g. 1 million PCs) is used in the 1st step (expensive) but in later steps an approximated spherical PCregion (cheap) is used during the QMsteps (e.g. a spherical 35 Å radius region) until step 50/100/150 etc. (if TruncatedPC_recalc_iter=50) where the full pointcharge field is recalculated. When the optimization converges, e.g step 80, a final energy evaluation is performed using the full PC field. For such an 80iteration job, the full PC gradient may be calculated only 3 times (instead of 80 times) that can result in considerable time savings.
Note that QM and QM/MM energies are approximate during the optimization steps where a truncated PC field is used. The final energy is always calculated using the full PC field. The error from the approximation depends on the TruncPCRadius parameter (smaller values than 30 not recommended) and TruncatedPC_recalc_iter (how often the full PC field is used). If TruncatedPC_recalc_iter=1 then no truncation is performed.
#QM/MM theory object defined with the truncated PC approximation
qmmm = QMMMTheory(qm_theory=qm, mm_theory=omm, fragment=frag, embedding="Elstat", qmatoms=qmatoms, printlevel=2,
TruncatedPC=True, TruncPCRadius=35, TruncatedPC_recalc_iter=50)
Example: QM/MM with ORCA and NonbondedTheory
Example for a H2OMeOH system where the MeOH is described by QM and H2O by MM. Here we read in a forcefieldfile containing a nonbonded forcefield (see MM Interfaces). The files for this example are available in the examples/QMMMexamples/QMMMORCAnonbondedtheory directory of the ASH repository.
from ash import *
#H2O...MeOH fragment defined. Reading XYZ file
H2O_MeOH = Fragment(xyzfile="h2o_MeOH.xyz")
# Specifying the QM atoms (38) by atom indices (MeOH). The other atoms (0,1,2) is the H2O and MM.
#IMPORTANT: atom indices begin at 0.
qmatoms=[3,4,5,6,7,8]
# Charge definitions for whole fragment. Charges for the QM atoms are not important (ASH will always set QM atoms to zero)
atomcharges=[0.8, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
#Defining atomtypes for whole system
atomtypes=['OT','HT','HT','CX','HX', 'HX', 'HX', 'OT', 'HT']
#Read forcefield (here containing LJpart only) from file
MM_forcefield=MMforcefield_read('MeOH_H2Osigma.ff')
#QM and MM objects
ORCAQMpart = ORCATheory(orcasimpleinput="!BP86 def2SVP def2/J tightscf", orcablocks="")
MMpart = NonBondedTheory(charges = atomcharges, atomtypes=atomtypes, forcefield=MM_forcefield,
LJcombrule='geometric', codeversion="py")
QMMMobject = QMMMTheory(fragment=H2O_MeOH, qm_theory=ORCAQMpart, mm_theory=MMpart, qmatoms=qmatoms,
charges=atomcharges, embedding='Elstat')
#Singlepoint energy calculation of QM/MM object
result = Singlepoint(theory=QMMMobject, fragment=H2O_MeOH, charge=0, mult=1)
print("Singlepoint QM/MM energy:", result.energy)
#Geometry optimization of QM/MM object (this may not converge)
result2 = Optimizer(fragment=H2O_MeOH, theory=QMMMobject, coordsystem='tric', ActiveRegion=True, actatoms=[3,4,5,6,7,8], charge=0, mult=1)
print("Optimized QM/MM energy:", result2.energy)
Example: QM/MM with ORCA and OpenMMTheory
See also QM/MM on a protein.
The files for this example (DHFR protein) are available in the examples/QMMMexamples/QMMMCHARMMexample directory of the ASH repository.
from ash import *
numcores=1
#Defining path to dir containing forcefield files and coordinates
forcefielddir="./"
psffile=forcefielddir+"step3_pbcsetup.psf"
topfile=forcefielddir+"top_all36_prot.rtf"
prmfile=forcefielddir+"par_all36_prot.prm"
xyzfile=forcefielddir+"coordinates.xyz"
#Read coordinates from XYZfile
frag = Fragment(xyzfile=xyzfile)
#Creating OpenMM object
openmmobject = OpenMMTheory(psffile=psffile, CHARMMfiles=True, charmmtopfile=topfile,
charmmprmfile=prmfile, periodic=True, periodic_cell_dimensions=[80.0, 80.0, 80.0, 90.0, 90.0, 90.0],
do_energy_decomposition=True, autoconstraints=None, rigidwater=False)
#Creating ORCATheory object
ORCAinpline="! HF3c tightscf"
ORCAblocklines="""
%maxcore 2000
"""
#Create ORCA QM object. Attaching numcores so that ORCA runs in parallel
orcaobject = ORCATheory(orcasimpleinput=ORCAinpline,
orcablocks=ORCAblocklines, numcores=numcores)
#act and qmatoms lists. Defines QMregion (atoms described by QM) and Activeregion (atoms allowed to move)
#IMPORTANT: atom indices begin at 0.
#Here selecting the sidechain of threonine
qmatoms = [569,570,571,572,573,574,575,576]
actatoms = qmatoms #Same active region as QMregion here
# Create QM/MM OBJECT by combining QM and MM objects above
qmmmobject = QMMMTheory(qm_theory=orcaobject, mm_theory=openmmobject, printlevel=2,
fragment=frag, embedding="Elstat", qmatoms=qmatoms)
#Run geometry optimization using geomeTRIC optimizer and HDLC coordinates. Using active region.
Optimizer(theory=qmmmobject, fragment=frag, ActiveRegion=True, actatoms=actatoms,
maxiter=500, coordsystem='hdlc', charge=0,mult=1)