PROGRAM OSCILL c to test harmonic oscillaotr properties IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (nvar0=4,nh0=4100,nbas0=100) logical lper(nvar0),lh dimension ngr(nvar0),ngrh(nvar0),nbase(nvar0),nh(nvar0), 1nbas(nvar0),ic(nvar0),ice(nvar0),ich(nvar0), 3ip(nh0) real*8 max(nvar0),min(nvar0),kelvin common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0) c pi=4.0d0*atan(1.0d0) c c angles supposed in deg (o) c other coords in A c write(6,*)' HO test, based on peswf' C call getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh, 1ice,ic,ich,wcm,x0,ipl,ip,ikin,kinddia,kelvin, 2w,nbas0,nh0,lh,emax1,emax2) c maxb=1 do 25 ii=1,nvar if(nbas(ii).gt.maxb)maxb=nbas(ii) if(nh(ii).gt.maxb)maxb=nh(ii) 25 if(nbase(ii).gt.maxb)maxb=nbase(ii) 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 do 105 ii=1,nvar 105 ank(ii,1)=dsqrt(dsqrt(w(ii)/pi)) k=1 102 k=k+1 akfac=1.0d0 do 101 kk=1,k-1 101 akfac=akfac*dble(kk) do 1011 ii=1,nvar 1011 ank(ii,k)=dsqrt(dsqrt(w(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) p(k,k )=2.0d0*p(k-1,k-1) if(k.lt.maxb)goto 102 c 4 write(6,6000)1.0d0,w(1),maxb 6000 format(' mass: ',f15.6,/, 1 ' freq: ',f15.6,/, 1 ' maxN: ',i15,/,/, 1 ' Input N, M, c, d:',$) read(5,*)n,m,c,d write(6,700) 700 format('state E0 Vii ', 1'pt2 v2ii E0-') do 1 is=1,m e0=h0(is,is) pt2=0.0d0 do 3 ii=1,n 3 if(ii.ne.is) 1pt2=pt2+vkn(is,c,d,ii)*vkn(ii,c,d,is)/(e0-h0(ii,ii)) vii=vkn(is,c,d,is) v2ii=v2nk(is,c,d,is) eexc=e0-(v2ii-vii**2)/pt2 1 write(6,17)is,e0,vii,pt2,v2ii,E0-eexc 17 format(i4,5f12.3,/) goto 4 end c ------------------------------------------------------------ subroutine report(s) character*(*) s write(6,*)s stop end c ------------------------------------------------------------ function h0(n,k) c c c implicit none integer nh0,nvar0,ivar parameter (nh0=4100,nvar0=4) real*8 h0,w,p,ank,x0,ff,dw,ei,s1,s2 integer*4 n,k,nn,kk common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0) ivar=1 dw=dsqrt(w(ivar)) ff=0.0d0 do 1 kk=1,k do 1 nn=1,n s1=-dble((kk-1)*(kk-2))/2.0d0*dw**2*ei(kk+nn-4) s2=dble(2*(kk-1)+1)/2.0d0*dw**2*ei(kk+nn-2) 1 ff=ff+p(k,kk)*p(n,nn)*(s1+s2)/dw h0=ff*ank(ivar,k)*ank(ivar,n) return end c ------------------------------------------------------------ function v2nk(n,c,d,k) c c c implicit none integer nh0,nvar0,ivar parameter (nh0=4100,nvar0=4) real*8 v2nk,w,p,ank,x0,ff,dw,ei,s1,s2,s3,c,d integer*4 n,k,nn,kk common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0) ivar=1 dw=dsqrt(w(ivar)) ff=0.0d0 do 1 kk=1,k do 1 nn=1,n s1= c**2*ei(kk+nn-2+6) s2= d**2*ei(kk+nn-2+8) s3=2.0d0*d*c*ei(kk+nn-2+7) 1 ff=ff+p(k,kk)*p(n,nn)*(s1+s2+s3) v2nk=ff*ank(ivar,k)*ank(ivar,n)/dw return end c ------------------------------------------------------------ function vkn(n,c,d,k) c c c implicit none integer nh0,nvar0,ivar parameter (nh0=4100,nvar0=4) real*8 vkn,w,p,ank,x0,ff,dw,ei,s1,s2,c,d integer*4 n,k,nn,kk common/ho/p(nh0,nh0),ank(nvar0,nh0),w(nvar0),x0(nvar0) ivar=1 dw=dsqrt(w(ivar)) ff=0.0d0 do 1 kk=1,k do 1 nn=1,n s1=c*ei(kk+nn-2+3) s2=d*ei(kk+nn-2+4) 1 ff=ff+p(k,kk)*p(n,nn)*(s1+s2) vkn=ff*ank(ivar,k)*ank(ivar,n)/dw return end subroutine getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh, 1ice,ic,ich,wcm,x0,ipl,ip,ikin,kinddia,T,w,nbas0,nh0, 2lh,emax1,emax2) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) character*4 key c icon=0 ... reads par for surfaces c icon=1 ... reads par for peswf dimension ip(*),ngr(*),min(*),max(*),ngrh(*),nbase(*), 1nbas(*),nh(*),ice(*),ic(*),ich(*),x0(*),w(*),lper(*) logical lper,lh real*8 min,max ipl=0 lh=.false. ikin=1 kinddia=1 T=298.15d0 il=0 open(8,file='SURF.OPT') 1 read(8,800,end=88,err=88)key 800 format(a4) il=il+1 c c Obligatory Parameters: if(key.eq.'NVAR'.and.il.eq.1)then read(8,*)nvar if(nvar.gt.nvar0)call report('Too many variables') do 2 ii=1,nvar write(6,6005)ii 6005 format(' Variable',i2,':') read(8,*)lper(ii) read(8,*)ngr(ii) read(8,*)min(ii),max(ii) read(8,*)ngrh(ii) read(8,*)nbase(ii),nbas(ii),nh(ii) read(8,*)ice(ii),ic(ii),ich(ii) read(8,*)wcm,x0(ii),dumm if(ice(ii).eq.0)then write(6,6002)nbase(ii), ' plane waves','E' 6002 format(i4,a12,' for ',a1) else if(ice(ii).eq.1)then write(6,6002)nbase(ii),' harmonics','E' else write(6,*)'Quadratic fit for E' endif endif if(ic(ii).eq.0)then write(6,6002)nbas(ii), ' plane waves','G' else if(ic(ii).eq.1)then write(6,6002)nbas(ii),' harmonics','G' else write(6,*)'Quadratic fit for G' endif endif if(ich(ii).eq.0)then write(6,6002)nh(ii), ' plane waves','H' else if(ich(ii).eq.1)then write(6,6002)nh(ii),' harmonics','H' else call report('Unknown basis option for H') endif endif if(nbase(ii).gt.nbas0)call report('too many E basis functions') if(nbas(ii). gt.nbas0)call report('too many G basis functions') if(nh(ii). gt.nh0 )call report('too many H basis functions') if(ic(ii).eq.1.or.ich(ii).eq.1)write(6,5006)wcm,x0(ii),1.0d0 5006 format(' HO basis',f12.2,' cm-1, x0 = ',f12.6,' mass = ',f10.2) 2 w(ii)=wcm else if(il.eq.1)call report('NVAR must be defined at the first line.') endif c c Optional Parameters if(key.eq.'STAT')read(8,*)ipl,(ip(ii),ii=1,ipl) if(key.eq.'KINE')read(8,*)ikin if(key.eq.'DIAG')read(8,*)kinddia if(key.eq.'TEMP')read(8,*)T if(key.eq.'HOLE')then lh=.true. read(8,*)emax1,emax2 endif goto 1 88 close(8) return end c function ei(N) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c infinity c / 2 c | N -y c i= | y e dy c | c / c -infinity c nt=N/2 if(nt+nt.ne.N)then ei=0.0d0 else ho=-1.0d0 zl=1.0d0 do 1 i=1,nt ho=ho+2.0d0 1 zl=zl*ho*0.5d0 ei=zl*dsqrt(4.0d0*atan(1.0d0)) endif return end