QM-only Calculation

Objective

The objective of this tutorial is to use QMzyme to create a truncated active site composed of residues within 5 Å of the bound ligand. All C-alpha atoms are set to be frozen during geometry optimization.

This workflow allows you to:

• Acquire basic understanding behind QM-only input file generation

• Understand the foundation behind each

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

• 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.data import PDB

Initialize Model and Define Regions

To begin our example workflow, we need to start by initializing QMzymeModel using GenerateModel().When loading the PDB file, we quickly run into an issue: a nonconventional residue is found. This simply means that there is a non-AMBER protein residue with an unknown charge. This can be resolved using residue_charges.update(), where we can manually update the charge of an unknown residue. Afterwards, set_catalytic_center() is used to designate a catalytic center in QMzymeModel. This is utilized in the SelectionScheme class to acquire the desired QM region.

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

# Since our bound ligand, equilenin, is not a standard AMBER protein residue,
# we need to specify the chage.
QMzyme.data.residue_charges.update({'EQU': -1})

# This is used to set catalytic center.
model.set_catalytic_center(selection='resid 263')

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).

For this workflow, we will utilize the DistanceCutoff class, which selects residues based on their distance from a user-defined catalytic center. Specifically, the workflow identifies atoms within a specified cutoff distance of the catalytic center and determines which residues should be included in the resulting region.

# Import selection scheme 'DistanceCutoff' and use it to create 5 Å region around catalytic center.
from QMzyme.SelectionSchemes import DistanceCutoff
model.set_region(selection=DistanceCutoff, cutoff=5)
print(model.regions)

[<QMzymeRegion catalytic_center contains 37 atom(s) and 1 residue(s)>, <QMzymeRegion cutoff_5 contains 427 atom(s) and 33 residue(s)>]

Designate Coordinate Constraints

After acquiring the region, it is important to set a fixed atom. This is done by selecting a specific atom of choice and freezing them using set_fixed_atoms(). For our case, we will freeze C-alpha atoms for geometry optimization. To achieve, this, we will create a selection of atoms with the name CA, and freezing these atoms.

# 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)

print("Fixed atoms: ",model.cutoff_5.get_atoms(attribute='is_fixed', value=True))

Fixed atoms:  [<QMzymeAtom 235: CA of resname TYR, resid 16>, <QMzymeAtom 256: CA of resname ILE, resid 17>, <QMzymeAtom 309: CA of resname VAL, resid 20>, <QMzymeAtom 595: CA of resname ASP, resid 40>, <QMzymeAtom 615: CA of resname PRO, resid 41>, <QMzymeAtom 851: CA of resname TYR, resid 57>, <QMzymeAtom 896: CA of resname GLN, resid 59>, <QMzymeAtom 913: CA of resname GLY, resid 60>, <QMzymeAtom 920: CA of resname LEU, resid 61>, <QMzymeAtom 982: CA of resname VAL, resid 66>, <QMzymeAtom 1022: CA of resname ALA, resid 68>, <QMzymeAtom 1225: CA of resname MET, resid 84>, <QMzymeAtom 1256: CA of resname PHE, resid 86>, <QMzymeAtom 1300: CA of resname VAL, resid 88>, <QMzymeAtom 1331: CA of resname MET, resid 90>, <QMzymeAtom 1461: CA of resname LEU, resid 99>, <QMzymeAtom 1492: CA of resname VAL, resid 101>, <QMzymeAtom 1527: CA of resname ASH, resid 103>, <QMzymeAtom 1743: CA of resname MET, resid 116>, <QMzymeAtom 1777: CA of resname ALA, resid 118>, <QMzymeAtom 1808: CA of resname TRP, resid 120>, <QMzymeAtom 1888: CA of resname LEU, resid 125>]

Build the QM Method and Assign to Region

Ultimately, QMzymeModel is a collection of QMzymeRegion and a QM method for calculation. This can be achieved by utilizing QM_Method() and assigning QM method to the 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.

For this workflow, we’ll show how to write for both ORCA and Gaussian. When looking at two methods, you can see that the only difference between two methods are program, where one is ORCA and one is Gaussian. This attribute is used to specify the input file format.

# Use this for creating ORCA input file
qm_method = QMzyme.QM_Method(
    basis_set='6-31G*',
    functional='wB97X-D3',
    qm_input='OPT FREQ',
    program='orca'
)

# Use this for creating Gaussian input file
#qm_method = QMzyme.QM_Method(
#    basis_set='6-31G*',
#    functional='wB97XD',
#    qm_input='OPT FREQ',
#    program='gaussian'
#)

After setting the method, we have to assign the method to the region. This can be done with assign_to_region() with your region of choice! When doing so, QMzyme will estimate the charge based on the commbination of charge information provided by the user with residue_charge.update() and standard AMBER residue charges. Alternatively, you can specify the charge and multiplicity of the system with charge= and mult=.

# Since we are not specifying the charge in this method below, the method
# will estimate the charge based on residue naming conventions
qm_method.assign_to_region(region=model.cutoff_5)

# You would alternatively uncomment and run:
#qm_method.assign_to_region(region=model.cutoff_5, charge=-2, mult=1)

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.

model.truncate()
print("Truncated model: ", QMzyme.CalculateModel.calculation['QM'])
print("Method details: ", QMzyme.CalculateModel.calculation['QM'].method)

Truncated model has been created and saved to attribute 'truncated' and stored in QMzyme.CalculateModel.calculation under key QM. This model will be used to write the calculation input.
Truncated model:  <QMzymeRegion cutoff_5_truncated contains 391 atom(s) and 33 residue(s)>
Method details:  {'type': 'QM', 'qm_input': '6-31G* wB97X-D3 OPT FREQ', 'basis_set': '6-31G*', 'functional': 'wB97X-D3', 'qm_end': '', 'program': 'orca', 'freeze_atoms': [2, 23, 42, 58, 78, 84, 105, 122, 129, 148, 164, 174, 191, 211, 227, 244, 263, 279, 292, 309, 319, 343], 'mult': 1, 'charge': -2}

Write QM Input File

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.

model.write_input()