program goverlap implicit none integer*4 nb0,is,ia,nat0,iwr,i,ip,ix,n2,IERR,jx,ig,ic,ibt, 1ias,iae parameter (nb0=14000,nat0=500) real*8 rs(3,nat0,nb0) integer*4 nas(nb0),ind(nat0),nn,s1s(nat0,nb0),bt(7,nat0,nb0) real*8 A,XS,XB,RTOL,XT(nat0,3) integer*4 NAT,ITER COMMON/MATVAR/A(3,3),XS(nat0,3),XB(nat0,3),RTOL,NAT,ITER character*10 s10 c write(6,7000) 7000 format('Reads a set of (MD) geometries and ', 1'makes best overlap to the first.',/,/, 2' Input: GEO.LST list of geometry files ',/, 5' CCT.INP atom map',/,/) c ********************* ibt=0 ic=1 c ********************* write(6,6006) 6006 format(60(1h*)) if(ibt.eq.0)then write(6,*)'version without bondtable' else write(6,*)'version with bondtable' endif if(ic.eq.0)then write(6,*)'for Tinker' else write(6,*)'for MCM' endif write(6,6006) iwr=0 call RL('GEO.LST',is,s1s,rs,nas,nb0,nat0,ic,ig,bt,ibt) do 2 i=1,nat0 2 ind(i)=0 open(2,file='CCT.INP') READ(2,*) READ(2,*) ias=1 23 iae=min(ias+19,nas(1)) READ(2,*)(IND(ia),ia=ias,iae) ias=ias+20 if(ias.le.nas(1))goto 23 ias=1 24 iae=min(ias+19,nas(1)) READ(2,*)(IND(ia),ia=ias,iae) ias=ias+20 if(ias.le.nas(1))goto 24 close(2) c c subtract centers of mass: call SC(ig,rs,nas,nb0,nat0,ind) c c take first as BIG: n2=0 do 17 ia=1,nas(1) if(ind(ia).ne.0)then n2=n2+1 do 18 ix=1,3 18 XB(n2,ix)=rs(ix,ia,1) endif 17 continue c c overlap with all: c write one big file, if this was in the beginning (is<>ig), otherwise c separate files: open(120,file='AMAT.LST') write(120,*)'Transformation matrices:' do 19 ip=1,ig nn=0 do 171 ia=1,nas(ip) if(ind(ia).ne.0)then nn=nn+1 do 22 ix=1,3 22 XS(nn,ix)=rs(ix,ia,ip) endif 171 continue do 20 ix=1,3 do 20 jx=1,3 20 A(ix,jx)=0.0d0 NAT=nn call DOMATRIX(IERR) write(120,*)ip write(120,6001)A 6001 format(3(3F10.4,/)) c do 511 ia=1,nas(ip) XT(ia,1)=A(1,1)*rs(1,ia,ip)+A(1,2)*rs(2,ia,ip)+A(1,3)*rs(3,ia,ip) XT(ia,2)=A(2,1)*rs(1,ia,ip)+A(2,2)*rs(2,ia,ip)+A(2,3)*rs(3,ia,ip) 511 XT(ia,3)=A(3,1)*rs(1,ia,ip)+A(3,2)*rs(2,ia,ip)+A(3,3)*rs(3,ia,ip) c write(s10,1010)ip 1010 format(i10) if(is.ne.ig)then if(ip.eq.1)open(90,file='all.x') if(iwr.gt.0)write(6,*)'all.x' else open(90,file=s10//'.x') if(iwr.gt.0)write(6,*)s10//'.x' endif write(90,*) write(90,*)nas(ip) do 512 ia=1,nas(ip) 512 write(90,9000)s1s(ia,ip),(XT(ia,ix),ix=1,3), 1(bt(ix,ia,ip),ix=1,7) 9000 format(i4,3f15.6,7i4,' 0.0') if(is.ne.ig)then if(ip.eq.ig)close(90) else close(90) endif 19 continue close(120) stop end SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) PARAMETER (nat0=500) COMMON/MATVAR/A(3,3),XS(nat0,3),XB(nat0,3),RTOL,NAT,ITER IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=0.0000001d0 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (nat0=500) COMMON/MATVAR/A(3,3),XS(nat0,3),XB(nat0,3),RTOL,NAT,ITER DIMENSION ANG(3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 R=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 R=R+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=R RETURN END subroutine RL(f1,i,s,r,n,nb0,nat0,ic,ig,bt,ibt) c ic=0 ... tinker, 1 ... MCM implicit none integer*4 i,ia,ix,ic,nb0,nat0,ig,nnext,ibt integer*4 s(nat0,nb0),n(nb0),bt(7,nat0,nb0) character*1 s1 character*(*) f1 character*80 fn real*8 r(3,nat0,nb0) open(2,file=f1,status='OLD') i=0 ig=0 3 read(2,300,end=333,err=333)fn 300 format(a80) i=i+1 write(6,*)fn open(21,file=fn,status='OLD') if(ic.eq.0)then c tinker ig=ig+1 read(21,*)n(ig) if(ig.gt.nb0)call report('too many structures') if(n(ig).gt.nat0)call report('too many atoms') do 22 ia=1,n(ig) read(21,2121)s1,(r(ix,ia,ig),ix=1,3) 2121 format(6x,2x,a1,2x,3f12.6,i6) do 4 ix=1,7 4 bt(ix,ia,ig)=0 s(ia,ig)=0 if(s1.eq.'C')s(ia,ig)=6 if(s1.eq.'N')s(ia,ig)=7 if(s1.eq.'O')s(ia,ig)=8 22 if(s1.eq.'H')s(ia,ig)=1 if(s(ia,ig).eq.0)write(6,*)'Unknown atoms symbol ',s1 else 88 read(21,*,err=99,end=99) read(21,*,err=99,end=99)nnext ig=ig+1 n(ig)=nnext if(ig.gt.nb0)call report('too many structures') if(n(ig).gt.nat0)call report('too many atoms') do 221 ia=1,n(ig) if(ibt.eq.1)then read(21,*)s(ia,ig),(r(ix,ia,ig),ix=1,3), 1 (bt(ix,ia,ig),ix=1,7) else read(21,*)s(ia,ig),(r(ix,ia,ig),ix=1,3) do 41 ix=1,7 41 bt(ix,ia,ig)=0 endif 221 continue goto 88 endif 99 close(21) goto 3 333 close(2) c write(6,*)i,' files in ',f1 write(6,*)ig,' geometries ' return end subroutine SC(i,r,n,nb0,nat0,ind) implicit none integer*4 i,ii,ia,ix,ind(*),nb0,nat0,n(nb0),n2 real*8 r(3,nat0,nb0),cb(3) c c subtract centers: open(123,file='CM.LST') write(123,1231) 1231 format('Centers of mass:') do 18 ii=1,i do 15 ix=1,3 cb(ix)=0.0d0 n2=0 do 16 ia=1,n(ii) if(ind(ia).ne.0)then cb(ix)=cb(ix)+r(ix,ia,ii) n2=n2+1 endif 16 continue 15 cb(ix)=cb(ix)/dble(n2) write(123,1232)ii,cb(1),cb(2),cb(3) 1232 format(i5,3f15.6) do 18 ia=1,n(ii) do 18 ix=1,3 18 r(ix,ia,ii)=r(ix,ia,ii)-cb(ix) close(123) c write(6,*)i,' centers subtracted' return end SUBROUTINE REPORT(STRING) CHARACTER*(*) STRING WRITE(6,1000) 1000 FORMAT(80(1H*)) WRITE(6,*)STRING WRITE(6,1000) WRITE(6,2000) 2000 FORMAT(' Program terminated') STOP END