program gcomp c c reads a set of (MD) geometries and either compares to c a library set or generates basis by a factor analysis c implicit none integer*4 ib,nb0,is,ia,nat0,iwr,np0 parameter (nb0=1000,nat0=800,np0=4000) real*8 rb(3,nat0,nb0),rs(3,nat0,nb0),sum,d0,dist,cut, 1so(nb0,nb0),alpha,e(nb0,2*nb0),TOL,ai(nb0,nb0),vo(nb0), 1ao(nb0,nb0),sy(nb0,np0),sr(np0),sx(np0),x,y integer*4 s1b(nat0,nb0),s1s(nat0,nb0),IERR,i, 1iap,ic,ii,ip,ix,job,jx,n2,np(nb0) integer*4 nab(nb0),nas(nb0),ind(nat0),nn,indi(nat0) 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*5 s5 character*80 fn c write(6,7000) 7000 format('Reads a set of (MD) geometries and either compares to',/, 1'a library set or generates basis by a factor analysis.',/,/, 2' Input: BIG.LST list of geometry files [.tinker form]',/, 4' SMALL.LST list of library structures [.x]',/, 3' SPEC.LST list of spectra [.prn]',/, 3' GCOMP.OPT list of options',/, 5' CCT.INP atom map',/,/, 6' Jobs 1: make spectrum for each .tinker geom and everage',/, 7' 2: find basis in .tinker geoms',/,/) job=1 cut=2.6d0 alpha=1.0d1 iwr=0 open(8,file='GCOMP.OPT') 1001 read(8,800,end=888,err=888)s5 800 format(a5) if(s5(1:3).eq.'JOB')read(8,*)job if(s5(1:3).eq.'CUT')read(8,*)cut if(s5(1:3).eq.'IWR')read(8,*)iwr if(s5(1:3).eq.'ALP')read(8,*)alpha goto 1001 888 close(8) write(6,609)job,iwr,cut,alpha 609 format(' Job type ',i10,/, 1 ' Print level ',i10,/, 1 ' Cuttoff distance ',f10.3,/, 1 ' Alpha ',f10.3) c c load geometries: call RL('BIG.LST',ib,s1b,rb,nab,nb0,nat0,0) call RL('SMALL.LST',is,s1s,rs,nas,nb0,nat0,1) do 2 i=1,nat0 ind(i)=0 2 indi(i)=0 open(2,file='CCT.INP') READ(2,*) READ(2,*) READ(2,2001)(IND(ia),ia=1,nab(ib)) READ(2,2001)(IND(ia),ia=1,nab(ib)) do 21 i=1,nab(ib) 21 if(IND(i).ne.0)indi(IND(i))=IND(i) 2001 FORMAT(20I3) close(2) write(6,*)'Map big:' write(6,2001)(ia,ia=1,nab(ib)) write(6,2001)(IND(ia),ia=1,nab(ib)) write(6,*)'Map small:' write(6,2001)(ia,ia=1,nas(is)) write(6,2001)(INDI(ia),ia=1,nas(is)) c c delete distant solvent atoms: call CL(ib,s1b,rb,nab,nb0,nat0,ind,cut,iwr) call CL(is,s1s,rs,nas,nb0,nat0,indi,cut,iwr) c c subtract centers of mass: call SC(ib,rb,nab,nb0,nat0,ind) call SC(is,rs,nas,nb0,nat0,indi) open(2,file='SPEC.LST',status='OLD') i=0 300 read(2,3001,end=333,err=333)fn 3001 format(a80) i=i+1 open(22,file=fn,status='OLD') ix=0 100 read(22,*,end=200,err=200)x,y ix=ix+1 if(ix.gt.np0)call report('too many spectral points') sx(ix)=x sy(i,ix)=y np(i)=ix goto 100 200 close(22) write(6,*)np(i) goto 300 333 close(2) write(6,*)i,' spectra read' if(job.eq.1)then c c correlation between basis c center of transfered atoms only write(6,*)'Subgeometry overlaps:' do 19 ii=1,is n2=0 do 17 ia=1,nas(ii) if(indi(ia).ne.0)then n2=n2+1 do 18 ix=1,3 18 XB(n2,ix)=rs(ix,ia,ii) endif 17 continue do 19 ip=1,is nn=0 do 171 ia=1,nas(ip) if(indi(ia).ne.0)then nn=nn+1 do 22 ix=1,3 22 XS(nn,ix)=rs(ix,ia,ip) endif 171 continue if(nn.ne.n2)call report('nn<>n2') do 20 ix=1,3 do 20 jx=1,3 20 A(ix,jx)=0.0d0 NAT=nn call DOMATRIX(IERR) if(iwr.gt.0)write(6,6001)A 6001 format(' Transformation matrix:',/,3(3F10.4,/)) c c scalar product: sum=0.0d0 do 51 ia=1,nn XT(ia,1)=A(1,1)*XS(ia,1)+A(1,2)*XS(ia,2)+A(1,3)*XS(ia,3) XT(ia,2)=A(2,1)*XS(ia,1)+A(2,2)*XS(ia,2)+A(2,3)*XS(ia,3) XT(ia,3)=A(3,1)*XS(ia,1)+A(3,2)*XS(ia,2)+A(3,3)*XS(ia,3) do 512 ix=1,3 512 sum=sum+XT(ia,ix)*XB(ia,ix) 51 if(iwr.gt.0)write(6,6009)(XB(ia,ix),ix=1,3),(xt(ia,ix),ix=1,3) 6009 format(6f12.5) 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 c scalar product for non-overlapped atoms: do 191 ia=1,nas(ii) if(indi(ia).eq.0)then c find closest atom of the same kind,non-overlapped d0=0.0d0 ic=0 do 192 iap=1,nas(ip) if(indi(iap).eq.0.and.s1s(ia,ii).eq.s1s(iap,ip))then dist=0.0d0 do 1922 ix=1,3 1922 dist=dist+(rs(ix,ia,ii)-XT(iap,ix))**2 if(ic.eq.0)then d0=dist ic=iap else if(dist.lt.d0)then ic=iap d0=dist endif endif endif 192 continue if(ic.gt.0)then do 1921 ix=1,3 1921 sum=sum+rs(ix,ia,ii)*XT(ic,ix)*0.5d0 endif endif 191 continue 19 so(ip,ii)=sum do 1911 ii=1,is do 1911 ip=ii+1,is so(ip,ii)=0.50d0*(so(ip,ii)+so(ii,ip)) 1911 so(ii,ip)=so(ip,ii) write(6,*)'Sub-geoms correlation matrix' do 198 ii=1,is 198 write(6,6003)(so(ii,ip),ip=1,is) 6003 format(10f8.2) do 193 ii=1,is 193 so(ii,ii)=so(ii,ii)+alpha if(is.gt.nb0+1)call report('is>nb0+1') do 194 ii=1,is so(is+1,ii)=1.0d0 194 so(ii,is+1)=1.0d0 so(is+1,is+1)=0.0d0 TOL=1.0d-10 call INV(nb0,so,ai,is+1,TOL,e,IERR) if(IERR.ne.0)call report('unsuccessful inversion') do 1 i=1,ib write(6,*)'Decomposing geometry ',i c c record big coordinates nn=0 do 13 ia=1,nab(i) if(ind(ia).ne.0)then nn=nn+1 if(nn.gt.nat0)call report('nn > nat0!') do 14 ix=1,3 14 XB(nn,ix)=rb(ix,ia,i) endif 13 continue c c record small coordinates c center of transfered atoms only if(iwr.gt.0)write(6,*)'Subgeometry overlaps:' do 27 ii=1,is n2=0 do 24 ia=1,nas(ii) if(indi(ia).ne.0)then n2=n2+1 do 23 ix=1,3 23 XS(n2,ix)=rs(ix,ia,ii) endif 24 continue if(n2.ne.nn)call report('n2<>nn') do 25 ix=1,3 do 25 jx=1,3 25 A(ix,jx)=0.0d0 CALL DOMATRIX(IERR) if(iwr.gt.0)write(6,6001)A sum=0.0d0 do 26 ia=1,nn XT(ia,1)=A(1,1)*XS(ia,1)+A(1,2)*XS(ia,2)+A(1,3)*XS(ia,3) XT(ia,2)=A(2,1)*XS(ia,1)+A(2,2)*XS(ia,2)+A(2,3)*XS(ia,3) XT(ia,3)=A(3,1)*XS(ia,1)+A(3,2)*XS(ia,2)+A(3,3)*XS(ia,3) do 521 ix=1,3 521 sum=sum+XT(ia,ix)*XB(ia,ix) 26 if(iwr.gt.0)write(6,6009)(xb(ia,ix),ix=1,3),(xt(ia,ix),ix=1,3) c do 261 ia=1,nas(ii) XT(ia,1)=A(1,1)*rs(1,ia,ii)+A(1,2)*rs(2,ia,ii)+A(1,3)*rs(3,ia,ii) XT(ia,2)=A(2,1)*rs(1,ia,ii)+A(2,2)*rs(2,ia,ii)+A(2,3)*rs(3,ia,ii) 261 XT(ia,3)=A(3,1)*rs(1,ia,ii)+A(3,2)*rs(2,ia,ii)+A(3,3)*rs(3,ia,ii) c c scalar product for non-overlapped atoms: i-> ii do 195 ia=1,nab(i) if(ind(ia).eq.0)then c find closest atom of the same kind,non-overlapped d0=0.0d0 ic=0 do 199 iap=1,nas(ii) if(indi(iap).eq.0.and.s1b(ia,i).eq.s1s(iap,ii))then dist=0.0d0 do 1923 ix=1,3 1923 dist=dist+(rb(ix,ia,i)-XT(iap,ix))**2 if(ic.eq.0)then d0=dist ic=iap else if(dist.lt.d0)then ic=iap d0=dist endif endif endif 199 continue if(ic.gt.0)then do 1926 ix=1,3 1926 sum=sum+rb(ix,ia,i)*XT(ic,ix)*0.5d0 endif endif 195 continue c c scalar product for non-overlapped atoms: ii-> i do 196 ia=1,nas(ii) if(indi(ia).eq.0)then c find closest atom of the same kind,non-overlapped d0=0.0d0 ic=0 do 197 iap=1,nab(i) if(ind(iap).eq.0.and.s1s(ia,ii).eq.s1b(iap,i))then dist=0.0d0 do 1924 ix=1,3 1924 dist=dist+(rb(ix,iap,i)-XT(ia,ix))**2 if(ic.eq.0)then d0=dist ic=iap else if(dist.lt.d0)then ic=iap d0=dist endif endif endif 197 continue if(ic.gt.0)then do 1925 ix=1,3 1925 sum=sum+rb(ix,ic,i)*XT(ia,ix)*0.5d0 endif endif 196 continue 27 vo(ii)=sum+alpha/dble(nas(is)) vo(is+1)=1.0d0 do 28 ix=1,is ao(ix,i)=0.0d0 do 28 ii=1,is+1 28 ao(ix,i)=ao(ix,i)+ai(ix,ii)*vo(ii) write(6,*)' Coefficients:' write(6,678)(ix,ix=1,is) 678 format(8i10) write(6,679)(ao(ix,i),ix=1,is) 679 format(8f10.3) 1 continue do 101 ix=1,np(1) sr(ix)=0.0d0 do 101 i=1,ib do 101 ii=1,is 101 sr(ix)=sr(ix)+ao(ii,i)*sy(ii,ix)/dble(ib) open(22,file='AV.PRN') do 102 ix=1,np(1) 102 write(22,2222)sx(ix),sr(ix) 2222 format(2g20.5) close(22) write(6,*)' AV.PRN written' endif if(job.ne.1)call report('not implemented yet') 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=800) 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=800) 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) implicit none integer*4 i,ia,ix,ic,nb0,nat0 c ic=0 ... tinker, 1 ... MCM integer*4 s(nat0,nb0),n(nb0) character*1 s1 character*(*) f1 character*80 fn real*8 r(3,nat0,nb0) open(2,file=f1,status='OLD') i=0 3 read(2,300,end=333,err=333)fn 300 format(a80) i=i+1 if(i.gt.nb0)call report('too many structures') open(21,file=fn,status='OLD') if(ic.eq.0)then read(21,*)n(i) if(n(i).gt.nat0)call report('too many atoms') do 22 ia=1,n(i) read(21,2121)s1,(r(ix,ia,i),ix=1,3) 2121 format(6x,2x,a1,2x,3f12.6,i6) s(ia,i)=0 if(s1.eq.'C')s(ia,i)=6 if(s1.eq.'N')s(ia,i)=7 if(s1.eq.'O')s(ia,i)=8 22 if(s1.eq.'H')s(ia,i)=1 if(s(ia,i).eq.0)write(6,*)'Unknown atoms symbol ',s1 else read(21,*) read(21,*)n(i) if(n(i).gt.nat0)call report('too many atoms') do 221 ia=1,n(i) 221 read(21,*)s(ia,i),(r(ix,ia,i),ix=1,3) endif close(21) goto 3 333 close(2) c write(6,*)i,' geometries in ',f1 return end subroutine CL(ns,s,r,n,nb0,nat0,ind,cut,iwr) implicit none integer*4 i,ia,ix,ns,ic,nb0,nat0,ind(*),iap,iwr c ic=0 ... tinker, 1 ... MCM integer*4 s(nat0,nb0),n(nb0) real*8 r(3,nat0,nb0),cut,dist,distmin c do 1 i=1,ns do 2 ia=1,n(i) if(ind(ia).eq.0.and.s(ia,i).eq.8)then ic=0 c find minimal distance from the solute distmin=0.0d0 do 21 iap=1,n(i) if(ind(iap).ne.0)then c this is solute atom dist=0.0d0 do 22 ix=1,3 22 dist=dist+(r(ix,ia,i)-r(ix,iap,i))**2 ic=ic+1 if(ic.eq.1)then distmin=dist else if(dist.lt.distmin)distmin=dist endif endif 21 continue if(distmin.gt.cut**2)then c this oxygen is far from solute,delete it with bond hydrogens s(ia,i)=0 do 24 iap=1,n(i) if(s(iap,i).eq.1)then dist=0.0d0 do 25 ix=1,3 25 dist=dist+(r(ix,ia,i)-r(ix,iap,i))**2 if(dist.lt.1.2d0**2)s(iap,i)=0 endif 24 continue endif endif 2 continue ic=0 ia=0 3 ia=ia+1 if(ia.le.n(i))then if(s(ia,i).eq.0)then do 4 iap=ia,n(i)-1 s(iap,i)=s(iap+1,i) do 4 ix=1,3 4 r(ix,iap,i)=r(ix,iap+1,i) n(i)=n(i)-1 ic=ic+1 endif goto 3 endif if(iwr.gt.0)write(6,600)i,ic 600 format('cluster ',i6,', ',i6,' atoms deleted') 1 continue write(6,*)'far atoms deleted' 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: 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) do 18 ia=1,n(ii) do 18 ix=1,3 18 r(ix,ia,ii)=r(ix,ia,ii)-cb(ix) 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 SUBROUTINE INV(nr,a,ai,n,TOL,e,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(nr,nr),a(nr,nr),e(nr,2*nr) C 10000 IERR=0 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 IERR=1 write(6,*)ii write(6,*)'tol = ',tol tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 RETURN c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END