IRed Analysis of MD Simulations: Difference between revisions
From Powers Wiki
(Created page with "Calculating order parameters from MD Simulations The analysis utilizes iRed software provided by the Dr. Rafael Bruschweiler Group. All iRed calculations were completed on th...") |
No edit summary |
||
Line 11: | Line 11: | ||
! /usr/bin/env python | ! /usr/bin/env python | ||
from __future__ import division | from __future__ import division | ||
import numpy as np | import numpy as np | ||
from optparse import OptionParser | from optparse import OptionParser | ||
import sys | import sys | ||
def setupParserOptions(): | 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 | |||
Revision as of 16:48, 14 May 2021
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.
- Create a pdb file of the trajectory in the same manner as the .xtc file.
- 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
- awk ‘{if (($3 == “N”) && ($4 != “PRO”)) {print $6, $7, $8, $9:}} xxx.pdb > N-HCoor.dat
- awk ‘{if (($3 == “H”) && ($4 != “H1”)) {print $6, $7, $8, $9:}} xxx.pdb > H-NCoor.dat
- Once .dat files are created, use the python script below to run the iRed calculation
- 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)
- Python script will output a .out file that contains the order parameter values for each amino acid residue.