program gfix implicit none integer*4 nat,n,iargc real*8 hartreetoJ,p,bohr,pau real*8,allocatable::r(:),g(:) integer*4,allocatable::z(:) character*80 s80 hartreetoJ=4.3597482d-18 bohr=0.529177d0 write(6,600) 600 format( 1' gfix: reads G98.OUT and adds gradient for external pressure',/, 1' Input: G98.OUT',/, 1' Output: GM.OUT',/) if(iargc().ne.1)call report('Usage: gfix

') call getarg(1,s80) read(s80,*)p pau=p/hartreetoJ*1.0d-30*bohr**3*1.0d9 write(6,601)p,pau 601 format(' p = ',e12.4,' Gp = ',e12.4,' au') allocate(z(1),r(1)) call rg('G98.OUT',0,nat,z,r,.true.) deallocate(z,r) n=3*nat allocate(z(nat),r(n),g(n)) call rg('G98.OUT',1,nat,z,r,.true.) call mkg(n,g,pau,r) end subroutine rg(fo,ic,nat,z,r,lzmat) implicit none integer*4 ic,nat,z(*),I,l,ix real*8 r(*),q character*(*) fo character*80 s80 logical lzmat open(2,file=fo) 1 read(2,2000,end=99,err=99)s80 2000 format(a80) IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(20:37).EQ.'Input orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard orientation:')))THEN DO 2004 I=1,4 2004 READ(2,*) if(ic.eq.0)then l=0 2005 READ(2,2000)s80 IF(s80(2:4).NE.'---')THEN l=l+1 BACKSPACE 2 READ(2,*)q,q IF(q.EQ.-1)l=l-1 GOTO 2005 ENDIF nat=l else do 2 l=1,nat 2 read(2,*)z(l),z(l),r(1+3*(l-1)),(r(ix+3*(l-1)),ix=1,3) write(6,*)nat,' atoms read' endif close(2) return endif goto 1 99 nat=0 close(2) call report('geometry not found') end subroutine mkg(n,g,p,r) implicit none integer*4 n,np,i,j,isp,nat,ia,iz real*8 g(n),x,y,z,s,st,bohr,xij,yij,zij,d,r(n),p,rc real*8,allocatable::rs(:),go(:) character*80 s80 character*160 s160 nat=n/3 allocate(rs(nat),go(n)) bohr=0.529177d0 np=0 isp=0 st=0.0d0 g=0.0d0 open(17,file='G98.OUT',status='old') 101 read(17,80,end=77,err=77)s80 if(s80(38:59).eq.'Forces (Hartrees/Bohr)')then read(17,*) read(17,*) do 5 i=1,nat 5 read(17,*)(go(1+3*(i-1)),j=1,3),(go(j+3*(i-1)),j=2,3) call wg(' Mechanical force:',n,go) call pg(n,go,r) call wg(' Projected:',n,go) endif goto 101 77 close(17) open(17,file='G98.OUT',status='old') 1 read(17,80,end=88,err=88)s80 80 format(a80) if(s80(4:14).eq.'ISph Re'.or. 1 s80(4:13).eq.'ISph Re')then isp=nat do 3 i=1,isp 3 read(17,*)rs(i),rs(i) endif if(s80(2:26).eq.'GePol: Number of points ')read(s80(56:62),*)np if(s80(4:28).eq.'ITs Coordinates'.and.np.ne.0. 1and.isp.ne.0)then read(17,*) read(17,*) write(6,603)np 603 format(' Number of points:',i6) do 2 i=1,np read(17,*)x,x,y,z,s c which atom: do 4 j=1,nat rc=rs(j) xij=x-r(1+3*(j-1)) yij=y-r(2+3*(j-1)) zij=z-r(3+3*(j-1)) d=dsqrt(xij**2+yij**2+zij**2) 4 if(dabs(d-rc).lt.1.0d-3)goto 41 write(6,604)i,x,y,z 604 format(i6,3f12.6) call report('unassigned tile') 41 st=st+s g(1+3*(j-1))=g(1+3*(j-1))-xij*s*p/bohr**2/rc g(2+3*(j-1))=g(2+3*(j-1))-yij*s*p/bohr**2/rc 2 g(3+3*(j-1))=g(3+3*(j-1))-zij*s*p/bohr**2/rc write(6,600)st 600 format(' Cavity surface area (A^2):',f10.2) call wg(' Pressure force:',n,g) call pg(n,g,r) call wg(' Projected:',n,g) close(17) open(17,file='G98.OUT',status='old') open(18,file='GN.OUT') do 203 i=len(s160),1,-1 203 s160(i:i)=' ' 201 read(17,160,end=66,err=66)s160 160 format(a160) if(s160(38:59).eq.'Forces (Hartrees/Bohr)')then call ws(s160) read(17,160)s160 call ws(s160) read(17,160)s160 call ws(s160) do 6 i=1,nat read(17,*)ia,iz,(go(j+3*(i-1)),j=1,3) 6 write(18,183)ia,iz,(go(j+3*(i-1))+g(j+3*(i-1)),j=1,3) 183 format(i7,i9,7x,3f15.9) goto 201 endif call ws(s160) goto 201 66 close(17) close(18) write(6,*)'GN.OUT written' return endif goto 1 88 close(17) call report('Tesserae not found in G98.OUT (try iop(3/33=2)') end subroutine ws(s) character*(*) s integer*4 i,j do 202 i=len(s),1,-1 202 if(s(i:i).ne.' ')goto 204 204 do 205 j=1,i 205 write(18,181)s(j:j) 181 format(a1,$) write(18,*) return end subroutine wg(s,n,g) implicit none character*(*)s integer*4 n,i real*8 g(n) write(6,'(a)')s do 5 i=1,n/3 5 write(6,602)i,g(1+3*(i-1)),g(2+3*(i-1)),g(3+3*(i-1)) 602 format(i6,3f12.6) return end subroutine report(r) character*(*) r write(6,*)r stop end subroutine pg(n,g,r) c project translations/rotations from gradient implicit none integer*4 n,i,j real*8 g(n),r(n) real*8,allocatable::A(:,:),t(:) allocate(A(n,n),t(n)) call DOMA(A,n,r) t=0.0d0 do 1 i=1,N do 1 j=1,N 1 t(i)=t(i)+A(i,j)*g(j) g=t return end SUBROUTINE DOMA(A,N,r) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),r(N) real*8,allocatable::C(:,:),TEM(:,:) C allocate(C(3,N/3),TEM(N,N)) AMACH=0.00000000000001d0 NAT=N/3 do 11 I=1,NAT do 11 J=1,3 11 C(J,I)=r(J+3*(I-1)) DO 12 I=1,N DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) N6=6 CALL ORT(A,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(*,*)N6,' coordinates projected' DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END SUBROUTINE ORT(A,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(3*NAT,N6) AMACH=0.00000000000001d0 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END