      subroutine waters(itot,did,iclust,idim,iortho,cosi,iring,irf,irff,WATE001
     &    irfdid,ihand,rindid,ringsm,iri,cel,lo,hi,loang,hiang)         WATE002
c ----------------------------------------------------------------------WATE003
c     This program module is part of the author's X-ray Library         WATE004
c     and may not be sold or reproduced for commercial purposes.        WATE005
c     Although it has been tested carefully, there is no warranty       WATE006
c     regarding accuracy or freedom from side-effects.                  WATE007
c                   (C) 1986-1994 Bernhard Rupp                         WATE008
c ----------------------------------------------------------------------WATE009
c     WATERS searches for water polygons (pentagons, hexagons, etc)     WATE010
c ----------------------------------------------------------------------WATE011
                                                                        WATE012
      integer ipg(7),irff(iri)                                          WATE013
      real hi,lo,dd(7),hiang,loang,ang(7),rindid(iri*7,6),ringsm(iri)   WATE014
      dimension po(3),p1(3),p2(3),p3(3),p4(3),p5(3),p6(3),p7(3),xa(3)   WATE015
      dimension cel(6), cosi(3), did(iclust,idim)                       WATE016
                                                                        WATE017
c      print*,'-',ide(did(1,5)),'-'                                     WATE018
c      if (ide(did(1,5)).eq.'103 ') then                                WATE019
c         do 1000 i=1,itot                                              WATE020
c            write(7,*) (did(i,j),j=1,5)                                WATE021
c 1000    continue                                                      WATE022
c         pause                                                         WATE023
c      end if                                                           WATE024
                                                                        WATE025
                                                                        WATE026
                                                                        WATE027
                                                                        WATE028
c ---    reset the atom label storage --                                WATE029
         do 0100 j=1,7                                                  WATE030
            ipg(j)=0                                                    WATE031
 0100    continue                                                       WATE032
c ---    search did list for all atoms within water H-bond              WATE033
c        start atom is did(1) --                                        WATE034
         do 0011 j=1,3                                                  WATE035
          po(j)=did(1,j)                                                WATE036
 0011    continue                                                       WATE037
c ---    search for first occurance of H-bonded atoms                   WATE038
         do 0005 i1=2,itot                                              WATE039
          if(did(i1,4).gt.hi) goto 0006                                 WATE040
          if((did(i1,4).ge.lo).and.(did(i1,4).le.hi)) then              WATE041
            ipg(1)=i1                                                   WATE042
            dd(1)=did(i1,4)                                             WATE043
c           write(*,*)' Found suitable atom !',ifirs, ia                WATE044
c           write(*,'(8f8.4)') (did(i1,j),j=1,8)                        WATE045
c ---      make this atom new centre atom --                            WATE046
           do 0010 j=1,3                                                WATE047
            p1(j)=did(i1,j)                                             WATE048
 0010      continue                                                     WATE049
c ---      search did for atoms within water bond ---                   WATE050
           do 0008 i2=1,itot                                            WATE051
            do 0009 j=1,3                                               WATE052
             p2(j)=did(i2,j)                                            WATE053
             xa(j)=cel(j)*(p2(j)-p1(j))                                 WATE054
 0009       continue                                                    WATE055
            dd(2)=distan(iortho,xa,cosi)                                WATE056
c ---       check whether within H-bond --                              WATE057
            if((dd(2).ge.lo).and.(dd(2).le.hi)) then                    WATE058
             ipg(2)=i2                                                  WATE059
c ---        check if o.k.angle ---                                     WATE060
             call vecang(p1,po,p1,p2,ang(2))                            WATE061
             if((ang(2).ge.loang).and.(ang(2).le.hiang)) then           WATE062
c ---         and now for the third run..                               WATE063
              do 0012 i3=1,itot                                         WATE064
               do 0013 j=1,3                                            WATE065
                p3(j)=did(i3,j)                                         WATE066
                xa(j)=cel(j)*(p3(j)-p2(j))                              WATE067
 0013          continue                                                 WATE068
               dd(3)=distan(iortho,xa,cosi)                             WATE069
c ---          check whether within H-bond --                           WATE070
               if((dd(3).ge.lo).and.(dd(3).le.hi)) then                 WATE071
                ipg(3)=i3                                               WATE072
c ---           check if o.k.angle ---                                  WATE073
                call vecang(p2,p1,p2,p3,ang(3))                         WATE074
                if((ang(3).ge.loang).and.(ang(3).le.hiang))then         WATE075
c ---            and now for the fourth run ---                         WATE076
                 do 0014 i4=1,itot                                      WATE077
                  do 0015 j=1,3                                         WATE078
                   p4(j)=did(i4,j)                                      WATE079
                   xa(j)=cel(j)*(p4(j)-p3(j))                           WATE080
 0015             continue                                              WATE081
                  dd(4)=distan(iortho,xa,cosi)                          WATE082
c ---             check whether within H-bond --                        WATE083
                  if((dd(4).ge.lo).and.(dd(4).le.hi)) then              WATE084
                   ipg(4)=i4                                            WATE085
c ---              check if o.k.angle ---                               WATE086
                   call vecang(p3,p2,p3,p4,ang(4))                      WATE087
                   if((ang(4).ge.loang).and.(ang(4).le.hiang))then      WATE088
c ---               and now for the fifth run..                         WATE089
                    do 0016 i5=1,itot                                   WATE090
                     do 0017 j=1,3                                      WATE091
                      p5(j)=did(i5,j)                                   WATE092
                      xa(j)=cel(j)*(p5(j)-p4(j))                        WATE093
 0017                continue                                           WATE094
                     dd(5)=distan(iortho,xa,cosi)                       WATE095
c ---                check whether within H-bond --                     WATE096
                     if((dd(5).ge.lo).and.(dd(5).le.hi)) then           WATE097
                      ipg(5)=i5                                         WATE098
c ---                 check if o.k.angle ---                            WATE099
                      call vecang(p4,p3,p4,p5,ang(5))                   WATE100
                      if((ang(5).ge.loang).and.(ang(5).le.hiang))then   WATE101
c ---                  check if this point is starting point ---        WATE102
                       if((abs(did(i5,1)-did(1,1)).lt.0.0001).and.      WATE103
     &                    (abs(did(i5,2)-did(1,2)).lt.0.0001).and.      WATE104
     &                    (abs(did(i5,3)-did(1,3)).lt.0.0001)) then     WATE105
c ---                   calculate ringsm for this configuration ---     WATE106
                        ringsm(iring)=0                                 WATE107
                        do 0534 k=1,5                                   WATE108
                         ringsm(iring)=ringsm(iring)+dd(k)              WATE109
 0534                   continue                                        WATE110
c                        print*,'Ringsum=',ringsm                       WATE111
                        iprint=1                                        WATE112
                        do 0535 k=1,iring-1                             WATE113
                           if(abs(ringsm(iring)-ringsm(k)).lt.0.0001)   WATE114
     &                     iprint=0                                     WATE115
 0535                   continue                                        WATE116
c ---                   fill rindid with the atom id and cartesian XYZ -WATE117
                        if (iprint.eq.1) then                           WATE118
                         print*,' Found new pentagon'                   WATE119
c                         pause                                         WATE120
                         iring=iring+1                                  WATE121
                         irf=irf+1                                      WATE122
                         irff(irf)=5                                    WATE123
c ---                    at this point, we can also get first angle --  WATE124
                         call vecang(p5,p4,p5,p1,ang(1))                WATE125
c ---                    fill rindid with the atom id and cartesian XYZ WATE126
                         ik=5                                           WATE127
                         call rinfil(ik,rindid,iri,irfdid,ipg,did,      WATE128
     &                        iclust,idim,dd,ang)                       WATE129
                         else                                           WATE130
                          print*,' Pentagon already in list '           WATE131
c                          pause                                        WATE132
                         end if                                         WATE133
c ---                   skip the search for hexagons, new 5th atom --   WATE134
                        goto 0116                                       WATE135
                       end if                                           WATE136
c ---                  and now for the 6th run..                        WATE137
                       do 0018 i6=1,itot                                WATE138
                        do 0019 j=1,3                                   WATE139
                         p6(j)=did(i6,j)                                WATE140
                         xa(j)=cel(j)*(p6(j)-p5(j))                     WATE141
 0019                   continue                                        WATE142
                        dd(6)=distan(iortho,xa,cosi)                    WATE143
c ---                   check whether within H-bond --                  WATE144
                        if((dd(6).ge.lo).and.(dd(6).le.hi)) then        WATE145
                         ipg(6)=i6                                      WATE146
c ---                    check if o.k.angle ---                         WATE147
                         call vecang(p5,p4,p5,p6,ang(6))                WATE148
                         if((ang(6).ge.loang).and.(ang(6).le.hiang))thenWATE149
c ---                     check if this point is starting point ---     WATE150
                          if((abs(did(i6,1)-did(1,1)).lt.0.0001).and.   WATE151
     &                     (abs(did(i6,2)-did(1,2)).lt.0.0001).and.     WATE152
     &                     (abs(did(i6,3)-did(1,3)).lt.0.0001)) then    WATE153
c ---                      calculate ringsm for this configuration ---  WATE154
                           ringsm(iring)=0                              WATE155
                           do 0634 k=1,6                                WATE156
                            ringsm(iring)=ringsm(iring)+dd(k)           WATE157
 0634                      continue                                     WATE158
                           iprint=1                                     WATE159
                           do 0635 k=1,iring-1                          WATE160
                            if(abs(ringsm(iring)-ringsm(k)).lt.0.0001)  WATE161
     &                      iprint=0                                    WATE162
 0635                      continue                                     WATE163
c ---                      fill rindid with the atom id and cartesian XYWATE164
                           if (iprint.eq.1) then                        WATE165
                            print*,' Found new hexagon'                 WATE166
c                            pause                                      WATE167
                            iring=iring+1                               WATE168
                            irf=irf+1                                   WATE169
                            irff(irf)=6                                 WATE170
c ---                       at this point, we can also get first angle -WATE171
                            call vecang(p6,p5,p6,p1,ang(1))             WATE172
c ---                       fill rindid with the atom id and cartesian XWATE173
                            ik=6                                        WATE174
                            call rinfil(ik,rindid,iri,irfdid,ipg,did,   WATE175
     &                           iclust,idim,dd,ang)                    WATE176
                           else                                         WATE177
                            print*,' Hexagon already in list'           WATE178
c                            pause                                      WATE179
                           end if                                       WATE180
c ---                      skip the search for heptagons, new 6th atom -WATE181
                           goto 0118                                    WATE182
                          end if                                        WATE183
c ---                     and now for the 7th and last run ---          WATE184
                          do 0020 i7=1,itot                             WATE185
                           do 0021 j=1,3                                WATE186
                            p7(j)=did(i7,j)                             WATE187
                            xa(j)=cel(j)*(p7(j)-p6(j))                  WATE188
 0021                      continue                                     WATE189
                           dd(7)=distan(iortho,xa,cosi)                 WATE190
c ---                      check whether within H-bond --               WATE191
                           if((dd(7).ge.lo).and.(dd(7).le.hi)) then     WATE192
                            ipg(7)=i7                                   WATE193
c ---                       check if o.k.angle ---                      WATE194
                            call vecang(p6,p5,p6,p7,ang(7))             WATE195
                            if((ang(7).ge.loang).and.(ang(7).le.hiang)) WATE196
     &                       then                                       WATE197
c ---                        check if this point is starting point ---  WATE198
                             if((abs(did(i7,1)-did(1,1)).lt.0.0001).and.WATE199
     &                        (abs(did(i7,2)-did(1,2)).lt.0.0001).and.  WATE200
     &                        (abs(did(i7,3)-did(1,3)).lt.0.0001)) then WATE201
c ---                         calculate ringsm for this configuration --WATE202
                              ringsm(iring)=0                           WATE203
                              do 0734 k=1,7                             WATE204
                               ringsm(iring)=ringsm(iring)+dd(k)        WATE205
 0734                         continue                                  WATE206
                              iprint=1                                  WATE207
                              do 0735 k=1,iring-1                       WATE208
                              if(abs(ringsm(iring)-ringsm(k)).lt.0.0001)WATE209
     &                        iprint=0                                  WATE210
 0735                         continue                                  WATE211
c ---                         fill rindid with the atom id and cartesianWATE212
                              if (iprint.eq.1) then                     WATE213
                               iring=iring+1                            WATE214
                               irf=irf+1                                WATE215
                               irff(irf)=7                              WATE216
c ---                          at this point, we can also get first anglWATE217
                               call vecang(p7,p6,p7,p1,ang(1))          WATE218
                               if (ang(1).lt.0.001) ihand=ihand+1       WATE219
                               print*,' Found new seven membered ring'  WATE220
c ---                          fill rindid with the atom id and cartesiaWATE221
                               ik=7                                     WATE222
                               call rinfil(ik,rindid,iri,irfdid,ipg,did,WATE223
     &                             iclust,idim,dd,ang)                  WATE224
                              end if                                    WATE225
                              goto 0120                                 WATE226
                             end if                                     WATE227
                            end if                                      WATE228
                           end if                                       WATE229
 0120                      continue                                     WATE230
 0020                     continue                                      WATE231
                         end if                                         WATE232
                        end if                                          WATE233
 0118                   continue                                        WATE234
 0018                  continue                                         WATE235
                      end if                                            WATE236
                     end if                                             WATE237
 0116                continue                                           WATE238
 0016               continue                                            WATE239
                   end if                                               WATE240
                  end if                                                WATE241
 0014            continue                                               WATE242
                end if                                                  WATE243
               end if                                                   WATE244
 0012         continue                                                  WATE245
             end if                                                     WATE246
            end if                                                      WATE247
 0008      continue                                                     WATE248
          end if                                                        WATE249
 0005    continue                                                       WATE250
 0006    continue                                                       WATE251
c        end if water search stuff........                              WATE252
         return                                                         WATE253
         end                                                            WATE254
