PROGRAM NDFT_W c c vibrational harmonic wafefunction/density c for each nucleus, write the parameters in the analytical c expression on the disk c IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NAT0=900) DIMENSION x(3*NAT0),iz(NAT0),am(NAT0),w(3*NAT0) DIMENSION av(3),bv(3),cv(3),s0(3) character*80 denfile common/cx/s(3*NAT0,3*NAT0) c call loadx(x,nat,iz) call loadm(am,nat) call loads(NAT0,s,nat,w,nm,am) call inputcage(nat,na,nb,nc,av,bv,cv,s0,ic,denfile) call outputf(nat,na,nb,nc,av,bv,cv,s0,ic,denfile,x,iz,am,w,nm) stop end c ====================================================================== subroutine outputf(nat,na,nb,nc,av,bv,cv,s0,ic,denfile, 1x,iz,am,w,nm) PARAMETER (NAT0=900,nb0=1000,N0=3*NAT0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION av(3),bv(3),cv(3),s0(3),x(*),iz(*),r(3),am(*),w(*) DIMENSION ebuff(N0),aj(N0) character*80 denfile common/bb/a(N0,N0),aa(N0,N0),u(N0,N0),uu(N0,N0),au(3,N0), 1c(NAT0,3,3),p(NAT0) common/cx/s(N0,N0) common/bu/buf(nb0) CM=219470.0d0 amu=1822.888d0 idebug=1 N=3*nat open(30,file='NJR.SCR',form='unformatted') do 6 ia=1,nat do 6 ix=1,3 i=3*(ia-1)+ix do 6 ja=1,nat do 6 jx=1,3 j=3*(ja-1)+jx sum=0.0d0 do 7 ii=1,nm 7 sum=sum+s(i,ii)*am(ia)*w(ii)*am(ja)*s(j,ii) 6 a(i,j)=sum c ajacobian=1.0d0 do 23 i=1,nat 23 ajacobian=ajacobian*dsqrt(am(i))**3 write(30)nat,ajacobian c DO 1000 JB=1,nat na1=JB na2=JB c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU do 8 ja=na1, na2 if(idebug.eq.1)write(6,*)' atom ',ja irow=0 do 93 ia=1,nat do 93 ix=1,3 if(ia.ne.ja)then irow=irow+1 i=3*(ia-1)+ix icol=0 do 9 la=1,nat do 9 lx=1,3 if(la.ne.ja)then icol=icol+1 l=3*(la-1)+lx aa(irow,icol)=a(i,l) endif 9 continue endif 93 continue call TRED12(N0,N-3,aa,u,aj,2,IERR,ebuff) c u:second index - mode (new) if(idebug.eq.1)write(6,600)(aj(j)*CM/amu,j=1,N-3) 600 format(6f13.1) tol=1.0d-4 do 11 i=1,N-3 do 11 j=1,N-3 sum=0.0d0 do 12 ii=1,N-3 12 sum=sum+u(i,ii)*u(j,ii)/aj(ii) 11 uu(i,j)=sum do 13 jx=1,3 j=3*(ja-1)+jx do 13 i=1,N-3 sum=0.0d0 icol=0 do 14 ka=1,nat do 14 kx=1,3 if(ka.ne.ja)then icol=icol+1 sum=sum+a(3*(ka-1)+kx,j)*uu(icol,i) endif 14 continue 13 au(jx,i)=sum do 15 jx=1,3 j=3*(ja-1)+jx do 15 ix=1,3 i=3*(ja-1)+ix sum=0.0d0 icol=0 do 16 ka=1,nat do 16 kx=1,3 if(ka.ne.ja)then icol=icol+1 k=3*(ka-1)+kx sum=sum+au(jx,icol)*a(k,i) endif 16 continue 15 c(ja,jx,ix)=a(j,i)-sum prod=1.0d0 pi=3.141592653589d0 do 17 i=1,nm c write(6,*)i,w(i),prod 17 prod=prod*w(i)/pi do 18 k=1,N-3 c write(6,*)k,aj(k),prod 18 prod=prod*pi/aj(k) c write(6,*)dsqrt(prod)*dble(iz(ja)) 8 p(ja)=dsqrt(prod)*dble(iz(ja))*ajacobian c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU write(30)p(JB),iz(JB),((c(JB,i,j),i=1,3),j=1,3) 1000 write(30)(x(3*(JB-1)+i),i=1,3) close(30) return end c ====================================================================== subroutine loadm(am,nat) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION am(*) character*80 filename c amu=1822.888d0 filename='FTRY.INP' c open(10,file=filename) read(10,*,end=9999,err=9999) read(10,*)nsf do 2 i=1,nsf read(10,*) 2 read(10,*) read(10,444)(qdum,i=1,nat) 444 format(6f12.6) read(10,*) read(10,444)(am(i),i=1,nat) do 3 i=1,nat 3 am(i)=am(i)*amu close(10) call msg('Atomic masses loaded from '//filename) return 9999 call report(filename//' not found') end c ====================================================================== subroutine loadx(x,nat,iz) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension x(*),iz(*) open(7,file='FILE.X',status='old') read(7,*) read(7,*)nat do 1 i=1,nat 1 read(7,*)iz(i),(x(3*(i-1)+j),j=1,3) close(7) write(6,*)nat,' atoms in FILE.X' bohr=0.529177d0 do 2 i=1,3*nat 2 x(i)=x(i)/bohr return end c ====================================================================== function notemptyline(ts) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) LOGICAL*4 notemptyline CHARACTER*(*) ts c l=len(ts) notemptyline=.false. do 1 i=1,l 1 if(ts(i:i).ne.' ')notemptyline=.true. return end c ====================================================================== subroutine loads(NAT0,s,nat,w,nm,am) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION s(3*NAT0,3*NAT0),w(*),am(*) c CM=219470.0d0 amu=1822.888d0 open(10,file='F.INP') read(10,*,end=9999,err=9999)nm if(nm.ne.3*nat)call report('nm must ne equal to 3*nat') do 2 i=1,nat+1 2 read(10,*) do 3 ia=1,nat do 3 im=1,nm 3 read(10,*)idum,idum,(s((ia-1)*3+ix,im),ix=1,3) read(10,*)nm read(10,*)(w(im),im=1,nm) do 1 i=1,nm 1 w(i)=w(i)/CM close(10) write(6,*)nm,' modes in F.INP' c write(6,*)'How many zero modes?' read(5,*)nz arb=100.0d0 do 4 i=nm-nz+1,nm 4 w(i)=arb write(6,*)nz,' last (zero) modes frozen' 5 write(6,699)(w(i)*CM,i=1,nm) 699 format(6f13.1) write(6,*)' Mode to freeze (0 to continue):' read(5,*)io if(io.gt.0)then w(io)=arb goto 5 endif 6 write(6,699)(w(i)*CM,i=1,nm) write(6,*)' Mode to allow (0 to continue):' read(5,*)io if(io.gt.0)then do 7 i=1,nm 7 if(i.ne.io)w(i)=arb goto 6 endif c do 8 i=1,nm do 8 k=1,3*nat 8 s(k,i)=s(k,i)/dsqrt(amu) c c orthonormality test tol=1.0d-4 do 9 i=1,nm do 9 j=1,nm aij=0.0d0 do 10 ka=1,nat do 10 kx=1,3 k=3*(ka-1)+kx 10 aij=aij+s(k,i)*s(k,j)*am(ka) if(abs(aij-1.0d0).gt.tol.and.i.eq.j)write(6,605)i,j,aij 605 format(' St.S <> E: ',2i4,f20.10) 9 if(abs(aij).gt.tol.and.i.ne.j)write(6,605)i,j,aij return 9999 call report('F.INP not found') end c ====================================================================== subroutine msg(s) character*(*) s write(6,90)s 90 format(a80) return end c ====================================================================== subroutine report(s) character*(*) s write(6,90)s 90 format(a80) write(6,90)' **** Program Stopped ***' stop end c ====================================================================== SUBROUTINE TRED12(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) 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 (MX3,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(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) 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 subroutine inputcage(nat,na,nb,nc,av,bv,cv,s0,ic,denfile) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER*80 filename,ts,denfile DIMENSION av(3),bv(3),cv(3),s0(3) logical notemptyline c filename='CAGE' open(10,file=filename) c c header 2 read(10,1000,end=899,err=899)ts 1000 format(a80) if(notemptyline(ts))goto 2 c c name 21 read(10,1000,end=899,err=899)ts if(notemptyline(ts))goto 21 c c charge and multiplicity read(10,1000,end=899,err=899)ts c c Cartesian coordinates - replace by new guess do 3 ia=1,nat 3 read(10,1000)ts c read(10,1000)ts read(10,1000)denfile read(10,*)ic,(s0(i),i=1,3) read(10,*)na,(av(i),i=1,3) read(10,*)nb,(bv(i),i=1,3) read(10,*)nc,(cv(i),i=1,3) close(10) bohr=0.529177d0 do 1 i=1,3 s0(i)=s0(i)/bohr av(i)=av(i)/bohr bv(i)=bv(i)/bohr 1 cv(i)=cv(i)/bohr return 899 close(10) call report('cannot read cage') end