#Load topologyy files

parm ../../../01_6DDE_BU72/rep1/01_6DDE_BU72_nowater.prmtop [BU72]
parm ../02_6DDE_BU72_PAM_nowater.prmtop [BU72_PAM]

#Load trajectories, focus on the first 100 ns (with 100,000 frames)

trajin ../../../01_6DDE_BU72/rep1/01_6DDE_BU72_nowater.dcd 1 100000 parm [BU72]
trajin ../02_6DDE_BU72_PAM_nowater.dcd 1 100000 parm [BU72_PAM]

#fit with muOR residues, exclude non-muOR residues
#I will use all non-H atoms
rms first :972-1252&!@H=

#create average structure
average crdset muOR-average

createcrd muOR-trajectories

run 

crdaction muOR-trajectories rms ref muOR-average :972-1252&!@H=

#Calculate coordinate covariance matrix
crdaction muOR-trajectories matrix covar name muOR-covar :972-1252&!@H=
##
###Diagonalize coordinate covariance matrix
###Get first 3 eigenvectors
##
runanalysis diagmatrix muOR-covar out muOR-evecs.dat vecs 3 name myEvecs nmwiz nmwizvecs 3 nmwizfile muOR.nmd nmwizmask :972-1252&!@H=
##
#Create separate projections for each set of trajectories
crdaction muOR-trajectories projection bu72 modes myEvecs beg 1 end 3 :972-1252&!@H= crdframes 1,100000
crdaction muOR-trajectories projection bu72_pam modes myEvecs beg 1 end 3 :972-1252&!@H= crdframes 100001,last
##
###Create normalized histogram of the 3 calculated projections
##
hist bu72:1 bins 100 out bu72-hist.agr norm name bu72-1
hist bu72:2 bins 100 out bu72-hist.agr norm name bu72-2
hist bu72:3 bins 100 out bu72-hist.agr norm name bu72-3
##
hist bu72_pam:1 bins 100 out bu72-pam-hist.agr norm name bu72-pam-1
hist bu72_pam:2 bins 100 out bu72-pam-hist.agr norm name bu72-pam-2
hist bu72_pam:3 bins 100 out bu72-pam-hist.agr norm name bu72-pam-3

run

clear all

readdata muOR-evecs.dat name Evecs

parm ../../../01_6DDE_BU72/rep1/01_6DDE_BU72_nowater.prmtop
parmstrip !(:972-1252&!@H=)
parmwrite out muOR-modes.prmtop

#Create NetCDF trajectory file with PC 1
runanalysis mode name Evecs trajout muOR-mode1.nc pcmin -100 pcmax 100 tmode 1 trajout :972-1252&!@H= trajout netcdf