Get started

Defining molecule

The definition of the molecule is done creating an instance of PyQchem Structure class. See PyQchem manual for more details. Charge and multiplicity are ignored.

Example of ethylene:

from pyqchem.structure import Structure

structure = Structure(coordinates=[[ 0.6695, 0.0000, 0.0000],
                                   [-0.6695, 0.0000, 0.0000],
                                   [ 1.2321, 0.9289, 0.0000],
                                   [ 1.2321,-0.9289, 0.0000],
                                   [-1.2321, 0.9289, 0.0000],
                                   [-1.2321,-0.9289, 0.0000]],
                      symbols=['C', 'C', 'H', 'H', 'H', 'H'])

Defining GROMACS input parameters

GROMACS input parameters are defined in a dictionary. For example:

gmx_params = {
             'integrator': 'md-vv',     # Verlet integrator
             'nsteps': 10000,           # 0.001 * 10000 = 100 ps
             'dt': 0.001,               # time step, in ps
             # Temperature coupling is on
             'tcoupl': 'nose-hoover',    # Nose-Hoover thermostat
             'tau_t': 0.3,               # time constant, in ps
             'ref_t': 300,               # reference temperature, one for each group, in K
             # Bond parameters
             'gen_vel': 'yes',           # assign velocities from Maxwell distributio
             'gen_temp': 300,            # temperature for Maxwell distribution
             'gen_seed': -1,             # generate a random seed
             }

The defined parameters override/add to the default values set by gromorg. This way for most purposes it is not necessary to define most parameters. These default values are:

'integrator': 'md-vv',     # Verlet integrator
'nsteps': 5000,            # 0.001 * 5000 = 50 ps
'dt': 0.001,               # ps
# Output control
'nstxout': 1,              # save coordinates every 0.001 ps
'nstvout': 1,              # save velocities every 0.001 ps
'nstenergy': 1,            # save energies every 0.001 ps
'nstlog': 100,             # update log file every 0.1 ps
# Bond parameters
'continuation': 'no',       # first dynamics run
'cutoff-scheme': 'Verlet',  # Buffered neighbor searching
'verlet-buffer-tolerance': 3.3e-03,
# 'ns_type': 'grid',          # search neighboring grid cells
'nstlist': 10,              # 20 fs, largely irrelevant with Verlet
'rcoulomb': 1.0,            # short-range electrostatic cutoff (in nm)
'rvdw': 1.0,                # short-range van der Waals cutoff (in nm)
'DispCorr': 'EnerPres',     # account for cut-off vdW scheme
# Electrostatics
'coulombtype': 'PME',       # Particle Mesh Ewald for long-range electrostatics
'pme_order': 4,             # cubic interpolation
'fourierspacing': 0.16,     # grid spacing for FFT
# Temperature coupling is on
'tcoupl': 'nose-hoover',    # Nose-Hoover thermostat
'tc-grps': 'system',        # one coupling group
'tau_t': 0.3,               # time constant, in ps
'ref_t': 100,               # reference temperature, one for each group, in K
# Pressure coupling is off
'pcoupl': 'no',             # no pressure coupling in NVT
# Periodic boundary conditions
'pbc': 'xyz',               # 3-D PBC
# Velocity generation
'gen_vel': 'yes',           # assign velocities from Maxwell distributio
'gen_temp': 10,             # temperature for Maxwell distribution
'gen_seed': -1,             # generate a random seed

Setting up the calculation

Example of simple parallel(openMP) calculation using 4 threads:

calc = GromOrg(structure,
               params=gmx_params,        # MDP parms
               box=[10, 10, 10],         # unitcell a, b, c in angstrom
               angles=[90, 90, 90],      # unitcell alpha, beta, gamma in degree
               supercell=[3, 3, 3],      # size of supercell
               delete_scratch=True,      # if true delete temp files when finished (default: True)
               silent=False,             # if true print MD log info in screen (default: False)
               omp_num_threads=False,    # number of parallel threads used (default: 1)
               maxwarn=0,                # max number of GROMACS warnings (default: 0)
               )

The defined structure correspond to the unit cell and it is replicated according to the supercell argument. This means that while usually will contain a single molecule, it is possible to include several units particularly to define crystals. The number of MPI processes is automatically set according to the cores available and the omp_num_threads parameter. The calculation is run in a temporary folder in the current working directory. This directory is deleted by default when the calculation is finished. delete_scratch keyword can be set to change this behavior.

maxwarn keyword can be used to set the maximum number of warnings GROMACS can ignore before aborting the calculation. It is recommended to keep this value to 0 (default) and only change it once you have revised all warnings and made sure you can ignore them.

Run the calculation

Once the calculation is set up, it can be run with the following command:

trajectory, energy = calc.run_md(whole=True)  # unwraps the trajectory

The return of this function is the trajectory as a MDtraj object and the energy as a dictionary. MDtraj is a flexible format to store trajectory data. Check the documentation of MDtraj for more information. (https://www.mdtraj.org/1.9.5/load_functions.html). The only argument of this function is whole, this optional arguments controls the unwrapping of the trajectory. That is, if whole is True, the trajectory is unwrapped such that the molecules are shown as whole for each step of the trajectory.

A simple way to visualize the trajectory is to store it in the disk as a common format. This can be done using save method:

trajectory.save('trajectory.gro')

MDtraj supports different formats, such as GROMACS (gro), Protein Data Bank (pdb) and xyz.

Energy dictionary contains the total energy, the kinetic energy and the potential energy as lists. This can be plotted, for example, as:

import matplotlib.pyplot as plt
plt.plot(energy['potential'], label='potential')
plt.plot(energy['kinetic'], label='kinetic')
plt.plot(energy['total'], label='total')
plt.legend()
plt.show()

Extracting molecules from trajectory

One of the main applications of Gromorg is use the geometries of the molecules within the trajectory. For this purpose gromorg package provides some functions to do this. mdtraj_to_pyqchem function is used to extract a particular molecule from a trajectory.

from gromorg.tools import mdtraj_to_pyqchem

molecule = mdtraj_to_pyqchem(trajectory, ifame, ires, center=True)

where trajectory is the trajectory as a MDtraj object, ifame is the index of the frame and ires is the index of the residue. center is an optional argument that centers the molecule in the origin. Each residue correspond to a unit cell (structure in the first example), this means that in the case that the unit cell contains multiple molecules, the separation of them will need to be done manually.

The return of this function is a PyQChem Structure containing the molecule. This object can be used directly with PyQchem to perform quantum chemistry calculations. For example:

from gromorg.tools import mdtraj_to_pyqchem
from pyqchem import QchemInput, get_output_from_qchem

qc_input = QchemInput(molecule,
                      jobtype='sp',
                      exchange='hf',
                      basis='6-31G')

output = get_output_from_qchem(qc_input)

print(output)