IRed Analysis of MD Simulations - Revision history

Mjeppesen at 05:16, 20 January 2022

2022-01-20T05:16:08Z

← Older revision		Revision as of 05:16, 20 January 2022
Line 105:		Line 105:
	#Python script will output a .out file that contains the order parameter values for each amino acid residue.		#Python script will output a .out file that contains the order parameter values for each amino acid residue.

	~~[[category:Protocols]]~~
	~~[[category:Data_Processing_and_Analysis]]~~
	[[category:Computer_Simulation]]		[[category:Computer_Simulation]]

Mjeppesen: Added to Category

2021-09-11T03:50:35Z

Added to Category

← Older revision		Revision as of 03:50, 11 September 2021
Line 104:		Line 104:

	#Python script will output a .out file that contains the order parameter values for each amino acid residue.		#Python script will output a .out file that contains the order parameter values for each amino acid residue.

			[[category:Protocols]]
			[[category:Data_Processing_and_Analysis]]
			[[category:Computer_Simulation]]

Tessa17 at 16:52, 14 May 2021

2021-05-14T16:52:11Z

← Older revision		Revision as of 16:52, 14 May 2021
Line 72:		Line 72:
	eigVal = np.diag(dSort)		eigVal = np.diag(dSort)
	eigVec = v[:, Ind]		eigVec = v[:, Ind]

	# Calculate the squared order parameters		# Calculate the squared order parameters
	s2_block = np.zeros(numRes)		s2_block = np.zeros(numRes)

Tessa17 at 16:51, 14 May 2021

2021-05-14T16:51:56Z

← Older revision		Revision as of 16:51, 14 May 2021
Line 54:		Line 54:
	mMat = np.zeros([numRes, numRes])		mMat = np.zeros([numRes, numRes])
	#print numRes, numSnapshot, blockNum, blockSize		#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]

	~~for a in range(0, numRes*blockSize, numRes):~~		# Calculate the squared order parameters
	~~# N-H vectors for the current snapshot~~		s2_block = np.zeros(numRes)
	~~hnTemp = hnC[a:a+numRes, :]~~		for b in range(numRes):
	~~# normalize the vectors~~		sumOverModes = 0;
	~~for b in range(numRes):~~		for a in range(5,numRes):
	~~hnTemp[b,:] = hnTemp[b,:] / np.linalg.norm(hnTemp[b,:])~~		sumOverModes += eigVal[a,a] * (eigVec[b,a])**2
			s2_block[b] = 1 - sumOverModes
	~~# Construct Rank 2 Legendre polynomial~~		s2_sum += s2_block
	~~c = np.dot(hnTemp, hnTemp.T)~~		s2 = s2_sum / blockNum
	~~mMat = mMat + (3*np.multiply(c, c)-1)/2~~		# output
			out = np.concatenate(([resNum], s2), axis = 0)
	~~# Ensemble averaging~~		with open('ired_s2_1vec_%dblock.out' %(block), 'wb') as f:
	~~mMat = mMat / blockSize~~		np.savetxt(f, np.transpose(out), fmt='%i %.3f')
	~~#print np.shape(hnTemp), np.shape(hnTemp.T), np.shape(c), np.shape(mMat)~~		if __name__ == '__main__':
			parser = setupParserOptions()
	~~# Diagonalize matrix~~		if len(sys.argv) <2:
	~~dVals, v = np.linalg.eig(mMat)~~		parser.print_help()
	~~Ind = dVals.argsort()[::-1] # sort eigenvalues in descending order~~		sys.exit()
	~~dSort = dVals[Ind]~~		options, args=parser.parse_args()
	~~eigVal = np.diag(dSort)~~		print "**************************************************\nMD Block-averaged iRED S2 Calculation"
	~~eigVec = v[:, Ind]~~		print "Input coordinate file: %s"%(options.coor1)
			print "Input coordinate file: %s"%(options.coor2)
	# Calculate the squared order parameters		print "Number of MD block(s): %s"%(options.block)
	s2_block = np.zeros(numRes)		f1 = options.coor1
	for b in range(numRes):		f2 = options.coor2
	sumOverModes = 0;		try:
	for a in range(5,numRes):		bn = int(options.block)
	sumOverModes += eigVal[a,a] * (eigVec[b,a])**2		except ValueError:
	s2_block[b] = 1 - sumOverModes		print "Error: Number of block should be an integer."
	s2_sum += s2_block		ired(f1, f2, bn)

	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.		#Python script will output a .out file that contains the order parameter values for each amino acid residue.

Tessa17 at 16:50, 14 May 2021

2021-05-14T16:50:18Z

← Older revision		Revision as of 16:50, 14 May 2021
Line 12:		Line 12:

	! /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=OptionParser()
Line 24:		Line 22:
	parser.add_option('--block', dest='block', default = 1, help = 'Number of blocks')		parser.add_option('--block', dest='block', default = 1, help = 'Number of blocks')
	return parser		return parser
			def ired(hCoor_file, nCoor_file, block):
			try:
	def ired(hCoor_file, nCoor_file, block):		h_all = np.loadtxt(hCoor_file)
	try:		except IOError:
	h_all = np.loadtxt(hCoor_file)		print 'cannot open ', hCoor_file
	except IOError:		try:
	print 'cannot open ', hCoor_file		n_all = np.loadtxt(nCoor_file)
	try:		except IOError:
	n_all = np.loadtxt(nCoor_file)		print 'cannot open ', nCoor_file
	except IOError:		numLines = len(h_all)
	print 'cannot open ', nCoor_file		numRes = len(set(h_all[:,0]))
			numSnapshot = numLines // numRes
	numLines = len(h_all)		blockNum = int(block)
	numRes = len(set(h_all[:,0]))		blockSize = numSnapshot // blockNum
	numSnapshot = numLines // numRes		if numLines != len(n_all):
	blockNum = int(block)		print 'Error: The two coordinate files do not have same length.'
	blockSize = numSnapshot // blockNum		return
			if np.mod(numLines, numRes) != 0:
	if numLines != len(n_all):		print 'Error: Coordinate files do not formatted correctly.'
	print 'Error: The two coordinate files do not have same length.'		return
	return		# initialize s2
	if np.mod(numLines, numRes) != 0:		s2_sum = np.zeros([1, numRes])
	print 'Error: Coordinate files do not formatted correctly.'		for i in range(blockNum):
	return		hC = h_all[inumResblockSize:(i+1)numResblockSize, :]
			nC = n_all[inumResblockSize:(i+1)numResblockSize, :]
	# initialize s2		# Calculate the bond vector orientations over the trajectory
	s2_sum = np.zeros([1, numRes])		hnC = hC[:,1:] - nC[:,1:]
			# residue index
	for i in range(blockNum):		resNum = hC[0:numRes, 0]
	hC = h_all[inumResblockSize:(i+1)numResblockSize, :]		# Construct matrix m
	nC = n_all[inumResblockSize:(i+1)numResblockSize, :]		mMat = np.zeros([numRes, numRes])
			#print numRes, numSnapshot, blockNum, 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):		for a in range(0, numRes*blockSize, numRes):

Tessa17 at 16:48, 14 May 2021

2021-05-14T16:48:21Z

← Older revision		Revision as of 16:48, 14 May 2021
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=OptionParser()
	parser.add_option('--coor1', dest='coor1',default='H-NCoor.dat',help = "Input coordinate file 1")		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('--coor2', dest='coor2',default='N-HCoor.dat',help = "Input coordinate file 2")
	parser.add_option('--block', dest='block', default = 1, help = 'Number of blocks')		parser.add_option('--block', dest='block', default = 1, help = 'Number of blocks')
	return parser		return parser

Tessa17: 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..."

2021-05-14T16:46:21Z

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..."

New page

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.