program cffm implicit none integer*4 m,n,ii,jj,i1,idiff,jdiff,oi,oj,ndiff,imin,i,nmo real*8 li,mi,si,lj,mj,sj,rthrj,sgn,z0,one,two,ds,rs,emin,emax parameter (emax=30000.0d0) real*8,allocatable::a(:),ux(:),uy(:),uz(:),mx(:),my(:),mz(:), 1e(:),rcij(:),icij(:),rbuf(:),ibuf(:), 1uxmor(:),uxmoi(:),uymor(:),uymoi(:),uzmor(:),uzmoi(:), 1mxmor(:),mxmoi(:),mymor(:),mymoi(:),mzmor(:),mzmoi(:) integer*4,allocatable::match(:) logical lex data z0,one,two/0.0d0,1.0d0,2.0d0/ c write(6,600) 600 format(/,' Crystal field f- wavefunctions, dipole moments',/,/, 1 ' (for the lanthanide program)',/,/, 1 ' Input: pauli.inp ',/, 1 ' eigenfunctions.txt ',/, 1 ' energies ',/, 2 ' output: ux uy uz mx my mz.scr.txt',/, 2 ' wf.txt',/) c initialize single and double factorials: call InitFac c read number of basis functions n and electrons m call rd(0,a,n,m) allocate(a(n*m*3)) c read basis functions call rd(1,a,n,m) c dipole moments allocate(ux(n*n)) allocate(uy(n*n)) allocate(uz(n*n)) allocate(mx(n*n)) allocate(my(n*n)) allocate(mz(n*n)) allocate(match(m)) inquire(file='ux.scr.txt',exist=lex) if(.not.lex) then call vz(ux,n*n) call vz(uy,n*n) call vz(uz,n*n) call vz(mx,n*n) call vz(my,n*n) call vz(mz,n*n) do 1 ii=1,n c |ii> = |1 ... oi ... m> do 1 jj=1,n c |jj> = |1 ... oj ... m> i1=ii+n*(jj-1) idiff=0 jdiff=0 ndiff=0 do 101 oj=1,m 101 match(oj)=0 do 2 oi=1,m li=a(1+3*(oi-1)+3*m*(ii-1)) mi=a(2+3*(oi-1)+3*m*(ii-1)) si=a(3+3*(oi-1)+3*m*(ii-1)) c find this substate state = orbital in jj: do 21 oj=1,m lj=a(1+3*(oj-1)+3*m*(jj-1)) mj=a(2+3*(oj-1)+3*m*(jj-1)) sj=a(3+3*(oj-1)+3*m*(jj-1)) if(li.eq.lj.and.mi.eq.mj.and.si.ne.sj)then c orbital oi has counterpart in jj match(oj)=1 goto 2 endif 21 continue c orbital oi in ii does not exist in jj idiff=oi ndiff=ndiff+1 if(ndiff.gt.1)goto 1 2 continue if(ndiff.eq.0)then c diagonal matrix elements do 102 oi=1,m li=a(1+3*(oi-1)+3*m*(ii-1)) mi=a(2+3*(oi-1)+3*m*(ii-1)) c = sum (io) if(mod(NINT(mi),2).eq.0)then sgn= one else sgn=-one endif uz(i1)=uz(i1)+sgn*(two*li+one)* 1 rthrj(li,one,li,-mi,z0,mi)*rthrj(li,one,li,z0,z0,z0) c = sum (io) 102 mz(i1)=mz(i1)+mi else c difference in one orbital, suppose li=lj do 103 oj=1,m 103 if(match(oj).eq.0)jdiff=oj li=a(1+3*(idiff-1)+3*m*(ii-1)) mi=a(2+3*(idiff-1)+3*m*(ii-1)) si=a(3+3*(idiff-1)+3*m*(ii-1)) lj=a(1+3*(jdiff-1)+3*m*(ii-1)) mj=a(2+3*(jdiff-1)+3*m*(ii-1)) sj=a(3+3*(jdiff-1)+3*m*(ii-1)) if(si.eq.sj.and.li.eq.lj)then if(mod(NINT(mj),2).eq.0)then sgn= one else sgn=-one endif if(mi.eq.mj+one)then mx(i1)= 0.5d0*dsqrt(li*(li+1)-mj*(mj+1)) c my is actually imaginary: my(i1)=-0.5d0*dsqrt(li*(li+1)-mj*(mj+1)) ux(i1)= sgn*dsqrt((two*li+one)*(two*lj+one)/two)* 1 rthrj(li,one,lj,-mi,one,mj)*rthrj(li,one,lj,z0,z0,z0) c uy is actually imaginary: uy(i1)=-sgn*dsqrt((two*li+one)*(two*lj+one)/two)* 1 rthrj(li,one,lj,-mi,one,mj)*rthrj(li,one,lj,z0,z0,z0) endif if(mi.eq.mj-one)then mx(i1)= 0.5d0*dsqrt(li*(li+1)-mj*(mj-1)) my(i1)= 0.5d0*dsqrt(li*(li+1)-mj*(mj-1)) ux(i1)= sgn*dsqrt((two*li+one)*(two*lj+one)/two)* 1 rthrj(li,one,lj,-mi,-one,mj)*rthrj(li,one,lj,z0,z0,z0) uy(i1)= sgn*dsqrt((two*li+one)*(two*lj+one)/two)* 1 rthrj(li,one,lj,-mi,-one,mj)*rthrj(li,one,lj,z0,z0,z0) endif endif c si=sj endif c ndiff=1 1 continue open(9,file='ux.scr.txt') call wmtr39(9,ux,n,n,'ux',1) close(9) open(9,file='uy.scr.txt') call wmtr39(9,uy,n,n,'uy',1) close(9) open(9,file='uz.scr.txt') call wmtr39(9,uz,n,n,'uz',1) close(9) open(9,file='mx.scr.txt') call wmtr39(9,mx,n,n,'mx',1) close(9) open(9,file='my.scr.txt') call wmtr39(9,my,n,n,'my',1) close(9) open(9,file='mz.scr.txt') call wmtr39(9,mz,n,n,'mz',1) close(9) write(6,*)'Dipoles written' else open(9,file='ux.scr.txt') call rmtr39(9,ux,n,n,1) close(9) open(9,file='uy.scr.txt') call rmtr39(9,uy,n,n,1) close(9) open(9,file='uz.scr.txt') call rmtr39(9,uz,n,n,1) close(9) open(9,file='mx.scr.txt') call rmtr39(9,mx,n,n,1) close(9) open(9,file='my.scr.txt') call rmtr39(9,my,n,n,1) close(9) open(9,file='mz.scr.txt') call rmtr39(9,mz,n,n,1) close(9) write(6,*)'Dipoles read' endif allocate(e(n),rcij(n*n),icij(n*n),rbuf(n*n),ibuf(n*n), 1uxmor(n*n),uxmoi(n*n),uymor(n*n),uymoi(n*n),uzmor(n*n),uzmoi(n*n), 1mxmor(n*n),mxmoi(n*n),mymor(n*n),mymoi(n*n),mzmor(n*n),mzmoi(n*n)) call reade(e,rcij,icij,n) emin=e(1) imin=1 do 3 i=2,n if(e(i).lt.emin)then emin=e(i) imin=i endif 3 continue c cut-off energies at emax,nmo: do 31 nmo=n,1,-1 31 if(e(nmo)-emin.lt.emax)goto 311 311 write(6,*)nmo,' energies taken' write(6,*)'transforming ux' call taomoc(ux,rbuf,ibuf,uxmor,uxmoi,rcij,icij,nmo,n) write(6,*)'transforming uy' call taomoc(uy,rbuf,ibuf,uymor,uymoi,rcij,icij,nmo,n) write(6,*)'transforming uz' call taomoc(uz,rbuf,ibuf,uzmor,uzmoi,rcij,icij,nmo,n) write(6,*)'transforming mx' call taomoc(mx,rbuf,ibuf,mxmor,mxmoi,rcij,icij,nmo,n) write(6,*)'transforming my' call taomoc(my,rbuf,ibuf,mymor,mymoi,rcij,icij,nmo,n) write(6,*)'transforming mz' call taomoc(mz,rbuf,ibuf,mzmor,mzmoi,rcij,icij,nmo,n) open(9,file='uxr.scr.txt') call wmtr39(9,uxmor,n,n,' ',1) close(9) open(9,file='uyr.scr.txt') call wmtr39(9,uymor,n,n,' ',1) close(9) open(9,file='uzr.scr.txt') call wmtr39(9,uzmor,n,n,' ',1) close(9) open(9,file='mxr.scr.txt') call wmtr39(9,mxmor,n,n,' ',1) close(9) open(9,file='myr.scr.txt') call wmtr39(9,mymor,n,n,' ',1) close(9) open(9,file='mzr.scr.txt') call wmtr39(9,mzmor,n,n,' ',1) close(9) open(9,file='uxi.scr.txt') call wmtr39(9,uxmoi,n,n,' ',1) close(9) open(9,file='uyi.scr.txt') call wmtr39(9,uymoi,n,n,' ',1) close(9) open(9,file='uzi.scr.txt') call wmtr39(9,uzmoi,n,n,' ',1) close(9) open(9,file='mxi.scr.txt') call wmtr39(9,mxmoi,n,n,' ',1) close(9) open(9,file='myi.scr.txt') call wmtr39(9,mymoi,n,n,' ',1) close(9) open(9,file='mzi.scr.txt') call wmtr39(9,mzmoi,n,n,' ',1) close(9) write(6,*)'MO dipoles written' open(9,file='SPEC.TAB') write(9,911) 911 format('lanthanide',/, 1' n wavelength (nm) dipole strength (D^2)', 1' rotatory strength (cgs/10**-36)',/,80(1h-)) do 4 i=1,nmo ds=uxmor(imin+n*(i-1))*uxmor(imin+n*(i-1)) 1 +uymor(imin+n*(i-1))*uymor(imin+n*(i-1)) 1 +uzmor(imin+n*(i-1))*uzmor(imin+n*(i-1)) 1 -uxmoi(imin+n*(i-1))*uxmoi(imin+n*(i-1)) 1 -uymoi(imin+n*(i-1))*uymoi(imin+n*(i-1)) 1 +uzmoi(imin+n*(i-1))*uzmoi(imin+n*(i-1)) 1 +mxmor(imin+n*(i-1))*mxmor(imin+n*(i-1)) 1 +mymor(imin+n*(i-1))*mymor(imin+n*(i-1)) 1 +mzmor(imin+n*(i-1))*mzmor(imin+n*(i-1)) 1 -mxmoi(imin+n*(i-1))*mxmoi(imin+n*(i-1)) 1 -mymoi(imin+n*(i-1))*mymoi(imin+n*(i-1)) 1 +mzmoi(imin+n*(i-1))*mzmoi(imin+n*(i-1)) rs=uxmor(imin+n*(i-1))*mxmor(imin+n*(i-1)) 1 +uymor(imin+n*(i-1))*mymor(imin+n*(i-1)) 1 +uzmor(imin+n*(i-1))*mzmor(imin+n*(i-1)) 1 +uxmoi(imin+n*(i-1))*mxmoi(imin+n*(i-1)) 1 +uymoi(imin+n*(i-1))*mymoi(imin+n*(i-1)) 1 +uzmoi(imin+n*(i-1))*mzmoi(imin+n*(i-1)) 4 write(9,3001)i,e(i)-emin,ds,rs 3001 format(i4,f12.2,2g20.10) write(9,3002) 3002 format(80(1h-)) close(9) write(6,*)'SPEC.TAB written' end subroutine taomoc(ao,rt,it,or,oi,rc,ic,nmo,nao) c transformation with complex coeff, real or purely imaginary c = c*ia cjb implicit none real*8 ao(*),rt(*),it(*),or(*),oi(*),rc(*),ic(*),sumr,sumi,z0 integer*4 nmo,nao,j,i,jj,ii data z0/0.0d0/ c do 1 i=1,nmo do 1 jj=1,nao sumr=z0 sumi=z0 do 2 ii=1,nao sumr=sumr+rc(ii+nao*(i-1))*ao(ii+nao*(jj-1)) 2 sumi=sumi-ic(ii+nao*(i-1))*ao(ii+nao*(jj-1)) rt(i+nao*(jj-1))=sumr 1 it(i+nao*(jj-1))=sumi c do 3 i=1,nmo do 3 j=1,nmo sumr=z0 sumi=z0 do 4 jj=1,nao sumr=sumr+rc(jj+nao*(j-1))*rt(i+nao*(jj-1)) 1 -ic(jj+nao*(j-1))*it(i+nao*(jj-1)) 4 sumi=sumi+ic(jj+nao*(j-1))*rt(i+nao*(jj-1)) 1 +rc(jj+nao*(j-1))*it(i+nao*(jj-1)) or(i+nao*(j-1))=sumr 3 oi(i+nao*(j-1))=sumi return end subroutine reade(e,rc,ic,n) implicit none integer*4 n,i,j real*8 e(*),rc(*),ic(*) open(9,file='energies',status='old') do 1 i=1,n 1 read(9,*)e(i) close(9) write(6,*)'energies red from energies' open(9,file='eigenfunctions.txt',status='old') do 2 i=1,n read(9,*) do 2 j=1,n 2 read(9,*)rc(j+n*(i-1)),rc(j+n*(i-1)),ic(j+n*(i-1)) close(9) write(6,*)'cij red from eigenfunctions.txt' return end subroutine rd(ic,a,n,m) implicit none integer*4 ic,n,m,i,j,k real*8 a(*),ms character*196 s196 open(2,file='pauli.inp',status='old') if(ic.eq.0)then n=0 1 READ(2,*,END=1000,ERR=1000)ms n=n+1 goto 1 1000 write(6,*)n,' basis functions' if(n.eq.0)call report('pogram has nothing to do') rewind 2 read(2,2012)s196 2012 format(a196) m=0 do 5 j=1,len(s196) 5 if(s196(j:j).eq.'.')m=m+1 write(6,*)m,' electrons' if(m.eq.0)call report('pogram has nothing to do') else do 4 i=1,n 4 read(2,*)((a(k+3*(j-1)+3*m*(i-1)),k=1,3),j=1,m) endif close(2) return end subroutine report(s) character*(*) s write(6,*)s stop end Function Fac(n) implicit none integer*4 n,ii real*8 Fac,bb bb=1.0d0 do 1 ii=2 , n 1 bb=bb*dble(ii) fac=bb return end Function DFac(n) c double factorial = n!! implicit none integer*4 n,ii real*8 DFac,bb if(n.eq.-1.or.n.eq.0)then bb=1.0d0 else bb=dble(n) if(mod(n,2).ne.0)then do 1 ii=n-2,1,-2 1 bb=bb*dble(ii) else do 2 ii=n-2,2,-2 2 bb=bb*dble(ii) endif endif dfac=bb return end subroutine Initfac implicit none integer*4 nfac0,i parameter (nfac0=100) real*8 f0,factorial,fac common/npar/f0,factorial(nfac0) real*8 dm,d0,df,dfac common/dpar/dm,d0,df(nfac0) do 1 i=0 , nfac0 1 factorial(i)=fac(i) do 2 i=-1, nfac0 2 df(i)=dfac(i) return end Function rthrj(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none integer*4 nfac0 parameter (nfac0=100) real*8 f0,factorial common/npar/f0,factorial(nfac0) real*8 rthrj,x1,x2,x3,y1,y2,y3,ar,u,sum,msign integer*4 j1d,j2d,j3d,m1d,m2d,m3d,s,ni,i1,i2,i3,i4,i5,i6 j1d=NINT(x1*2.0d0) j2d=NINT(x2*2.0d0) j3d=NINT(x3*2.0d0) m1d=NINT(y1*2.0d0) m2d=NINT(y2*2.0d0) m3d=NINT(y3*2.0d0) if(abs(m1d).gt.j1d.or.abs(m2d).gt.j2d.or.abs(m3d).gt.j3d.or. 1m1d+m2d .ne.-m3d.or.j1d+j2d.lt.j3d.or.j2d+j3d.lt.j1d.or. 2j3d+j1d.lt.j2d .or.mod(j1d-m1d,2).ne.0 .or. 3mod(j2d-m2d,2) .ne.0.or.mod(j3d-m3d,2).ne. 0.or. 4mod(j1d+j2d+j3d,2).ne.0) then ar=0.0d0 else if (mod((j1d-j2d-m3d)/2,2).eq.0) then msign=1.0d0 else msign=-1.0d0 endif s=NINT(x1+x2+x3) ar=factorial(s-j3d)*factorial(s-j2d)*factorial(s-j1d) 1 /factorial(s+1) ar=ar*factorial((j1d+m1d) / 2)*factorial((j1d-m1d) / 2) ar=ar*factorial((j2d+m2d) / 2)*factorial((j2d-m2d) / 2) ar=ar*factorial((j3d+m3d) / 2)*factorial((j3d-m3d) / 2) ar=msign*dsqrt(ar) sum=0.0d0 do 1 ni=0 , s i1=ni i2=(j1d+j2d-j3d) / 2-ni i3=(j1d-m1d) / 2-ni i4=(j2d+m2d) / 2-ni i5=(j3d-j2d+m1d) / 2+ni i6=(j3d-j1d-m2d) / 2+ni if( mod(ni, 2) .ne. 0) then msign=-1.0d0 else msign=1.0d0 endif if(((((i2.ge.0).and.(i3.ge.0)).and.(i4.ge.0)).and.(i5.ge.0)) 1 .and.(i6.ge.0)) then u=factorial(i1)*factorial(i2) u=u*factorial(i3)*factorial(i4) u=u*factorial(i5)*factorial(i6) sum=sum+msign/u endif 1 continue ar=ar*sum if(abs(ar).lt.1.0d-38) ar=0.0d0 endif rthrj=ar return end subroutine vz(v,n) implicit none real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end subroutine wmtr39(io,A,n0,n,st,ic) c c io number of the file c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(*) character*(*) st write(io,*)st write(io,*) N1=1 1 N3=MIN(N1+4,N) WRITE(io,18)(I,I=N1,N3) 18 FORMAT(4X,5I14) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 WRITE(io,17)LN,(A(LN+n0*(J-1)),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(io,*) return end subroutine rmtr39(io,A,n0,n,ic) c c io number of the file c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(*) read(io,*) read(io,*) N1=1 1 N3=MIN(N1+4,N) read(io,*) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 read(io,*)LNdum,(A(LN+n0*(J-1)),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 read(io,*) return end