PROGRAM NDFT_M4 c c Displays vibrational anharmonic wafefunction/density c IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NT0=27) DIMENSION x(3*NT0),iz(NT0),am(NT0) DIMENSION av(3),bv(3),cv(3),s0(3) character*80 denfile c call loadx(x,nat,iz) call loadm(am,nat) call inputcage(nat,na,nb,nc,av,bv,cv,s0,ic,denfile) call loadanharm call init_hermite(am) call outputcage(nat,na,nb,nc,av,bv,cv,s0,ic,denfile,x,iz) stop end c ====================================================================== subroutine loadanharm IMPLICIT none integer*4 NT0,ME0,NM0,NT3,NGEO,N7,N,LEXCL,NHBIG,MV,KG,ISB,IS, 1ICC,IN,ISI,I,J,K,i0,ig,NDIM,IDUMM PARAMETER (NT0=27,ME0=20,NM0=99999,NT3=3*NT0,NGEO=24) real*8 EAH,W,CFI,WDUMM common/anharm/N7,N,LEXCL,NHBIG,MV,KG,ISB,IS(NM0),ICC(NM0,ME0), 1IN(NM0,ME0),EAH(NM0),W(NT3),ISI(NGEO,4),NDIM,i0,ig,CFI(NM0) OPEN(4,FILE='EIGEN.SCR',FORM='UNFORMATTED') READ(4)N7,N,LEXCL,NHBIG,MV,KG,ISB write(6,6000)N7,N,LEXCL,NHBIG,MV,KG,ISB 6000 format('N7:',I3,' N:',I3,' LEXCL:',I3,' NHBIG:',i9,' MV:',I10, 1' KG:',i4,' ISB:',i2) if(NHBIG.GT.NM0)call report('NHBIG>NM0') if(LEXCL.GT.ME0)call report('LEXCL>ME0') if(N.GT.NT3)call report('N>NT3') if(ISB.GT.NGEO)call report('ISB>NGEO') if(KG.GT.NGEO)call report('KG>NGEO') READ(4)(IDUMM,I=1,NHBIG) ,(IS(I),I=1,NHBIG), 1((ICC(I,J),I=1,NHBIG),J=1,LEXCL) ,((IN(I,J),I=1,NHBIG),J=1,LEXCL), 2(EAH(I),I=1,MV),(W(I),I=N7,N),(WDUMM,I=N7,N), 3((IDUMM,K=1,2),J=1,ISB),((ISI(J,K),K=1,4),J=1,KG) CLOSE(4) write(6,*)' Symmetry and state:' read(5,*)i0,ig NDIM=ISI(i0,2)-ISI(i0,1) write(6,*)' NDIM =',NDIM call loadc(ig,CFI,NDIM) return end c ====================================================================== subroutine outputcage(nat,na,nb,nc,av,bv,cv,s0,ic,denfile, 1x,iz) PARAMETER (NT0=27,nb0=1000) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION av(3),bv(3),cv(3),s0(3),x(*),iz(*) character*80 denfile common/bu/buf(nb0) integer*4 ME0,NM0,NT3,NGEO,N7,N,LEXCL,NHBIG,MV,KG,ISB,IS, 1ICC,IN,ISI,I,i0,ig,NDIM PARAMETER (ME0=20,NM0=99999,NT3=3*NT0,NGEO=24) real*8 EAH,W,whm(NT3),bu(NT0,3,3),ahm(NT3) common/anharm/N7,N,LEXCL,NHBIG,MV,KG,ISB,IS(NM0),ICC(NM0,ME0), 1IN(NM0,ME0),EAH(NM0),W(NT3),ISI(NGEO,4),NDIM,i0,ig,CFI(NM0) integer*4 si(NT3),sip(NT3) real*8 p(ME0,ME0),ank(NT3,ME0),wh(NT3),am0(NT3) common/ho/p,ank,wh,am0 N=3*nat write(6,*)' Charge density or probability (0,1):' read(5,*)icp write(6,*)' Interval of atoms ("0 0" for all):' read(5,*)na1,na2 if(na1.eq.0)then na1=1 na2=nat endif open(2,file='BU.SCR',form='unformatted') read(2)(((bu(ia,ix,jx),ia=1,nat),ix=1,3),jx=1,3) close(2) do 3 i=1,N ahm(i)=dsqrt(am0(i)) 3 whm(i)=dsqrt(wh(i)) write(6,5000)ig,EAH(ig)*219470.0d0 5000 format(' state ',i10,' ',F10.2,' cm-1') an=0.0d0 do i=1,NDIM an=an+CFI(i)**2 enddo write(6,*)' norm =',an open(11,file=denfile) write(11,1101) 1101 format(' Nuclear density') write(11,1102) 1102 format(' SCF Total Density') write(11,1103)nat,(s0(i),i=1,3) 1103 format(i5,4f12.6) write(11,1103)na,(av(i),i=1,3) write(11,1103)nb,(bv(i),i=1,3) write(11,1103)nc,(cv(i),i=1,3) cvx=cv(1) cvy=cv(2) cvz=cv(3) if(nc.gt.nb0)call report('nc>nb0') do 1 ia=1,nat 1 write(11,1103)iz(ia),dble(iz(ia)),(x(3*(ia-1)+i),i=1,3) c ainteg=0.0d0 do 2 i1=0,na-1 write(6,*)i1,ainteg do 2 ib=0,nb-1 rx=s0(1)+dble(i1)*av(1)+dble(ib)*bv(1) ry=s0(2)+dble(i1)*av(2)+dble(ib)*bv(2) rz=s0(3)+dble(i1)*av(3)+dble(ib)*bv(3) do 41 ic=1,nc 41 buf(ic)=0.0d0 c call fillbuf(buf,nc,NDIM,N,si,sip,IS,ICC,IN,NM0,ME0,CFI,na1,na2, 1nat,ank,NT3,icp,iz,cvx,cvy,cvz,x,rx,ry,rz,ahm,whm,bu,NT0,p) do 4 ic=1,nc 4 ainteg=ainteg+buf(ic) write(11,1104)(buf(ic),ic=1,nc) 2 continue 1104 format(6G13.5) close(11) write(6,607)ainteg*av(1)*bv(2)*cv(3) 607 format(' Integral over cage: ',f12.4) c c planes open(66,file='xy.prn') do 21 i1=0,na-1 do 22 ib=0,nb-1 22 buf(ib+1)=0.0d0 cvx=bv(1) cvy=bv(2) cvz=0.0d0 rx=s0(1)+dble(i1)*av(1) ry=s0(2)+dble(i1)*av(2) rz=0.0d0 call fillbuf(buf,nb,NDIM,N,si,sip,IS,ICC,IN,NM0,ME0,CFI,na1,na2, 1nat,ank,NT3,icp,iz,cvx,cvy,cvz,x,rx,ry,rz,ahm,whm,bu,NT0,p) do 21 ib=1,nb 21 write(66,88)rx+cvx*dble(ib-1),ry+cvy*dble(ib-1),buf(ib) 88 format(2f15.6,g15.6) close(66) open(66,file='xz.prn') do 23 i1=0,na-1 do 24 ic=0,nc-1 24 buf(ic+1)=0.0d0 cvx=cv(1) cvy=0.0d0 cvz=cv(3) rx=s0(1)+dble(i1)*av(1) ry=0.0d0 rz=s0(3)+dble(i1)*av(3) call fillbuf(buf,nc,NDIM,N,si,sip,IS,ICC,IN,NM0,ME0,CFI,na1,na2, 1nat,ank,NT3,icp,iz,cvx,cvy,cvz,x,rx,ry,rz,ahm,whm,bu,NT0,p) do 23 ic=1,nc 23 write(66,88)rx+cvx*dble(ic-1),rz+cvz*dble(ic-1),buf(ic) close(66) open(66,file='yz.prn') do 25 ib=0,nb-1 do 26 ic=0,nc-1 26 buf(ic+1)=0.0d0 cvx=0.0d0 cvy=cv(2) cvz=cv(3) rx=0.0d0 ry=s0(2)+dble(ib)*bv(2) rz=s0(3)+dble(ib)*bv(3) call fillbuf(buf,nc,NDIM,N,si,sip,IS,ICC,IN,NM0,ME0,CFI,na1,na2, 1nat,ank,NT3,icp,iz,cvx,cvy,cvz,x,rx,ry,rz,ahm,whm,bu,NT0,p) do 25 ic=1,nc 25 write(66,88)ry+cvy*dble(ic-1),rz+cvz*dble(ic-1),buf(ic) close(66) 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 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 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 LOADC(J,CFJ,N) INTEGER*4 J,I,N REAL*8 CFJ(*) OPEN(45,FILE='CFI.SCR',FORM='UNFORMATTED',STATUS='OLD') DO 1 K=1,J-1 1 READ(45) READ(45,err=999)(CFJ(I),I=1,N) CLOSE(45) RETURN 999 write(6,*)J,N,I call report('reading error') END c ========================================================= subroutine init_hermite(am) implicit none integer*4 ME0,NT0,NT3,maxb,ii,jj,k,kk,ize parameter (ME0=20,NT0=27,NT3=3*NT0) real*8 p(ME0,ME0),ank(NT3,ME0),wh(NT3),am0(NT3),am(*) common/ho/p,ank,wh,am0 integer*4 NM0,NGEO,N7,N,LEXCL,NHBIG,MV,KG,ISB,IS, 1ICC,IN,ISI,J,i0,ig,NDIM PARAMETER (NM0=99999,NGEO=24) real*8 EAH,W,akfac,pi,CFI common/anharm/N7,N,LEXCL,NHBIG,MV,KG,ISB,IS(NM0),ICC(NM0,ME0), 1IN(NM0,ME0),EAH(NM0),W(NT3),ISI(NGEO,4),NDIM,i0,ig,CFI(NM0) c pi=4.0d0*atan(1.0d0) maxb=ME0 write(6,*)'Initializing HO basis, kmax = ',maxb do 103 ii=1,maxb do 103 jj=1,maxb 103 p(jj,ii)=0.0d0 p(1,1)=1.0d0 ize=0 do 105 ii=1,N c restrict zero modes arbitrarily: if(W(ii)*219470.0d0.lt.11.0d0)then wh(ii)=20000.0d0/219470.0d0 W(ii)=20000.0d0/219470.0d0 ize=ize+1 else wh(ii)=W(ii) endif am0(ii)=am((ii+2)/3) 105 ank(ii,1)=dsqrt(dsqrt(am0(ii)*wh(ii)/pi)) write(6,*)ize,' zero modes' do 102 k=2,maxb akfac=1.0d0 do 101 kk=1,k-1 101 akfac=akfac*dble(kk) do 106 ii=1,N 106 ank(ii,k)=dsqrt(dsqrt(am0(ii)*wh(ii)/pi)/2.0d0**(k-1)/akfac) p(k,1 )=-p(k-1,2) do 104 j=2,k-2 104 p(k,j )=(2.0d0*p(k-1,j-1)-dble(j)*p(k-1,j+1)) if(k-2.gt.0)p(k,k-1)=2.0d0*p(k-1,k-2) 102 p(k,k )=2.0d0*p(k-1,k-1) return end