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
 
(5 intermediate revisions by 2 users not shown)
Line 11: Line 11:




! /usr/bin/env python
        ! /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)


from __future__ import division
#Python script will output a .out file that contains the order parameter values for each amino acid residue.
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)
[[category:Computer_Simulation]]
 
#Python script will output a .out file that contains the order parameter values for each amino acid residue.

Latest revision as of 05:16, 20 January 2022

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.