program makegrad implicit none integer*4 N,i,il,im,ih,iargc,is,ix,j,nat,io,ioo,oo,ia,ibest, 1ierr,ii,jj,ja,joo,jx,ifound,i1 logical letter character*80 s80 real*8 d0,x,y,z,q,di,pbest,e,Esum logical lff,lwr character*1 ,allocatable::le(:) integer*4 ,allocatable::ll(:),lm(:),lh(:),no(:),lo(:,:) real*8 ,allocatable::g(:),r(:),cm(:,:),p(:,:),pp(:,:,:), 1f(:,:),fo(:,:) write(6,6009) 6009 format(' makegrad: joining ONIOM gradients/FFs ',/,/, 1 ' Input: FILE.X ',/, 1 ' */G.OUT Gaussian outputs ',/, 1 ' MAKEGRAD.OPT options ',/,/, 1 ' Output: FILE.GR gradient ',/) if(iargc().ne.1)call report(' Usage: makegrad N') call getarg(1,s80) read(s80,*)N call readopt(lff,lwr) open(8,file='FILE.X',status='old') read(8,*) read(8,*)nat close(8) allocate(le(nat),ll(nat),lm(nat),lh(nat),g(3*nat),no(nat), 1r(3*nat),lo(nat,N),cm(N,3),p(nat,N)) c for each atom, get no geometries that contain it at high level c lo ... list of the geometry umbers c mass centers of high level clusters c g ... gradient no=0 cm=0.0d0 g=0.0d0 d0=0.5d0 p=0.0d0 c first reading .... coordinates do 1 io=1,N write(s80,80)io 80 format(i80) do 2 is=1,len(s80) 2 if(s80(is:is).ne.' ')goto 3 3 open(8,file=s80(is:len(s80))//'/G.OUT') write(6,*)s80(is:len(s80))//'/G.OUT' 6 read(8,800,end=1,err=1)s80 800 format(a80) if(s80(2:19).eq.'Symbolic Z-matrix:')then write(6,*)s80(2:18) 111 read(8,800,end=1,err=1)s80 if(s80(2:9).eq.'Charge =')goto 111 backspace 8 ih=0 im=0 il=0 do 4 i=1,nat read(8,800)s80 c read geometry only for io=1, it is the same for all: if(io.eq.1)read(s80(6:len(s80)),*)(r(ix+3*(i-1)),ix=1,3) le(i)='X' do 5 j=1,len(s80) if(letter(s80(j:j)))then le(i)=s80(j:j) goto 51 endif 5 continue 51 if(le(i).ne.'L'.and.le(i).ne.'H'.and.le(i).ne.'M')then write(6,601)i 601 format(' Atom',i6,' type not defined, counted as high') ih=ih+1 lh(ih)=i no(i)=no(i)+1 lo(i,no(i))=io endif if(le(i).eq.'L')then il=il+1 ll(il)=i endif if(le(i).eq.'M')then im=im+1 lm(im)=i endif if(le(i).eq.'H')then ih=ih+1 lh(ih)=i no(i)=no(i)+1 lo(i,no(i))=io endif 4 continue if(io.eq.1)then if(il.ne.0)then write(6,602)il,' low-level atoms ' 602 format(i6,A20) if(lwr)write(6,603)(ll(i),i=1,il) endif if(im.ne.0)then write(6,602)im,' medium-level atoms ' if(lwr)write(6,603)(lm(i),i=1,im) endif endif if(ih.ne.0)then if(io.eq.1)write(6,602)ih,' high-level atoms: ' if(lwr.and.io.eq.1)write(6,603)(lh(i),i=1,ih) 603 format(12i6) do 7 ix=1,3 do 71 i=1,ih 71 cm(io,ix)=cm(io,ix)+r(ix+3*(lh(i)-1)) 7 cm(io,ix)=cm(io,ix)/dble(ih) write(6,605)(cm(io,ix),ix=1,3) 605 format(' Mass center:',3f12.4,' A') endif if(im+il+ih.ne.nat)call report('Atoms missed') goto 1 endif goto 6 1 close(8) c determine probabilities of each atom in each overlap geometry if(lwr)write(6,608) 608 format(/,' Atom best ov. overlap/probability(%):') ifound=0 do 9 ia=1,nat if(lwr)write(6,607)ia 607 format(i6,$) x=r(1+3*(ia-1)) y=r(2+3*(ia-1)) z=r(3+3*(ia-1)) q=0.0d0 do 10 ioo=1,no(ia) oo=lo(ia,ioo) c distance from the mass center di=dsqrt((x-cm(oo,1))**2+(y-cm(oo,2))**2+(z-cm(oo,3))**2) p(ia,oo)=1.0d0/(d0+di) 10 q=q+p(ia,oo) ibest=0 pbest=0.0d0 do 11 ioo=1,no(ia) oo=lo(ia,ioo) p(ia,oo)=p(ia,oo)/q if(p(ia,oo).gt.pbest)then pbest=p(ia,oo) ibest=oo endif 11 continue if(no(ia).eq.0)then if(lwr)write(6,617) 617 format(' not found') else ifound=ifound+1 if(lwr)then write(6,606)ibest 606 format(i6,$) do 12 ioo=1,no(ia) oo=lo(ia,ioo) 12 write(6,609)oo,p(ia,oo)*100.0d0 609 format(i4,'/',f4.0,$) write(6,*) endif endif 9 continue write(6,618)ifound,'atoms',nat 618 format(i10,' ',a5,' of ',i10,' found') c second reading ... energy and gradient Esum=0.0d0 do 101 io=1,N write(s80,80)io do 201 is=1,len(s80) 201 if(s80(is:is).ne.' ')goto 301 301 open(8,file=s80(is:len(s80))//'/G.OUT') do 14 i=is,len(s80) 14 write(6,610)s80(i:i) 610 format(A1,$) write(6,611)'/G.OUT' 611 format(A6,$) call addgrad(g,p,nat,N,io,Esum) 101 close(8) call svgr(g,nat,Esum) c determine probabilities of each atom pair in each overlap geometry if(lff)then allocate(pp(nat,nat,N),f(3*nat,3*nat),fo(3*nat,3*nat)) pp=0.0d0 f=0.0d0 if(lwr)write(6,612) 612 format(/,' Pair best ov. overlap/probability(%):',/, 1 ' (only found pair listed)') ifound=0 do 15 ia=1,nat do 15 ja=ia,nat x=(r(1+3*(ia-1))+r(1+3*(ja-1)))/2.0d0 y=(r(2+3*(ia-1))+r(2+3*(ja-1)))/2.0d0 z=(r(3+3*(ia-1))+r(3+3*(ja-1)))/2.0d0 q=0.0d0 do 16 ioo=1,no(ia) oo=lo(ia,ioo) c does this overlap contain also ja? do 16 joo=1,no(ja) if(lo(ja,joo).eq.oo)then c distance of ij from the mass center di=dsqrt((x-cm(oo,1))**2+(y-cm(oo,2))**2+(z-cm(oo,3))**2) pp(ia,ja,oo)=1.0d0/(d0+di) q=q+pp(ia,ja,oo) endif 16 continue ibest=0 pbest=0.0d0 do 19 ioo=1,no(ia) oo=lo(ia,ioo) do 19 joo=1,no(ja) if(lo(ja,joo).eq.oo)then pp(ia,ja,oo)=pp(ia,ja,oo)/q if(pp(ia,ja,oo).gt.pbest)then pbest=pp(ia,ja,oo) ibest=oo endif endif 19 continue if(ibest.ne.0)then ifound=ifound+1 if(lwr)then write(6,613)ia,ja 613 format(2i6,$) write(6,606)ibest do 20 ioo=1,no(ia) oo=lo(ia,ioo) do 20 joo=1,no(ja) 20 if(lo(ja,joo).eq.oo)write(6,609)oo,pp(ia,ja,oo)*100.0d0 write(6,*) endif endif 15 continue write(6,618)ifound,'pairs',(nat*(nat+1))/2 c third reading ... force field do 501 io=1,N write(s80,80)io do 401 is=1,len(s80) 401 if(s80(is:is).ne.' ')goto 661 661 open(7,file=s80(is:len(s80))//'/G.OUT') do 21 i=is,len(s80) 21 write(6,610)s80(i:i) write(6,611)'/G.OUT' call gar9f(3*nat,fo,ierr) if(ierr.ne.0)call report('FF not found') do 22 i=1,nat do 22 j=i,nat if(pp(i,j,io).gt.1.0d-10)then do 23 ix=1,3 ii=ix+3*(i-1) if(i.eq.j)then i1=ix else i1=1 endif do 23 jx=i1,3 jj=jx+3*(j-1) f(ii,jj)=f(ii,jj)+pp(i,j,io)*fo(ii,jj) 23 f(jj,ii)=f(ii,jj) else c atom pair not in high, take average of all: do 24 ix=1,3 ii=ix+3*(i-1) if(i.eq.j)then i1=ix else i1=1 endif do 24 jx=i1,3 jj=jx+3*(j-1) f(ii,jj)=f(ii,jj)+fo(ii,jj)/dble(N) 24 f(jj,ii)=f(ii,jj) endif 22 continue 501 close(7) CALL WRITEFF(3*nat,3*nat,f) endif end subroutine svgr(g,n,Esum) implicit none real*8 g(*),gn,Esum integer*4 n,i,ix open(70,file='FILE.E') 2 read(70,*,end=99,err=99)gn goto 2 99 backspace 70 write(70,701)Esum 701 format(f24.12) write(6,600)Esum 600 format(' Energy ',f24.12,' appended to FILE.E') close(70) open(70,file='FILE.GR') gn=0.0d0 do 1 i=1,n gn=gn+g(3*i-2)**2+g(3*i-1)**2+g(3*i)**2 1 write(70,700)(g(ix+3*(i-1)),ix=1,3) 700 format(3f15.9) close(70) write(6,*) write(6,*)'Gradient written to FILE.GR' gn=dsqrt(gn/dble(3*n)) write(6,604)gn 604 format(' Gnorm:',f12.6) return end function letter(c) character*1 c logical letter c capital letters A...Z: letter=ichar(c).ge.65.and.ichar(c).le.90 return end subroutine report(s) character*(*) s write(6,*) s stop end subroutine readopt(lff,lwr) logical lff,lwr,lex character*3 s3 c make also force field: lff=.false. c large output: lwr=.false. inquire(file='MAKEGRAD.OPT',exist=lex) if(lex)then open(9,file='MAKEGRAD.OPT') 1 read(9,3,end=99,err=99)s3 3 format(A3) if(s3.eq.'LFF')read(9,*)lff if(s3.eq.'LWR')read(9,*)lwr goto 1 99 close(9) endif return end subroutine gar9f(N,F,IER) c extracts FF from Gaussian archive IMPLICIT none CHARACTER*1 BUFF(140),L,NUMBER(80) CHARACTER*5 SSTR CHARACTER*8 s8 real*8 F(N,N),FU integer*4 I,J,K,NB,IX,IY,N,IN,IER IER=1 NB=70 c c pre-read for faster run: 1 read(7,71,END=8888)s8 71 format(1x,a8) if(s8.eq.'Test job')goto 5444 goto 1 5444 READ(7,7000,END=8888)(BUFF(I),I=1,NB) READ(7,7000,END=8888)(BUFF(I),I=NB+1,2*NB) 7000 FORMAT(1X,70A1) BACKSPACE 7 DO 51 I=1,2*NB-5 DO 52 J=1,5 52 SSTR(J:J)=BUFF(I+J-1) IF(SSTR.EQ.'NImag'.OR.SSTR.EQ.'NIMAG')THEN WRITE(6,*)' Second derivatives found' GOTO 5555 ENDIF 51 CONTINUE GOTO 5444 C C Force Field: 5555 I=I+1 C I .. index of the letter to be gotten IF(I.GT.70)THEN DO 4 K=1,70 4 BUFF(K)=BUFF(K+70) I=I-70 READ(7,*) ENDIF C NB .. line length 5556 I=I+1 CALL GETLET(L,BUFF,NB,I) IF(L.NE.char(92).and.L.NE.char(124))GOTO 5556 5557 I=I+1 CALL GETLET(L,BUFF,NB,I) IF(L.NE.char(92).and.L.NE.char(124))GOTO 5557 DO 23 IX=1,N DO 23 IY=1,IX C C READ F(IX,IY): DO 5 J=1,20 5 NUMBER(J)=' ' IN=0 5558 IN=IN+1 I=I+1 CALL GETLET(L,BUFF,NB,I) NUMBER(IN)=L c IF(L.NE.','.AND.L.NE.'\\'.and.L.NE.'\')GOTO 5558 c ichar('|')=124 IF(ichar(L).NE.44.AND.ichar(L).NE.92.and.ichar(L).ne.124)GOTO 5558 NUMBER(IN)=' ' F(IX,IY)=FU(NUMBER) 23 F(IY,IX)=FU(NUMBER) IER=0 8888 return end SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FCAR(MX3,N) C CONST=4.359828/0.5291772**2 CONST=1.0d0 OPEN(20,FILE='FILE.FC') DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) CLOSE(20) WRITE(*,*)' FF written into FILE.FC' RETURN END SUBROUTINE GETLET(L,BUFF,NB,I) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*1 BUFF(140) CHARACTER*1 L IF(I.GT.NB)THEN IF(I.GT.NB+1)THEN WRITE(*,*)' Cannot read letter ',I c STOP ENDIF I=1 READ(7,7000)(BUFF(J),J=1,70) 7000 FORMAT(1X,70A1) ENDIF L=BUFF(I) RETURN END function FU(NUMBER) implicit none real*8 FU,x CHARACTER*1 NUMBER(*) CHARACTER*20 s20 integer*4 i do 1 i=1,20 1 s20(i:i)=NUMBER(i) read(s20,*,err=99)x FU=x return 99 write(6,*)s20 write(6,*)'format error' stop end subroutine addgrad(g,p,nat,N,io,Esum) implicit none integer*4 N,nat,i,io real*8 g(*),p(nat,N),gn,gx,gy,gz,Esum,e character*80 s80 e=0.0d0 66 read(8,800,end=101,err=101)s80 800 format(a80) if(s80(2:10).eq.'SCF Done:')then do 1 i=11,len(s80) 1 if(s80(i:i).eq.'=')goto 2 2 read(s80(i+1:i+21),*)e Esum=Esum+e endif if(s80(2:29).eq.'ONIOM: extrapolated energy =')then read(s80(30:len(s80)),*)e Esum=Esum+e endif if(s80(32:64).eq.'Integrated Forces (Hartrees/Bohr)'.or. 1 s80(38:59).eq.'Forces (Hartrees/Bohr)')then READ(8,*) READ(8,*) gn=0.0d0 do 13 i=1,nat read(8,*)gx,gx,gx,gy,gz if(p(i,io).gt.0.0d0)then g(1+3*(i-1))=g(1+3*(i-1))+p(i,io)*gx g(2+3*(i-1))=g(2+3*(i-1))+p(i,io)*gy g(3+3*(i-1))=g(3+3*(i-1))+p(i,io)*gz endif c if(i.eq.32)then c write(6,6100)io,gx,gy,gz,p(i,io) c100 format(' 610. atom:',i3,3f10.2,2x,f10.3) c endif 13 gn=gn+gx**2+gy**2+gz**2 gn=dsqrt(gn/dble(3*nat)) write(6,604)gn,e 604 format(' Gnorm:',f12.6,' E:',f24.12) return endif goto 66 101 call report('gradient not found') end