program stick implicit none character*3 key integer*4 nat,nats,i,j,nm,ib,ii,kk,k,jj,ik,ij,io,natd,i0,icomb, 1ib1,ib2 integer*4, allocatable::t(:),ts(:),m(:),icc(:),over(:),p(:), 1isc(:),td(:),item(:) real*8,allocatable::r(:),rs(:) real*8 cut,dist,dd,rms,rmsbest LOGICAL MTC,EVEN i0=692 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)) 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),icc(nat),over(nats),p(nats), 1isc(nats),item(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.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') do 6 i=1,nat if(ts(m(1)).eq.t(i))then c define initial permutation in order: DO 12 kk = 2, nm 12 p(kk-1) = kk-1 MTC = .FALSE. c try if atoms would be close enough: do 10 ii=1,nat 10 icc(ii)=0 do 13 ii=2,nm nn(ii)=0 13 isc(ii)=0 c nn .. number of close atoms at place ii ib=1 over(1)=i isc(1)=i icc(i)=m(1) do 8 k=1,nat do 8 ik=1,nm-1 kk=p(ik)+1 if(k.ne.i.and.t(k).eq.ts(m(kk)). 1 and.icc(k).eq.0.and.isc(kk).eq.0)then if(dabs(dist(r,i,k)-dist(rs,m(1),m(kk))).lt.cut)then ib=ib+1 if(i.eq.i0) 1 write(6,*)'ib=',ib, 1 'd=',dabs(dist(r,i,k)-dist(rs,m(1),m(kk))), 1 ' k kk:',k,kk if(ib.gt.nm)goto 6 over(kk)=k icc(k)=m(kk) isc(kk)=k endif endif 8 continue c if all atoms can be assigned, investigate all combinations: rmsbest=1.0d9 if(ib.eq.nm)then icomb=0 20 CONTINUE icomb=icomb+1 if(i.eq.i0)write(6,*)' Combination ',icomb CALL NEXPER (nm-1, p, MTC, EVEN) if(i.eq.i0)then WRITE (*,"(I6), ",advance="no") 0 DO 30 kk = 2, nm if(p(kk-1).eq.0)goto 6 30 WRITE (*,"(I6), ",advance="no") p(kk-1) PRINT * WRITE (*,"(I6), ",advance="no") i DO 31 kk = 2, nm 31 WRITE (*,"(I6), ",advance="no") isc(kk) PRINT * WRITE (*,"(I6), ",advance="no") m(1) DO 32 kk = 2, nm 32 WRITE (*,"(I6), ",advance="no") m(p(kk-1)+1) PRINT * endif c check that all distances match big atoms isc(1)....isc(nm-1): rms=0.0d0 do 11 ij=2,nm jj=p(ij-1)+1 ib1=isc(ij) do 11 ik=ij+1,nm ib2=isc(ik) kk=p(ik-1)+1 if(t(ib1).ne.ts(m(jj)).or.t(ib2).ne.ts(m(kk)))then if(i.eq.i0)then write(6,*)ib1,t(ib1),ib2,t(ib2),ts(m(jj)),ts(m(kk)) write(6,*)'impossible combination ' endif goto 61 endif dd=dabs(dist(r,ib1,ib2)-dist(rs,m(jj),m(kk))) if(dd.gt.cut)then if(i.eq.i0)then write(6,*)'combination failed' write(6,*)ib1,ib2,m(jj),m(kk),dd endif goto 61 endif rms=rms+dd 11 continue if(i.eq.i0)write(6,*)' deviation:',rms if(rms.lt.rmsbest)then rmsbest=rms do 33 ij=2,nm 33 item(ij)=p(ij-1) endif c if possible, try a next permutation: 61 IF (MTC) GOTO 20 endif if(rmsbest.lt.1.0d4)then write(6,'(a8)',advance="no")' Match: ' io=io+1 write(6,6003,advance="no")io 6003 format(i6) do 35 ii=1,nat 35 icc(ii)=0 icc(i)=m(1) do 34 ij=2,nm icc(isc(item(ij)+1))=m(item(ij)+1) 34 over(item(ij)+1)=isc(ij-1) do 9 j=1,ib td(over(j))=0 9 write(6,6003,advance="no")over(j) write(6,*) write(9,*)io, 'overlap' write(9,900)(ii,ii=1,nat) 900 format(20i6) write(9,900)(icc(ii),ii=1,nat) goto 6 endif endif 6 continue close(9) natd=0 do 14 i=1,nat 14 if(td(i).ne.0)natd=natd+1 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) 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 c ============================================================ SUBROUTINE NEXPER(N, A, MTC, EVEN) implicit none INTEGER*4 N,A(*),I, S, D, NM3, IA, I1, L, J, M LOGICAL MTC, EVEN IF (MTC) GOTO 10 NM3 = N-3 DO 1 I=1,N 1 A(I)=I MTC=.TRUE. EVEN=.TRUE. IF(N .EQ. 1) GOTO 8 6 IF(A(N) .NE. 1 .OR. A(1) .NE. 2+MOD(N,2)) RETURN IF(N .LE. 3) GOTO 8 DO 7 I=1,NM3 IF(A(I+1) .NE. A(I)+1) RETURN 7 CONTINUE 8 MTC=.FALSE. RETURN 10 IF(N .EQ. 1) GOTO 27 IF(.NOT. EVEN) GOTO 20 IA=A(1) A(1)=A(2) A(2)=IA EVEN=.FALSE. GOTO 6 20 S=0 DO 26 I1=2,N IA=A(I1) I=I1-1 D=0 DO 30 J=1,I 30 IF(A(J) .GT. IA) D=D+1 S=D+S IF(D .NE. I*MOD(S,2)) GOTO 35 26 CONTINUE 27 A(1)=0 GOTO 8 35 M=MOD(S+1,2)*(N+1) DO 40 J=1,I IF(ISIGN(1,A(J)-IA) .EQ. ISIGN(1,A(J)-M)) GOTO 40 M=A(J) L=J 40 CONTINUE A(L)=IA A(I1)=M EVEN=.TRUE. RETURN 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