program pfix implicit none integer*4 nat,n,i,ix,iargc,IERR,iq real*8 hartreetoJ,p,bohr,pau real*8,allocatable::f(:,:),r(:),g(:),s(:,:),rn(:),e(:),gq(:) integer*4,allocatable::z(:) character*80 s80,fn hartreetoJ=4.3597482d-18 bohr=0.529177d0 write(6,600) 600 format(' pfix: fix molecular geometry for external pressure',/, 1 ' Input: FILE.X',/, 1 ' FILE.FC',/,/, 1 ' G98.OUT',/,/, 1 ' Output: NEWFILE.X',/) if(iargc().ne.1)call report('Usage: pfix

') 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') OPEN(17,FILE='FILE.X',status='old') READ(17,178)fn 178 format(a80) READ(17,*)nat n=3*nat CLOSE(17) allocate(f(n,n),r(n),z(nat),g(n),s(n,n),rn(n),e(n),gq(n)) OPEN(17,FILE='FILE.X') READ(17,*) READ(17,*) do 1 i=1,nat read(17,*)z(i),(r(ix+3*(i-1)),ix=1,3) do 1 ix=1,3 1 r(ix+3*(i-1))=r(ix+3*(i-1)) CLOSE(17) CALL READFF(n,f) call mkg(n,g,pau,r) IERR=0 CALL TRED12(n,f,s,e,2,IERR) if(IERR.ne.0)call report('Diagonalization error') write(6,610) 610 format(' Eigenvalues:') write(6,611)e 611 format(6e12.4) do 6 iq=1,n gq(iq)=0.0d0 do 6 i=1,n 6 gq(iq)=gq(iq)+s(i,iq)*g(i) do 2 i=1,n rn(i)=r(i) do 2 iq=1,n 2 if(e(iq).gt.1.0d-3)rn(i)=rn(i)-s(i,iq)*gq(iq)/e(iq) open(17,file='NEWFILE.X') write(17,178)fn write(17,170)nat do 3 i=1,nat 3 write(17,170)z(i),(rn(ix+3*(i-1)),ix=1,3) 170 format(i6,3f12.6) close(17) 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 gradient:',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 gradient:',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 READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) OPEN(20,FILE='FILE.FC',STATUS='OLD') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) CLOSE(20) WRITE(*,*)' Cartesian FF read in.' RETURN END subroutine report(s) character*(*) s write(6,*)s stop end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END