program stick2 implicit none character*3 key integer*4 nat,nats,i,nm,ii,kk,k,io,natd,i0,icomb, 1ib1,ib2,mk,m1,ic,iex,LEXCL,N7,NEXC,NSL integer*4, allocatable::t(:),ts(:),m(:),icc(:),over(:), 1isc(:),td(:),inn(:,:),STATE(:),si(:),nn(:) real*8,allocatable::r(:),rs(:) real*8 cut,dist,dd,rms,dik i0=0 write(6,600) 600 format(' Automatic chromophore recognition',/,/, 1 ' Input: Big.X',/, 1 ' SMALL.X ("chromophore")',/, 1 ' STICK.OPT (options)',/) call rdx('BIG.X',nat,r,t,0) allocate(t(nat),r(3*nat),td(nat),icc(nat)) call rdx('BIG.X',nat,r,t,1) do 16 i=1,nat 16 td(i)=t(i) call rdx('SMALL.X',nats,rs,ts,0) allocate(ts(nat),rs(3*nats),m(nats),over(nats), 1isc(nats),inn(nats,nats),STATE(nats),nn(nats)) call rdx('SMALL.X',nats,rs,ts,1) nm=0 cut=0.5d0 open(9,file='STICK.OPT') 1 read(9,901,end=99,err=99)key 901 format(a3) if(key.eq.'CUT')read(9,*)cut if(key(1:2).eq.'I0')read(9,*)i0 if(key.eq.'MAP')then read(9,*)nm if(nm.gt.nats)call report('nm > nats') read(9,*)(m(i),i=1,nm) endif goto 1 99 close(9) if(nm.eq.0)then nm=nats do 3 i=1,nm 3 m(i)=i endif write(6,601)nat,nats,nm 601 format(i10,' atoms in BIG.X',/, 1 i10,' atoms in SMALL.X',/, 1 i10,' atoms in SMALL.X will be taken as chromophore:') do 7 i=1,nm 7 write(6,602)m(i) 602 format(i5,$) write(6,605)cut 605 format(/,' Distance cutoff:',f10.5,' A',/) io=0 open(9,file='CCT.TEM') m1=m(1) do 6 i=1,nat if(ts(m1).eq.t(i))then do 13 ii=2,nm 13 nn(ii)=0 c find all atoms close enough to each from small: do 8 k=1,nat if(k.ne.i)then dik=dist(r,i,k) do 81 kk=2,nm mk=m(kk) if(t(k).eq.ts(mk).and.(dabs(dik-dist(rs,m1,mk))).lt.cut)then if(i.eq.i0)write(6,*)'i k m1 mk ',i,k,m1,mk if(i.eq.i0)write(6,*)'d=',dabs(dik-dist(rs,m1,mk)) nn(kk)=nn(kk)+1 inn(kk,nn(kk))=k endif 81 continue endif 8 continue c check if all centers occupied do 36 kk=2,nm 36 if(nn(kk).eq.0)goto 6 if(i.eq.i0)write(6,*)'possible candidate',m1,i c extra atoms on each center, "excitations" LEXCL=0 do 37 kk=2,nm 37 LEXCL=LEXCL+nn(kk)-1 if(i.eq.i0)write(6,*)'LEXCL ',LEXCL allocate(si(LEXCL)) c generate all combinations c states Nexc "excited"=occupied icomb=0 N7=2 NSL=nm c comb comb comb comb comb comb comb comb comb comb do 103 Nexc=1,LEXCL c distribute Nexc excitations upon centers N7..NSL=2..nm c initial indices - all to the first center: do 141 iex=1,Nexc 141 si(iex)=N7 105 DO 1966 kk=N7,NSL 1966 STATE(kk)=0 do 104 iex=1,Nexc 104 state(si(iex))=state(si(iex))+1 c check that the state is allowed: do 111 kk=N7,NSL 111 if(state(kk).gt.nn(kk)-1)goto 112 do 40 ii=1,nat 40 icc(ii)=0 do 39 kk=2,nm over(kk)=0 39 isc(kk)=0 isc(1)=i over(1)=i icc(i)=m(1) do 38 kk=2,nm icc(inn(kk,state(kk)+1))=m(kk) over(kk)=inn(kk,state(kk)+1) 38 isc(kk)=inn(kk,state(kk)+1) c check that all atoms unique do 41 ii=1,nm do 41 kk=ii+1,nm 41 if(over(kk).eq.over(ii))goto 112 c c investigate this combination icomb=icomb+1 if(i.eq.i0)then write(6,*)' Combination ',icomb DO 31 kk = 1, nm 31 WRITE (*,"(I6) ",advance="no") isc(kk) PRINT * DO 32 kk = 1, nm 32 WRITE (*,"(I6) ",advance="no") m(kk) PRINT * endif c check that all distances match big atoms isc(1)....isc(nm-1): rms=0.0d0 do 11 ii=2,nm ib1=isc(ii) do 11 kk=ii+1,nm ib2=isc(kk) dd=dabs(dist(r,ib1,ib2)-dist(rs,m(ii),m(kk))) if(dd.gt.cut)then if(i.eq.i0)then write(6,*)'combination failed' write(6,*)ib1,ib2,m(ii),m(kk),dd endif goto 112 endif 11 rms=rms+dd if(i.eq.i0)write(6,*)' deviation:',rms write(6,'(a8)',advance="no")' Match: ' io=io+1 write(6,6003,advance="no")io 6003 format(i6) do 9 kk=1,nm td(over(kk))=0 9 write(6,6003,advance="no")over(kk) write(6,*) write(9,*)io, 'overlap' write(9,900)(ii,ii=1,nat) 900 format(20i6) write(9,900)(icc(ii),ii=1,nat) deallocate(si) goto 6 c c find index to be changed 112 do 108 ic=Nexc,1,-1 108 if(si(ic).lt.NSL)goto 109 goto 103 109 do 110 ii=ic+1,Nexc 110 si(ii)=si(ic)+1 si(ic)=si(ic)+1 c goto 105 c 103 continue c comb comb comb comb comb comb comb comb comb comb deallocate(si) endif 6 continue close(9) write(6,*)'CCT.TEM written' natd=0 do 14 i=1,nat 14 if(td(i).ne.0)natd=natd+1 write(6,*)nat,natd open(9,file='t.x') write(9,*)natd write(9,*)natd do 15 i=1,nat 15 if(td(i).ne.0)write(9,9009)td(i),(r(ii+3*(i-1)),ii=1,3) 9009 format(i6,3f12.6) close(9) write(6,*)'t.x control sritten' stop end c ============================================================ subroutine rdx(f,n,r,t,ic) implicit none character*(*) f integer*4 n,t(*),i,ix,ic real*8 r(*) open(9,file=f,status='old') read(9,*) read(9,*)n if(ic.ne.0)then do 1 i=1,n 1 read(9,*)t(i),(r(ix+3*(i-1)),ix=1,3) endif close(9) return end c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end function dist(r,i,k) implicit none real*8 r(*),dist integer*4 i,k dist=dsqrt((r(3*i-2)-r(3*k-2))**2 1+ (r(3*i-1)-r(3*k-1))**2+(r(3*i)-r(3*k))**2) return end