# Copyright (c) 2002-08 Peter Guntert. All rights reserved.
## 7MACROS: cashiftaco - CYANA macro
##
## Parameters: offset=<real>                  (default: 0.0)
##             prooffset=<real>               (default: 2.0)
##
## Generates restraints for the backbone dihedral angles f and y in proteins
## by comparing the Ca chemical shifts with the corresponding random coil
## values of Spera & Bax (1991). Angle restraints are derived according to
## the rules of Luginbuehl et al. (1995). The Ca random coil shifts are
## relative to internal TSP. Optionally, an offset Dw is added to the
## chemical shifts in the proton list. A warning is printed for Ca chemical
## shifts that deviate by more than 15 ppm from their random coil value.

var echo aco rnames crc i j c r rn s ir jr n

syntax $macro:offset=@r=0.0 prooffset=@r=2.0

echo:=off
rnames:= ALA  ARG  ASN  ASP  GLN  GLU  ILE  LEU \
         LYS  MET  PHE  SER  THR  TRP  TYR  VAL
rnames='$rnames'
#Random coil shifts from Spera & Bax:
crc:=    52.3,56.1,52.8,54.0,56.1,56.4,61.3,55.1,\
         56.5,55.3,58.0,58.2,62.1,57.7,58.1,62.3
#Random coil shifts from Richarz & W"uthrich:
#crc:=   52.3,56.1,53.0,54.2,55.6,56.9,61.1,55.3,\
#        56.1,55.5,57.7,58.1,61.6,57.2,57.8,62.2

aco:=d$getpid.ca
print "\# Angle restraints from CA chemical shifts" >$aco
n=0
do j 1 na
  if (anam(j).eq.'CA') then
    if (shift(j,1).lt.990.0) then
      s=shift(j,1); jr=iar(j); r=rnam(jr)
      i=(index(rnames,r(1:3))+3)/4
      if (i.ge.1) then
        c=s+offset-crc(i); ir=rnum(jr)
        if (jr.lt.nr) then
          rn=rnam(jr+1)
          if (rn(1:3).eq.'PRO') c=c+prooffset
        end if
        if (abs(c).gt.15.0) \
          print "\# *** WARNING: Unusual CA shift for $r $ir: $s(F5.1) ppm."
        if (abs(c).gt.1.5) n=n+2
        if (c.gt.2.0) then
          print "$ir(I4) $r(1:4)  PHI    -120.0   -20.0" >>
          print "$ir(I4) $r(1:4)  PSI    -100.0     0.0" >>
        else if (c.gt.1.5) then
          print "$ir(I4) $r(1:4)  PHI    -120.0    80.0" >>
          print "$ir(I4) $r(1:4)  PSI    -100.0    60.0" >>
        else if (c.lt.-1.5) then
          print "$ir(I4) $r(1:4)  PHI    -200.0   -80.0" >>
          print "$ir(I4) $r(1:4)  PSI      40.0   220.0" >>
        end if
      end if
    end if
  end if
end do
print >>.
read aco $aco append info=none
remove $aco
print "    $n angle restraints derived from CA chemical shifts."