program xshell
      implicit none
      integer*4 ti,tj,i,iargc,ic,ie,ih,ih2,ihe,ii,in,inwmax,io,is,
     1ist,it,ite,iuu,iw,ix,iy,iz,j,jj,jjk,k,n22,nat,natt,natu,nw,nwmax,
     1natw,ntx,natuo
      real*8 a1,a2,cov,cov2,dd,dij,dist,distmin,doh,dx,dy,px,py,pz,dz,
     1rb,rmax,rmin,rt,x1,xi,xii,xj,xjj,xk,xl,y1,yi,yii,yj,yjj,yk,yl,
     1z1,zi,zii,zj,zjj,zk,zl
      real*4 bonding(88)
c      derived from atomic radii, some made longer for MD snapshots:
c       Du    H   He   Li   Be    B    C    N    O    F
c       Ne   Na   Mg   Al   Si    P    S   Cl   Ar    K   Ca   Sc   Ti
c        V   Cr   Mn   Fe   Co   Ni   Cu   Zn   Ga   Ge   As 
      data bonding/
     10.50,0.42,0.98,1.28,0.95,0.87,0.92,0.90,0.88,0.77,
     10.76,1.59,1.41,1.23,1.16,1.11,1.20,1.04,1.03,2.08,1.79,1.49,1.37,
     11.27,1.23,1.22,1.22,1.21,1.20,1.22,1.30,1.31,1.27,1.25,1.21,1.19,
     11.17,2.21,1.96,1.67,1.50,1.39,1.35,1.32,1.30,1.30,1.33,1.39,1.53,
     11.49,1.46,1.45,1.41,1.38,1.36,2.40,2.03,1.74,1.70,1.70,1.69,1.68,
     11.67,1.90,1.66,1.64,1.64,1.63,1.62,1.61,1.75,1.61,1.49,1.39,1.35,
     11.33,1.31,1.32,1.35,1.39,1.54,1.53,1.52,1.51,1.51,1.50,1.49,0.10/
      character*50 st
      character*80 s80
      logical llayer,lsolulst,lper,lexcl,lhyd,lrev,ldw,dwb
      integer*4,allocatable:: as(:),qt(:),isu(:),i22(:),kt(:),ity(:),
     1jo(:),jh(:),jt(:),io1(:),io2(:),item(:),isw(:),isuo(:),asx(:)
      real*8   ,allocatable::  r(:,:),rtem(:,:),rx(:,:)
      logical,  allocatable::lsa(:)

      write(6,600)
600   format(' Extracts hydration shell from FILE.X',/,/,
     1' XSHELL.OPT  ....   list of options',/,
     2' line arguments: rmax rmin (for rmax<0 only)',/)

      call readopt(cov,lsolulst,lper,lrev,lexcl,nwmax,lhyd,ldw,px,py,pz,
     1dwb)

      call setatmax(ntx)
      allocate(as(ntx),qt(ntx),isu(ntx),i22(ntx),kt(ntx),ity(ntx),
     1jo(ntx),jh(ntx),jt(ntx),io1(ntx),io2(ntx),r(ntx,3),rtem(ntx,3),
     1lsa(ntx),item(ntx),isw(ntx),isuo(ntx),rx(ntx,3),asx(ntx))
c
      if(lsolulst)then
c      solute defined in SOLU.LST:
       call ral('SOLU.LST',natuo,isuo)
      else
       write(6,*)'Water as solvent supposed'
      endif

c     "dead wood" atoms in DW.LST:
      if(ldw)call ral('DW.LST',natw,isw)

      rmax=0.0d0
      rmin=0.10d0
      if(nwmax.ne.0)write(6,*)' Maximum number of solvent molecules :',
     1nwmax
      if(lhyd)then
       write(6,*)' Hydrogen bond solvent only (close to polar atoms) '
      else
       write(6,*)' All solvent'
      endif
      write(6,800)cov
800   format(' covalent limit ',f10.2,' A')
      cov2=cov**2

      if(iargc().eq.0)then
       write(6,*)' Cut-off distance (A):'
       read(5,*)rmax
       llayer=rmax.lt.0.0d0
       if(llayer)then
        rmax=abs(rmax)
        write(6,*)' rmin:'
        read(5,*)rmin
       endif
      else
       call getarg(1,st)
       read(st,*)rmax
       write(6,602)rmax
602    format(' Cut-off distance:',f8.2,' A')
       llayer=rmax.lt.0.0d0
       if(llayer)then
        rmax=abs(rmax)
        call getarg(2,st)
        read(st,*)rmin
        write(6,*)' rmin ',rmin
       endif
      endif

      is=0
      inwmax=0
      ix=0
      open(3,file='FILE.X',status='old')
      open(33,file='NEWFILE.X')
1     read(3,300,end=999,err=999)s80
300   format(a80)
      read(3,*,end=999,err=999)nat
      is=is+1
      do 2 i=1,nat
2     read(3,*)as(i),(r(i,k),k=1,3)

c     remember for later:
      rx=r
      asx=as

      natu=natuo
      isu=isuo
c     reduce number of atoms by the read wood list
      if(ldw)call as2(ntx,nat,natu,natw,isu,isw,lsolulst,as,r)

c     assign solute atoms:
      if(lsolulst)call as1(ntx,nat,natu,lsa,isu)


c     iw ... total number of solvent atoms:
      iw=0
c     numbers of H(-N) and O(=C):
      do 3 i=1,nat
c     33333333333333333333333333333333333333333333333333333333333333333
      jt(i)=0
c     jt ... indicate, if this atom already taken
      xi=r(i,1)
      yi=r(i,2)
      zi=r(i,3)
c
c     io, in, ih,,it
c      ... how many covalently bond Os,Ns,Hs,electronegative and any atoms:
      io=0
      in=0
      ih=0
      ihe=0
      it=0

      do 4 j=1,nat
      if(j.ne.i)then
       xj=r(j,1)
       yj=r(j,2)
       zj=r(j,3)
       dij=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
       if(dij.lt.1.44d0.and.as(j).eq.8)then
        io=io+1
        jo(io)=j
       endif
       if(dij.lt.1.96d0.and.as(j).eq.7)in=in+1
       if(dij.lt.1.44d0.and.as(j).eq.1)then
        ih=ih+1
        jh(ih)=j
       endif
       if(dij.lt.1.96d0             )it=it+1
      endif
4     continue
c
      ity(i)=0
c     ity : 0 -  general (aliphatic)
c     ity : 1 -  water atom
c     ity : 2 -  (non-water) acidic hydrogen
c     ity : 3 -  (non-water) electronegative atom
c     ity : 5 -  coordinating metal

c     water oxygen, remember the other two hydrogens
      if(as(i).eq.8.and.ih.eq.2)then
       ity(i)=1
       iw=iw+1
       io1(i)=jh(1)
       io2(i)=jh(2)
      endif

      if(as(i).eq.1.and.io.eq.1)then
c      check water hydrogen:
       xk=r(jo(1),1)
       yk=r(jo(1),2)
       zk=r(jo(1),3)
       ih2=0
       jjk=0
       do 42 jj=1,nat
       xl=r(jj,1)
       yl=r(jj,2)
       zl=r(jj,3)
       doh=(xk-xl)**2+(yk-yl)**2+(zk-zl)**2
       if(jj.ne.i.and.jj.ne.jo(1).and.doh.lt.1.44d0)then
c       1.44=(1.2A)**2
        ih2=ih2+1
        jjk=jj
       endif
42     continue
       if(ih2.eq.1)then
        ity(i)=1
        iw=iw+1
c       to each water atom, define the other two
        io1(i)=jo(1)
        io2(i)=jjk
       endif
      endif

      if(lsolulst.or.ity(i).eq.0)then
c      hydrogen bonded to electronegative atom
       if(as(i).eq.1)then
        do 41 j=1,nat
        if((as(i).eq.8.or.as(i).eq.7.or.as(i).eq.16.
     1  or.as(i).eq.34.or.as(i).eq.9.or.as(i).eq.17.
     2  or.as(i).eq.35.or.as(i).eq.53).and.j.ne.i)then
         xj=r(j,1)
         yj=r(j,2)
         zj=r(j,3)
         dij=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
         if(dij.lt.1.44d0.and.as(j).eq. 7) ihe=ihe+1
         if(dij.lt.1.44d0.and.as(j).eq. 8) ihe=ihe+1
         if(dij.lt.1.44d0.and.as(j).eq. 9) ihe=ihe+1
         if(dij.lt.1.94d0.and.as(j).eq.17) ihe=ihe+1
         if(dij.lt.1.94d0.and.as(j).eq.16) ihe=ihe+1
         if(dij.lt.2.24d0.and.as(j).eq.34) ihe=ihe+1
         if(dij.lt.2.24d0.and.as(j).eq.35) ihe=ihe+1
        endif
41      continue
        if(ihe.eq.1)ity(i)=2
       endif
c      electronegative atom
c      O N S Se F Cl Br I:
       if(as(i).eq.8.or.as(i).eq.7.or.as(i).eq.16.
     1 or.as(i).eq.34.or.as(i).eq.9.or.as(i).eq.17.
     2 or.as(i).eq.35.or.as(i).eq.53)ity(i)=3
c      Sc-Zn:
       if(as(i).le.30.and.as(i).ge.21)ity(i)=5
      endif
3     continue
c     33333333333333333333333333333333333333333333333333333333333333333
c
c     number of solute atoms:
      iuu=0
c     number of atoms in the chosen cluster:
      natt=0
c     number of waters in the chosen cluster:
      nw=0
c 
c     31313131313131313131313131313131313131313131313131313131313131313131
      do 31 i=1,nat
c     atoms not taken yet:
      if(jt(i).eq.0)then
       ic=0
c      take solute atoms:
       if(lsolulst)then
        do 311 ii=1,natu
        if(isu(ii).eq.i)then
         ic=1
         goto 312
        endif
311     continue
       else
c       take non-water atoms if solute not defined:
        if(ity(i).ne.1)ic=1
       endif
312    if(ic.eq.1)then
        iuu=iuu+1
        natt=natt+1
        do 32 k=1,3
32      rtem(natt,k)=r(i,k)
        item(natt)=i
        qt(natt)=as(i)
        jt(i)=1
       else
        x1=r(i,1)
        y1=r(i,2)
        z1=r(i,3)
c
c       look for covalently-bond atoms:
        n22=1
        i22(1)=i
        if(lsolulst)then
         do 72 ii=1,nat
72       kt(ii)=0
         kt(i)=1
71       do 7 jj=1,n22
         x1=r(i22(jj),1)
         y1=r(i22(jj),2)
         z1=r(i22(jj),3)
         tj=as(i22(jj))+1
         rb=dble(bonding(tj))
         do 7 ii=1,nat
         if(.not.lsa(ii).and.kt(ii).eq.0)then
          a1=(r(ii,1)-x1)**2
          ti=as(ii)+1
          if(ti.ne.2.or.tj.ne.2)then
           if(tj.lt.1.or.tj.gt.88)then
            rt=cov2
           else
            rt=(rb+dble(bonding(ti)))**2
           endif
           if(a1.lt.rt)then
            a2=a1+(r(ii,2)-y1)**2
            if(a2.lt.rt)then
             if(a2+(r(ii,3)-z1)**2.lt.rt)then
              kt(ii)=1
              n22=n22+1
              i22(n22)=ii
              goto 71
             endif
            endif
           endif
          endif
         endif
7        continue
        else
c        the other two water atoms:
         n22=3
         i22(2)=io1(i)
         i22(3)=io2(i)
        endif

c       clone the geometry for the periodic box:
        if(lper)then
         ist=-1
         ie=1
        else
         ist=0
         ie=0
        endif
c       clonecloneclonecloneclonecloneclonecloneclonecloneclone
        do 1010 ix=ist,ie
        dx=px*dble(ix)
        do 1010 iy=ist,ie
        dy=py*dble(iy)
        do 1010 iz=ist,ie
        dz=pz*dble(iz)
c
c       how far are molecule's atoms from the solute:
        distmin=9.09d9
        it=0
c       loop over atoms of molecule of atom i:
        do 6 jj=1,n22
        xjj=r(i22(jj),1)+dx
        yjj=r(i22(jj),2)+dy
        zjj=r(i22(jj),3)+dz

        dist=9.0d9
        do 61 ii=1,nat
        ic=0
        if(lsolulst)then
         if(lsa(ii))then
          if(lhyd)then
c          polar solute atom:
           if(ity(ii).eq.2.or.ity(ii).eq.3.or.ity(ii).eq.5)ic=1
          else
c          any solute atom:
           ic=1
          endif
         endif
        else
         if(lhyd)then
c         polar solute atom:
          if(ity(ii).eq.2.or.ity(ii).eq.3)ic=1
c         coordinating metal solute atom:
          if(ity(ii).eq.5)ic=1
         else
c         any solute atom:
          if(ity(ii).ne.1)ic=1
         endif
        endif
        if(ic.eq.1)then
         xii=r(ii,1)
         yii=r(ii,2)
         zii=r(ii,3)
         dd=(xii-xjj)**2+(yii-yjj)**2+(zii-zjj)**2
         if(dd.lt.dist)dist=dd
        endif
61      continue
        dist=sqrt(dist)
c
        if(dist.lt.distmin)distmin=dist
        if(dist.gt.rmax.or.(llayer.and.(dist.le.rmin)))it=1
6       continue


        dist=distmin
        ite=0
        if(dist.gt.rmax.or.(llayer.and.(dist.le.rmin)))ite=1
c
c       ite ... for exclusive shells, minimum distance is in the interval
c       it ... all atoms are in the interval
c

        if((lexcl.and.ite.eq.0).or.(it.eq.0))then
         do 8 ii=1,n22
8        jt(i22(ii))=1
         if((nw+1.eq.nwmax).and.nwmax.ne.0)then
          write(6,*)' Maximum number of waters overflow'
         else
c         all atoms close, take them:
          do 9 ii=1,n22
          natt=natt+1
          qt(natt)=as(i22(ii))
          item(natt)=i22(ii)
          rtem(natt,1)=r(i22(ii),1)+dx
          rtem(natt,2)=r(i22(ii),2)+dy
9         rtem(natt,3)=r(i22(ii),3)+dz
          nw=nw+1
         endif
        endif
1010    continue
c       clonecloneclonecloneclonecloneclonecloneclonecloneclone
       endif

      endif
31    continue
c     31313131313131313131313131313131313131313131313131313131313131313131
c      put dead wood atoms back if they are solute:
      if(ldw.and.dwb.and.lsolulst)call pbw(natw,isw,natuo,isuo,natu,
     1isu,natt,qt,rtem,asx,rx,ntx)
      call wgne(ntx,s80,natt,natu,item,isu,qt,rtem,lrev,lsolulst)
      if(lsolulst)iw=nat-iuu
      write(6,6001)is,iw,nw,iuu
6001  format(i4,': ',i6,' solvent atoms,',i4,' selected molecules.',i4,
     1' solute atoms')
      if(nw.gt.inwmax)then
       inwmax=nw
       ix=is
      endif

      goto 1

999   close(3)
      close(33)

      write(6,6003)is,inwmax,ix
6003  format(/,i4,' structures, maximum of ',i3,' sols in number',i4)
      end

      subroutine readopt(cov,lsolulst,lper,lrev,lexcl,nwmax,lhyd,ldw,
     1px,py,pz,dwb)
      implicit none
      real*8 cov
      logical lsolulst,lhyd,lper,lrev,lexcl,lopt,ldw,dwb
      integer*4 nwmax
      real*8 px,py,pz
      character*50 st
c     covalent bond limit:
      cov=1.8d0
      lsolulst=.false.
      lper=.false.
      ldw=.false.
      dwb=.true.
      lrev=.false.
      lexcl=.true.
      nwmax=0
      lhyd=.false.
      px=0.0d0
      py=0.0d0
      pz=0.0d0
      inquire(file='XSHELL.OPT',exist=lopt)
      if(lopt)then
       write(6,*)'XSHELL.OPT found'
       open(44,file='XSHELL.OPT')
101    read(44,440,end=102,err=102)st
440    format(a50)
       if(st(1:3).eq.'LST')read(44,*)lsolulst
       if(st(1:3).eq.'EXC')read(44,*)lexcl
       if(st(1:3).eq.'COV')read(44,*)cov
       if(st(1:3).eq.'HYD')read(44,*)lhyd
       if(st(1:3).eq.'NWM')read(44,*)nwmax
       if(st(1:3).eq.'REV')read(44,*)lrev
       if(st(1:3).eq.'LDW')read(44,*)ldw
       if(st(1:3).eq.'DWB')read(44,*)dwb
       if(st(1:3).eq.'PER')then
        read(44,*)lper
        read(44,*)px,py,pz
       endif
       if(st(1:3).eq.'END')goto 102
       goto 101
102    close(44)
      endif

      if(lrev.and.lsolulst)
     1write(6,*)'Solvent forced to beginning of NEWFILE.X'
      if(lper)then
       write(6,*)'Periodic boundary conditions on'
      else
       write(6,*)'Periodic boundary conditions off'
      endif
      write(6,*)

      return
      end

      subroutine wgne(ntx,s80,natt,natu,item,isu,qt,rtem,lrev,
     1lsolulst)
      implicit none
      integer*4 ntx
      real*8 rtem(ntx,3)
      integer*4 natt,natu,ic,i,j,k,item(*),isu(*),qt(*)
      logical lrev,lsolulst
      character*80 s80
      write(33,300)s80
300   format(a80)
      write(33,*)natt
c     if required, gather solute atoms first:
      if(lrev.and.lsolulst)then
       do 325 i=1,natt
       do 325 j=1,natu
325    if(item(i).eq.isu(j))write(33,3333)qt(i),(rtem(i,k),k=1,3)
       do 326 i=1,natt
       ic=0
       do 327 j=1,natu
327    if(item(i).eq.isu(j))ic=1
326    if(ic.eq.0)write(33,3333)qt(i),(rtem(i,k),k=1,3)
      else
       do 324 i=1,natt
324    write(33,3333)qt(i),(rtem(i,k),k=1,3)
3333   format(i3,3f12.6)
      endif
      return
      end

      subroutine ral(f,n,isu)
      implicit none
      integer*4 n,isu(*),i
      character*(*) f
      open(9,file=f,status='old')
      read(9,*)n
      read(9,*)(isu(i),i=1,n)
      close(9)
      write(6,*)n,' atoms from '//f
      return
      end

      subroutine setatmax(natmax)
      implicit none
      integer*4 natmax,nat,i
      natmax=0
      open(9,file='FILE.X',status='old')
1     read(9,300,end=999,err=999)
300   format(a80)
      read(9,*,end=999,err=999)nat
      if(nat.gt.natmax)natmax=nat
      do 2 i=1,nat
2     read(9,*)
      goto 1
999   close(9)
      write(6,601)natmax
601   format(i6,': maximal number of atoms in FILE.X')
      return
      end

      subroutine as1(ntx,nat,natu,lsa,isu)
      implicit none
      integer*4 nat,i,ii,ntx,natu,isu(ntx)
      logical lsa(ntx)
      do 21 i=1,nat
      lsa(i)=.false.
      do 211 ii=1,natu
      if(isu(ii).eq.i)then
       lsa(i)=.true.
       goto 21
      endif
211   continue
21    continue
      return
      end

      subroutine as2(ntx,nat,natu,natw,isu,isw,lsolulst,as,r)
      implicit none
      integer*4 nat,i,ii,ntx,natu,isu(ntx),isw(ntx),natw,as(ntx),ix,
     1natnew,natun
      real*8 r(ntx,3)
      logical lsolulst
      integer*4,allocatable::asn(:),jnd(:),isun(:)
      real*8,allocatable::rn(:,:)
      allocate(rn(nat,3),asn(nat),jnd(nat),isun(nat))
      jnd=0
      natnew=0
      do 1 i=1,nat
      do 2 ii=1,natw
2     if(isw(ii).eq.i)goto 1
      natnew=natnew+1
      jnd(i)=natnew
      asn(natnew)=as(i)
      do 3 ix=1,3
3     rn(natnew,ix)=r(i,ix)
1     continue     
      write(6,*)natnew,' atoms left from',nat

      if(lsolulst)then

       natun=0
       do 5 i=1,natu
       if(jnd(isu(i)).ne.0)then
        natun=natun+1
        isun(natun)=jnd(isu(i))
       endif
5      continue

       natu=natun
       do 7 i=1,natu
 7     isu(i)=isun(i)

      endif

      do 6 i=1,natnew
      as(i)=asn(i)
      do 6 ix=1,3
6     r(i,ix)=rn(i,ix)

      nat=natnew


      return
      end

      subroutine pbw(natw,isw,natuo,isuo,natu,isu,natt,qt,rtem,asx,
     1rx,ntx)
      implicit none
      integer*4 natw,isw(*),natuo,isuo(*),natu,isu(*),natt,qt(*),
     1ntx,i,jj,ii,asx(*),ix,j
      real*8 rtem(ntx,3),rx(ntx,3)
      do 10 i=1,natw
      jj=isw(i)
      do 11 ii=1,natuo
      if(isuo(ii).eq.jj)then
       natu=natu+1
       isu(natu)=jj
       natt=natt+1
       do 12 j=natt+1,2,-1
       qt(j)=qt(j-1)
       do 12 ix=1,3
12     rtem(j,ix)=rtem(j-1,ix)
       qt(1)=asx(jj)
       do 13 ix=1,3
13     rtem(1,ix)=rx(jj,ix)
      endif
 11   continue
 10   continue
      return
      end