-
Notifications
You must be signed in to change notification settings - Fork 3
Using CodeEntropy as a module
This example demonstrates the use of CodeEntropy as a python package to compute the vibrational entropy of a molecule by Multiscale-cell-correlation method using its API features. The example shows entropy calculation at 3 different levels of hierarchy - molecule (M) level, nucleotide (N) / residue (R) level and united-atom (UA) level.
The trajectory was obtained using GROMACS. All that is needed is the TPR file and a TRR file with both, coordinates and force information. Its crucial that the number of atoms in the TPR and TRR files are identical.
Import neccesary packages.
import os, sys
import numpy as nmp
from CodeEntropy.Reader import GromacsReader
from CodeEntropy.ClassCollection import Atomselection as SEL
from CodeEntropy.ClassCollection import BeadClasses as BC
from CodeEntropy.ClassCollection import ModeClasses as MC
from CodeEntropy.FunctionCollection import Utils
from CodeEntropy.FunctionCollection import UnitsAndConversions as UAC
from CodeEntropy.FunctionCollection import CustomFunctions as CF
from CodeEntropy.FunctionCollection import EntropyFunctions as EF
Assign tprfile and trrfile with the appropriate file names. Former is the binary format topology/structure file and latter the binary-format trajectory file with coordinates and forces. Also assign a file name to outfile on which the output will be written.
Use the reader function to read the data into an object of class BaseMolecule and simultaneously into an object of Class DataContainer.
mainMolecule, dataContainer = GromacsReader.read_gromacs_input(arg_tprFile=tprfile, \
arg_beginTime= -1, \
arg_stride=-1,\
arg_endTime=-1, \
arg_trajFile=trrfile,\
arg_verbose= 3, \
arg_outFile=outfile)
The snippets below work only if dataContainer (object of class DataContainer) and mainMolecule (object of class BaseMolecule) are properly defined and populated.
Entropy at this level is obtained by treating all the atoms in the molecule as one cell which has three translational and 3 rotational degrees of freedom. The script and the purpose they serve (in comments) are shown below:
Read the nuber of frames read (optional) into a variable for use later in the program.
numFrames = len(dataContainer.trajSnapshots)
Utils.printflush(f'Total number of frame = {numFrames}')
A level is characetrised or rather depicted in the program as a collection of Cell(s) or Bead(s). Create an instance of BeadCollection by providing a name for it as well as connecting with the dataContainer variable. The create a single bead representing the whole dna molecule. This is an object of Class Bead which requires, at the outset, the indices of atoms that will comprise it along with other descriptors. Add it to the BeadCollection.
# Define a bead collection for this level
mlevel = BC.BeadCollection("mlevel", dataContainer)
allSel = SEL.Atomselection(mainMolecule, "ALL")
# create a bead
wholeDNABead = BC.Bead(arg_atomList=allSel.get_indices(), \
arg_numFrames=numFrames, \
arg_hostDataContainer = dataContainer,\
arg_beadName = "WMOL",
arg_beadResi = 0,
arg_beadResn = "WMOL",
arg_beadChid = "X")
# add the bead to the bead colleciton
mlevel.listOfBeads = [wholeDNABead]
Utils.printflush(f"Total number of beads at M-level = {len(mlevel.listOfBeads)}")
Before starting computations, it is useful to reset weighted vectors for each bead and and assign a representative position (would be used for visualization of modes) using some visualiation program like VMD. These weighted vectors will be later used to store the mass-weighted forces and inertia-weighted torques of each cell or bead making up the level.
for iBead in mlevel.listOfBeads:
iBead.reset_totalWeightedVectors( (numFrames,3) )
iBead.position = iBead.get_center_of_mass_lab(arg_frame = 0)
It is also recommended that the force and torque covariance matrices be reset using an inbuilt function that has the right dimensions when the BeadCollection was initialzed.
mlevel.reinitialize_matrices()
Setup translational and rotational axes for all the atoms at this level. Setting up these axes systems must be done at any level of hierarchy. Essentially, to each atom across all beads at that level, a translational and a rotational axes are assigned per frame.
At the M-level, prinicpal axes of the whole molecule centered at its center of mass is used. User has the liberty to chose any axes system they find appropriate. However, an axes system is a matrix of shape 4 x 3 with the first three rows containing orthonormal vectors (in a right-handed order) and the 4th row the origin of the axes system.
Inbuilt functions are written in the code that wrap the process of assigning the translational and rotational axes to specific atoms at desired frames.
# USE whole molecule's principal axes system (which changes every frame) for each atom
allAtoms = allSel.get_indices()
for iFrame in range(numFrames):
selMOI, selAxes = mainMolecule\
.get_principal_axes(arg_atomList = allAtoms,arg_frame = iFrame, arg_sorted=False)
selCOM = mainMolecule\
.get_center_of_mass(arg_atomList = allAtoms, arg_frame = iFrame)
dataContainer.update_translationAxesArray_at(arg_frame = iFrame, arg_atomList = allAtoms, arg_pAxes = selAxes, arg_orig = selCOM)
dataContainer.update_rotationAxesArray_at(arg_frame = iFrame, arg_atomList = allAtoms, arg_pAxes = selAxes, arg_orig = selCOM)
After the translational and rotational axes of the desired atoms (all here) at desired frames (all here) have been determined and stored in the dataContainer, they are used to determine the local coordinates and forces acting on each atom. Local coordinates are essential to compute the torques within the framework of the rotational axes system assigned to each atom. Inhouse wrappers are written to simplify the process as shown below. Essentially, the coordinates from the lab-frame are used to obtain the coordinates in the local-frame.
The snippet below demands that local coordinates be obtained with respect to the rotational axes using the token "R". At M-level, the translational and rotational axes are the molecular principal axes, this could have easily been replaced by token "T" (for translational axes).
# update local coordinates.
# using the rotational axes.
dataContainer.update_localCoords_of_all_atoms(arg_type="R")
The next step is calculating the local forces and torques acting on each atom using the translational and rotational axes respectively. For this, the value of forces in the lab-frame is transformed into force in the local translation frame and torques are computed by transforming the same forces in the local rotational frame instead of the translational frame and using the local coordinates (
dataContainer.update_localForces_of_all_atoms(arg_type="T")
The above snippet determines and populates the array storing local forces using a wrapper. The snippet below computes the local torques without using any wrapper function but ne can be written for it.
for iFrame in range(numFrames):
for iAtom in allSel.get_indices():
coords_i = dataContainer.localCoords[iFrame, iAtom]
forces_i = dataContainer.localForces[iFrame, iAtom]
dataContainer.localTorques[iFrame,iAtom] = CF.cross_product(coords_i,forces_i)
After obtaining the local forces and torques, they are vectorially summed and then weighted by square root of the total mass and moment-of-inertia respectively for each bead. Weighting is done to normalize and homogenise the units of forces and torques.
for iBead in mlevel.listOfBeads:
# mass weighting the forces for each bead (iBead) in each direction (j)
# inertia weighting the torques for each bead (iBead) in each direction (j)
for iFrame in range(numFrames):
# define local basis as the rotationalAxes of the first atom in the atomList of iBead
# doesnt matter because they all have the same R and T axes
iLocalBasis = dataContainer.rotationAxesArray[iFrame][iBead.atomList[0]]
#get the moment of inertia tensor for ibead in thid local basis
beadMOITensor = iBead.get_moment_of_inertia_tensor_local(arg_localBasis = iLocalBasis, arg_frame = iFrame)
# get total force and torque in each direction and weight them
for j in range(3):
iBead.totalWeightedForces[iFrame,j] = nmp.sum(nmp.asarray([dataContainer.localForces[iFrame,aid,j] for aid in iBead.atomList]))
iBead.totalWeightedForces[iFrame,j] /= nmp.sqrt(iBead.get_total_mass())
iBead.totalWeightedTorques[iFrame,j] = nmp.sum(nmp.asarray([dataContainer.localTorques[iFrame,aid,j] for aid in iBead.atomList]))
iBead.totalWeightedTorques[iFrame,j] /= nmp.sqrt(beadMOITensor[j,j])
The steps from this point onward are common to all levels. It starts with compiling the force and torque covariance matrices.
# now fill in the matrices
mlevel.update_subMatrix(arg_pairString="FF", arg_verbose=3)
mlevel.update_subMatrix(arg_pairString="TT", arg_verbose=3)
# make quadrant from subMatrices
# FF and TT quadrants must be symmetric
ffQuadrant = mlevel.generate_quadrant(arg_pairString="FF",arg_filterZeros=1)
ttQuadrant = mlevel.generate_quadrant(arg_pairString="TT",arg_filterZeros=1)
Use of in-house functions to compile and build matrices is apparent in the snippet above. These functions are member functions of class BeadCollection and are therefore generally appilcable to any level. Essentially, the weighted force and torque arrays are used to determine the elements of the covariance matrices.
One other thing to note is the use of arg_filterZeros=1 argument in the generate_quadrant() function. When set to 1, rows (and columns becaause of the symmetric nature of the matrices) which are all zeroes are eliminated. Such can happen when a bead has zero moment of inertia. If set to 0, zero rows (and columns) are retained.
To allow halving of forces and torques, one can also do the following:
# scale forces/torques of these quadrants
fScale = 0.5
tScale = 0.5
ffQuadrant = nmp.multiply(fScale**2, ffQuadrant)
ttQuadrant = nmp.multiply(tScale**2, ttQuadrant)
The scaling terms are squared because each element in the matrix is a product of two force or two torque components.
After the matrices are created, they should be diagonalized. An inhouse wrapper function exists for that purpose alone as part of the Utils module.
#diagnolaize
lambdasFF, eigVectorsFF = Utils.diagonalize(ffQuadrant)
lambdasTT, eigVectorsTT = Utils.diagonalize(ttQuadrant)
We recommend that the units of the eigen values from the above operation be changed to the SI units but it is optional. If one chooses to do so, however, they can use the inhouse change_lambda_unit() function that converts units strictly from
lambdasFF = UAC.change_lambda_units(lambdasFF)
lambdasTT = UAC.change_lambda_units(lambdasTT)
Diagonalization yields eigen values and eigen modes that collectively produce the vibrational motion. These modes and their frequencies make up a spectrum. On occassions it is required that certain modes of the spectrum be eliminated when computing entropy to avoid double-counting, e.g. removing the 6 smallest modes in terms of their frequencies. That requires ordering these modes.
Objects of class Mode store these information, can be used to deal with specific regions of the spectrum and sort them. Every mode is assigned a unique object of this class where it's frequency, eigen value and eigen mode is stored. Vibrational modes from translational (forces) as well as rotational (torques) can have their separate spectra.
For translational vibrational modes the following can be written. Temperature (in K) is required for these computations.
temper = 300 #K
modeSpectraFF = []
for midx, mcombo in enumerate(zip(lambdasFF, eigVectorsFF)):
fflmb, evec = mcombo
# compute mode frequencies
# nu = sqrt(lambda/kT)*(1/2pi)
# Units: 1/s
mfreq = EF.compute_frequency_from_lambda(fflmb, temper)
newMode = MC.Mode(arg_modeIdx = midx + 1, \
arg_modeEval = fflmb, \
arg_modeEvec = evec, \
arg_modeFreq = mfreq)
newMode.modeAmpl = EF.compute_ampfac_from_lambda(fflmb, temper)
modeSpectraFF.append(newMode)
# ... same for rotational modes
Setting up modes in form of a python list can prove useful when one is interested in specific modes only.
Sorting of modes comes in handy when certain modes are to be eliminated based on their frequency values. The following snippet shows how that can be done (for the translational modes).
modeSpectraFF = MC.sort_modes(modeSpectraFF)
Entropy is computed from the eigen modes (their frequencies) using the QM harmonic oscillator model. This is implemented via calculate_entropy_per_dof() function in EntropyFunctions module of FunctionCollection module.
The code below shows how all the degrees of freedom are used to obtain the total translational componenent of vibrational entropy. If the conversion of units shown far above has been done, the entropy is rendered in
s = 0
for m in modeSpectraFF:
s += EF.calculate_entropy_per_dof(m.modeFreq, temper)
Utils.printflush(f"{'FF Entropy':<40s} : {s:.4f} J/mol/K")