Working with Trajectories

Objective

The objective of this tutorial is to show how to create QM input files using trajectory files. In this example you will see QMzyme can be used to generate multiple QMzymeModels from different trajectory snapshots.

It is recommended to first visit the QM-only calculation and QM/xTB calculation workflows, since both workflows are utilized. You can find it here for QM-only and here for QM/xTB.

This workflow allows you to:

Be able to extract specific frames from the trajectory (PQR and DCD) files.
Create QM input files from trajectory files

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:

Required Files

To start, you will need:

A fully prepped and protonated PDB.

[5]:

# Here are the necesary imports for this tutorial!

import QMzyme
from QMzyme.SelectionSchemes import DistanceCutoff
from QMzyme.data import PQR, DCD

Initialize Model and Define Regions

To begin our example workflow, we need to start by initializing QMzymeModel using GenerateModel().When loading PQR and DCD files, QMzyme will detect the number of frames

[6]:

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

# We can use this print statement to show how many frames are present.
print("Number of Frames: ", model.universe.trajectory.n_frames)

Number of Frames:  1

QM-only: Loop over frames to generate QMzymeModels

In this first scenario, you will generate five QM-only models using the DistanceCutoff selection scheme with a cutoff of 3 Angstroms. To achieve this, we need to first set the QM_Method(). From there, we can use the for loop to choose our frame range. In each of the snapshot, we can choose our QMzymeModel to be the structure snapshot at the given frame, which is used for set_catalytic_center and set_region to get our QM region. Afterwards, we can freeze C-alpha atoms with set_fixed_atoms, truncate(), and write_input.

[ ]:

# We first decide on the QM_Method for input file.
qm_method = QMzyme.QM_Method(
    basis_set='6-31G*',
    functional='wB97X-D3',
    qm_input='OPT FREQ',
    program='orca'
)

# We use for loop to extract snapshops every 10 frames across 50 frames.
for frame in range(0, 50, 10):
    print('\n====================================')
    print(f'             Frame {frame}')
    print('====================================')

    # For each 10 frames, we extract structure information and set it as QMzymeModel.
    m = QMzyme.GenerateModel(PQR, DCD, frame=frame)

    # We set bound ligand as the catalytic center.
    m.set_catalytic_center('resid 263')

    # DistanceCutoff class is used to get 3 Angstrom region.
    m.set_region(selection=DistanceCutoff, cutoff=3)

    # We rename the region.
    m.cutoff_3.rename(f'{m.cutoff_3.name}_frame_{frame}')

    # We set C-alpha atoms as fixed atoms.
    c_alpha_atoms = m.regions[-1].get_atoms(attribute='name', value='CA')
    m.regions[-1].set_fixed_atoms(atoms=c_alpha_atoms)

    # We assign QM_Method to the region.
    qm_method.assign_to_region(region=m.regions[-1])

    # We truncate the region using TerminalAlphaCarbon class.
    m.truncate()

    # And at last, we write the input!
    m.write_input(filename=f"{m.name}_cutoff3_frame{frame}_qm")

QM/xTB: Loop over frames to generate QMzymeModels

In this second scenario, you will generate five QM/xTB models with a QM region consisting of the ligand, EQU and the catalytic ASP103 residue and a large xTB region containing residues within 8 Angstroms of EQU. To achieve this, we need to first set QM_Method() and XTB_Method(). From there, we can use the for loop to choose our frame range. In each of the snapshot, we can choose our QMzymeModel to be the structure snapshot at the given frame, which is used for set_catalytic_center and set_region to get our QM region. Afterwards, we can truncate() and write_input.

[ ]:

# We first decide on the QM_Method for input file.
qm_method = QMzyme.QM_Method(
    basis_set='6-31G*',
    functional='wB97X-D3',
    qm_input='OPT FREQ',
    program='orca'
)
# xTB method does not need any arguments.
xtb_method = QMzyme.XTB_Method()

# We use for loop to extract snapshops every 10 frames across 50 frames.
for frame in range(0, 50, 10):
    print('\n====================================')
    print(f'             Frame {frame}')
    print('====================================')

    # For each 10 frames, we extract structure information and set it as QMzymeModel.
    m = QMzyme.GenerateModel(PQR, DCD, frame=frame)

    # We set bound ligand as the catalytic center.
    m.set_catalytic_center('resid 263')

    # We specify which regions will be treated with QM_Method
    m.set_region(name='qm_region', selection='resid 263 or resid 103')

    # DistanceCutoff class is used to get 3 Angstrom region.
    m.set_region(selection=DistanceCutoff, cutoff=8)

    # We rename the region.
    m.cutoff_8.rename(f'{m.cutoff_8.name}_frame_{frame}')

    # We assign QM_Method and XTB_Method to respective regions.
    qm_method.assign_to_region(region=m.qm_region)
    xtb_method.assign_to_region(region=m.cutoff_8)

    # We truncate the region using TerminalAlphaCarbon class.
    m.truncate()

    # And at last, we write the input!
    m.write_input(filename=f"{m.name}_cutoff8_frame{frame}_qmxtb")