{+ file: delta_phipsi.inp +}
{+ directory: general +}
{+ description: Calculate Phi/Psi difference between 2 structures +}
{+ comment: 
            calculates overall phi/psi difference between the structure
            and lists phi/psi difference as a function of residue number +}
{+ authors: Axel T. Brunger, and Paul D. Adams +}
{+ copyright: Yale University +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) are not quoted
   - do not remove any evaluate statements from the file -}

{- begin block parameter definition -} define(

{======================= molecular structure =========================}

{* structure file *}
{===>} structure_infile="amy.mtf";

{* reference coordinate file *}
{===>} comp_coordinate_infile="amy.pdb";

{* coordinate file *}
{===>} main_coordinate_infile="amy_anneal.pdb";

{========================== atom selection ===========================}

{* select atoms to be included in calculation *}
{===>} atom_select=(known and not hydrogen);

{=========================== output files ============================}

{* output listing file *}
{===>} list_outfile="delta_phipsi.list";

{===========================================================================}
{        things below this line do not normally need to be changed          }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates disp=comp @&comp_coordinate_infile

 coordinates @&main_coordinate_infile

 set display=&list_outfile end

 display >>>> Phi/psi difference calculated between 2 structures:
 display >>>>   structure file= &STRIP%structure_infile
 display >>>>   main coordinate set= &STRIP%main_coordinate_infile
 display >>>>   reference coordinate set= &STRIP%comp_coordinate_infile
 display
 display >>>> the difference is defined as:
 display >>>>   difference = sqrt((phi_1 - phi_2)^2 + (psi_1 - psi_2)^2)
 display
 display  segid  resid   delta phi/psi
 display

 set display=OUTPUT end

 evaluate ($overall_phipsi=0)
 evaluate ($nphipsi=0)

 for $id in id ( &atom_select and name ca ) loop dram

   show (segid) (id $id)
   evaluate ($segid=$result)
   show (resid) (id $id)
   evaluate ($resid=$result)

   evaluate ($phi1=-9999)
   evaluate ($psi1=-9999)
   evaluate ($phi2=-9999)
   evaluate ($psi2=-9999)

   identity (store1) (name c and bondedto
            (name n and resid $resid and segid $segid ))

   if ( $select = 1 ) then

     show element (store1) (attribute store1 > 0)
     evaluate ($id_prev=$result)

     show (segid) (id $id_prev)
     evaluate ($segid_prev=$result)
     show (resid) (id $id_prev)
     evaluate ($resid_prev=$result)

     pick dihedral
       (name  c and segid $segid_prev and resid $resid_prev)
       (name  n and segid $segid and resid $resid)
       (name ca and segid $segid and resid $resid)
       (name  c and segid $segid and resid $resid)
       geometry

     evaluate ($phi1=mod($result+360,360))

     coord swap end

     pick dihedral
       (name  c and segid $segid_prev and resid $resid_prev)
       (name  n and segid $segid and resid $resid)
       (name ca and segid $segid and resid $resid)
       (name  c and segid $segid and resid $resid)
       geometry

     evaluate ($phi2=mod($result+360,360))

     coord swap end

   end if

   identity (store1) (name n and bondedto
            (name c and resid $resid and segid $segid ))

   if ( $select = 1 ) then

     show element (store1) (attribute store1 > 0)
     evaluate ($id_next=$result)

     show (segid) (id $id_next)
     evaluate ($segid_next=$result)
     show (resid) (id $id_next)
     evaluate ($resid_next=$result)

     pick dihedral
       (name  n and segid $segid and resid $resid)
       (name ca and segid $segid and resid $resid)
       (name  c and segid $segid and resid $resid)
       (name  n and segid $segid_next and resid $resid_next)
       geometry

     evaluate ($psi1=mod($result+360,360))

     coord swap end

     pick dihedral
       (name  n and segid $segid and resid $resid)
       (name ca and segid $segid and resid $resid)
       (name  c and segid $segid and resid $resid)
       (name  n and segid $segid_next and resid $resid_next)
       geometry

     evaluate ($psi2=mod($result+360,360))

     coord swap end

   end if

   if ( $phi1 # -9999 ) then
     if ( $psi1 # -9999 ) then
       evaluate ($dphi=abs($phi1-$phi2))
       if ( $dphi > 180 ) then
         evaluate ($dphi=360-$dphi)
       end if
       evaluate ($dpsi=abs($psi1-$psi2))
       if ( $dpsi > 180 ) then
         evaluate ($dpsi=360-$dpsi)
       end if
       evaluate ($distance=sqrt($dphi^2 + $dpsi^2))
       set display=&list_outfile end
       display  $segid[a4]    $resid[a4]    $distance[f8.4]
       set display=OUTPUT end
       evaluate ($overall_phipsi=$overall_phipsi+$distance)
       evaluate ($nphipsi=$nphipsi+1)
     end if
   end if
   
 end loop dram

 evaluate ($overall_phipsi=$overall_phipsi/$nphipsi)
 set display=&list_outfile end
 display
 display >>>> overall difference= $overall_phipsi[f8.4] degrees
 display
 set display=OUTPUT end

 stop