QM/xTB Calculation

Objective

The objective of this tutorial is to create a multi-level ORCA input file using QM/xTB method. In this example you will define a catalytic center as the bound ligand, residue 263 (EQU), and a second region containing atoms within 5 Å of the catalytic center using the Selection Scheme DistanceCutoff. All alpha carbons will be constrained for optimization using set_fixed_atoms(), a method of the QMzymeRegion class. You will set a high-level QM method for the catalytic center region, and a lower-lever xTB (XTB) method to treat the larger system. After assigning the calculation methods you will truncated the model using the default Truncation Scheme, AlphaCarbonTerminal.

It is recommended to first visit the QM-only calculation workflow since it provides a more basic understanding of the QMzyme workflow. You can find it here.

This workflow allows you to:

• Create QM/xTB input file using QMzyme.

In this specific example, we are using ketosteroid isomerase (KSI) as the model system. The structure for KSI is obtained from the PDB 1OH0 and MM-minimized prior to this tutorial. We will utilize it to create an input file for ORCA for geometry optimization and frequency calculation.

Classes used in this example

• GenerateModel

• QM_Method

• XTB_Method

• CA_terminal TruncationScheme

• DistanceCutoff SelectionScheme

Required Files

To start, you will need:

• A fully prepped and protonated PDB.

---

# Here are the necesary imports for this tutorial!

import QMzyme
from QMzyme.SelectionSchemes import DistanceCutoff
from QMzyme.data import PDB

Initialize Model and Define Regions

For doing QM/xtB calculation, we need a QMzymeModel with two regions, one QM method, and xTB method. Thus, the first step is to create 2 distinct QMzymeRegion. For our enzyme, we will utilize the bound ligand (residue 263) and the 5 Å distance cutoff region as our two regions for QM/QM2 calculation. To initialize and define our regions, we will first load our structure file and set regions using set_catalytic_center() and set_region(). You can find more information about DistanceCutoff class here.

# We first initilize QMzymeModel by loading structure file.
# The initial pdb file should be prepared prior (hydrogens must be present).
model = QMzyme.GenerateModel(PDB)

# This adds unknown residue charge.
QMzyme.data.residue_charges.update({'EQU': -1})

# We select first region using set_catalytic_center and update the charge of the region.
model.set_catalytic_center(selection='resid 263')
model.catalytic_center.set_charge(-1)

# We use DistanceCutoff class to create our second region.
model.set_region(selection=DistanceCutoff, cutoff=5)

# Selecting and freeznig CA atoms.
c_alpha_atoms = model.cutoff_5.get_atoms(attribute='name', value='CA')
model.cutoff_5.set_fixed_atoms(atoms=c_alpha_atoms)

Charge information not present. QMzyme will try to guess region charges based on residue names consistent with AMBER naming conventions (i.e., aspartate: ASP --> Charge: -1, aspartic acid: ASH --> Charge: 0.). See QMzyme.data.residue_charges for the full set.

        Nonconventional Residues Found
        ------------------------------
        EQU --> Charge: UNK, defaulting to 0

You can update charge information for nonconventional residues by running
        >>>QMzyme.data.residue_charges.update({'3LETTER_RESNAME':INTEGER_CHARGE}).
Note your changes will not be stored after you exit your session. It is recommended to only alter the residue_charges dictionary. If you alter the protein_residues dictionary instead that could cause unintended bugs in other modules (TruncationSchemes).

Build the QM Method and Assign to Region

For the QM region, you need to define QM method and designate the method to region. There are multiple attributes that can be given to the method, which will impact the input file! Some of the essential options are:

• basis_set: defines the basis set to use for calculation.

• functional: defines the functional to use for calculation.

• qm_input: keywords to include in the input file route line to declare any details.

• program: defines the program that will be used for QM calculation.

You can find more information about QM_Method here.

# We first create two QM method, each with desired level of theory.
qm_method = QMzyme.QM_Method(basis_set='6-31G*',
               functional='wB97X-D3',
               qm_input='OPT FREQ',
               program='orca')

# Afterwards, we assign two methods to their own respective region.
qm_method.assign_to_region(region=model.catalytic_center)

Initialize the xTB method and Assign to Region

For the xTB region, we can use XTB_Method(). In QMzyme, XTB_Method() does not require any arguments to function. The only thing that is needed is to assign XTB_Method() it to a region using assign_to_region().

# no arguments are passed to initialize XTB_METHOD.
QMzyme.XTB_Method().assign_to_region(region=model.cutoff_5)

QMzymeRegion cutoff_5 has an estimated charge of -2.

Truncate Model

The last step before creating an input file is to use model.truncate(). This takes the full QMzymeModel (which, if you recall, is a region with a method), truncates it, and prepares it for input file generation. If not specified, model.truncate() will undergo TerminalAlphaCarbon truncation class, which will preserves peptide bonds between consecutive residues and remove terminal backbone amine and carboxylate groups. You can find more information about it here.

# We truncate the model with TerminalAlphaCarbon class.
model.truncate()

QMzymeRegion catalytic_center_cutoff_5_combined has an estimated charge of -2.

Truncated model has been created and saved to attribute 'truncated' and stored in QMzyme.CalculateModel.calculation under key QMXTB. This model will be used to write the calculation input.

Write QM Input File for ORCA

The final step now is to create the calculation input file. This can be done by using model.write_input() method, which will create a directory named QCALC and save the input file there. When this method is called, the store_pickle() method is triggered, and you will have the corresponding .pkl file for your QMzymeModel. This is especially useful because after your calculation is completed, you can then load the calculation results back into the corresponding QMzymeRegion using the QMzymeRegion method store_calculation_results(), passing it the calculation file.

Be mindful that xTB method is exclusive to ORCA, thus you can only generate QM input file for ORCA.

model.write_input()