remarks Generate dihedral series in a form suitable
remarks for a dials plot with DialX.
remarks this file: DialX.inp
remarks Thomas Simonson, Yale University, Sept. 1991

set message=off end
set echo=off end

!---------------------------
! user files and parameters
!---------------------------
evaluate ($coor="il8_template.pdb")
evaluate ($struc="dials.psf")
evaluate ($traj="tad.crd")
evaluate ($begin=1)
evaluate ($stop=400)
evaluate ($skip=1)
evaluate ($outp="dials.data")


!--------------------------------------------------
! steer the relevant output to a display/print file
!--------------------------------------------------
!set display=$outp end
!set print=$outp end


!-----------------------------
! patches for making dihedrals
!-----------------------------
topology @patches.pro end


!---------------------------------------------------------
! read existing structure file, edited to remove dihedrals
!---------------------------------------------------------
structure @$struc end


!----------------------------------------------------------
! read coordinates to check connectivity for phi,psi setup
!----------------------------------------------------------
coordinates @$coor
coor disp=comp @$coor
!vector do (segid="1") (segid PROT)
!vector do (segid="M13M") (segid 1)


!----------------------------------------------------------------
! now generate the dihedrals in the proper order for a dials plot
! If a subset of the structure is to be plotted, change the
! atom selection for the loop.
! NB: don't use vector stores to define the atom selection!!!
! they will get erased!!
!----------------------------------------------------------------
for $i in id (name CA) loop li
vector show (segid) (id $i)
evaluate ($seg=$result)
vector show (resnam) (id $i)
evaluate ($name=$result)
vector show (resid) (id $i)
evaluate ($iden=$result)
! phi angle: look for residue with its C close to current N
patch PPHI 
   refe=1=(byres (name C and ((atom $seg $iden N) around 1.7)) ) 
   refe=2=(byres id $i)
end
! psi,omega angles: look for residue with its N close to current C
patch PPSI 
   refe=2=(byres id $i)
   refe=3=(byres (name N and ((atom $seg $iden C) around 1.7)) ) 
end
! chi angles: use special patches
evaluate ($rname="P"+$name)
patch $rname refe=nil=(byres id $i) end 
end loop li


!--------------------------------------
! write structure file for future plots
!--------------------------------------
!write structure output=dials.psf end


!------------------
! dummy parameters
!------------------
parameters 
  BOND  (all) (all) 0. 0.
  ANGL  (all) (all) (all) 0. 0.
  DIHE  (all) (all) (all) (all)  0. 3  0.
  IMPR  (all) (all) (all) (all)  0. 3  0.      
end


!-----------------------------------------------------------
! display list of amino acid types with number of chi angles
!-----------------------------------------------------------
display amino acid ALA 0
display amino acid ARG 5
display amino acid ASN 2
display amino acid ASP 2
display amino acid CYS 1
display amino acid GLN 3
display amino acid GLU 3
display amino acid GLY 0
display amino acid HIP 3
display amino acid HIS 3
display amino acid ILE 2
display amino acid LEU 2
display amino acid LYS 4
display amino acid MET 3
display amino acid PHE 4
display amino acid PRO 2
display amino acid SER 1
display amino acid THR 1
display amino acid TRP 2
display amino acid TYR 4
display amino acid VAL 1
display end amino acid list


!--------------------------------------------------------------------
! display backbone dihedral values; set unknown dihedrals to 9999.
! If a subset of the structure is to be plotted, change the
! atom selection for the loop.
! NB: don't use vector stores to define the atom selection!!!
! they will get erased!!
!--------------------------------------------------------------------
for $i in id (name CA) loop li
vector show (segid) (id $i)
evaluate ($seg=$result)
vector show (resnam) (id $i)
evaluate ($name=$result)
vector show (resid) (id $i)
evaluate ($iden=$result)
vector show (x) (id $i)
! phi
evaluate ($result=9999.)
pick dihedral
   (name C and ((atom $seg $iden N) around 1.7) ) 
   (atom $seg $iden N)
   (id $i)
   (atom $seg $iden C)
   geom
evaluate ($phi=$result)
! psi
evaluate ($result=9999.)
pick dihedral
   (atom $seg $iden N)
   (id $i)
   (atom $seg $iden C)
   (name N and ((atom $seg $iden C) around 1.7) ) 
   geom
evaluate ($psi=$result)
! omega
evaluate ($result=9999.)
pick dihedral
   (id $i)
   (atom $seg $iden C)
   (name N and ((atom $seg $iden C) around 1.7) ) 
   (name CA and segid $seg
            and byres (name N and ((atom $seg $iden C) around 1.7)) ) 
   geom
evaluate ($omega=$result)
! display results
display residue list:: $name $iden $phi $psi $omega
end loop li

!--------------------------------------------
! number of static structures to be displayed
!--------------------------------------------
display 2 static structures

display selected angles PHI PSI OMEGA CHI1 CHI2

!---------------- 
! print dihedrals
!---------------- 
vector do (segid="") (all)   ! shortens output a bit
vector do (x=999.) (not known)
vector do (y=999.) (not known)
vector do (z=999.) (not known)
print dihedrals


!-------------------------
! second static structure
coor swap end
vector do (x=999.) (not known)
vector do (y=999.) (not known)
vector do (z=999.) (not known)
print dihedrals


!----------------------------------------
! read trajectory and print out dihedrals
!----------------------------------------
evaluate ($nframe=($stop-$begin+$skip)/$skip +0.5)
read dynamics
  ascii=false
  input=$traj
  begin=$begin skip=$skip stop=$stop
end
print dihedrals
evaluate ($i=2)
while ($i < $nframe) loop li
read dynamics next end 
print dihedrals
evaluate ($i=$i+1)
end loop li


stop

.