IRed Analysis of MD Simulations

From Powers Wiki
Revision as of 16:48, 14 May 2021 by Tessa17 (talk | contribs)

Calculating order parameters from MD Simulations

The analysis utilizes iRed software provided by the Dr. Rafael Bruschweiler Group. All iRed calculations were completed on the HCC cluster.

  1. Create a pdb file of the trajectory in the same manner as the .xtc file.
  2. After the pdb file is created (will be very large), use the following scripts to extract the N and H coordinates for each atom at each time point in the trajectory
    1. awk ‘{if (($3 == “N”) && ($4 != “PRO”)) {print $6, $7, $8, $9:}} xxx.pdb > N-HCoor.dat
    2. awk ‘{if (($3 == “H”) && ($4 != “H1”)) {print $6, $7, $8, $9:}} xxx.pdb > H-NCoor.dat
  3. Once .dat files are created, use the python script below to run the iRed calculation
    1. In command line, run “ired_1vec_block.py --coor1 N-HCoor.dat --coor2 H-NCoor.dat --block 100


       ! /usr/bin/env python
       from __future__ import division
       import numpy as np
       from optparse import OptionParser
       import sys
       def setupParserOptions():
               parser=OptionParser()
               parser.add_option('--coor1', dest='coor1',default='H-NCoor.dat',help = "Input coordinate file 1")
               parser.add_option('--coor2', dest='coor2',default='N-HCoor.dat',help = "Input coordinate file 2")
               parser.add_option('--block', dest='block', default = 1, help = 'Number of blocks')
               return parser


def ired(hCoor_file, nCoor_file, block):

   try:
       h_all = np.loadtxt(hCoor_file)
   except IOError:
       print 'cannot open ', hCoor_file
   try:
       n_all = np.loadtxt(nCoor_file)
   except IOError:
       print 'cannot open ', nCoor_file
   numLines = len(h_all)
   numRes = len(set(h_all[:,0]))
   numSnapshot = numLines // numRes
   blockNum = int(block)
   blockSize = numSnapshot // blockNum
   if numLines != len(n_all):
       print 'Error: The two coordinate files do not have same length.'
       return
   if np.mod(numLines, numRes) != 0:
       print 'Error: Coordinate files do not formatted correctly.'
       return
   # initialize s2
   s2_sum = np.zeros([1, numRes])
   for i in range(blockNum):
       hC = h_all[i*numRes*blockSize:(i+1)*numRes*blockSize, :]
       nC = n_all[i*numRes*blockSize:(i+1)*numRes*blockSize, :]
       # Calculate the bond vector orientations over the trajectory
       hnC = hC[:,1:] - nC[:,1:]
       # residue index
       resNum = hC[0:numRes, 0]
       # Construct matrix m
       mMat = np.zeros([numRes, numRes])
       #print numRes, numSnapshot, blockNum, blockSize
       for a in range(0, numRes*blockSize, numRes):
           # N-H vectors for the current snapshot
           hnTemp = hnC[a:a+numRes, :]
           # normalize the vectors
           for b in range(numRes):
               hnTemp[b,:] = hnTemp[b,:] / np.linalg.norm(hnTemp[b,:])
           # Construct Rank 2 Legendre polynomial
           c = np.dot(hnTemp, hnTemp.T)
           mMat = mMat + (3*np.multiply(c, c)-1)/2
       # Ensemble averaging
       mMat = mMat / blockSize
       #print np.shape(hnTemp), np.shape(hnTemp.T), np.shape(c), np.shape(mMat)
       # Diagonalize matrix
       dVals, v = np.linalg.eig(mMat)
       Ind = dVals.argsort()[::-1] # sort eigenvalues in descending order
       dSort = dVals[Ind]
       eigVal = np.diag(dSort)
       eigVec = v[:, Ind]
       # Calculate the squared order parameters
       s2_block = np.zeros(numRes)
       for b in range(numRes):
           sumOverModes = 0;
           for a in range(5,numRes):
               sumOverModes += eigVal[a,a] * (eigVec[b,a])**2
           s2_block[b] = 1 - sumOverModes
       s2_sum += s2_block
   s2 = s2_sum / blockNum
   # output
   out = np.concatenate(([resNum], s2), axis = 0)
   with open('ired_s2_1vec_%dblock.out' %(block), 'wb') as f:
       np.savetxt(f, np.transpose(out), fmt='%i %.3f')


if __name__ == '__main__':

   parser = setupParserOptions()
   if len(sys.argv) <2:
       parser.print_help()
       sys.exit()
   options, args=parser.parse_args()
   print "**************************************************\nMD Block-averaged iRED S2 Calculation"
   print "Input coordinate file:  %s"%(options.coor1)
   print "Input coordinate file:  %s"%(options.coor2)
   print "Number of MD block(s):  %s"%(options.block)
   f1 = options.coor1
   f2 = options.coor2
   try:
       bn = int(options.block)
   except ValueError:
       print "Error: Number of block should be an integer."
   ired(f1, f2, bn)
  1. Python script will output a .out file that contains the order parameter values for each amino acid residue.