# Copyright (c) 2002-08 Peter Guntert. All rights reserved.
## 7MACROS: stereoassign - CYANA macro
##
## Parameters: selection=<string>             (default: <blank>)
##             tfcut=<real>                   (default: 0.0)
##             conformations=<string>         (default: <blank>)
##
## Tries to find stereospecific assignments by systematic analysis of the
## local conformation of a molecular fragment with grid searches. If there
## are n pairs of (stereospecificly unassigned) diastereotopic substituents
## within the molecular fragment, 2^n grid searches will be performed, one for
## each possible combination of stereospecific assignments. If a (connected)
## angle selection is specified, then it defines the molecular fragment that
## will be analyzed; otherwise, the fragment set in the preceding grid
## fragment command will be used. The parameter tfcut has the same meaning as
## in the grid search command. Optionally, the total number of allowed
## conformations can be stored in a variable with the given name. Grid
## searches are restricted to values of the torsion angles given in the
## standard grid memory, A, on input. On output, the allowed values of the
## torsion angles are again stored in grid memory A.

var echo info i j l ip jp np ninp nrev mconf str swapp swapped

syntax $macro:selection=*= tfcut=@r=0.0 conformations=*=

echo:=off; info:=none

# ----- if a selection is given on the command line, define a new fragment -----

if ('$selection'.ne.' ') grid fragment "$selection"

# ----- find pairs of diasteretopic atoms -----

do np 1 15
  j=diastereotopic(np)
  if (j.eq.0) break
  ip(np)=j
  swapped(np)=.false.
  ninp(np)=0
  nrev(np)=0
end do
np=np-1
do i 1 np
  atom stereo "$atom(ip(i))"
end do

# ----- loop over all permutations of stereospecific assignments -----

grid memory input=A
grid memory combi=.false.
mconf=0
do jp 1 2**np
  l=jp-1
  do i 1 np
    swapp=mod(l,2).eq.1
    l=l/2
    if (swapp.neqv.swapped(i)) grid swap "$atom(ip(i))"
    swapped(i)=swapp
  end do
  grid memory A=input
  grid search tfcut=$tfcut
  if (nconf.lt.0) then
    do i 1 np
      atom stereo "$atom(ip(i))" delete
    end do
    grid memory A=input
    mconf=-1
    go to end
  else if (nconf.gt.0) then
    grid memory combi=combi.or.A
    mconf=mconf+nconf
    l=jp-1
    do i 1 np
      if (mod(l,2).eq.0) then
        ninp(i)=ninp(i)+nconf
      else
        nrev(i)=nrev(i)+nconf
      end if
      l=l/2
    end do
  end if
end do

# ----- set stereospecific assignments and angles ranges found -----

do i 1 np
  if (ninp(i).eq.0 .eqv. nrev(i).eq.0) atom stereo "$atom(ip(i))" delete
end do
if (mconf.eq.0) grid memory combi=input
grid memory A=combi

# ----- write output -----

unset info
do i 1 np
  if (ninp(i).gt.0 .and. nrev(i).eq.0) then
    str:=assigned as input
  else if (ninp(i).eq.0 .and. nrev(i).gt.0) then
    str:=assigned reversed, swapped
    atom swap "$atom(ip(i))" info=none
  else if (ninp(i).eq.0 .and. nrev(i).eq.0) then
    str:=inconsistent data
  else if (ninp(i).lt.0 .or. nrev(i).lt.0) then
    str:=not finished
  else
    str:=unassigned
  end if
  print "    $anam(ip(i))/$atomlong(abs(stereopartner(ip(i)))): $ninp(i)/$nrev(i) conformations, $str."
end do
if (np.eq.0) then
  if (nconf.gt.0) then
    print "    No stereo pair: $nconf conformations."
  else
    print "    No stereo pair: $nconf conformations, inconsistent data."
  end if
end if

end:
grid memory input=
grid memory combi=
grid fragment none
if ('$conformations'.ne.' ') external $conformations=mconf