diff --git a/README b/README
new file mode 100644
index 0000000..fb4ad63
--- /dev/null
+++ b/README
@@ -0,0 +1,123 @@
+After extracting both path.py and func.py to the same folder, the PATH program
+can be run using the command
+
+path.py -h
+
+This gives a list of all the flags for customizing the simulation.
+
+The program can be used in two different modes - "path" and "rock" mode. "path"
+mode calculates the most probable pathways between two equilibrium states of a
+protein. "rock" mode calculates the trajectory of a single structure along its
+nth vibrational mode. The default mode is "path" which can be changed to "rock"
+using the '-ty' flag. Further customization options with each mode is described
+below.
+
+#####################
+    "path" mode:
+#####################
+
+The most probable pathway between two equilibrium states of a protein can be
+calculated used the following command.
+
+path.py -f1 <start structure> -f2 <end structure>
+
+The two structures are PDB files. To successfully calculate the conformational
+change trajectories using the PATH program, the number of atoms in both
+the structures must be the same and so should the order of atoms in both the
+files.
+
+The program outputs the transition state (trans.pdb) and the trajectory (traj.pdb)
+connecting the two end states. Both these files are in the PDB format. It also
+outputs the energy (engy) of each structure along the trajectory relative to the
+two end states.
+
+The program also generates an output parameter file (path_log) which contains
+the RMSD between the two states of the protein, the difference in energy between
+the two minima, based on the "path" potential (described below), the time to
+the transition state from the starting structure (tbar left), time to the final
+state from the transition state (tbar right) and the energy of the transition
+state relative to the final state (Utrans).
+
+-n flag:
+
+The number of structures in the trajectory that the program outputs can be
+specified using the '-n' flag in the above command in the following way
+
+path.py -f1 <start structure> -f2 <end structure> -n 11
+
+This commands outputs a trajectory with 11 structures, which includes the two
+end states as well.
+
+-ca flag:
+
+By default, the above commands performs an all atom simulation. But, it is also
+possible to run an all atom simulation or a C-alpha only simulation using the
+'-ca' flag.
+
+path.py -f1 <start structure> -f2 <end structure> -ca 1
+
+This command performs a C-alpha only simulation.
+
+As PATH builds a Hessian matrix for both the structures which is later
+diagonalized for calculating the trajectory, the size of the protein determines
+how long it takes it calculate the trajectory. An all atom simulation of
+proteins which have less than 1000 atoms can be run on Desktop computers (8 GB
+RAM with i5 processor). For large systems, it is recommended to run PATH
+calculations on computers with larger memory. If only the overall dynamics of
+the protein is of interest then an alternative method to simulate large proteins
+is to run C-alpha only simulations. We have observed that the C-alpha only
+simulation can reproduce most of the essential dynamic information about the
+trajectory that an all atom simulation generates.
+
+
+#####################
+    "rock" mode:
+#####################
+
+This mode of the program calculates the trajectory along the nth vibrational
+mode of a single structure. It can be run using the command
+
+path.py -f1 <structure> -ty rock
+
+Apart from the '-n' and '-ca' flags that behave in the same manner as in the
+"path" mode, there are several additions flags that could be used with the
+"rock" mode.
+
+-m flag:
+
+This flag can be used to specify the nth vibrational mode along which the
+rocking trajectory is calculated. By default the first vibrational mode is used.
+
+-exag flag:
+
+Often the magnitude of displacement along the nth vibrational mode is small and
+may require an increase in the magnitude to make the displacement more
+perceptible. This flag can use used to specify a factor that is then multiplied
+to the displacement. The default value is 10.
+
+-c flag:
+
+Constraints can be applied between CA atoms of pairs of aminoacids by inputting
+a constraints file using the '-c' flag. The constraints file consists three
+columns - the first column specifies the magnitude of the constraints relative
+to the regular force constants, and the other two columns specify the aminoacid
+pairs that are constrained. For example the constraint files can consist of the
+following lines
+
+10 A37 A43
+100 A57 B74
+
+The first line constrains the CA atoms of aminoacids 37 and 43 in chain A with a
+force constant of magnitude 10 times greater than the regular force constants.
+The second line constraints CA atoms of aminoacid 57 in chain A to aminoacid 74
+in chain B with a force constant 100 times greater than the regular force
+constants.
+
+#####################
+  "path" potential:
+#####################
+
+For all atom simulations, the program uses the ANM potential energy function
+described in the Chandrasekaran et al. paper (doi: 10.1063/1.4941599). For CA
+only simulations, the program uses a mass weighted version of the empirical
+potential described in the Hinsen et al. paper (doi: 10.1016/S0301-0104(00)00222-6)
diff --git a/func.py b/func.py
new file mode 100644
index 0000000..dadf0b6
--- /dev/null
+++ b/func.py
@@ -0,0 +1,440 @@
+import numpy as np
+import Bio.PDB
+import sys
+from Bio.SVDSuperimposer import SVDSuperimposer
+from math import *
+import os
+
+###########################
+####Read a PDB file####
+###########################
+
+def pdb_read(pdb,atoms):
+	if atoms == "all":
+		sub_atm = []
+		sub_aa = []
+		sub_chain = []
+		sub_aano = []
+		sub_coord = []
+		ca_arr = []
+		j = 0
+		for i in range(len(pdb)):
+			if pdb[i][0:4] == 'ATOM' or pdb[i][0:6] == 'HETATM':
+				sub_atm.append(pdb[i][12:16])
+				sub_aa.append(pdb[i][17:20])
+				sub_chain.append(pdb[i][21:22])
+				sub_aano.append(pdb[i][23:26])
+				sub_coord.append(float(pdb[i][30:38]))
+				sub_coord.append(float(pdb[i][38:46]))
+				sub_coord.append(float(pdb[i][46:54]))
+				if (pdb[i][12:16].rstrip()).lstrip() == 'CA':
+					ca_arr.append([])
+					ca_arr[j].append(pdb[i][21:22])
+					ca_arr[j].append(pdb[i][23:26])
+					ca_arr[j].append(i)
+					j += 1
+
+		return sub_atm,sub_aa,sub_chain,sub_aano,sub_coord,ca_arr
+	elif atoms == "CA":
+		sub_atm = []
+		sub_aa = []
+		sub_chain = []
+		sub_aano = []
+		sub_coord = []
+		ca_arr = []
+		j = 0
+		for i in range(len(pdb)):
+			if pdb[i][0:4] == 'ATOM' and (pdb[i][12:16].rstrip()).lstrip() == 'CA' or pdb[i][0:6] == 'HETATM':
+				sub_atm.append(pdb[i][12:16])
+				sub_aa.append(pdb[i][17:20])
+				sub_chain.append(pdb[i][21:22])
+				sub_aano.append(pdb[i][23:26])
+				sub_coord.append(float(pdb[i][30:38]))
+				sub_coord.append(float(pdb[i][38:46]))
+				sub_coord.append(float(pdb[i][46:54]))
+				ca_arr.append([])
+				ca_arr[j].append(pdb[i][21:22])
+				ca_arr[j].append(pdb[i][23:26])
+				ca_arr[j].append(j)
+				j += 1
+
+		return sub_atm,sub_aa,sub_chain,sub_aano,sub_coord,ca_arr
+
+########################
+####Write a PDB file####
+########################
+
+def print_pdb(sub_atm,sub_aa,sub_chain,sub_aano,sub_coord,file1):
+	for i in range(len(sub_atm)):
+		file1.write("ATOM   %4d %s %s %s%4d    %8.3f%8.3f%8.3f\n" % (i+1,sub_atm[i],sub_aa[i],sub_chain[i],int(sub_aano[i]),sub_coord[3*i],sub_coord[3*i+1],sub_coord[3*i+2]))
+	file1.write("TER\n")
+	file1.write("ENDMDL\n")
+	file1.close()
+
+########################
+####Print trajectory####
+########################
+
+def print_traj_pdb(sub_atm,sub_aa,sub_chain,sub_aano,sub_coord,file1):
+	for i in range(0,len(sub_coord)):
+		for j in range(0,len(sub_coord[0])/3):
+			file1.write("ATOM   %4d %s %s %s%4d    %8.3f%8.3f%8.3f\n" % (j+1,sub_atm[j],sub_aa[j],sub_chain[j],int(sub_aano[j]),sub_coord[i][3*j],sub_coord[i][3*j+1],sub_coord[i][3*j+2]))
+		file1.write("TER\n")
+		file1.write("ENDMDL\n")
+	file1.close()
+
+###############################
+####Building Hessian Matrix####
+###############################
+
+def build_hessian(pdb,fconst,natm,ndim,cutoff):
+	hess = [[0 for x in xrange(ndim*natm)] for x in xrange(ndim*natm)]
+	for i in range(0,natm):
+		for j in range(0,natm):
+			if i == j:
+				continue
+			else:
+				distsq = (pdb[3*i]-pdb[3*j])**2 + (pdb[(3*i)+1]-pdb[(3*j)+1])**2 + (pdb[(3*i)+2]-pdb[(3*j)+2])**2
+				if distsq <= cutoff**2:
+					for k in range(0,ndim):
+						for l in range(0,ndim):
+							hess[(3*i)+k][(3*j)+l] = -(fconst*(pdb[(3*i)+k]-pdb[(3*j)+k])*(pdb[(3*i)+l]-pdb[(3*j)+l]))/distsq
+							hess[(3*i)+k][(3*i)+l] += -hess[(3*i)+k][(3*j)+l]
+	return hess
+
+#######################################
+####Hinsen Hessian Matrix with mass####
+#######################################
+
+def build_hessian_hinsen_mass(pdb,natm,ndim,aa):
+	hess = [[0 for x in xrange(ndim*natm)] for x in xrange(ndim*natm)]
+	for i in range(0,natm):
+		for j in range(0,natm):
+			if i == j:
+				continue
+			else:
+				distsq = (pdb[3*i]-pdb[3*j])**2 + (pdb[(3*i)+1]-pdb[(3*j)+1])**2 + (pdb[(3*i)+2]-pdb[(3*j)+2])**2
+				if distsq < 16:
+					fconst = 8.6e2*np.sqrt(distsq)-2.39e3
+				else:
+					fconst = float(128e4)/float(distsq**3)
+				for k in range(0,ndim):
+					for l in range(0,ndim):
+						hess[(3*i)+k][(3*j)+l] = float(-(fconst*(pdb[(3*i)+k]-pdb[(3*j)+k])*(pdb[(3*i)+l]-pdb[(3*j)+l]))/distsq)/float(np.sqrt(aa_mass(aa[i]))*np.sqrt(aa_mass(aa[j])))
+						hess[(3*i)+k][(3*i)+l] += -hess[(3*i)+k][(3*j)+l]
+	return hess
+
+#############################
+####Constaints processing####
+#############################
+
+def constr_process(constr,ca):
+	arr = []
+	for i in range(0,len(constr)):
+		temp = (constr[i].rstrip()).split()
+		arr.append([])
+		arr[i].append(float(temp[0]))
+		for j in range(0,len(ca)):
+			if temp[1][0:1] == ca[j][0] and temp[1][1:] == ca[j][1]:
+				arr[i].append(ca[j][2])
+				continue
+			if temp[2][0:1] == ca[j][0] and temp[2][1:] == ca[j][1]:
+				arr[i].append(ca[j][2])
+				break
+	return arr
+
+##########################
+####Adding Constraints####
+##########################
+
+def add_constr(hess,constr,ca,pdb,opt,cutoff,aa):
+	if constr == None:
+		return hess
+	else:
+		arr = constr_process(constr,ca)
+		print arr
+		for idx in range(0,len(arr)):
+			mat = [[ 0 for x in xrange(3)] for x in xrange(3)]
+			mag = arr[idx][0]
+			i = arr[idx][1]
+			j = arr[idx][2]
+			distsq = (pdb[3*i]-pdb[3*j])**2 + (pdb[(3*i)+1]-pdb[(3*j)+1])**2 + (pdb[(3*i)+2]-pdb[(3*j)+2])**2
+			if opt  == 'CA':
+				if distsq < 16:
+					fconst = mag*(8.6e2*np.sqrt(distsq)-2.39e3)
+				else:
+					fconst = mag*(float(128e4)/float(distsq**3))
+				for k in range(0,3):
+					for l in range(0,3):
+						hess[(3*i)+k][(3*j)+l] += float(-(fconst*(pdb[(3*i)+k]-pdb[(3*j)+k])*(pdb[(3*i)+l]-pdb[(3*j)+l]))/distsq)/float(np.sqrt(aa_mass(aa[i]))*np.sqrt(aa_mass(aa[j])))
+						hess[(3*j)+k][(3*i)+l] += hess[(3*i)+k][(3*j)+l]
+						hess[(3*i)+k][(3*i)+l] += -hess[(3*i)+k][(3*j)+l]
+						hess[(3*j)+k][(3*j)+l] += -hess[(3*i)+k][(3*j)+l]
+			elif opt == 'all':
+				fconst = mag*0.01
+				for k in range(0,3):
+					for l in range(0,3):
+						hess[(3*i)+k][(3*j)+l] += -(fconst*(pdb[(3*i)+k]-pdb[(3*j)+k])*(pdb[(3*i)+l]-pdb[(3*j)+l]))/distsq
+						hess[(3*j)+k][(3*k)+l] += hess[(3*i)+k][(3*j)+l]
+						hess[(3*i)+k][(3*i)+l] += -hess[(3*i)+k][(3*j)+l]
+						hess[(3*j)+k][(3*j)+l] += -hess[(3*i)+k][(3*j)+l]
+		return hess
+
+##########################
+####Aligning molecules####
+##########################
+
+def align(mol1,mol2,natm,ndim):
+	ref_atoms = []
+	alt_atoms = []
+	moving2 = []
+
+	for i in range(0,ndim*natm):
+		ref_atoms.append(mol1[i])
+		alt_atoms.append(mol2[i])
+
+	fixed = np.reshape(ref_atoms,(natm,3))
+	moving = np.reshape(alt_atoms,(natm,3))
+	sup=SVDSuperimposer()
+	sup.set(fixed,moving)
+	sup.run()
+	rot,tran = sup.get_rotran()
+	rot=rot.astype('f')
+	tran=tran.astype('f')
+	for j in range(0,natm):
+		moving2.append(np.dot(moving[j],rot)+tran)
+
+	j=0
+
+	mol2_a = [0 for x in xrange(0,natm*ndim)]
+
+	for i in range(0,len(moving2)):
+		mol2_a[j]=moving2[i][0]
+		mol2_a[j+1]=moving2[i][1]
+		mol2_a[j+2]=moving2[i][2]
+		j = j + 3
+
+	return ref_atoms,mol2_a,sup.get_rms()
+
+##################
+####Read files####
+##################
+
+def readfile(fname):
+	fin = open(fname,"r")
+	arr = fin.readlines()
+	fin.close()
+	return arr
+
+###############################
+####Work Matrix Calculation####
+###############################
+
+def work(array1,array2):
+	sub_work = []
+	for i in range(0,len(array1)):
+		sub_work.append(array1[i]-array2[i])
+	return sub_work
+
+##########################
+####Energy Calculation####
+##########################
+
+def ener_calc(hess,work):
+	ener = 0.5*np.dot(work.T,np.dot(np.asarray(hess),work))
+	return ener
+
+####################################
+####Velocity continuity equation####
+####################################
+
+def calc_xbar(eva1,eva2,evec1,evec2,mol1_d,mol2_d,tl,tr,tf,natm,ndim):
+	bl = [[0 for x in xrange(natm*ndim)] for x in xrange(natm*ndim)]
+	ar = [[0 for x in xrange(natm*ndim)] for x in xrange(natm*ndim)]
+
+	for i in range(0,len(eva1)):
+		if eva1[i]<0.000001:
+			bl[i][i] = 1/tl
+			eva2[i]=0
+		else:
+			if eva1[i]*tl > 707:
+				bl[i][i] = eva1[i]
+			else:
+				bl[i][i] = (eva1[i]*cosh(eva1[i]*tl))/(sinh(eva1[i]*tl))
+		if eva2[i]<0.000001:
+			ar[i][i] = 1/(tl-tf)
+			eva2[i]=0
+		else:
+			if eva2[i]*(tf-tl) > 707:
+				ar[i][i] = -eva2[i]
+			else:
+				ar[i][i] = (eva2[i]*cosh(eva2[i]*(tl-tf)))/(sinh(eva2[i]*(tl-tf)))
+
+	den1 = np.dot(evec1,np.dot(bl,evec1.T))
+	den2 = np.dot(evec2,np.dot(ar,evec2.T))
+
+	num1 = np.dot(den1,mol1_d)
+	num2 = np.dot(den2,mol2_d)
+
+	num = num1 - num2
+	den = den1 - den2
+
+	xbar = np.dot(num,np.linalg.inv(den))
+
+	return xbar
+
+########################
+####Calculating tbar####
+########################
+
+def calc_tbar(eval):
+	if len(eval) == 6:
+		skip = 5
+	else:
+		skip = 6
+
+	reci_eval = []
+
+	for i in range(skip,len(eval)):
+		reci_eval.append(float(1)/float(eval[i]))
+
+	sum_eval = sum(reci_eval)
+
+	incr = 0
+	i = 0
+
+	while incr <= 0.95*sum_eval:
+		incr += reci_eval[i]
+		i += 1
+
+	if skip != 5:
+		incr += reci_eval[i]
+		i += 1
+
+	avg_k = 1/(float(incr)/float(i))
+
+	return (float(7)/float(avg_k)),avg_k
+
+#################################
+####Calculate path trajectory####
+#################################
+
+def path_traj(nconf,hess_mol1,hess_mol2,mol1a,mol2a,eva1,eva2,evec1,evec2,kl,kr,tl,tr,tf,workl,workr,left_e,right_e,fengy,natm,ndim):
+	dt = tf/(nconf-1)
+	t = 0.0
+
+	xmat = []
+	xlr = [[0 for x in xrange(natm*ndim)] for x in xrange(natm*ndim)]
+
+	t_arr = []
+	t_arr.append(0)
+
+	h = (nconf-1)/2
+	de = float(1)/float(h)
+
+	for i in range(0,h):
+		t = float(7+np.log((i+1)*de))/float(kl)
+		t_arr.append(t)
+
+	for i in range(0,h-1):
+		t = float(7+np.log(1-((i+1)*de)))/float(kr)
+		t_arr.append(tf-t)
+
+	t_arr.append(tf)
+
+	for i in range(0,nconf):
+		xmat.append([])
+		for j in range(0,len(eva1)):
+			if t_arr[i] <= tl:
+				if eva1[j]<0.000001:
+					xlr[j][j] = t_arr[i]/tl
+				else:
+					if eva1[j]*tl > 707:
+						xlr[j][j] = 1
+					else:
+						xlr[j][j] = (sinh(eva1[j]*t_arr[i]))/(sinh(eva1[j]*tl))
+			elif t_arr[i] > tl:
+				if eva2[j]<0.000001:
+					xlr[j][j] = (tf-t_arr[i])/(tf-tl)
+				else:
+					if eva2[j]*(tf-tl) > 707:
+						xlr[j][j] = -1
+					else:
+						xlr[j][j] = -(sinh(eva2[j]*(t_arr[i]-tf)))/(sinh(eva2[j]*(tf-tl)))
+		if t_arr[i] <= tl:
+			xt = np.dot(np.dot(evec1,np.dot(xlr,evec1.T)),workl) + mol1a
+			fengy.write("%2.3f\t%+2.3f\n" % (t_arr[i],ener_calc(hess_mol1,np.asarray(work(xt,mol1a))) + float(right_e-left_e)))
+		elif t_arr[i] > tl:
+			xt = np.dot(np.dot(evec2,np.dot(xlr,evec2.T)),workr) + mol2a
+			fengy.write("%2.3f\t%+2.3f\n" % (t_arr[i],ener_calc(hess_mol2,np.asarray(work(xt,mol2a)))))
+		for j in range(len(xt)):
+			xmat[i].append(xt[j])
+
+	fengy.close()
+
+	return xmat
+
+####################################
+####Calculate rocking trajectory####
+####################################
+
+def rock_traj(nconf,hess_mol,mol,eva,evec,natm,ndim,mode,exag):
+	if len(eva) == 6:
+		if mode == 1:
+			m = 5
+		else:
+			print 'For a two atom system there is only 1 vibrational mode'
+			sys.exit(1)
+	else:
+		m = 5 + mode
+
+	kval = eva[m]
+
+	tf = float(1)/float(kval)
+
+	t = 0.0
+	dt = float(tf)/float(nconf-1)
+
+	xmat = []
+
+	for i in range(0,nconf):
+		xmat.append([])
+		temp = [0 for x in xrange(len(mol))]
+		for j in range(len(eva)):
+			for k in range(m,m+1):
+				temp[j] += (exag*evec[j][k]*sin(2*pi*eva[k]*t))
+		for j in range(len(temp)):
+			xmat[i].append(temp[j]+mol[j])
+		t = t + dt
+
+	return xmat
+
+#######################
+#####Aminoacid Mass####
+#######################
+
+def aa_mass(aa):
+	mass_aa = {'ALA':71.0788,
+				'CYS':103.1388,
+				'ASP':115.0886,
+				'GLU':129.1155,
+				'PHE':147.1766,
+				'GLY':57.0519,
+				'HIS':137.1411,
+				'ILE':113.1594,
+				'LYS':128.1741,
+				'LEU':113.1594,
+				'MET':131.1926,
+				'ASN':114.1038,
+				'PRO':97.1167,
+				'GLN':128.1307,
+				'ARG':156.1875,
+				'SER':87.0782,
+				'THR':101.1051,
+				'VAL':99.1326,
+				'TRP':186.2132,
+				'TYR':163.1760,
+				'TRX':186.2132}
+	for mass in mass_aa:
+		if aa == mass:
+			return mass_aa[mass]
diff --git a/path.py b/path.py
new file mode 100644
index 0000000..6fecf65
--- /dev/null
+++ b/path.py
@@ -0,0 +1,268 @@
+#!/usr/bin/env python
+
+"""PATH v3.0
+Author: Srinivas Niranj Chandrasekaran
+Affliation: Carter Lab, UNC-Chapel Hill; Bailey Lab, UMass-Worcester
+Email: Srinivas.Chandrasekaran@umassmed.edu
+Date: August 23rd 2016"""
+
+from math import *
+import sys
+import os.path
+import numpy as np
+import time
+import argparse
+import func
+import scipy.linalg as sp
+
+start_time = time.time()
+
+parser = argparse.ArgumentParser(description = 'Computes the most probable pathway connecting two equilibrium states of a protein or rock a structure along the normal modes',\
+	epilog='For more information check the README file')
+parser.add_argument('-ty',metavar='',type=str,default="path",help='use "rock" for rocking the structures [default:"path"]')
+parser.add_argument('-m',metavar='',type=int,default=1,help='The nth mode along which the structure should be rocked [default:1] [optional]')
+parser.add_argument('-f1',metavar='',required=True,help='Initial pdb [required]')
+parser.add_argument('-f2',metavar='',help='Final pdb [required only for "path"]')
+parser.add_argument('-n',metavar='',type=int,default='3',help='Number of conformations [default:3] [optional]')
+parser.add_argument('-ca',metavar='',type = int,default=0,help='1 for CA only [default:0 all atom] [optional]')
+parser.add_argument('-c',metavar='',default=None,help='Constraints between CA atoms of aminoacids [can be used only with "rock"] [optional]')
+parser.add_argument('-exag',metavar='',default=10,type=float,help='The motion is exaggerated such that the displacement from equilibrium is perceptible [default:10] [optional]')
+
+args = parser.parse_args()
+
+#############################
+####Functions and Classes####
+#############################
+
+class Input(object):
+	def __init__(self,args):
+		self.pdb1 = func.readfile(args.f1)
+		self.pdb1name = args.f1
+		self.mode = args.m
+		self.type = args.ty
+		self.nconf = args.n
+		self.exag = args.exag
+
+		if self.type == 'path':
+			if args.f2:
+				self.pdb2 = func.readfile(args.f2)
+				self.pdb2name = args.f2
+			else:
+				print '\nThe second structure id required for calculating the most probable path'
+				sys.exit(1)
+
+		if args.ca == 0:
+			self.opt = 'all'
+			self.space = ''
+		elif args.ca == 1:
+			self.opt = 'CA'
+			self.space = ' CA'
+
+		if args.c != None and self.type == 'rock':
+			self.c = func.readfile(args.c)
+		else:
+			self.c = None
+
+def open_file(opt):
+	if opt == 'log':
+		fname = 'path_log'
+	if opt == 'trans':
+		fname = 'trans.pdb'
+	if opt == 'traj':
+		fname = 'traj.pdb'
+	elif opt == 'engy':
+		fname  = 'engy'
+
+	if os.path.isfile(fname):
+		os.remove(fname)
+
+	out = open(fname,'w')
+	return out
+
+class Const(object):
+	def __init__(self):
+		self.dim = 3
+		self.cutoff = 8.0
+		self.k = 0.01
+
+def len_check(arr1,arr2,space):
+	if len(arr1) != len(arr2):
+		print "\n@> The two structures don't have the same number of atoms. Please check!\n"
+		sys.exit(1)
+	else:
+		natm = len(arr1)/3
+
+	if natm > 1:
+		print "\n@> There are %d%s atoms in your molecule\n" % (natm,space)
+	else:
+		print "\n@> There is %d atom in your molecule. PATH requires at least two atoms\n" % natm
+		sys.exit(1)
+
+	return natm
+
+#####################
+####Main Function####
+#####################
+
+if __name__ == '__main__':
+
+	inp = Input(args)
+
+	ftraj = open_file('traj')
+
+	const = Const()
+
+	###############
+	####Rocking####
+	###############
+
+	if inp.type == 'rock':
+
+		###########################
+		####Reading Coordinates####
+		###########################
+
+		atm,aa,chain,aano,mol,ca = func.pdb_read(inp.pdb1,inp.opt)
+
+		natm = len(atm)
+
+		#############################
+		####Building Hessian matrices
+		#############################
+
+		if inp.opt == 'all':
+			hess = func.build_hessian(mol,const.k,natm,const.dim,const.cutoff)
+		elif inp.opt == 'CA':
+			hess = func.build_hessian_hinsen_mass(mol,natm,const.dim,aa)
+
+		##########################
+		####Adding Constraints####
+		##########################
+
+		hess_mol = func.add_constr(hess,inp.c,ca,mol,inp.opt,const.cutoff,aa)
+
+		print "\n@> %2.3fs - The Hessian matrix has been built.\n" %(time.time()-start_time)
+
+		####################################
+		####Eigenvalues and Eigenvectors####
+		####################################
+
+		eva, evec = sp.eigh(hess_mol,eigvals=(0,(natm*const.dim)-1))
+
+		print "@> %2.3fs - The Eigenvalues and the Eigenvectors have been calculated.\n" % (time.time()-start_time)
+
+		##############################
+		####Calculating trajectory####
+		##############################
+
+		xmat = func.rock_traj(inp.nconf,hess_mol,mol,eva,evec,natm,const.dim,inp.mode,inp.exag)
+
+		func.print_traj_pdb(atm,aa,chain,aano,xmat,ftraj)
+
+		print "@> %2.3fs - The trajectory has been computed.\n" % (time.time()-start_time)
+
+	#############################
+	####Most probable pathway####
+	#############################
+
+	elif inp.type == 'path':
+
+		flog = open_file('log')
+
+		###########################
+		####Reading Coordinates####
+		###########################
+
+		atm1,aa1,chain1,aano1,mol1,ca1 = func.pdb_read(inp.pdb1,inp.opt)
+		atm2,aa2,chain2,aano2,mol2,ca2 = func.pdb_read(inp.pdb2,inp.opt)
+
+		natm = len_check(mol1,mol2,inp.space)
+
+		mol1a,mol2a,rmsd = func.align(mol1,mol2,natm,const.dim)
+
+		flog.write('Initial State = %s\n\n' % inp.pdb1name)
+		flog.write('Final State = %s\n\n' % inp.pdb2name)
+
+		flog.write("The rmsd between the two structures is %f Angstroms\n\n" % rmsd)
+
+		#############################
+		####Building Hessian matrices
+		#############################
+
+		if inp.opt == 'all':
+			hess1 = func.build_hessian(mol1a,const.k,natm,const.dim,const.cutoff)
+			hess2 = func.build_hessian(mol2a,const.k,natm,const.dim,const.cutoff)
+		elif inp.opt == 'CA':
+			hess1 = func.build_hessian_hinsen_mass(mol1a,natm,const.dim,aa1)
+			hess2 = func.build_hessian_hinsen_mass(mol2a,natm,const.dim,aa2)
+
+		##########################
+		####Adding Constraints####
+		##########################
+
+		hess_mol1 = func.add_constr(hess1,inp.c,ca1,mol1a,inp.opt,const.cutoff,aa1)
+		hess_mol2 = func.add_constr(hess2,inp.c,ca2,mol2a,inp.opt,const.cutoff,aa2)
+
+		print "@> %2.3fs - The Hessian matrices have been built.\n" %(time.time()-start_time)
+
+		####################################
+		####Eigenvalues and Eigenvectors####
+		####################################
+
+		eva1, evec1 = sp.eigh(hess_mol1,eigvals=(0,(natm*const.dim)-1))
+		eva2, evec2 = sp.eigh(hess_mol2,eigvals=(0,(natm*const.dim)-1))
+
+		print "@> %2.3fs - The Eigenvalues and the Eigenvectors have been calculated.\n" % (time.time()-start_time)
+
+		########################
+		####Calculating tbar####
+		########################
+
+		tl,kl = func.calc_tbar(eva1)
+		tr,kr = func.calc_tbar(eva2)
+
+		tf = tl + tr
+
+		########################
+		####Identifying xbar####
+		########################
+
+		xbar = func.calc_xbar(eva1,eva2,evec1,evec2,mol1a,mol2a,tl,tr,tf,natm,const.dim)
+
+		ftrans = open_file('trans')
+		func.print_pdb(atm1,aa1,chain1,aano1,xbar,ftrans)
+
+		workl = np.asarray(func.work(xbar,mol1a))
+		workr = np.asarray(func.work(xbar,mol2a))
+
+		###########################################
+		####Calculating transition state energy####
+		###########################################
+
+		left_e = func.ener_calc(hess_mol1,workl)
+		right_e = func.ener_calc(hess_mol2,workr)
+
+		flog.write("The difference in energy between the two wells is : %2.3f\n" % (float(right_e-left_e)))
+		flog.write("\ntbar left  = %.3e\n" % tl)
+		flog.write("\ntbar right = %.3e\n" % tr)
+		flog.write("\nUtrans  = %2.3f\n" % right_e)
+
+		flog.close()
+
+		##############################
+		####Calculating trajectory####
+		##############################
+
+		fengy = open_file('engy')
+
+		xmat = func.path_traj(inp.nconf,hess_mol1,hess_mol2,mol1a,mol2a,eva1,eva2,evec1,evec2,kl,kr,tl,tr,tf,workl,workr,left_e,right_e,fengy,natm,const.dim)
+
+		func.print_traj_pdb(atm1,aa1,chain1,aano1,xmat,ftraj)
+
+		print "@> %2.3fs - The trajectory has been computed.\n" % (time.time()-start_time)
+
+	#################################
+	####Time taken by the program####
+	#################################
+
+	print "Total time taken %2.3fs\n" % (time.time()-start_time)