! Module file: mtfautogenerate
!
! CNS module
! **********
!
! Authors: Axel Brunger
!
! copyright Yale University
!
! Function:
!
!    Automatically generates the molecular topology (mtf) from coordinates. 
!      Important notes:
!      1. Each distinct protein, DNA, or RNA chain must have a separate segid 
!         (or chainid if the chainid is non-blank and the "convert" 
!          option is set to true). 
!      2. Each contiguous protein, RNA, or RNA chain must not be disrupted by 
!         other types of residues or ligands.  Rather, these other residues 
!         should be listed after protein, RNA/DNA chains. 
!
! Requirements:
!    Must be called from the main level. 
!    Topology files must be read prior to executing this module.
!    The linkage files must be defined in the calling task file using
!    the declaraion symbols &link_infile_1,&link_infile_2,..


module {mtfautogenerate}

(
  &coordinate_infile; ! coordinate file 
  &convert=true;    ! if true, converts chain identifier to segid if the chain identifier is non-blank
  &separate=true;   ! if true, starts a new chain if the chain identifier or segid changes. 

  &atom_delete=(none); ! Deletes specified atoms after molecular topology generation and 
                       ! reading of coordinates. 

  &hydrogen_flag=false; ! Hydrogen flag - if false, hydrogens will be deleted from molecular 
                        ! topology regardless if they are present in the coordinate input file
                        ! or the topology file

  &break_cutoff=2.5; ! If greater than zero, automatically breaks polypeptide and 
                     ! polynucleotide chains based on distance.The default of 2.5 A 
                     ! should be reasonable for most cases. If the input structure 
                     ! has bad geometry it may be necessary to increase this distance 

  &disulphide_dist=3.0;  ! If greater than zero, automatically determines disulphides based on distance.
                         ! The default of 3.0 A should be reasonable for most cases. If the input
                         ! structure has bad geometry it may be necessary to increase this distance 
                         ! Note: the atom positions for the sulfurs have to be known.

  &carbo_dist=2.75;  ! If greater than zero, automatically determines links between N-linked carbohydrates.
                     ! The default of 2.75 A should be reasonable for most cases. It is not
                     ! recommended to increase this distance to prevent incorrect linkages.
                     ! Note: the atom positions for the linked atoms have to be known. 

  &patch_infile=""; ! For custom patches (to be applied after linear chain generation and before reading coordinates)

  &O5_becomes="O" ! if specified, renames the internal O5 in NDG and A2G carbohydrate residues to O (current PDB convention)
)

 checkversion 1.3

 set message ? end 
 evaluate ($message_old_mtfauto=$result)
 set echo ? end 
 evaluate ($echo_old_mtfauto=$result)
 
 set message=off echo=off end

 display 
 display  automatically generating molecular topology 
 display 

 set message=off echo=off end

 {- generating molecular topology -}
 
 segment
   chain

     if ( &convert = true ) then
       convert=true
     end if
     if ( &separate = true ) then
       separate=true
     end if

     evaluate ($counter=1)
     evaluate ($done=false)
     while ( $done = false ) loop read
        if ( &exist_link_infile_$counter = true ) then
           if ( &BLANK%link_infile_$counter = false ) then
              @@&link_infile_$counter
           end if
        else
           evaluate ($done=true)
        end if
        evaluate ($counter=$counter+1)
     end loop read
     set message=on end
     coordinates @@&coordinate_infile
     set message=off end
   end
 end
 
 {- execute custom patches -}
 if ( &BLANK%patch_infile = false ) then
    display
    display  executing external patch file
    display
    inline @&patch_infile
 end if

 {- rename atoms if requested -}

!!! ILE CD, and CTER OT1, OT2 have been renamed in the topology files
!!! if ( &BLANK%ile_CD_becomes = false ) then
!!!   do (name=&ile_CD_becomes) (resname ILE and name CD)
!!! end if
!!! if ( &BLANK%OT1_becomes = false ) then
!!!   do (name=&OT1_becomes) (name OT1)
!!! end if
!!! if ( &BLANK%OT2_becomes = false ) then
!!!   do (name=&OT2_becomes) (name OT2)
!!! end if

 if ( &BLANK%O5_becomes = false ) then
   do (name=&O5_becomes) (name O5 and ( resname NDG or resname A2G ))
 end if

 {- Done with all patches that potentially create new atoms and done with renaming atoms. -}
 {- We're ready to read the coordinates at this point. -}
 coordinates 
     if ( &convert = true ) then
       convert=true
     else
       convert=false
     end if
     @@&coordinate_infile

 {- we need to change O back to O5 in NDG and A2G since the links (see below) require it -}
 if ( &BLANK%O5_becomes = false ) then
   do (name=O5) (name &O5_becomes and ( resname NDG or resname A2G ))
 end if

!!! not needed anymore since there are now deoxy residues in the topology files
!!! {- patching of RNA to DNA -}
!!! evaluate ($counter=0)
!!! for $id in id ( tag and (&dna_sele) ) loop dna
!!!   evaluate ($counter=$counter+1)
!!!   show (segid) (id $id)
!!!   evaluate ($dna.segid.$counter=$result)
!!!   show (resid) (id $id)
!!!   evaluate ($dna.resid.$counter=$result)
!!! end loop dna
!!! evaluate ($dna.num=$counter)

!!! evaluate ($counter=0)
!!! while ($counter < $dna.num) loop dnap
!!!   evaluate ($counter=$counter+1)
!!!   patch deox reference=nil=(segid $dna.segid.$counter and
!!!                             resid $dna.resid.$counter) end
!!!   display  deox patch applied to $dna.segid.$counter[a4] $dna.resid.$counter[a4] 
!!! end loop dnap

 {- auto breakage based on distance -}
 if ( &break_cutoff > 0 ) then

   {- This is OK here after reading coordinates since breaking bonds does        -}
   {- does not change atom numbers or atom mappings.                             -} 

   {- polypeptides -}
   
   identity (store1) ( name C and bondedto(name CA) and bondedto(name O)
                  and bondedto(name N))
    show element (store1) (attribute store1 > 0 and not known)
    if ($select > 0) then
       eval ($id1=$result)
       show (segid) (id $id1)
       evaluate ($segid1=$result)
       show (resid) (id $id1)
       evaluate ($resid1=$result)
       show (resname) (id $id1)
       evaluate ($resname1=$result)
       display  unknown coordinates for segid $segid1 resname $resname1 resid $resid1 name C (cannot test for chain break)
    end if

    identity (store2) (name N and bondedto store1)
    show element (store2) (attribute store2 > 0 and not known)
    if ($select > 0) then
      eval ($id2=$result)
      show (segid) (id $id2)
      evaluate ($segid2=$result)
      show (resid) (id $id2)
      evaluate ($resid2=$result)
      show (resname) (id $id2)
      evaluate ($resname2=$result)
      display  unknown coordinates for segid $segid2 resname $resname2 resid $resid2 name N (cannot test for chain break)
    end if

    evaluate ($break=0)

    for $id1 in id (store1 and known) loop break

       show elem (store2) (store2 and bondedto id $id1 and known)

       if ( $select = 1 ) then
          evaluate ($id2=$result)

          show (x) ( id $id1 ) evaluate ($x1=$result)
          show (y) ( id $id1 ) evaluate ($y1=$result)
          show (z) ( id $id1 ) evaluate ($z1=$result)
          show (x) ( id $id2 ) evaluate ($x2=$result)
          show (y) ( id $id2 ) evaluate ($y2=$result)
          show (z) ( id $id2 ) evaluate ($z2=$result)
          evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

          if ( $distance > &break_cutoff ) then

             evaluate ($break=$break+1)
             show (segid) (id $id1)
             evaluate ($seg1.$break=$result)
             show (resid) (id $id1)
             evaluate ($res1.$break=$result)
             show (segid) (id $id2)
             evaluate ($seg2.$break=$result)
             show (resid) (id $id2)
             evaluate ($res2.$break=$result)
             show (resname) (id $id2)
             if ( $result = PRO ) then
                evaluate ($patch.$break=DPPP)
             elseif ( $result = CPR ) then
                evaluate ($patch.$break=DPPP)
             else
               evaluate ($patch.$break=DPEP)
             end if
           end if
         end if

    end loop break

   {- polynucleotides -}
   
   identity (store1) ( name O3' and bondedto(name C3') 
                  and bondedto(name P))
    show element (store1) (attribute store1 > 0 and not known)
    if ($select > 0) then
       eval ($id1=$result)
       show (segid) (id $id1)
       evaluate ($segid1=$result)
       show (resid) (id $id1)
       evaluate ($resid1=$result)
       show (resname) (id $id1)
       evaluate ($resname1=$result)
       display  unknown coordinates for segid $segid1 resname $resname1 resid $resid1 name O3' (cannot test for chain break)
    end if

    identity (store2) (name P and bondedto store1)
    show element (store2) (attribute store2 > 0 and not known)
    if ($select > 0) then
      eval ($id2=$result)
      show (segid) (id $id2)
      evaluate ($segid2=$result)
      show (resid) (id $id2)
      evaluate ($resid2=$result)
      show (resname) (id $id2)
      evaluate ($resname2=$result)
      display  unknown coordinates for segid $segid2 resname $resname2 resid $resid2 name P (cannot test for chain break)
    end if

    for $id1 in id (store1 and known) loop break

       show elem (store2) (store2 and bondedto id $id1 and known)

       if ( $select = 1 ) then
          evaluate ($id2=$result)

          show (x) ( id $id1 ) evaluate ($x1=$result)
          show (y) ( id $id1 ) evaluate ($y1=$result)
          show (z) ( id $id1 ) evaluate ($z1=$result)
          show (x) ( id $id2 ) evaluate ($x2=$result)
          show (y) ( id $id2 ) evaluate ($y2=$result)
          show (z) ( id $id2 ) evaluate ($z2=$result)
          evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

          if ( $distance > &break_cutoff ) then
             evaluate ($break=$break+1)
             show (segid) (id $id1)
             evaluate ($seg1.$break=$result)
             show (resid) (id $id1)
             evaluate ($res1.$break=$result)
             show (segid) (id $id2)
             evaluate ($seg2.$break=$result)
             show (resid) (id $id2)
             evaluate ($res2.$break=$result)
             show (resname) (id $id2)
             evaluate ($patch.$break=DNUC)
           end if
         end if

    end loop break

    evaluate ($counter=1)

    while ($counter <= $break) loop delete
      patch $patch.$counter
        reference=-=(segid $seg1.$counter and resid $res1.$counter)
        reference=+=(segid $seg2.$counter and resid $res2.$counter)
      end
      display  link removed (applied $patch.$counter): from \
$seg1.$counter[a4] $res1.$counter[a4] to $seg2.$counter[a4] $res2.$counter[a4]
      evaluate ($counter=$counter+1)
    end loop delete
 end if

 {- auto disulphide bond formation based on distance -}
 if (&disulphide_dist > 0) then

    {- this is OK here after reading the coordinates since the DISUL patch -}
    {- does not create or rename atoms -}
    evaluate ($disu=0)
 
    for $id1 in id ( resname CYS and name SG and known ) loop dis1

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)

      identity (store1) (all)

      for $id2 in id ( resname CYS and name SG and 
                  ( attr store1 > $id1 ) and known ) loop dis2

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        show (x) ( id $id1 ) evaluate ($x1=$result)
        show (y) ( id $id1 ) evaluate ($y1=$result)
        show (z) ( id $id1 ) evaluate ($z1=$result)
        show (x) ( id $id2 ) evaluate ($x2=$result)
        show (y) ( id $id2 ) evaluate ($y2=$result)
        show (z) ( id $id2 ) evaluate ($z2=$result)
        evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

        if ( $distance <= &disulphide_dist ) then

          evaluate ($disu=$disu+1)
          evaluate ($seg1.$disu=$segid1)
          evaluate ($seg2.$disu=$segid2)
          evaluate ($res1.$disu=$resid1)
          evaluate ($res2.$disu=$resid2)
        end if

      end loop dis2

    end loop dis1

    evaluate ($counter=1)
    while ( $counter <= $disu ) loop disu
       patch disu
          reference=1=(segid $seg1.$counter and resid $res1.$counter)
          reference=2=(segid $seg2.$counter and resid $res2.$counter)
       end

       display  disulphide added: from \
$seg1.$counter[a4] $res1.$counter[a4] to $seg2.$counter[a4] $res2.$counter[a4]
       evaluate ($counter=$counter+1)
    end loop disu

 end if

 {- auto linkage of n-linked carbohydrates -}
 if (&carbo_dist > 0) then

    {- this is OK here after reading the coordinates since the various -}
    {= patches used here do not create or rename atoms -}

    {- first we build a data structure for the patches, and then apply the patches later -}
    {- the reason we do it this way is that patch may delete or add atoms which will currently 
       reset the STORE and loop arrays -}

    eval ($npatch=0)

    {- NAG-C2 (or NDG-C2) to ASN-ND2 -}
    for $id1 in id ( (resname NAG or resname NDG) and name C1 and known ) loop B1Na

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)
      show (resname) (id $id1)
      evaluate ($resname1=$result)

      for $id2 in id ( resname ASN and name ND2 and known ) loop B1Nb

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        show (x) ( id $id1 ) evaluate ($x1=$result)
        show (y) ( id $id1 ) evaluate ($y1=$result)
        show (z) ( id $id1 ) evaluate ($z1=$result)
        show (x) ( id $id2 ) evaluate ($x2=$result)
        show (y) ( id $id2 ) evaluate ($y2=$result)
        show (z) ( id $id2 ) evaluate ($z2=$result)
        evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

        if ( $distance <= &carbo_dist ) then
           evaluate ($npatch=$npatch+1)
           evaluate ($seg1.$npatch=$segid1)
           evaluate ($seg2.$npatch=$segid2)
           evaluate ($res1.$npatch=$resid1)
           evaluate ($res2.$npatch=$resid2)
           if ($resname1="NDG") then
              evaluate ($pres.$npatch="B1N")  ! force to beta anomer, ATB 02/03/2010
              display   warning: unusual NDG to ASN patch for NDG segid= $segid1 resid= $resid1 - forced to beta anomer
           else
              evaluate ($pres.$npatch="B1N")
           end if
       end if

      end loop B1Nb

    end loop B1Na

    {- A23: GAL to SIA -}
    for $id1 in id ( resname GAL and name C2 and known ) loop A23a

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)

      for $id2 in id ( resname SIA and name O3 and known ) loop A23b

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        show (x) ( id $id1 ) evaluate ($x1=$result)
        show (y) ( id $id1 ) evaluate ($y1=$result)
        show (z) ( id $id1 ) evaluate ($z1=$result)
        show (x) ( id $id2 ) evaluate ($x2=$result)
        show (y) ( id $id2 ) evaluate ($y2=$result)
        show (z) ( id $id2 ) evaluate ($z2=$result)
        evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

        if ( $distance <= &carbo_dist ) then
           evaluate ($npatch=$npatch+1)
           evaluate ($seg1.$npatch=$segid1)
           evaluate ($seg2.$npatch=$segid2)
           evaluate ($res1.$npatch=$resid1)
           evaluate ($res2.$npatch=$resid2)
           evaluate ($pres.$npatch="A23")
       end if

      end loop A23b

    end loop A23a

    {- A26: GAL to SIA -}
    for $id1 in id ( resname GAL and name C2 and known ) loop A26a

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)

      for $id2 in id ( resname SIA and name O6 and known ) loop A26b

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        show (x) ( id $id1 ) evaluate ($x1=$result)
        show (y) ( id $id1 ) evaluate ($y1=$result)
        show (z) ( id $id1 ) evaluate ($z1=$result)
        show (x) ( id $id2 ) evaluate ($x2=$result)
        show (y) ( id $id2 ) evaluate ($y2=$result)
        show (z) ( id $id2 ) evaluate ($z2=$result)
        evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

        if ( $distance <= &carbo_dist ) then
           evaluate ($npatch=$npatch+1)
           evaluate ($seg1.$npatch=$segid1)
           evaluate ($seg2.$npatch=$segid2)
           evaluate ($res1.$npatch=$resid1)
           evaluate ($res2.$npatch=$resid2)
           evaluate ($pres.$npatch="A26")
       end if

      end loop A26b

    end loop A26a
 
    define 
    (
        {- define all alpha anomer carbo residues for the links below -}
      alpha=   ( resname GLC or resname GLA or resname MAN or resname NDG or resname A2G or resname FUC or resname FCA );
        {- define all beta anomer carbo residues for the links below -}
      beta=    ( resname BGC or resname GAL or resname BMA or resname NAG or resname FUL or resname FCB or resname XYP )
    )

    {- now the beta-anomer to any links -}
    {- for B12: beta(1,2)  link from beta-anomer to any -}
    ident (store2) ( ( &alpha or &beta) and name O2 and known ) 
    {- for B14: beta(1,4)  link from beta-anomer to any -}
    ident (store4) ( ( &alpha or &beta) and name O4 and known ) 

    for $id1 in id ( &beta and name C1 and known ) loop B1

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)
      show (resname) (id $id1)
      evaluate ($resname1=$result)

      for $id2 in id ( recall2 ) loop B12

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then

           show (x) ( id $id1 ) evaluate ($x1=$result)
           show (y) ( id $id1 ) evaluate ($y1=$result)
           show (z) ( id $id1 ) evaluate ($z1=$result)
           show (x) ( id $id2 ) evaluate ($x2=$result)
           show (y) ( id $id2 ) evaluate ($y2=$result)
           show (z) ( id $id2 ) evaluate ($z2=$result)
           evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

           if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUL") then
                 evaluate ($pres.$npatch="B12L")
              else
                 evaluate ($pres.$npatch="B12")
              end if
          end if
        end if

      end loop B12

      for $id2 in id ( recall4 ) loop B14

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)
        
        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then

          show (x) ( id $id1 ) evaluate ($x1=$result)
          show (y) ( id $id1 ) evaluate ($y1=$result)
          show (z) ( id $id1 ) evaluate ($z1=$result)
          show (x) ( id $id2 ) evaluate ($x2=$result)
          show (y) ( id $id2 ) evaluate ($y2=$result)
          show (z) ( id $id2 ) evaluate ($z2=$result)
          evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

          if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUL") then
                 evaluate ($pres.$npatch="B14L")
              else
                 evaluate ($pres.$npatch="B14")
              end if
          end if
        end if

      end loop B14

    end loop B1

    {- now the alpha anomer to any links -}
    {- for A12: alpha(1,2)  link from alpha-anomer to any -}
    ident (store2) ( ( &alpha or &beta) and name O2 and known ) 
    {- for A13: alpha(1,3)  link from alpha-anomer to any -}
    ident (store3) ( ( &alpha or &beta) and name O3 and known ) 
    {- for A14: alpha(1,4)  link from alpha-anomer to any -}
    ident (store4) ( ( &alpha or &beta) and name O4 and known ) 
    {- for A16: alpha(1,6)  link from alpha-anomer to any -}
    ident (store6) ( ( &alpha or &beta) and name O6 and known ) 

    for $id1 in id ( &alpha and name C1 and known ) loop A1

      show (segid) (id $id1)
      evaluate ($segid1=$result)
      show (resid) (id $id1)
      evaluate ($resid1=$result)
      show (resname) (id $id1)
      evaluate ($resname1=$result)

      for $id2 in id ( recall2 ) loop A12

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then

           show (x) ( id $id1 ) evaluate ($x1=$result)
           show (y) ( id $id1 ) evaluate ($y1=$result)
           show (z) ( id $id1 ) evaluate ($z1=$result)
           show (x) ( id $id2 ) evaluate ($x2=$result)
           show (y) ( id $id2 ) evaluate ($y2=$result)
           show (z) ( id $id2 ) evaluate ($z2=$result)
           evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

           if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUC") then
                 evaluate ($pres.$npatch="A12L")
              else
                 evaluate ($pres.$npatch="A12")
              end if
          end if
        end if

      end loop A12

      for $id2 in id ( recall3 ) loop A13

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then

           show (x) ( id $id1 ) evaluate ($x1=$result)
           show (y) ( id $id1 ) evaluate ($y1=$result)
           show (z) ( id $id1 ) evaluate ($z1=$result)
           show (x) ( id $id2 ) evaluate ($x2=$result)
           show (y) ( id $id2 ) evaluate ($y2=$result)
           show (z) ( id $id2 ) evaluate ($z2=$result)
           evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

           if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUC") then
                 evaluate ($pres.$npatch="A13L")
              else
                 evaluate ($pres.$npatch="A13")
              end if
          end if
        end if

      end loop A13

      for $id2 in id ( recall4 ) loop A14

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)
        
        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then
 
           show (x) ( id $id1 ) evaluate ($x1=$result)
           show (y) ( id $id1 ) evaluate ($y1=$result)
           show (z) ( id $id1 ) evaluate ($z1=$result)
           show (x) ( id $id2 ) evaluate ($x2=$result)
           show (y) ( id $id2 ) evaluate ($y2=$result)
           show (z) ( id $id2 ) evaluate ($z2=$result)
           evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

           if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUC") then
                 evaluate ($pres.$npatch="A14L")
              else
                 evaluate ($pres.$npatch="A14")
              end if
          end if
        end if

      end loop A14

      for $id2 in id ( recall6 ) loop A16

        show (segid) (id $id2)
        evaluate ($segid2=$result)
        show (resid) (id $id2)
        evaluate ($resid2=$result)

        evaluate ($identical=false)
        if ($resid1 = $resid2) then
            if ($segid1 = $segid2) then
               evaluate ($identical=true)
            end if
        end if

        if ($identical=false) then
           show (x) ( id $id1 ) evaluate ($x1=$result)
           show (y) ( id $id1 ) evaluate ($y1=$result)
           show (z) ( id $id1 ) evaluate ($z1=$result)
           show (x) ( id $id2 ) evaluate ($x2=$result)
           show (y) ( id $id2 ) evaluate ($y2=$result)
           show (z) ( id $id2 ) evaluate ($z2=$result)
           evaluate ($distance=sqrt( ($x1-$x2)^2 + ($y1-$y2)^2 + ($z1-$z2)^2 ))

           if ( $distance <= &carbo_dist ) then
              evaluate ($npatch=$npatch+1)
              evaluate ($seg1.$npatch=$segid1)
              evaluate ($seg2.$npatch=$segid2)
              evaluate ($res1.$npatch=$resid1)
              evaluate ($res2.$npatch=$resid2)
              if ($resname1="FUC") then
                 evaluate ($pres.$npatch="A16L")
              else
                 evaluate ($pres.$npatch="A16")
              end if
          end if
        end if

      end loop A16

    end loop A1


    {- now we apply all the patches -}
    evaluate ($counter=1)
    while ( $counter <= $npatch ) loop patc

       patch $pres.$counter
          reference=-=(segid $seg1.$counter and resid $res1.$counter)
          reference=+=(segid $seg2.$counter and resid $res2.$counter)
       end

       display  $pres.$counter link added: $seg1.$counter[a4] $res1.$counter[a4] to $seg2.$counter[a4] $res2.$counter[a4]
       evaluate ($counter=$counter+1)
    end loop patc

 end if


 {- list all non-hydrogen atoms that have unknown coordinates -}
 show sum(1) (not (known) and not hydrogen)
 if ( $result > 0) then
  if ( $result < 1000 ) then
   for $id in id (not(known) and not hydrogen) loop print
     show (segid) (id $id)
     evaluate ($segid=$result)
     show (resname) (id $id)
     evaluate ($resname=$result)
     show (resid) (id $id)
     evaluate ($resid=$result)
     show (name) (id $id)
     evaluate ($name=$result)
     display  unknown coordinates for atom: $segid[a4] $resname[a4] $resid[a4] $name[a4]
   end loop print
  else
   display  unknown coordinates for more than 100 atoms
  end if
 end if

 if (&hydrogen_flag=false) then
   delete selection=( hydrogen ) end
 end if
 
 delete selection=&atom_delete end

 {- reapply the renaming of NDG, A2G O5 to O -}
 if ( &BLANK%O5_becomes = false ) then
   do (name=&O5_becomes) (name O5 and ( resname NDG or resname A2G ))
 end if

 display 

 set message=$message_old_mtfauto echo=$echo_old_mtfauto end