# system: 7mer peptide in a box of
# periodic water (1577 Waters)


# Read in the trajectory file 
# starting at 1
# 100 total snapshots
# consider each one 
#

trajin trj.5 1 100 1


#  first calculate some things that don't depend on
#  the position of the solute in the box

# do dihedral angles

dihedral phi_1 :1@C  :2@N  :2@CA :2@C out phi_1
dihedral psi_1   :1@N  :1@CA :1@C  :2@N out psi_1
dihedral omega_1   :1@CA :1@C  :2@N  :2@CA out omega_1

dihedral phi_2 :2@C  :3@N  :3@CA :3@C out phi_2 
dihedral psi_2   :2@N  :2@CA :2@C  :3@N out psi_2
dihedral omega_2   :2@CA :2@C  :3@N  :3@CA out omega_2

dihedral phi_3 :3@C  :4@N  :4@CA :4@C out phi_3 
dihedral psi_3   :3@N  :3@CA :3@C  :4@N out psi_3
dihedral omega_3   :3@CA :3@C  :4@N  :4@CA out omega_3

dihedral phi_4 :4@C  :5@N  :5@CA :5@C out phi_4
dihedral psi_4   :4@N  :4@CA :4@C  :5@N out psi_4
dihedral omega_4   :4@CA :4@C  :5@N  :5@CA out omega_4

dihedral phi_5 :5@C  :6@N  :6@CA :6@C out phi_5
dihedral psi_5   :5@N  :5@CA :5@C  :6@N out psi_5
dihedral omega_5   :5@CA :5@C  :6@N  :6@CA out omega_5

dihedral phi_6 :6@C  :7@N  :7@CA :7@C out phi_6
dihedral psi_6   :6@N  :6@CA :6@C  :7@N out psi_6
dihedral omega_6   :6@CA :6@C  :7@N  :7@CA out omega_6

# calculate the distance between two atoms 
# over the 10 snapshots
# output to the file "dist.list"

distance end_to_end :1@N :7@N out  dist_end_to_end.list

##calculate an "average" structure for the 10 
## snapshots and output to file "avgb4.pdb"

average avgb4.pdb pdb 


#get the rms values for the whole system 
#referenced to the first snapshot

rms first :1-1584


#now deal with drift

# the peptide (residues 1->7) has drifted around 
# the box during MD.  We'd like to view the system
# with a "constant" solute.

# note that this changes the coordinates you have in memory

# first center the solute by itself

center origin :1-7 

# now image the whole system about the centered origin

image origin center 

# now calculate a pseudo occupancy for the Waters about
# the solute

# in this case 9*3.0 = boxx
#             13*3.0 = boxy
#              9*3.0 = boxz
# where the integer is the number of bins and the real
# is the size of the bin (in Angstroms).
#
# max is the "minimum" occupancy "time" the site
# must be occupied to be output to grid.list
# (which is  pseudo PDB file in the "new" reference frame)

grid grid.list 9 3.0 13 3.0 9 3.0 :WAT max 0.8  

# get an average pdbfile from the "new" reference
# frame to use to with grid.lis from the last step
# to display 
 
average avg.pdb pdb 

# save the final "new" snapshot as a PDB file
# for alternate veiwing

trajout test.pdb PDB 

# now do some H-Bond analysis


# specify the electron pair DONOR
#
donor mask :1-7@O

# specify the ACCEPTOR(s)

acceptor WAT O H1 
acceptor WAT O H2 

# do the Hbond search/output

hbond solventacceptor O H1 solventacceptor O H2 solventneighbor 2 series hbond

# calculate the waters in the first and second solvation shells
# (0-3.5A and 3.5-5.0A) and output to watershell.list 
 
watershell :1-7 watershell.list  
 
# save only the peptide for easy graphics viewing
 
strip :WAT
 
trajout trj.peptide nobox