program analcfp implicit none integer*4 n0,nmax,n14 parameter (n0=3503,n14=14) real*8 lmax,ll,smax,SM,SI,pexc(4) integer*4 LM,nn,kev,ts,nb,iargc,timemy,maxns real*8,allocatable:: x(:,:),F(:,:),DR(:,:,:,:),DI(:,:,:,:), 1CI(:),cfp(:,:,:), 1x6m(:,:,:,:),ss(:,:),CR(:) integer*4,allocatable::mu(:,:),ls(:,:),ps(:,:),ns(:),bl(:,:), 1ml(:,:) character*20 s20 character*30 date character*1 ssy(17) data ssy/'s','p','d','f','g','h','i','k','l','m','n','o','q','r', 1 't','u','v'/ logical lex,lpolar,lopt(10) write(6,800) 800 format(/, 1' Analyzing lanthanide wavefunction',/,/, 1' Petr Bour, 2015',/,/, 1' Input:ANAL.OPT options, optional',/, 1' energies from the lanthanide.c program',/, 1' eigenfunctions - " -',/, 1' cfp.inp crystal field potential (charges)',/, 1' cfpp.inp crystal field potential', 1 ' (polarizabilities, optional)',/, 1' ',/, 1' Output: STATES.LST list of states and labels',/, 1' spectrum.tab magnetic dipole absorption',/, 1' spectrum2.tab as spectrum.tab, via CFP',/, 1' ecd.tab ECD spectrum for CF',/, 1' mcd.tab MCD spectrum for CF',/) call Initfac call ctime(timemy(),date) write(6,4444)date 4444 format(A30) allocate(cfp(n14,119,119),mu(n14,119),ls(n14,119),ps(n14,119), 1ns(n14)) c read coefficients of fractional parentage: call read_cfp(cfp,mu,ls,ps,ns,maxns) c read basis set: c (the only "fixed" dimension fields on the next line:) allocate(bl(n0,n14),ml(n0,n14),ss(n0,n14)) call readbs(n0,ml,ss,nn,nb,bl,n14) call readopt(lex,kev,nmax,pexc,lpolar,lopt) if(.not.lex)then if(iargc().eq.2)then call getarg(1,s20) read(s20,*)kev call getarg(2,s20) read(s20,*)nmax write(6,6013)kev,nmax 6013 format(' Eigenvalue of interest: ',i5,/, 1 ' Maximal excitation : ',i5) else write(6,6009) 6009 format(/,' Eigenvalue of interest number: ',$) read(5,*)kev write(6,6011) 6011 format(' Maximal excitation: ',$) read(5,*) nmax endif endif if(nmax.gt.nb)call report('Maximum exceeded') if(kev.gt.nmax)call report('Chosen eigenvalue outside range') allocate(CR(nb*nmax),CI(nb*nmax)) c read the wavefunction(s): call readwf(CR,CI,n0,ss,ml,nmax,kev,nb,nn,n14) c analyze basis set occurence within 1...nmax: c call occbf(CR,CI,nmax,nb,n0,n14,bl,ml,ss,nn) c transform basis set to alpha S L wavefunctions: ll=3.0d0 LM=NINT(lmax(ll,nn)) SM=smax(ll,nn) SI=dble(mod(nn,2))*0.5d0 write(6,6010)ssy(NINT(ll)+1),nn,LM,SM,SI 6010 format(/,1x,a1,'-shell, ',i2,' electrons:', 1' LMAX = ',i3,' SMAX = ',f4.1,' SMIN = ',f4.1,/) ts=NINT(2.0d0*SM)+1 allocate(x(nn,maxns),F(nb,maxns), 1DR(nb,maxns,2*LM+1,NINT(2.0d0*SM)+1), 1DI(nb,maxns,2*LM+1,NINT(2.0d0*SM)+1), 1x6m(nb,maxns,2*LM+1,ts)) call trwf1(x6m,n0,119,LM,SM,ml,ss,x,ns,mu,ls,cfp,F,ts,ll,nn,nb, 1n14,maxns) call ctime(timemy(),date) write(6,4444)date c transform wf to alpha S L wavefunctions: call trwf2(kev,ns,LM,DR,DI,CR,CI,x6m,ts,mu,ps,ls,nmax,SM,nn,nb, 1maxns) deallocate(x6m) call ctime(timemy(),date) write(6,4444)date c estimate J: call Jestim(kev,LM,SM,ns,ls,mu,DR,DI,ts,nn,nb,maxns) call ctime(timemy(),date) write(6,4444)date c calculate LS matrices if(lopt(1))then call LSs(n0,ss,ml,ll,CR,CI,nmax,kev,nn,nb,bl,n14) call ctime(timemy(),date) write(6,4444)date endif c calculate magnetic moments if(lopt(2))then call mmoments(n0,ss,ml,ll,CR,CI,nmax,nn,nb,n14) call ctime(timemy(),date) write(6,4444)date endif c calculate electric moments if(lopt(3))then call emoments(n0,ss,ml,ll,CR,CI,nmax,nn,nb,n14,pexc,lpolar) call ctime(timemy(),date) write(6,4444)date endif c calculate magnetic moments using CFP if(lopt(4))then call mmomentsp(ll,LM,ns,SM,nn,cfp,DR,DI,ls,mu,ts,nmax,nb,maxns) call ctime(timemy(),date) write(6,4444)date endif end function lmax(ll,nn) c maximum L in ll-shell with nn electrons implicit none real*8 ll,lc,lmax integer*4 nn,m,ie,ic ie=0 lc=0.0d0 do 1 m=NINT(ll),-NINT(ll),-1 do 1 ic=1,2 if(ie.lt.nn)then ie=ie+1 lc=lc+dble(m) else goto 2 endif 1 continue 2 lmax=lc return end function smax(ll,nn) c maximum S in ll-shell with nn electrons implicit none real*8 ll,smax,s integer*4 nn,il il=NINT(ll) s=dble(min(nn,2*il+1))*0.50d0 if(nn.gt.2*il+1)s=2.0d0*ll+1.0d0-dble(nn)*0.5d0 smax=s return end Function rcg(j1,m1,j2,m2,j3,m3) c {Clebsh-Gordon coefficients} implicit none real*8 rcg,j1,m1,j2,m2,j3,m3,rthrj if (mod(NINT(j1-j2+m3) , 2).eq. 0) then rcg= dsqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) else rcg=-dsqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) endif return end Function rthrj(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none integer*4 nfac0 parameter (nfac0=170) 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 Initfac implicit none integer*4 nfac0,i parameter (nfac0=170) real*8 f0,factorial,fac common/npar/f0,factorial(nfac0) do 1 i=0 , nfac0 1 factorial(i)=fac(i) return 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 subroutine tz(t,n0,m0,n,m) implicit none integer*4 n0,m0,n,m,i,j real*8 t(n0,m0) do 1 i=1,n do 1 j=1,m 1 t(i,j)=0.0d0 return end subroutine read_cfp(cfp,mu,ls,ps,ns,maxns) implicit none integer*4 i,j,k,nn,ne,mu(14,119),ls(14,119),ps(14,119), 1ns(14),ip,jp,fc,snu,ol,os,op,it,maxns real*8 cfp(14,119,119) character*80 s80 logical lex nn=0 it=0 do 1 i=1,14 do 1 j=1,119 do 1 k=1,119 1 cfp(i,j,k)=0.0d0 inquire(file='CFP.LST',exist=lex) if(.not.lex)call report('CFP.LST not found') open(9,file='CFP.LST') do 2 i=1,4 2 read(9,*) 3 read(9,900,end=99,err=99)s80 900 format(a80) if(s80(7:7).ne.' ')then read(s80(6:7),*)ne if(ne.gt.nn)then ip=0 nn=ne endif ip=ip+1 ns(nn)=ip read(s80(10:11),*)mu(nn,ip) ls(nn,ip)=fc(s80(13:13)) read(s80(14:15),*)ps(nn,ip) if(nn.gt.1)then read(s80(37:38),*)os ol=fc(s80(40:40)) read(s80(41:42),*)op jp=snu(os,ol,op,mu,ls,ps,nn-1,ns) read(s80(60:79),*)cfp(nn,ip,jp) it=it+1 4 read(9,900,end=99,err=99)s80 if(s80(70:71).ne.'--')then read(s80(37:38),*)os ol=fc(s80(40:40)) read(s80(41:42),*)op jp=snu(os,ol,op,mu,ls,ps,nn-1,ns) read(s80(60:79),*)cfp(nn,ip,jp) it=it+1 goto 4 endif endif endif goto 3 99 close(9) maxns=ns(1) do 5 i=1,nn if(ns(i).gt.maxns)maxns=ns(i) 5 write(6,600)i,ns(i) 600 format(i4,' electrons ',i3,' states') write(6,*)it,' CFPs, CFP.LST read' return end function snu(os,ol,op,mu,ls,ps,nn,ns) c returns number of a state with nn electrons, os multiplicity, c ol angular momentum, op P number in the list implicit none integer*4 snu,os,ol,op,mu(14,119),ls(14,119),ps(14,119), 1ns(14),nn,j do 1 j=1,ns(nn) if(mu(nn,j).eq.os.and.ls(nn,j).eq.ol.and.ps(nn,j).eq.op)then snu=j return endif 1 continue write(6,*)'nn os ol op:',nn,os,ol,op call report('state not found') end function fc(i) implicit none character*1 i integer*4 fc,c c=-1 if(i.eq.'S')c= 0 if(i.eq.'P')c= 1 if(i.eq.'D')c= 2 if(i.eq.'F')c= 3 if(i.eq.'G')c= 4 if(i.eq.'H')c= 5 if(i.eq.'I')c= 6 if(i.eq.'K')c= 7 if(i.eq.'L')c= 8 if(i.eq.'M')c= 9 if(i.eq.'N')c=10 if(i.eq.'O')c=11 if(i.eq.'Q')c=12 if(i.eq.'R')c=13 if(i.eq.'T')c=14 if(i.eq.'U')c=15 if(i.eq.'V')c=16 if(c.eq.-1)then write(6,*)i call report('Unknown symbol') endif fc=c return end c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ subroutine readbs(n0,ml,ss,nn,nb,bl,n14) implicit none integer*4 n14 integer*4 n0,ml(n0,n14),i,k,ie,nn,nb,bl(n0,n14), 1istart,iend character*400 s400 real*8 ss(n0,n14) real*8,allocatable::t(:) c determine number of electrons nn=0 open(9,file='pauli.inp') read(9,90)s400 90 format(a400) close(9) do 1001 istart=1,len(s400) 1001 if(s400(istart:istart).ne.' ')goto 1002 1002 do 1003 iend=len(s400),1,-1 1003 if(s400(iend:iend).ne.' ')goto 1004 1004 do 1 i=istart,iend-1 1 if(s400(i:i).eq.' '.and.s400(i+1:i+1).ne.' ')nn=nn+1 nn=(nn+1)/3 write(6,*)nn,' electrons' if(nn.gt.n14)call report('Too many electrons') c determine number of basis functions allocate(t(nn*3)) nb=0 open(9,file='pauli.inp') 2 read(9,*,err=3,end=3)t(1),t(2),t(3) nb=nb+1 goto 2 3 close(9) write(6,*)nb,' basis functions' if(nb.gt.n0)call report('Too many basis functions') open(9,file='pauli.inp') do 4 i=1,nb read(9,*)(t(k),k=1,3*nn) k=0 do 4 ie=1,nn k=k+1 bl(i,ie)=nint(t(k)) k=k+1 ml(i,ie)=nint(t(k)) k=k+1 4 ss(i,ie)=t(k) close(9) write(6,70) 70 format(' Basis set read from pauli.inp.') return end c ============================================================ subroutine readwf(CR,CI,n0,ss,ml,nmax,kev,nb,nn,n14) implicit none integer*4 n10,n14 parameter (n10=5) integer*4 n0,nmax,nn,ii,ia(n10),j,kk,BM,ie,ml(n0,n14),i,kev,nb real*8 CR(*),CI(*),aa(n10),st,ss(n0,n14) logical lex inquire(file='eigenfunctions',exist=lex) if(lex)then open(9,file='eigenfunctions') do 1 i=1,nmax do 1 j=1,nb 1 read(9,*)CR(j+nb*(i-1)),CI(j+nb*(i-1)) close(9) write(6,71)nmax 71 format(i6,' eigenstates read from eigenfunctions.') else inquire(file='eigenfunctions.txt',exist=lex) if(lex)then open(9,file='eigenfunctions.txt') do 31 i=1,nmax read(9,*) do 31 j=1,nb 31 read(9,*)CR(j+nb*(i-1)),CR(j+nb*(i-1)),CI(j+nb*(i-1)) close(9) write(6,70)nmax 70 format(i6,' eigenstates read from eigenfunctions.txt.') else call report('eigenfunctions not found') endif endif do 101 ii=1,n10 aa(ii)=0.0d0 ia(ii)=0 do 101 j=1,nb do 102 kk=1,ii-1 102 if(ia(kk).eq.j)goto 101 if(CR(j+nb*(kev-1))**2+CI(j+nb*(kev-1))**2.gt.aa(ii))then aa(ii)=CR(j+nb*(kev-1))**2+CI(j+nb*(kev-1))**2 ia(ii)=j endif 101 continue c write the largest contributions: write(6,*)'Most contributing basis functions to eigenvalue',kev, 1':' do 3 kk=1,n10 st=0.0d0 BM=0 do 5 ie=1,nn st=st+ss(kk,ie) 5 BM=BM+ml(kk,ie) write(6,6001)ia(kk) 6001 format(i6,$) write(6,6005)' ' do 4 ie=1,nn 4 write(6,6002)ss(kk,ie) 6002 format(f4.1,$) write(6,6005)' ' 6005 format(a1,$) do 6 ie=1,nn 6 write(6,6003)ml(kk,ie) 6003 format(i2,$) 3 write(6,6004)st,BM,aa(kk)*100.d0 6004 format(' S = ',f4.1,' M = ',i3,f10.4,' %') return end c ============================================================ subroutine trwf1(x6m,n0,n119,LM,SM,ml,ss,x,ns,mu,ls, 1cfp,F,ts,ll,nn,nb,n14,maxns) implicit none integer*4 n14,k,n0,n119,LM,i,j,kk,ii,ml(n0,n14),ns(*), 1mu(14,n119),ls(14,n119),jj,ts,nn,nb,maxns real*8 SM,x6m(nb,maxns,2*LM+1,ts),x(nn,maxns), 1sn,so,jn,jo,cfp(14,119,119),rcg,F(nb,maxns), 2sh,ss(n0,n14),ll real*8,allocatable::m(:),z(:),ms(:),zs(:) logical lex integer*4 nfac0 parameter (nfac0=170) real*8 f0,factorial common/npar/f0,factorial(nfac0) c wavefunction P = sum(i) Ci Xi c Xi = l1m1s1z1 l2m2 s2z2 ... lnmnsnzn c = sum(LS) D^SL |LM>|SZ> c Xi = sum(LS) EiLS |LM>|SZ> c => D^SL = sum(i) Ci EiLS c c Xi=|L2M2>|S2Z2>|l3m3s3z3l4m4s4z4...... c = |l4m4s4z4... c allocate(m(nn),z(nn),ms(nn),zs(nn)) sh=0.50d0 inquire(file='F.SCR',exist=lex) if(lex)then open(9,file='F.SCR',form='unformatted') do 804 ii=1,ts do 804 k=1,2*LM+1 804 read(9)((x6m(i,j,k,ii),i=1,nb),j=1,ns(nn)) close(9) write(6,*)'F-coefficients read from F.SCR' else do 802 i=1,nb do 802 j=1,ns(nn) do 802 kk=1,2*LM+1 do 802 ii=1,ts 802 x6m(i,j,kk,ii)=0.0d0 c bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb do 6 i=1,nb do 801 kk=1,nn m(kk)=dble(ml(i,kk)) 801 z(kk)=dble(ss(i,kk)) do 800 kk=1,nn ms(kk)=0.0d0 zs(kk)=0.0d0 do 800 ii=1,kk ms(kk)=ms(kk)+m(ii) 800 zs(kk)=zs(kk)+z(ii) c {s1l1 s2l2 s3l3 .... snln} .. antisymmetrized WF c transcript to sum(S2L2) X_S2L2 |S2L2> {s3l3 ... snln} c p2=ps(2,ii) c j2=L2 sp=S2 c x2: cfp(l^2 L2S2P2]|sl) c s3=sum(P2L2S2) c x2: cfp(l^2 L2S2P2]|sl) c x cfp(l^3 L3S3P3]|l^2S3L2P2sl) c etc call tz(x,nn,maxns,nn,maxns) x(1,1)=1.0d0 do 803 kk=2,nn do 803 ii=1,ns(kk) sn=dble(mu(kk,ii)-1)/2.0d0 jn=dble(ls(kk,ii)) do 803 jj=1,ns(kk-1) so=dble(mu(kk-1,jj)-1)/2.0d0 jo=dble(ls(kk-1,jj)) 803 x(kk,ii)=x(kk,ii)+x(kk-1,jj)*rcg(jo,ms(kk-1),ll,m(kk),jn,ms(kk)) 1 * rcg(so,zs(kk-1),sh,z(kk),sn,zs(kk)) 1 * cfp(kk,ii,jj) do 13 ii=1,ns(nn) x(nn,ii)=x(nn,ii)*dsqrt(factorial(nn)) x6m(i,ii,NINT(ms(nn))+LM+1,NINT(zs(nn)+SM)+1)=x(nn,ii) 13 F(i,ii)=x(nn,ii) 6 continue c bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb open(9,file='F.SCR',form='unformatted') do 805 ii=1,ts do 805 k=1,2*LM+1 805 write(9)((x6m(i,j,k,ii),i=1,nb),j=1,ns(nn)) close(9) write(6,*)'F-coefficients written to F.SCR' endif return end c ============================================================ subroutine trwf2(kev,ns,LM,DR,DI,CR,CI,x6m,ts,mu,ps,ls,nmax,SM,nn, 1nb,maxns) implicit none integer*4 n10,maxns parameter (n10=5) integer*4 ns(*),kev,ev,LM,kk,i,ii,k,i3,nn,nb, 1ia(n10),is(n10),ip(n10),nns(n10),ts,mu(14,119),ls(14,119), 1ps(14,119),jj,cind,io,nmax,nmax_w,jmax,j,ij(n10) real*8 DR(nb,maxns,2*LM+1,ts),CR(*),aa(n10),anorm,SM,rcg, 1DI(nb,maxns,2*LM+1,ts), 1x6m(nb,maxns,2*LM+1,ts),amax,al,ml,as,ms,mj,CI(*) character*1 creturn,achar character*1 sy(17) data sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R', 1 'T','U','V'/ logical lex creturn=achar(13) inquire(file='D.SCR',exist=lex) if(lex)then open(9,file='D.SCR',form='unformatted') read(9)nmax_w close(9) endif if(lex.and.nmax_w.eq.nmax)then open(9,file='D.SCR',form='unformatted') read(9)nmax_w do 171 j=1,ns(nn) do 171 kk=1,2*LM+1 do 171 i=1,ts read(9)(DR(ev,j,kk,i),ev=1,nmax) 171 read(9)(DI(ev,j,kk,i),ev=1,nmax) close(9) write(6,*)'D.SCR for',nmax,' eigenvalues read' else io=-1 do 17 ev=1,nmax if(ev*100/nmax.gt.io)then io=ev*100/nmax write(6,777,advance='NO')creturn,io 777 format(a1,' Contracting F-coefficients ',i4,' %') endif cind=(ev-1)*nb do 17 j=1,ns(nn) do 17 kk=1,2*LM+1 do 17 i=1,ts DR(ev,j,kk,i)=0.0d0 DI(ev,j,kk,i)=0.0d0 do 17 ii=1,nb DR(ev,j,kk,i)=DR(ev,j,kk,i)+CR(ii+cind)*x6m(ii,j,kk,i) 17 DI(ev,j,kk,i)=DI(ev,j,kk,i)+CI(ii+cind)*x6m(ii,j,kk,i) open(9,file='D.SCR',form='unformatted') write(9)nmax do 172 j=1,ns(nn) do 172 kk=1,2*LM+1 do 172 i=1,ts write(9)(DR(ev,j,kk,i),ev=1,nmax) 172 write(9)(DI(ev,j,kk,i),ev=1,nmax) close(9) write(6,*)'D.SCR written' endif do 19 ii=1,n10 aa(ii)=0.0d0 ia(ii)=0 is(ii)=0 ip(ii)=0 nns(ii)=0 i3=0 do 19 jj=1,ns(nn) do 19 k=1,2*LM+1 do 19 i=1,ts i3=i3+1 do 20 kk=1,ii-1 20 if(nns(kk).eq.i3)goto 19 if(DR(kev,jj,k,i)**2+DI(kev,jj,k,i)**2.gt.aa(ii))then aa(ii)=DR(kev,jj,k,i)**2+DI(kev,jj,k,i)**2 ia(ii)=ls(nn,jj) is(ii)=mu(nn,jj) ip(ii)=ps(nn,jj) nns(ii)=i3 c get J jmax=0 amax=0.0d0 al=dble(ia(ii)) ml=dble(-LM +k-1) as=dble(is(ii)-1)/2.0d0 ms= -SM+dble(i-1) mj=ml+ms do 22 j=abs((2*ia(ii)-is(ii)+1)/2),(2*ia(ii)+is(ii)-1)/2 if((rcg(al,ml,as,ms,dble(j),mj))**2.gt.amax)then amax=rcg(al,ml,as,ms,dble(j),mj)**2 Jmax=j endif 22 continue ij(ii)=Jmax endif 19 continue anorm=0.0d0 do 21 j=1,ns(nn) do 21 k=1,2*LM+1 do 21 i=1,ts 21 anorm=anorm+DR(kev,j,k,i)**2+DI(kev,j,k,i)**2 write(6,6002)kev, 1 is(1), 1 is(2), 1 is(3), 1nns(1),ip(1), sy(ia(1)+1),aa(1)*100.d0, 1nns(2),ip(2), sy(ia(2)+1),aa(2)*100.d0, 1nns(3),ip(3), sy(ia(3)+1),aa(3)*100.d0,anorm, 1 ij(1), 1 ij(2), 1 ij(3) 6002 format(/,i5,' eigen state largest contributions:',/,/, 14x,3(5x,' ',2x,' ',i2,1x,1x, 7x,' '),/, 14x,3(i5,'(',i2,')',2x,A1,1x,f7.2,' %'),f8.4,/, 14x,3(5x,' ',2x,' ',2x, i2, 7x,' '),/) return end c ============================================================ subroutine Jestim(kev,LM,SM,ns,ls,mu,DR,DI,ts,nn,nb,maxns) implicit none integer*4 n10,maxns parameter (n10=5) integer*4 LM,d1,d2,kev,ns(*),mu(14,119),ls(14,119), 1i,jj,k,ii,i3,j,kk,ts,nn,nb,nns(n10),Jind,Mind,lj,sj,ijmax real*8 jmin,jmax,SM,rcg,DR(nb,maxns,2*LM+1,ts),ia(n10), 1DI(nb,maxns,2*LM+1,ts), 1al,so,aj,am,az,mj,aa(n10),anorm,is(n10) real*8,allocatable::GR(:,:),GI(:,:) jmin=0.0d0 jmax=DABS(DBLE(LM)+SM) ijmax=NINT(2.0d0*jmax) if(mod(ijmax,2).ne.0)jmin=0.5d0 d1=NINT(jmax-jmin)+1 d2=NINT(2.0d0*jmax)+1 write(6,6004)LM,SM,jmin,jmax 6004 format(' Mmax = ',i5 ,' Smax = ',f5.1,/, 1 ' Jmin = ',f5.1,' Jmax = ',f5.1) allocate(GR(d1,d2),GI(d1,d2)) do 1 i=1,d1 do 1 j=1,d2 GR(i,j)=0.0d0 1 GI(i,j)=0.0d0 c suppose that psi_kev=sum D_jj |aSjjZjjLjjMLjj> = sum G_JM |JM> c then G_JM = sum D_jj do 22 Jind=NINT(jmin*2.0d0),NINT(jmax*2.0d0),2 i=(Jind-NINT(2.0d0*jmin)+2)/2 aj=dble(Jind)/2.0d0 do 22 Mind=-Jind,Jind,2 j=(Mind+NINT(2.0d0*Jmax)+2)/2 mj=dble(Mind)/2.0d0 do 22 jj=1,ns(nn) lj=ls(nn,jj) sj=mu(nn,jj)-1 if(abs(2*lj-sj).le.Jind.and.abs(2*lj+sj).ge.Jind)then al=dble(lj) so=dble(sj)/2.0d0 do 23 k=1,2*LM+1 am=dble(-LM+k-1) do 23 ii=1,ts az=-SM+dble(ii-1) if(2*(-LM+k-1)-NINT(2.0d0*SM)+2*ii-2.eq.Mind)then GR(i,j)=GR(i,j)+DR(kev,jj,k,ii)*rcg(al,am,so,az,aj,mj) GI(i,j)=GI(i,j)+DI(kev,jj,k,ii)*rcg(al,am,so,az,aj,mj) endif 23 continue endif 22 continue do 24 ii=1,n10 aa(ii)=0.0d0 ia(ii)=0.0d0 is(ii)=0.0d0 nns(ii)=0 i3=0 do 24 i=1,d1 do 24 j=1,d2 i3=i3+1 do 26 kk=1,ii-1 26 if(nns(kk).eq.i3)goto 24 if(GR(i,j)**2+GI(i,j)**2.gt.aa(ii))then aa(ii)=GR(i,j)**2+GI(i,j)**2 c J: ia(ii)=jmin+dble(i-1) c MJ: is(ii)=dble(j-1)-jmax nns(ii)=i3 endif 24 continue anorm=0.0d0 do 25 i=1,d1 do 25 j=1,d2 25 anorm=anorm+GR(i,j)**2+GI(i,j)**2 write(6,6003) nns(1),ia(1),is(1),aa(1)*100.d0, 1 nns(2),ia(2),is(2),aa(2)*100.d0, 1 nns(3),ia(3),is(3),aa(3)*100.d0, 1 anorm 6003 format(3(i4,' J=',f3.1,' MJ=',f4.1,f7.2,'%'),g9.3) return end c ============================================================ function sxf(l1,s1,m1,z1,l2,s2,m2,z2) c TAS (1b) p221 (4^8) implicit none real*8 sxf,l1,s1,m1,z1,l2,s2,m2,z2,t,two,z0,half integer*4 is1,is2,im1,im2,iz1,iz2,il1,il2 data z0,half,two/0.0d0,0.5d0,2.0d0/ t=z0 il1=NINT(l1) il2=NINT(l2) is1=NINT(two*s1) is2=NINT(two*s2) im1=NINT(m1) im2=NINT(m2) iz1=NINT(two*z1) iz2=NINT(two*z2) if((il1.ne.il2).or.(is1.ne.is2).or.(im1.ne.im2))goto 99 if(iz1.eq.iz2+2.or.iz1.eq.iz2-2)t=half 99 sxf=t return end c ============================================================ function lsf(l1,s1,m1,z1,l2,s2,m2,z2) c TAS (1c) p221 (4^8) implicit none real*8 lsf,l1,s1,m1,z1,l2,s2,m2,ma,mb,z2,t,two,z0,half integer*4 is1,is2,im1,im2,ima,imb,iz1,iz2,il1,il2 data z0,half,two/0.0d0,0.5d0,2.0d0/ t=z0 il1=NINT(l1) il2=NINT(l2) is1=NINT(two*s1) is2=NINT(two*s2) im1=NINT(m1) im2=NINT(m2) iz1=NINT(two*z1) iz2=NINT(two*z2) ma=m1+z1 mb=m2+z2 ima=NINT(two*ma) imb=NINT(two*mb) if((il1.ne.il2).or.(is1.ne.is2).or.(ima.ne.imb))goto 99 if(im1.eq.im2.and.iz1.eq.iz2)t=z1*m1 if(iz1.eq.iz2+2.or.iz1.eq.iz2-2) 1t=t+half*dsqrt( (l1-ma+half)*(l1+ma+half)) 99 lsf=t return end c ============================================================ function lxf(l1,s1,m1,z1,l2,s2,m2,z2) c TAS (1a) p221 (4^8) implicit none real*8 lxf,l1,s1,m1,z1,l2,s2,m2,z2,t,two,z0,one,half integer*4 is1,is2,im1,im2,iz1,iz2,il1,il2 data z0,half,one,two/0.0d0,0.5d0,1.0d0,2.0d0/ t=z0 il1=NINT(l1) il2=NINT(l2) is1=NINT(two*s1) is2=NINT(two*s2) im1=NINT(m1) im2=NINT(m2) iz1=NINT(two*z1) iz2=NINT(two*z2) if((il1.ne.il2).or.(is1.ne.is2).or.(iz1.ne.iz2))goto 99 if(im1.eq.im2+1.or.im1.eq.im2-1)t=half*dsqrt(l1*(l1+one)-m1*m2) 99 lxf=t return end c ============================================================ subroutine LSs(n0,ss,ml,ll,CR,CI,nmax,kev,nn,nb,bl,n14) implicit none integer*4 n14,n0,ml(n0,n14),il,i,j,ie,ij,idiff,nmax,iep, 2kev,NP,nn,a(2),ap(2),nb,ix,b,bl(n0,n14),S_s,L_s,J_s real*8 ss(n0,n14),ll,sf,CR(*),as2,als,al2,EL,sh,MS,MLP,x,aj2, 1MSP,M,MP,ma,map,za,zap,mb,mbp,zb,zbp,lxf,sxf,lsf,CI(*),la,e_s, 1e_snm real*8,allocatable::maor(:),maoi(:),e(:) logical lex integer*4,allocatable::mi(:),mj(:),li(:),lj(:),lli(:) character*1 sy(17) data sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R', 1 'T','U','V'/ allocate(mi(nn),mj(nn),li(nn),lj(nn),lli(nn)) sh=0.5d0 il=NINT(ll) c AO LS matrices: allocate(maor(3*nb**2),maoi(3*nb**2)) inquire(file='LSAO.SCR',exist=lex) if(lex)then open(9,file='LSAO.SCR',form='unformatted') do 101 i=1,nb read(9)((maor(ix+3*((i-1)+nb*(j-1))),ix=1,3),j=1,nb) 101 read(9)((maoi(ix+3*((i-1)+nb*(j-1))),ix=1,3),j=1,nb) close(9) write(6,*)'LS AO integrals read from LSAO.SCR' else write(6,*)'calculating LS AO integrals' do 1 i=1,nb c spin-orbital numbers: EL=0.0d0 MS=0.0d0 do 2 ie=1,nn lli(ie)=bl(i,ie) mi(ie) =ml(i,ie) EL=EL+mi(ie) MS=MS+ss(i,ie) 2 li(ie)=mi(ie)+il+1+(2*il+1)*((NINT(2.0d0*ss(i,ie))+1)/2) M=EL+MS c number of pairs of individual sets with the same n l ml: NP=0 do 21 ie=1,nn do 21 iep=ie+1,nn 21 if(mi(ie).eq.mi(iep))NP=NP+1 do 1 j=1,nb MLP=0.0d0 MSP=0.0d0 do 3 ie=1,nn mj(ie)=ml(j,ie) MLP=MLP+mj(ie) MSP=MSP+ss(j,ie) 3 lj(ie)=mj(ie)+il+1+(2*il+1)*((NINT(2.0d0*ss(j,ie))+1)/2) MP=MLP+MSP als=0.0d0 as2=0.0d0 al2=0.0d0 c try to order lj according to li call odiff(li,lj,2,idiff,sf,nn,a,ap) if(idiff.eq.0)then c diagonal: al2=EL**2 do 4 ie=1,nn la=dble(lli(ie)) ma=dble(mi(ie)) za=ss(i,ie) al2=al2+la*(la+1.0d0)-ma**2 do 4 iep=ie+1,nn mb=dble(mi(iep)) zb=ss(i,iep) 4 al2=al2-4.0d0*(lxf(ll,sh,ma,za,ll,sh,mb,zb))**2 as2=as2+MS**2+dble(nn)/2.0d0-dble(NP) als=als+MS*EL else if(idiff.le.2)then c off-diagonal za=dble(2*((a(1)-1)/(2*il+1))-1)/2.0d0 ma=dble(-il+mod(a(1)-1,2*il+1)) zap=dble(2*((ap(1)-1)/(2*il+1))-1)/2.0d0 map=dble(-il+mod(ap(1)-1,2*il+1)) if(idiff.eq.1)then if(NINT(2.0d0*M).eq.NINT(2.0d0*MP))then c TAS p 222 (5b) (4^8): als=als+sf*lsf(ll,sh,ma,za,ll,sh,map,zap) do 22 ie=1,nn if(li(ie).ne.a(1))then b=li(ie) zb=dble(2*((b-1)/(2*il+1))-1)/2.0d0 mb=dble(-il+mod(b-1,2*il+1)) c TAS p 222 (5b) (4^8): als=als-2.0d0*sf*( 1 lxf(ll,sh,ma,za,ll,sh,mb,zb)*sxf(ll,sh,mb,zb,ll,sh,map,zap)+ 1 lxf(ll,sh,map,zap,ll,sh,mb,zb)*sxf(ll,sh,mb,zb,ll,sh,ma,za)) endif 22 continue endif else zb=dble(2*((a(2)-1)/(2*il+1))-1)/2.0d0 mb=dble(-il+mod(a(2)-1,2*il+1)) zbp=dble(2*((ap(2)-1)/(2*il+1))-1)/2.0d0 mbp=dble(-il+mod(ap(2)-1,2*il+1)) if(NINT(2.0d0*M).eq.NINT(2.0d0*MP))then c TAS p 222 (5a) (4^8): als=als+sf*2.0d0*( 1 lxf(ll,sh,ma,za,ll,sh,map,zap)*sxf(ll,sh,mb,zb,ll,sh,mbp,zbp)+ 1 lxf(ll,sh,mb,zb,ll,sh,mbp,zbp)*sxf(ll,sh,ma,za,ll,sh,map,zap)- 1 lxf(ll,sh,ma,za,ll,sh,mbp,zbp)*sxf(ll,sh,mb,zb,ll,sh,map,zap)- 1 lxf(ll,sh,mb,zb,ll,sh,map,zap)*sxf(ll,sh,ma,za,ll,sh,mbp,zbp)) c TAS p 221 (4a) (4^8): as2=as2+sf*4.0d0*( 1 sxf(ll,sh,ma,za,ll,sh,map,zap)*sxf(ll,sh,mb,zb,ll,sh,mbp,zbp)- 1 sxf(ll,sh,ma,za,ll,sh,mbp,zbp)*sxf(ll,sh,mb,zb,ll,sh,map,zap)) endif if(NINT(2.0d0*M).eq.NINT(2.0d0*MP))then c TAS p 221 (3a) (4^8): al2=sf*4.0d0*( 1 lxf(ll,sh,ma,za,ll,sh,map,zap)*lxf(ll,sh,mb,zb,ll,sh,mbp,zbp)- 1 lxf(ll,sh,ma,za,ll,sh,mbp,zbp)*lxf(ll,sh,mb,zb,ll,sh,map,zap)) endif endif endif endif ij=3*((i-1)+nb*(j-1)) maor(1+ij)=als maor(2+ij)=as2 maor(3+ij)=al2 maoi(1+ij)=0.0d0 maoi(2+ij)=0.0d0 maoi(3+ij)=0.0d0 1 continue open(9,file='LSAO.SCR',form='unformatted') do 102 i=1,nb write(9)((maor(ix+3*((i-1)+nb*(j-1))),ix=1,3),j=1,nb) 102 write(9)((maoi(ix+3*((i-1)+nb*(j-1))),ix=1,3),j=1,nb) close(9) write(6,*)'LS AO integrals written to LSAO.SCR' endif c transform AO to MO: inquire(file='LSMO.SCR',exist=lex) if(lex)then open(9,file='LSMO.SCR',form='unformatted') do 103 i=1,nmax read(9)((maor(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) 103 read(9)((maoi(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) close(9) write(6,*)'MO integrals read from LSMO.SCR' else call trafo3(nb,maor,maoi,nmax,CR,CI) open(9,file='LSMO.SCR',form='unformatted') do 104 i=1,nmax write(9)((maor(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) 104 write(9)((maoi(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) close(9) write(6,*)' MO integrals written to LSMO.SCR' endif c L^2+2L.S+S^2: als=maor(1+3*((kev-1)+nmax*(kev-1))) as2=maor(2+3*((kev-1)+nmax*(kev-1))) al2=maor(3+3*((kev-1)+nmax*(kev-1))) aj2=as2+2.0d0*als+al2 x=aj2 write(6,609)kev,as2,al2,als,dsqrt(x+0.25d0)-0.5d0 609 format(' State',i4,': S^2 = ',f10.4,/, 1 ' ',4x,' L^2 = ',f10.4,/, 1 ' ',4x,' L.S = ',f10.4,/, 1 ' ',4x,' J = ',f10.4) allocate(e(nmax)) call re(e,nmax) open(9,file='STATE.LST') do 5 i=1,nmax as2=maor(2+3*((i-1)+nmax*(i-1))) x=as2 c 2S: S_s=NINT(2.0d0*(dsqrt(x+0.25d0)-0.5d0)) al2=maor(3+3*((i-1)+nmax*(i-1))) x=al2 c L: L_s=NINT(dsqrt(x+0.25d0)-0.5d0) als=maor(1+3*((i-1)+nmax*(i-1))) aj2=as2+2.0d0*als+al2 x=aj2 c 2J: J_s=NINT(2.0d0*(dsqrt(x+0.25d0)-0.5d0)) e_s=e(i)-e(1) if(e_s.gt.0.01d0)then e_snm=min(1.0d7/e_s,99999.9999d0) else e_snm=0.0d0 endif if(L_s+1.le.17)then if(mod(J_s,2).eq.0)then write(9,909)i,e_s,e_snm,S_s+1,sy(L_s+1),J_s/2 else write(9,919)i,e_s,e_snm,S_s+1,sy(L_s+1),J_s,2 endif else write(9,919)i,e_s,e_snm,S_s+1,'?',J_s/2 endif 5 continue 909 format(i6,f13.3,' cm-1 ',f10.4,' nm',i5,A1,i2) 919 format(i6,f13.3,' cm-1 ',f10.4,' nm',i5,A1,i2,'/',i1) close(9) write(6,*)'State labels written into STATE.LST' return end c ============================================================ subroutine emoments(n0,ss,ml,ll,CR,CI,nmax,nn,nb,n14,pc,lpolar) implicit none integer*4 n14,n0,ml(n0,n14),nb,z1,z3,il,i,j,ie,ij,idiff,nmax, 2nn,a(2),ap(2),ikdiff,kjdiff,api(2), 1apj(2),ai(2),aj(2),l1,l3,m1,m3,ie_free,mlp,is,lko integer*4 md,mo,zo real*8 ss(n0,n14),ll,sf,CR(*),xr,yr,zr,so,sfi,sfj,tsp,CI(*), 1almi(64),almr(64),pi,xi,yi,zim,sikr,siki,skji,skjr, 1Y1zik,Y1pik,Y1mik,Y1zkj,Y1pkj,Y1mkj,pc(*) real*8,allocatable::maor(:),maoi(:),mgr(:),mgi(:),st(:), 1polr(:),poli(:),c(:) integer*4,allocatable::li(:),lj(:),zi(:),mi(:),zj(:),mj(:), 1mt(:),zt(:),lk(:),lid(:),ljd(:),nic(:), 1kic(:),mic(:) logical lex,lpolar character*1 achar write(6,6003)pc(1),pc(2),pc(3) 6003 format(' Electric moment contributions:',/, 1' f - > d ',f12.6,/, 1' d - > f ',f12.6,/, 1' f - > g ',f12.6,/) c read the crystal field Atp parameters: open(9,file='cfp.inp') do 6 i=1,64 6 read(9,*)almr(i),almi(i),almr(i),almi(i) close(9) write(6,*)'Crystal filed parameters read from cfp.inp' if(lpolar)then allocate(polr(120),poli(120),c(120),nic(120),kic(120), 1 mic(120)) open(9,file='cfpp.inp') do 20 i=1,120 20 read(9,*)nic(i),nic(i),kic(i),kic(i),mic(i),c(i),c(i), 1 c(i),polr(i),poli(i) close(9) write(6,*)'CP polarizability parameters read from cfpp.inp' endif il=NINT(ll) pi=4.0d0*atan(1.0d0) tsp=dsqrt(2.0d0*pi/3.0d0) allocate(zi(nn),mi(nn),zj(nn),mj(nn),li(nn),lj(nn), 1st(nn),zt(14),lk(nn+1),lid(nn+1),ljd(nn+1),mt(14)) c AO magnetic moment matrix: allocate(maor(3*nb**2),maoi(3*nb**2)) inquire(file='EAO.SCR',exist=lex) if(lex)then call rwi('r','EAO.SCR',nb,'Electric AO integrals',maor,maoi) else write(6,*)'calculating electric AO integrals' do 1 i=1,nb write(6,800,advance='NO')achar(13),i,nb 800 format(a1,i5,' of ',i5) c number spin-orbitals for all electrons: do 2 ie=1,nn mi(ie)=ml(i,ie) zi(ie)=NINT(2.0d0*ss(i,ie)) li(ie)=mi(ie)+il+1+(2*il+1)*((zi(ie)+1)/2) c remember this as a basis for excitations: lid(ie)=li(ie) 2 lk(ie)=li(ie) do 1 j=1,nb do 3 ie=1,nn mj(ie)=ml(j,ie) zj(ie)=NINT(2.0d0*ss(j,ie)) lj(ie)=mj(ie)+il+1+(2*il+1)*((zj(ie)+1)/2) 3 ljd(ie)=lj(ie) xr=0.0d0 yr=0.0d0 zr=0.0d0 xi=0.0d0 yi=0.0d0 zim=0.0d0 call odiff(li,lj,2,idiff,sf,nn,a,ap) if(idiff.le.2)then IF(pc(1).ne.0.0d0)THEN c generate |f^(n-1) d^1> states: c take the l^n i configuration as a basis: do 111 ie=1,nn mt(ie)=mi(ie) st(ie)=ss(i,ie) 111 zt(ie)=zi(ie) c one after nother, move electron ie to a d orbital md: do 11 ie=1,nn do 11 md=-2,2 c do not change spin: c remember old configuration: mo=mt(ie) so=st(ie) zo=zt(ie) lko=lk(ie) c define new by putting electron ie to d-shell: mt(ie)=md lk(ie)=14+md+3+5*((zt(ie)+1)/2) call odiff(li,lk,2,ikdiff,sfi,nn,ai,api) call odiff(lk,lj,2,kjdiff,sfj,nn,aj,apj) if(ikdiff.eq.1.and.kjdiff.eq.1)then c calculate the matrix elements , etc. c (diagonal (should not happen) sikr=0.0d0 siki=0.0d0 Y1zik=0.0d0 Y1pik=0.0d0 Y1mik=0.0d0 c off diagonal l1=3 z1=2*((ai(1)-1)/7)-1 m1=-il+mod(ai(1)-1,7) if(api(1).le.14)call report('api<=14') l3=2 z3=2*((api(1)-14-1)/5)-1 m3=-2+mod(api(1)-14-1,5) call uvp(Y1zik,Y1mik,Y1pik,sikr,siki,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) skjr=0.0d0 skji=0.0d0 Y1zkj=0.0d0 Y1pkj=0.0d0 Y1mkj=0.0d0 c off diagonal if(aj(1).le.14)call report('aj <=14') l1=2 z1=2*((aj(1)-14-1)/5)-1 m1=-2+mod(aj(1)-14-1,5) l3=3 z3=2*((apj(1)-1)/7)-1 m3=-il+mod(apj(1)-1,7) if(z1.eq.z3)then call uvp(Y1zkj,Y1mkj,Y1pkj,skjr,skji,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) endif xr=xr+ tsp*(sikr*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skjr)*pc(1) xi=xi+ tsp*(siki*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skji)*pc(1) yr=yr+ tsp*(siki*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skji)*pc(1) yi=yi- tsp*(sikr*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skjr)*pc(1) zr=zr +2.0d0*tsp*(sikr* Y1zkj + Y1zik *skjr)*pc(1) zim=zim+2.0d0*tsp*(siki* Y1zkj + Y1zik *skji)*pc(1) endif c move the electron ie back: lk(ie)=lko zt(ie)=zo mt(ie)=mo 11 st(ie)=so ENDIF c fN -> FN-1 d IF(pc(2).ne.0.0d0)THEN c write(6,6565)i,(li(is),is=1,nn),j,(lj(is),is=1,nn) c565 format(i4,6i3,'->',i4,6i3) c generate |d^9 f^(n+1) > states: c take FREE places in the l^n i configuration as a basis: ie_free=0 do 211 is=-1,1,2 do 211 mlp=-3,3 c orbital |MZ> do 2111 ie=1,nn c occupied-do nothing: 2111 if(mi(ie).eq.mlp.and.zi(ie).eq.is)goto 211 c free-record: ie_free=ie_free+1 mt(ie_free)=mlp zt(ie_free)=is 211 continue if(ie_free.ne.14-nn)call report('ie_free <> 14-nn') c one after another, move electron from a d orbital md to a free place c (the spins must be the same): do 12 md=-2,2 do 12 ie=1,8 c old configuration: d^1 f^6 lid(nn+1)=14+md+3+5*((zt(ie)+1)/2) ljd(nn+1)=14+md+3+5*((zt(ie)+1)/2) c define new f^7 lk(nn+1)=mt(ie)+4+7*((zt(ie)+1)/2) call odiff(lid,lk,2,ikdiff,sfi,nn+1,ai,api) call odiff(lk,ljd,2,kjdiff,sfj,nn+1,aj,apj) if(ikdiff.eq.1.and.kjdiff.eq.1)then c calculate the matrix elements , etc. c (diagonal (should not happen) if(ai(1).le.14)call report('ai<=14') if(api(1).gt.14)call report('api>14') if(aj(1).gt.14)call report('aj>14') if(apj(1).le.14)call report('aj<=14') sikr=0.0d0 siki=0.0d0 Y1zik=0.0d0 Y1mik=0.0d0 Y1pik=0.0d0 c off diagonal l1=2 z1=2*((ai(1)-14-1)/5)-1 m1=-2+mod(ai(1)-14-1,5) l3=3 z3=2*((api(1)-1)/7)-1 m3=-3+mod(api(1)-1,7) call uvp(Y1zik,Y1mik,Y1pik,sikr,siki,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) c write(6,6566)md,ie,l1,m1,z1,l3,m3,z3,Y1zik,sikr c566 format(2i2,1x,3i2,1x,3i2,2g10.4) skjr=0.0d0 skji=0.0d0 Y1zkj=0.0d0 Y1mkj=0.0d0 Y1pkj=0.0d0 c off diagonal l1=3 z1=2*((aj(1)-1)/7)-1 m1=-3+mod(aj(1)-1,7) l3=2 z3=2*((apj(1)-14-1)/5)-1 m3=-2+mod(apj(1)-14-1,5) if(z1.eq.z3)then call uvp(Y1zkj,Y1mkj,Y1pkj,skjr,skji,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) c write(6,6566)md,ie,l1,m1,z1,l3,m3,z3,Y1zik,sikr endif xr=xr+ tsp*(sikr*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skjr)*pc(2) xi=xi+ tsp*(siki*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skji)*pc(2) yr=yr+ tsp*(siki*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skji)*pc(2) yi=yi- tsp*(sikr*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skjr)*pc(2) zr=zr +2.0d0*tsp*(sikr* Y1zkj + Y1zik *skjr)*pc(2) zim=zim+2.0d0*tsp*(siki* Y1zkj + Y1zik *skji)*pc(2) endif 12 continue ENDIF c d^10 f^6 -> d^9 f^7 IF(pc(3).ne.0.0d0)THEN c generate |f^(n-1) g^1> states: c take the l^n i configuration as a basis: do 311 ie=1,nn mt(ie)=mi(ie) st(ie)=ss(i,ie) 311 zt(ie)=zi(ie) c one after nother, move electron ie to a g orbital md: do 13 ie=1,nn do 13 md=-4,4 c do not change spin: c remember old configuration: mo=mt(ie) so=st(ie) zo=zt(ie) lko=lk(ie) c define new by putting electron ie to g-shell: mt(ie)=md lk(ie)=14+md+5+9*((zt(ie)+1)/2) call odiff(li,lk,2,ikdiff,sfi,nn,ai,api) call odiff(lk,lj,2,kjdiff,sfj,nn,aj,apj) if(ikdiff.eq.1.and.kjdiff.eq.1)then c calculate the matrix elements , etc. c (diagonal (should not happen) sikr=0.0d0 siki=0.0d0 Y1zik=0.0d0 Y1mik=0.0d0 Y1pik=0.0d0 c off diagonal ~ l1=3 z1=2*((ai(1)-1)/7)-1 m1=-il+mod(ai(1)-1,7) if(api(1).le.14)call report('api g <= 14') l3=4 z3=2*((api(1)-14-1)/9)-1 m3=-4+mod(api(1)-14-1,9) call uvp(Y1zik,Y1mik,Y1pik,sikr,siki,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) skjr=0.0d0 skji=0.0d0 Y1zkj=0.0d0 Y1mkj=0.0d0 Y1pkj=0.0d0 c off diagonal ~ if(aj(1).le.14)call report('aj g <= 14') l1=4 z1=2*((aj(1)-14-1)/9)-1 m1=-4+mod(aj(1)-14-1,9) l3=3 z3=2*((apj(1)-1)/7)-1 m3=-il+mod(apj(1)-1,7) if(z1.eq.z3)then call uvp(Y1zkj,Y1mkj,Y1pkj,skjr,skji,l1,m1,l3,m3,almr,almi, 1 polr,poli,c,nic,kic,mic,lpolar) endif xr=xr+ tsp*(sikr*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skjr)*pc(3) xi=xi+ tsp*(siki*(Y1pkj+Y1mkj)+(Y1pik+Y1mik)*skji)*pc(3) yr=yr+ tsp*(siki*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skji)*pc(3) yi=yi- tsp*(sikr*(Y1pkj-Y1mkj)+(Y1pik-Y1mik)*skjr)*pc(3) zr=zr +2.0d0*tsp*(sikr* Y1zkj + Y1zik *skjr)*pc(3) zim=zim+2.0d0*tsp*(siki* Y1zkj + Y1zik *skji)*pc(3) endif c move the electron ie back: lk(ie)=lko zt(ie)=zo mt(ie)=mo 13 st(ie)=so ENDIF c fN -> FN-1 g endif ij=3*((i-1)+nb*(j-1)) maor(1+ij)=xr maor(2+ij)=yr maor(3+ij)=zr maoi(1+ij)=xi maoi(2+ij)=yi maoi(3+ij)=zim 1 continue call rwi('w','EAO.SCR',nb,'Electric AO integrals',maor,maoi) endif c transform AO to MO: inquire(file='EMO.SCR',exist=lex) if(lex)then call rwi('r','EMO.SCR',nmax,'Electric MO integrals',maor,maoi) else call trafo3(nb,maor,maoi,nmax,CR,CI) call rwi('w','EMO.SCR',nmax,'Electric MO integrals',maor,maoi) endif allocate(mgr(3*nmax**2),mgi(3*nmax**2)) call rwi('r','MMO.SCR',nmax,'Magnetic MO integrals',mgr,mgi) call dospecg(maor,maoi,mgr,mgi,nmax,'ecd.tab') call dospecf(maor,maoi,mgr,mgi,nmax,'fluor.tab') call dospece(maor,maoi,mgr,mgi,nmax,'mcd.tab') call dospecm(maor,maoi,mgr,mgi,nmax,'mcpl.tab') return end c ============================================================ subroutine rwi(sr,fn,nmax,s,ir,ii) implicit none integer*4 nmax,i,j,ix real*8 ir(*),ii(*) character*(*) fn,s,sr open(9,file=fn,form='unformatted') if(sr.eq.'r')then do 102 i=1,nmax read(9)((ir(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) 102 read(9)((ii(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) write(6,*)s//' read from '//fn else do 103 i=1,nmax write(9)((ir(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) 103 write(9)((ii(ix+3*((i-1)+nmax*(j-1))),ix=1,3),j=1,nmax) write(6,*)s//' written to '//fn endif close(9) return end c ============================================================ subroutine mmoments(n0,ss,ml,ll,CR,CI,nmax,nn,nb,n14) implicit none integer*4 n14,n0,ml(n0,n14),nb,il,i,j,ie,ij,idiff,z_j,m_j,z_i, 1m_i,l2,nmax,nn,a(1),ap(1) real*8 ss(n0,n14),ll,sf,CR(*),CI(*),tmx,tmy,tmz real*8,allocatable::maor(:),maoi(:) integer*4,allocatable::li(:),lj(:),zi(:),mi(:),zj(:),mj(:) logical lex allocate(zi(nn),mi(nn),zj(nn),mj(nn),li(nn),lj(nn)) il=NINT(ll) l2=il*(il+1) c AO magnetic moment matrix: allocate(maor(3*nb**2),maoi(3*nb**2)) inquire(file='MAO.SCR',exist=lex) if(lex)then call rwi('r','MAO.SCR',nb,'Magnetic AO integrals',maor,maoi) else write(6,*)'calculating magnetic AO integrals' do 1 i=1,nb c number spin-orbitals for all electrons: do 2 ie=1,nn mi(ie)=ml(i,ie) zi(ie)=NINT(2.0d0*ss(i,ie)) 2 li(ie)=mi(ie)+il+1+(2*il+1)*((zi(ie)+1)/2) do 1 j=1,nb do 3 ie=1,nn mj(ie)=ml(j,ie) zj(ie)=NINT(2.0d0*ss(j,ie)) 3 lj(ie)=mj(ie)+il+1+(2*il+1)*((zj(ie)+1)/2) tmx=0.0d0 tmy=0.0d0 tmz=0.0d0 call odiff(li,lj,2,idiff,sf,nn,a,ap) if(idiff.eq.0)then c Mz = L + 2S diagonal: do 4 ie=1,nn 4 tmz=tmz+dble(mi(ie)+zi(ie)) else c off-diagonal: if(idiff.eq.1)then z_i=2*((a(1)-1)/(2*il+1))-1 m_i=-il+mod(a(1)-1,2*il+1) z_j=2*((ap(1)-1)/(2*il+1))-1 m_j=-il+mod(ap(1)-1,2*il+1) if(z_j.eq.z_i)then if(m_i.eq.m_j+1)then tmx=tmx+0.5d0*dsqrt(dble(l2-m_j*m_i))*sf tmy=tmy-0.5d0*dsqrt(dble(l2-m_j*m_i))*sf endif if(m_i.eq.m_j-1)then tmx=tmx+0.5d0*dsqrt(dble(l2-m_j*m_i))*sf tmy=tmy+0.5d0*dsqrt(dble(l2-m_j*m_i))*sf endif endif if(m_i.eq.m_j)then if(z_i.eq.z_j+2)then tmx=tmx+sf tmy=tmy-sf endif if(z_i.eq.z_j-2)then tmx=tmx+sf tmy=tmy+sf endif endif endif endif ij=3*((i-1)+nb*(j-1)) maor(1+ij)=tmx maor(2+ij)=0.0d0 maor(3+ij)=tmz maoi(1+ij)=0.0d0 maoi(2+ij)=tmy maoi(3+ij)=0.0d0 1 continue call rwi('w','MAO.SCR',nb,'Magnetic AO integrals',maor,maoi) endif c transform AO to MO: inquire(file='MMO.SCR',exist=lex) if(lex)then call rwi('r','MMO.SCR',nmax,'Magnetic MO integrals',maor,maoi) else call trafo3(nb,maor,maoi,nmax,CR,CI) call rwi('w','MMO.SCR',nmax,'Magnetic MO integrals',maor,maoi) endif call dospec(maor,nb,nmax,'spectrum.tab') return end c ============================================================ function fiu(zj,mj,j,zi,mi,i,SM,ts,LM,tl,n14) implicit none integer*4 fiu,zj,mj,j,zi,mi,i,SM,ts,LM,tl,n14 fiu=3*(((zj+SM)/2)+ts*(mj+LM+tl*(j-1+ 1n14 *(((zi+SM)/2)+ts*(mi+LM+tl*(i-1)))))) return end c ============================================================ subroutine mmomentsp(ll,LM,ns,MS,nn,cfp,DR,DI,ls,mu,ts,nmax,nb, 1maxns) implicit none integer*4 ii,ix,iy,iz,ns(*),il,l2,tl,i,nb, 1li,lj,si,sj,mi,mj,SM,ts,k,j,LM,iu,iw, 1mnp,mk,mn,znp,nn,iix,iiy,iiz,tsl,zn, 1nmax,m_i,m_j,z_i,z_j,zi,zj,zk,lk,sk, 1mu(14,119),ls(14,119),fiu,maxns real*8 mz,ll,sur,sui,msx,msy,msz, 1mmx,mmy,mmz,om,os,MS,coe,cfp(14,119,119),rcg,co2, 1alk,amn,amnp,ali,ami,amj,ask,azk,sh,asi,azi,aznp, 1asj,azj,azn,DI(nb,maxns,2*LM+1,ts),alj,mx,my,amk,co1, 1DR(nb,maxns,2*LM+1,ts) real*8,allocatable::mao(:),tr(:),mstr(:),msti(:),ti(:) character*2 mxyz(3) data mxyz/'mx','my','mz'/ character*1 creturn,achar logical lex il=NINT(ll) l2=il*(il+1) tsl=2*il+1 sh=0.5d0 c AO magnetic moment matrix < l mi s zi|m_a| l mj s zj>: allocate( mao(3*(tsl*2)**2)) ii=0 do 1 m_i=-il,il do 1 z_i=-1,1,2 do 1 m_j=-il,il do 1 z_j=-1,1,2 ix=ii+1 iy=ii+2 iz=ii+3 ii=ii+3 mmx=0.0d0 mmy=0.0d0 mmz=0.0d0 msx=0.0d0 msy=0.0d0 msz=0.0d0 os=0.0d0 om=0.0d0 if(z_j.eq.z_i )then msz=msz+dble(z_i) os=1.0d0 endif if(z_i.eq.z_j+2)then msx=msx+0.5d0*dsqrt(dble(3-z_j*z_i)) msy=msy-0.5d0*dsqrt(dble(3-z_j*z_i)) endif if(z_i.eq.z_j-2)then msx=msx+0.5d0*dsqrt(dble(3-z_j*z_i)) msy=msy+0.5d0*dsqrt(dble(3-z_j*z_i)) endif if(m_i.eq.m_j)then mmz=mmz+dble(m_i) om=1.0d0 endif if(m_i.eq.m_j+1)then mmx=mmx+0.5d0*dsqrt(dble(l2-m_j*m_i)) mmy=mmy-0.5d0*dsqrt(dble(l2-m_j*m_i)) endif if(m_i.eq.m_j-1)then mmx=mmx+0.5d0*dsqrt(dble(l2-m_j*m_i)) mmy=mmy+0.5d0*dsqrt(dble(l2-m_j*m_i)) endif mao(ix)=mmx*os+msx*om mao(iy)=mmy*os+msy*om mao(iz)=mmz*os+msz*om 1 continue write(6,*)'One-electron integrals made' tl=2*LM+1 SM=NINT(2.0d0*MS) ts=SM+1 write(6,601)3*(ns(nn)*tl*ts)**2*8/1024000 601 format(' Allocating',i6,' MBytes') allocate(mstr(3*(ns(nn)*tl*ts)**2),msti(3*(ns(nn)*tl*ts)**2)) write(6,*)'Allocated' inquire(file='MM.SCR',exist=lex) if(lex)then open(9,file='MM.SCR',form='unformatted') do 22 i=1,ns(nn) li=ls(nn,i) si=mu(nn,i)-1 do 22 mi=-li,li do 22 zi=-si,si,2 do 22 j=1,ns(nn) lj=ls(nn,j) sj=mu(nn,j)-1 do 22 mj=-lj,lj do 22 zj=-sj,sj,2 iu=fiu(zj,mj,j,zi,mi,i,SM,ts,LM,tl,ns(nn)) read(9) mstr(iu+1),mstr(iu+2),mstr(iu+3) 22 read(9) msti(iu+1),msti(iu+2),msti(iu+3) close(9) write(6,*)'State integrals read from MM.SCR' else c make the matrix creturn=achar(13) do 2 i=1,ns(nn) write(6,777,advance='NO')creturn,i,ns(nn) 777 format(a1,' state ',i4,' of ',i4) c alphai=ps(6,i) li= ls(nn,i) ali=dble(li) si= mu(nn,i)-1 asi=dble(si)/2.0d0 do 2 mi=-li,li ami= dble(mi) do 2 zi=-si,si,2 azi= dble(zi)/2.0d0 do 2 j=1,ns(nn) c alphaj=ps(6,j) lj= ls(nn,j) alj=dble(lj) sj= mu(nn,j)-1 asj=dble(sj)/2.0d0 do 2 mj=-lj,lj amj= dble(mj) do 2 zj=-sj,sj,2 azj= dble(zj)/2.0d0 mx=0.0d0 my=0.0d0 mz=0.0d0 do 33 k=1,ns(nn-1) c alphak=ps(5,k) lk= ls(nn-1,k) if(lk.gt.li+il.or.lk.lt.abs(li-il))goto 33 if(lk.gt.lj+il.or.lk.lt.abs(lj-il))goto 33 alk=dble(lk) sk= mu(nn-1,k)-1 if(sk.gt.si+1.or.sk.lt.abs(si-1))goto 33 if(sk.gt.sj+1.or.sk.lt.abs(sj-1))goto 33 ask=dble(sk)/2.0d0 co1=dble(nn)*cfp(nn,i,k)*cfp(nn,j,k) if(co1.eq.0.0d0)goto 33 do 32 znp=-1,1,2 zk=zj-znp zn=zi-zk if(abs(zn).ne.1)goto 32 aznp=dble(znp)/2.0d0 azk =dble(zk) /2.0d0 azn =dble(zn) /2.0d0 co2=co1*rcg(ask,azk,sh, azn,asi,azi)*rcg(ask,azk,sh,aznp,asj,azj) do 31 mnp=-il,il mk=mj-mnp mn=mi-mk if(abs(mn).gt.il)goto 31 amnp=dble(mnp) amk=dble(mk) amn=dble(mn) coe=co2* 1 rcg(alk,amk,ll,amn ,ali,ami)*rcg(alk,amk,ll,amnp,alj,amj) iix=1+3*((znp+1)/2+2*((mnp+il)+tsl*(((zn+1)/2)+2*(mn+il)))) iiy=iix+1 iiz=iiy+1 mx=mx+coe*mao(iix) my=my+coe*mao(iiy) mz=mz+coe*mao(iiz) 31 continue 32 continue 33 continue iu=fiu(zj,mj,j,zi,mi,i,SM,ts,LM,tl,ns(nn)) mstr(iu+1)=mx mstr(iu+2)=0.0d0 mstr(iu+3)=mz msti(iu+1)=0.0d0 msti(iu+2)=my 2 msti(iu+3)=0.0d0 open(9,file='MM.SCR',form='unformatted') do 23 i=1,ns(nn) li=ls(nn,i) si=mu(nn,i)-1 do 23 mi=-li,li do 23 zi=-si,si,2 do 23 j=1,ns(nn) lj=ls(nn,j) sj=mu(nn,j)-1 do 23 mj=-lj,lj do 23 zj=-sj,sj,2 iu=fiu(zj,mj,j,zi,mi,i,SM,ts,LM,tl,ns(nn)) write(9) mstr(iu+1),mstr(iu+2),mstr(iu+3) 23 write(9) msti(iu+1),msti(iu+2),msti(iu+3) close(9) write(6,*)'State integrals written to MM.SCR' write(6,*)'State dipoles done' endif inquire(file='MP.SCR',exist=lex) if(lex)then call rwi('r','MP.SCR',nmax,'Magnetic MO-P integrals',mstr,msti) else c trafo with D write(6,600)ns(nn)*ts*tl,ns(nn)*ts*tl,nb,nb 600 format(' Space of SL states :',i8,' x ',i8,/, 1 ' Space of eigen vectors:',i8,' x ',i8) allocate (tr((max(nb,ns(nn)*ts*tl))**2)) allocate (ti((max(nb,ns(nn)*ts*tl))**2)) write(6,6003) (max(nb,ns(nn)*ts*tl))**2/1024000 6003 format(' Allocated ',i8,' MBytes') do 10 ix=1,3 write(6,*)'transforming ',mxyz(ix) write(6,*)'first index' c transcript only: do 9 i=1,ns(nn) li=ls(nn,i) si=mu(nn,i)-1 do 9 mi=-li,li do 9 zi=-si,si,2 do 9 j=1,ns(nn) lj=ls(nn,j) sj=mu(nn,j)-1 do 9 mj=-lj,lj do 9 zj=-sj,sj,2 iu=fiu(zj,mj,j,zi,mi,i,SM,ts,LM,tl,ns(nn)) tr(iu/3+1)=mstr(ix+iu) 9 ti(iu/3+1)=msti(ix+iu) write(6,*)'rewritten' c trafo: do 8 I=1,nmax do 8 j=1,ns(nn) lj=ls(nn,j) sj=mu(nn,j)-1 do 8 mj=-lj,lj do 8 zj=-sj,sj,2 sur=0.0d0 sui=0.0d0 do 81 k=1,ns(nn) lk=ls(nn,k) sk=mu(nn,k)-1 do 81 mk=-lk,lk do 81 zk=-sk,sk,2 iu=fiu(zj,mj,j,zk,mk,k,SM,ts,LM,tl,ns(nn)) sur=sur+tr(iu/3+1)*DR(I,k,mk+LM+1,(zk+SM)/2+1) 1 -ti(iu/3+1)*DI(I,k,mk+LM+1,(zk+SM)/2+1) 81 sui=sui+tr(iu/3+1)*DI(I,k,mk+LM+1,(zk+SM)/2+1) 1 +ti(iu/3+1)*DR(I,k,mk+LM+1,(zk+SM)/2+1) iw=3*((zj+SM)/2+ts*(mj+LM+tl*(j-1+ns(nn)*(i-1)))) mstr(ix+iw)=sur 8 msti(ix+iw)=sui write(6,*)'transformed' write(6,*)'second index' c transcript only do 11 I=1,nmax do 11 j=1,ns(nn) lj=ls(nn,j) sj=mu(nn,j)-1 do 11 mj=-lj,lj do 11 zj=-sj,sj,2 iw=3*((zj+SM)/2+ts*(mj+LM+tl*(j-1+ns(nn)*(i-1)))) tr(iw/3+1)=mstr(ix+iw) 11 ti(iw/3+1)=msti(ix+iw) write(6,*)'rewritten' c trafo: do 10 I=1,nmax do 10 J=1,nmax sur=0.0d0 sui=0.0d0 do 101 k=1,ns(nn) lk=ls(nn,k) sk=mu(nn,k)-1 do 101 mk=-lk,lk do 101 zk=-sk,sk,2 iw=3*((zk+SM)/2+ts*(mk+LM+tl*(k-1+ns(nn)*(i-1)))) sur=sur+tr(iw/3+1)*DR(J,k,mk+LM+1,(zk+SM)/2+1) 1 -ti(iw/3+1)*DI(J,k,mk+LM+1,(zk+SM)/2+1) 101 sui=sui+ti(iw/3+1)*DR(J,k,mk+LM+1,(zk+SM)/2+1) 1 +tr(iw/3+1)*DI(J,k,mk+LM+1,(zk+SM)/2+1) mstr(ix+3*((I-1)+nmax*(J-1)))=sur 10 msti(ix+3*((I-1)+nmax*(J-1)))=sui call rwi('w','MP.SCR',nmax,'Magnetic MO-P integrals',mstr,msti) endif call dospec(mstr,nb,nmax,'spectrum2.tab') return end c ============================================================ subroutine dospec(mao,nb,nmax,fn) implicit none integer*4 nb,i,nmax,ig real*8 mao(*),emin,rr,mx,my,mz,qpt,q,t,p,tcm,plim,de real*8,allocatable::e(:) character*(*) fn common/copt/t,plim allocate(e(nb)) call re(e,nb) tcm=t*0.695d0 emin=0.0d0 q=qpt(e,nb,emin) open(9,file=fn) write(9,900) 900 format(/,/) do 1 ig=1,nmax if((e(ig)-emin)/tcm.gt.40.0d0)then p=0.0d0 else p=exp(-(e(ig)-emin)/tcm)/q endif if(p.gt.plim)then do 3 i=1,nmax if(i.ne.ig)then mx=mao(1+3*((ig-1)+nmax*(i-1))) my=mao(2+3*((ig-1)+nmax*(i-1))) mz=mao(3+3*((ig-1)+nmax*(i-1))) rr=mx*mx+my*my+mz*mz if(dabs(p*rr).gt.1.0d-10)then de=e(i)-e(ig) if(de.lt.0.0d0)then de=-de rr=-rr endif if(dabs(de).gt.1.0d-10)write(9,901)i,de,p*rr,mx,my,mz,ig,i 901 format(i5,f12.2,4g12.3,i6,' ->',i5) endif endif 3 continue endif 1 continue write(9,902) 902 format(60(1h-)) close(9) write(6,*)fn//' written' return end c ============================================================ subroutine dospecf(egr,egi,mgr,mgi,nmax,fn) c absorption and CD implicit none integer*4 i,nmax,nt,n1,n2,j,nl real*8 egr(*),egi(*),mgr(*),mgi(*),rr,dd,est, 1mxr,myr,mzr,exr,eyr,ezr,si, 1mxi,myi,mzi,exi,eyi,ezi,de real*8,allocatable::e(:) character*(*) fn logical lex integer*4,allocatable::list(:) inquire(file='FLUOR.LST',exist=lex) if(lex)then write(6,*)'FLUOR.LST found' inquire(file='EST.TXT',exist=lex) if(lex)then open(9,file='EST.TXT') read(9,*)est si=-1.0d0 close(9) write(6,*)'Standard energy read from EST.TXT' else est=0.0d0 si=1.0d0 write(6,*)'Standard energy set to zero' endif allocate(e(nmax)) call re(e,nmax) open(9,file='FLUOR.LST') nt=0 1 read(9,*,end=99,err=99)i nt=nt+1 goto 1 99 close(9) write(6,*)nt,' transition requested' if(nt.gt.0)then open(10,file='FLUOR.LST') open(9,file=fn) write(9,900)est 900 format(' fluorescence, Est = ',f12.3,'cm-1',/,/) do 3 i=1,nt read(10,*)n1,nl allocate(list(nl)) backspace 10 read(10,*)n1,nl,(list(j),j=1,nl) do 31 j=1,nl n2=list(j) exr=egr(1+3*(n2-1+nmax*(n1-1))) eyr=egr(2+3*(n2-1+nmax*(n1-1))) ezr=egr(3+3*(n2-1+nmax*(n1-1))) mxr=mgr(1+3*(n2-1+nmax*(n1-1))) myr=mgr(2+3*(n2-1+nmax*(n1-1))) mzr=mgr(3+3*(n2-1+nmax*(n1-1))) exi=egi(1+3*(n2-1+nmax*(n1-1))) eyi=egi(2+3*(n2-1+nmax*(n1-1))) ezi=egi(3+3*(n2-1+nmax*(n1-1))) mxi=mgi(1+3*(n2-1+nmax*(n1-1))) myi=mgi(2+3*(n2-1+nmax*(n1-1))) mzi=mgi(3+3*(n2-1+nmax*(n1-1))) dd=exr*exr+eyr*eyr+ezr*ezr+exi*exi+eyi*eyi+ezi*ezi rr=exr*mxr+eyr*myr+ezr*mzr+exi*mxi+eyi*myi+ezi*mzi if(dabs(rr).lt.1.0d-10)rr=0.0d0 de=est+si*(e(n1)-e(n2)) 31 if(dabs(de).gt.1.0d-10)write(9,901)i,de,dd,rr,n1,n2 901 format(i5,f12.2,2g12.3,i5,' ->',i5) deallocate(list) 3 continue write(9,902) 902 format(60(1h-)) close(9) close(10) write(6,*)fn//' written' endif endif return end c ============================================================ subroutine dospecg(egr,egi,mgr,mgi,nmax,fn) c absorption and CD implicit none integer*4 i,nmax,ig real*8 egr(*),egi(*),mgr(*),mgi(*),emin,rr,dd, 1mxr,myr,mzr,exr,eyr,ezr, 1mxi,myi,mzi,exi,eyi,ezi,p real*8,allocatable::e(:) character*(*) fn real*8 t,plim,q,qpt,tcm,de common/copt/t,plim allocate(e(nmax)) call re(e,nmax) tcm=t*0.695d0 emin=0.0d0 q=qpt(e,nmax,emin) open(9,file=fn) write(9,900) 900 format(/,/) do 2 ig=1,nmax if((e(ig)-emin)/tcm.gt.40.0d0)then p=0.0d0 else p=exp(-(e(ig)-emin)/tcm)/q endif if(p.gt.plim)then do 3 i=1,nmax if(i.ne.ig)then exr=egr(1+3*((ig-1)+nmax*(i-1))) eyr=egr(2+3*((ig-1)+nmax*(i-1))) ezr=egr(3+3*((ig-1)+nmax*(i-1))) mxr=mgr(1+3*((ig-1)+nmax*(i-1))) myr=mgr(2+3*((ig-1)+nmax*(i-1))) mzr=mgr(3+3*((ig-1)+nmax*(i-1))) exi=egi(1+3*((ig-1)+nmax*(i-1))) eyi=egi(2+3*((ig-1)+nmax*(i-1))) ezi=egi(3+3*((ig-1)+nmax*(i-1))) mxi=mgi(1+3*((ig-1)+nmax*(i-1))) myi=mgi(2+3*((ig-1)+nmax*(i-1))) mzi=mgi(3+3*((ig-1)+nmax*(i-1))) dd=exr*exr+eyr*eyr+ezr*ezr+exi*exi+eyi*eyi+ezi*ezi rr=exr*mxr+eyr*myr+ezr*mzr+exi*mxi+eyi*myi+ezi*mzi if(dabs(rr).lt.1.0d-10)rr=0.0d0 de=e(i)-e(ig) if(de.lt.0.0d0)then de=-de rr=-rr dd=-dd endif if(dabs(de).gt.1.0d-2)write(9,901)i,de,p*dd,p*rr,ig,i 901 format(i5,f12.2,2g12.3,i10,' -> ',i5) endif 3 continue endif 2 continue write(9,902) 902 format(60(1h-)) close(9) write(6,*)fn//' written' return end c ============================================================ subroutine dospece(egr,egi,mgr,mgi,nmax,fn) c absorption and MCD implicit none integer*4 nmms parameter(nmms=5) integer*4 i,nmax,ie,immg(nmms),imme(nmms),iep,j,ig real*8 egr(*),egi(*),mgr(*),mgi(*),emin,gammaau,fjg, 1dsr,rsr1,rsr2,com1,com2,com3,ekg,ekj,fkj,fkg,dgjr(9),dgji(9), 1djkr(9),djki(9),rsr3,u1kgr,u1kgi,u2kgr,u2kgi,u3kgr,u3kgi,p0,p1, 1dif1r,dif2r,dif3r,dif1i,dif2i,dif3i,htcm,mjgr(3),mjgi(3), 1u1kjr,u1kji,u2kjr,u2kji,u3kjr,u3kji,ejg,cmmg(nmms),cmme(nmms), 1dkgr(9),dkgi(9),de1,de2,de3,dg1,dg2,dg3,ddu(3),rsr character*(*) fn real*8,allocatable::e(:) real*8 t,plim,tcm,p,q,qpt,de common/copt/t,plim gammaau=0.0001d0 htcm=219474.0d0 allocate(e(nmax)) call re(e,nmax) emin=0.0d0 q=qpt(e,nmax,emin) do 2 i=1,nmax 2 e(i)=e(i)/htcm emin=emin/htcm c temperature in hartrees: tcm=t*0.695d0/htcm call inispec(30,fn,'MCD spectrum in length formalism') do 1 ig=1,nmax if((e(ig)-emin)/tcm.gt.40.0d0)then p=0.0d0 else p=exp(-(e(ig)-emin)/tcm)/q endif c p plim p plim p plim p plim p plim p plim if(p.gt.plim)then do 304 ie=1,nmax c ie <> ig ie <> ig ie <> ig ie <> ig ie <> ig if(ie.ne.ig)then ejg=e(ie)-e(ig) if(dabs(ejg).gt.1.0d-10)then fjg=ejg/(ejg**2+gammaau**2) dsr=0.0d0 rsr1=0.0d0 rsr2=0.0d0 com1=0.0d0 com2=0.0d0 com3=0.0d0 c maximum contributions to ground and excited MCD: do 336 i=1,nmms immg(i)=0 imme(i)=0 cmmg(i)=0.0d0 336 cmme(i)=0.0d0 c transcript dipoles do 3191 j=1,3 dgjr(j)=egr(j+3*((ig-1)+nmax*(ie-1))) dgji(j)=egi(j+3*((ig-1)+nmax*(ie-1))) mjgr(j)=mgr(j+3*((ie-1)+nmax*(ig-1))) 3191 mjgi(j)=mgi(j+3*((ie-1)+nmax*(ig-1))) do 307 iep=1,nmax if(iep.ne.ie.and.iep.ne.ig)then ekg=e(iep) ekj=ekg-ejg fkj=ekj/(ekj**2+gammaau**2) fkg=ekg/(ekg**2+gammaau**2) c transcript dipoles do 319 j=1,3 dkgr(j )=egr(j+3*((iep-1)+nmax*(ig-1))) dkgi(j )=egi(j+3*((iep-1)+nmax*(ig-1))) dkgr(6+j)=mgr(j+3*((iep-1)+nmax*(ig-1))) dkgi(6+j)=mgi(j+3*((iep-1)+nmax*(ig-1))) djkr(j )=egr(j+3*((ie-1)+nmax*(iep-1))) djki(j )=egi(j+3*((ie-1)+nmax*(iep-1))) djkr(6+j)=mgr(j+3*((ie-1)+nmax*(iep-1))) 319 djki(6+j)=mgi(j+3*((ie-1)+nmax*(iep-1))) c ukg: x (=*) u1kgr= dgjr(2)*djkr(3)-dgjr(3)*djkr(2) 1 -(dgji(2)*djki(3)-dgji(3)*djki(2)) u1kgi= dgji(2)*djkr(3)-dgji(3)*djkr(2) 1 +dgjr(2)*djki(3)-dgjr(3)*djki(2) u2kgr= dgjr(3)*djkr(1)-dgjr(1)*djkr(3) 1 -(dgji(3)*djki(1)-dgji(1)*djki(3)) u2kgi= dgji(3)*djkr(1)-dgji(1)*djkr(3) 1 +dgjr(3)*djki(1)-dgjr(1)*djki(3) u3kgr= dgjr(1)*djkr(2)-dgjr(2)*djkr(1) 1 -(dgji(1)*djki(2)-dgji(2)*djki(1)) u3kgi= dgji(1)*djkr(2)-dgji(2)*djkr(1) 1 +dgjr(1)*djki(2)-dgjr(2)*djki(1) c ukj: u1kjr= dgjr(2)*dkgr(3)-dgjr(3)*dkgr(2) 1 -(dgji(2)*dkgi(3)-dgji(3)*dkgi(2)) u1kji= dgji(2)*dkgr(3)-dgji(3)*dkgr(2) 1 +dgjr(2)*dkgi(3)-dgjr(3)*dkgi(2) u2kjr=-(dgjr(3)*dkgr(1)-dgjr(1)*dkgr(3)) 1 +dgji(3)*dkgi(1)-dgji(1)*dkgi(3) u2kji= dgji(3)*dkgr(1)-dgji(1)*dkgr(3) 1 +dgjr(3)*dkgi(1)-dgjr(1)*dkgi(3) u3kjr= dgjr(1)*dkgr(2)-dgjr(2)*dkgr(1) 1 -(dgji(1)*dkgi(2)-dgji(2)*dkgi(1)) u3kji= dgji(1)*dkgr(2)-dgji(2)*dkgr(1) 1 +dgjr(1)*dkgi(2)-dgjr(2)*dkgi(1) c magnetic moment and dipole: p0=(dkgr(6+1)*u1kgr+dkgr(6+2)*u2kgr+dkgr(6+3)*u3kgr 1 +dkgi(6+1)*u1kgi+dkgi(6+2)*u2kgi+dkgi(6+3)*u3kgi)*fkg rsr1=rsr1+p0 com1=com1+p0 p1=(djkr(6+1)*u1kjr+djkr(6+2)*u2kjr+djkr(6+3)*u3kjr 1 +djki(6+1)*u1kji+djki(6+2)*u2kji+djki(6+3)*u3kji)*fkj rsr2=rsr2+p1 com2=com2+p1 c store the largest contributions: call dtm(iep,p0,immg,cmmg,nmms) call dtm(iep,p1,imme,cmme,nmms) endif c if(iep.ne.ie)then 307 continue c add parts dependent on permanent dipole moments: de1=egr(1+3*((ie-1)+nmax*(ie-1))) de2=egr(2+3*((ie-1)+nmax*(ie-1))) de3=egr(3+3*((ie-1)+nmax*(ie-1))) dg1=egr(1+3*((ig-1)+nmax*(ig-1))) dg2=egr(2+3*((ig-1)+nmax*(ig-1))) dg3=egr(3+3*((ig-1)+nmax*(ig-1))) c -: ddu(1)=de1-dg1 ddu(2)=de2-dg2 ddu(3)=de3-dg3 dif1r=dgjr(2)*ddu(3)-dgjr(3)*ddu(2) dif2r=dgjr(3)*ddu(1)-dgjr(1)*ddu(3) dif3r=dgjr(1)*ddu(2)-dgjr(2)*ddu(1) dif1i=dgji(2)*ddu(3)-dgji(3)*ddu(2) dif2i=dgji(3)*ddu(1)-dgji(1)*ddu(3) dif3i=dgji(1)*ddu(2)-dgji(2)*ddu(1) c .x(-): com3=(mjgr(1)*dif1r+mjgr(2)*dif2r+mjgr(3)*dif3r 1 +mjgi(1)*dif1i+mjgi(2)*dif2i+mjgi(3)*dif3i)*fjg rsr3=com3 dsr= dgjr(1)**2+dgjr(2)**2+dgjr(3)**2 1 +dgji(1)**2+dgji(2)**2+dgji(3)**2 rsr=(rsr1+rsr2+rsr3)/2.0d0 de=e(ie)-e(ig) if(de.lt.0.0d0)then de=-de dsr=-dsr rsr=-rsr endif write(30,3008)ie,htcm*de,p*dsr,p*rsr,ig,ie 3008 format(i4,f12.4,2g20.8,i10,' -> ',i5) endif endif c e-emin<>0 c ie <> ig ie <> ig ie <> ig ie <> ig ie <> ig 304 continue endif c p plim p plim p plim p plim p plim p plim 1 continue write(30,3002) 3002 format(80(1h-)) close(30) write(6,*)fn//' written' return end c ============================================================ subroutine dospecm(egr,egi,mgr,mgi,nmax,fn) c magnetic CPL (= MCD from excitedt states) implicit none integer*4 nmms parameter(nmms=5) integer*4 i,nmax,ie,immg(nmms),imme(nmms),iep,j,ig,nt,nl,ii, 1jj real*8 egr(*),egi(*),mgr(*),mgi(*),emin,gammaau,fjg, 1dsr,rsr1,rsr2,com1,com2,com3,ekg,ekj,fkj,fkg,dgjr(9),dgji(9), 1djkr(9),djki(9),rsr3,u1kgr,u1kgi,u2kgr,u2kgi,u3kgr,u3kgi,p0,p1, 1dif1r,dif2r,dif3r,dif1i,dif2i,dif3i,htcm,mjgr(3),mjgi(3), 1u1kjr,u1kji,u2kjr,u2kji,u3kjr,u3kji,ejg,cmmg(nmms),cmme(nmms), 1dkgr(9),dkgi(9),de1,de2,de3,dg1,dg2,dg3,ddu(3),rsr, 1si,est,de character*(*) fn real*8,allocatable::e(:) integer*4,allocatable::list(:) logical lex gammaau=0.0001d0 htcm=219474.0d0 inquire(file='FLUOR.LST',exist=lex) if(lex)then write(6,*)'FLUOR.LST found' inquire(file='EST.TXT',exist=lex) if(lex)then open(9,file='EST.TXT') read(9,*)est si=-1.0d0 close(9) write(6,*)'Standard energy read from EST.TXT' else est=0.0d0 si=1.0d0 write(6,*)'Standard energy set to zero' endif else return endif allocate(e(nmax)) call re(e,nmax) emin=0.0d0 write(6,*)'Calculating magnetic CPL' open(9,file='FLUOR.LST') nt=0 111 read(9,*,end=99,err=99)i nt=nt+1 goto 111 99 close(9) write(6,*)nt,' transition requested' if(nt.lt.1)return c all energies to hartree: do 2 i=1,nmax 2 e(i)=e(i)/htcm emin=emin/htcm call inispec(30,fn,'Magnetic CPL spectrum in length formalism') open(10,file='FLUOR.LST') do 1 ii=1,nt read(10,*)ie,nl allocate(list(nl)) backspace 10 read(10,*)ie,nl,(list(j),j=1,nl) do 304 jj=1,nl ig=list(jj) ejg=e(ie)-e(ig) fjg=ejg/(ejg**2+gammaau**2) dsr=0.0d0 rsr1=0.0d0 rsr2=0.0d0 com1=0.0d0 com2=0.0d0 com3=0.0d0 c maximum contributions to ground and excited MCD: do 336 i=1,nmms immg(i)=0 imme(i)=0 cmmg(i)=0.0d0 336 cmme(i)=0.0d0 c transcript dipoles do 3191 j=1,3 dgjr(j)=egr(j+3*((ig-1)+nmax*(ie-1))) dgji(j)=egi(j+3*((ig-1)+nmax*(ie-1))) mjgr(j)=mgr(j+3*((ie-1)+nmax*(ig-1))) 3191 mjgi(j)=mgi(j+3*((ie-1)+nmax*(ig-1))) do 307 iep=1,nmax if(iep.ne.ie.and.iep.ne.ig)then ekg=e(iep) ekj=ekg-ejg fkj=ekj/(ekj**2+gammaau**2) fkg=ekg/(ekg**2+gammaau**2) c transcript dipoles do 319 j=1,3 dkgr(j )=egr(j+3*((iep-1)+nmax*(ig-1))) dkgi(j )=egi(j+3*((iep-1)+nmax*(ig-1))) dkgr(6+j)=mgr(j+3*((iep-1)+nmax*(ig-1))) dkgi(6+j)=mgi(j+3*((iep-1)+nmax*(ig-1))) djkr(j )=egr(j+3*((ie-1)+nmax*(iep-1))) djki(j )=egi(j+3*((ie-1)+nmax*(iep-1))) djkr(6+j)=mgr(j+3*((ie-1)+nmax*(iep-1))) 319 djki(6+j)=mgi(j+3*((ie-1)+nmax*(iep-1))) c ukg: x (=*) u1kgr= dgjr(2)*djkr(3)-dgjr(3)*djkr(2) 1 -(dgji(2)*djki(3)-dgji(3)*djki(2)) u1kgi= dgji(2)*djkr(3)-dgji(3)*djkr(2) 1 +dgjr(2)*djki(3)-dgjr(3)*djki(2) u2kgr= dgjr(3)*djkr(1)-dgjr(1)*djkr(3) 1 -(dgji(3)*djki(1)-dgji(1)*djki(3)) u2kgi= dgji(3)*djkr(1)-dgji(1)*djkr(3) 1 +dgjr(3)*djki(1)-dgjr(1)*djki(3) u3kgr= dgjr(1)*djkr(2)-dgjr(2)*djkr(1) 1 -(dgji(1)*djki(2)-dgji(2)*djki(1)) u3kgi= dgji(1)*djkr(2)-dgji(2)*djkr(1) 1 +dgjr(1)*djki(2)-dgjr(2)*djki(1) c ukj: u1kjr= dgjr(2)*dkgr(3)-dgjr(3)*dkgr(2) 1 -(dgji(2)*dkgi(3)-dgji(3)*dkgi(2)) u1kji= dgji(2)*dkgr(3)-dgji(3)*dkgr(2) 1 +dgjr(2)*dkgi(3)-dgjr(3)*dkgi(2) u2kjr=-(dgjr(3)*dkgr(1)-dgjr(1)*dkgr(3)) 1 +dgji(3)*dkgi(1)-dgji(1)*dkgi(3) u2kji= dgji(3)*dkgr(1)-dgji(1)*dkgr(3) 1 +dgjr(3)*dkgi(1)-dgjr(1)*dkgi(3) u3kjr= dgjr(1)*dkgr(2)-dgjr(2)*dkgr(1) 1 -(dgji(1)*dkgi(2)-dgji(2)*dkgi(1)) u3kji= dgji(1)*dkgr(2)-dgji(2)*dkgr(1) 1 +dgjr(1)*dkgi(2)-dgjr(2)*dkgi(1) c magnetic moment and dipole: p0=(dkgr(6+1)*u1kgr+dkgr(6+2)*u2kgr+dkgr(6+3)*u3kgr 1 +dkgi(6+1)*u1kgi+dkgi(6+2)*u2kgi+dkgi(6+3)*u3kgi)*fkg rsr1=rsr1+p0 com1=com1+p0 p1=(djkr(6+1)*u1kjr+djkr(6+2)*u2kjr+djkr(6+3)*u3kjr 1 +djki(6+1)*u1kji+djki(6+2)*u2kji+djki(6+3)*u3kji)*fkj rsr2=rsr2+p1 com2=com2+p1 c store the largest contributions: call dtm(iep,p0,immg,cmmg,nmms) call dtm(iep,p1,imme,cmme,nmms) endif c if(iep.ne.ie)then 307 continue c add parts dependent on permanent dipole moments: de1=egr(1+3*((ie-1)+nmax*(ie-1))) de2=egr(2+3*((ie-1)+nmax*(ie-1))) de3=egr(3+3*((ie-1)+nmax*(ie-1))) dg1=egr(1+3*((ig-1)+nmax*(ig-1))) dg2=egr(2+3*((ig-1)+nmax*(ig-1))) dg3=egr(3+3*((ig-1)+nmax*(ig-1))) c -: ddu(1)=de1-dg1 ddu(2)=de2-dg2 ddu(3)=de3-dg3 dif1r=dgjr(2)*ddu(3)-dgjr(3)*ddu(2) dif2r=dgjr(3)*ddu(1)-dgjr(1)*ddu(3) dif3r=dgjr(1)*ddu(2)-dgjr(2)*ddu(1) dif1i=dgji(2)*ddu(3)-dgji(3)*ddu(2) dif2i=dgji(3)*ddu(1)-dgji(1)*ddu(3) dif3i=dgji(1)*ddu(2)-dgji(2)*ddu(1) c .x(-): com3=(mjgr(1)*dif1r+mjgr(2)*dif2r+mjgr(3)*dif3r 1 +mjgi(1)*dif1i+mjgi(2)*dif2i+mjgi(3)*dif3i)*fjg rsr3=com3 dsr= dgjr(1)**2+dgjr(2)**2+dgjr(3)**2 1 +dgji(1)**2+dgji(2)**2+dgji(3)**2 rsr=(rsr1+rsr2+rsr3)/2.0d0 de=e(ie)-e(ig) if(de.lt.0.0d0)then de=-de dsr=-dsr rsr=-rsr endif 304 write(30,3008)ie,est+si*htcm*de,dsr,rsr,ig,ie 3008 format(i4,f12.4,2g20.8,i10,' -> ',i5) deallocate(list) 1 continue write(30,3002) 3002 format(80(1h-)) close(30) write(6,*)fn//' written' return end c ============================================================ subroutine dtm(e,a,i,c,m) real*8 a,c(*) integer*4 e,i(*),m do 1 j=1,m if(dabs(a).gt.dabs(c(j)))then do 2 k=m,j+1,-1 i(k)=i(k-1) 2 c(k)=c(k-1) i(j)=e c(j)=a return endif 1 continue return end c ============================================================ subroutine trafo3(nb,maor,maoi,nmax,CR,CI) c = CIi* CJj implicit none integer*4 ix,nb,i,j,ii,nmax real*8 maor(*),maoi(*),sur,sui,CI(*),CR(*) character*2 mxyz(3) data mxyz/' x',' y',' z'/ real*8,allocatable::tr(:),ti(:) allocate(tr(nb**2),ti(nb**2)) do 10 ix=1,3 write(6,*)'transforming ',mxyz(ix) write(6,*)'first index' do 9 j=1,nb do 9 i=1,nb tr(i+(j-1)*nb)=maor(ix+3*(i-1+nb*(j-1))) 9 ti(i+(j-1)*nb)=maoi(ix+3*(i-1+nb*(j-1))) do 8 I=1,nmax do 8 j=1,nb sur=0.0d0 sui=0.0d0 do 81 ii=1,nb c = CIi* sur=sur+tr(ii+(j-1)*nb)*CR(ii+(I-1)*nb) 1 +ti(ii+(j-1)*nb)*CI(ii+(I-1)*nb) 81 sui=sui-tr(ii+(j-1)*nb)*CI(ii+(I-1)*nb) 1 +ti(ii+(j-1)*nb)*CR(ii+(I-1)*nb) maor(ix+3*(I-1+nb*(j-1)))=sur 8 maoi(ix+3*(I-1+nb*(j-1)))=sui write(6,*)'second index' do 11 I=1,nmax do 11 j=1,nb tr(I+(j-1)*nb)=maor(ix+3*(I-1+nb*(j-1))) 11 ti(I+(j-1)*nb)=maoi(ix+3*(I-1+nb*(j-1))) do 12 I=1,nmax do 12 J=1,nmax sur=0.0d0 sui=0.0d0 do 101 ii=1,nb c = CJj sur=sur+tr(I+nb*(ii-1))*CR(ii+(J-1)*nb) 1 -ti(I+nb*(ii-1))*CI(ii+(J-1)*nb) 101 sui=sui+tr(I+nb*(ii-1))*CI(ii+(J-1)*nb) 1 +ti(I+nb*(ii-1))*CR(ii+(J-1)*nb) maor(ix+3*(I-1+nmax*(J-1)))=sur 12 maoi(ix+3*(I-1+nmax*(J-1)))=sui 10 continue return end c ============================================================ subroutine odiff(li,ljo,d,idiff,sf,nn,a,ap) c find idiff = number of different orbital in the lists li and lj c d is maximal difference to be investigated c sf is the anti-symmetric phase sign c a,ap .. list of the different orbitals implicit none integer*4 it,li(*),ljo(*),d,idiff,ie,nn,iep,ic,id,iie,ii, 1a(*),ap(*) real*8 sf integer*4 ,allocatable::lj(:) c try to order lj according to li idiff=0 sf=1.0d0 allocate(lj(nn)) do 1 ie=1,nn 1 lj(ie)=ljo(ie) do 15 ie=1,nn c find orbital li(ie) in lj ic=0 do 16 iep=ie,nn 16 if(lj(iep).eq.li(ie))ic=iep if(ic.eq.0)then c this i-spinorbital is unique idiff=idiff+1 if(idiff.gt.d)return c put on its place ie in j also a unique one: id=0 do 23 iep=ie,nn do 24 iie=1,nn 24 if(lj(iep).eq.li(iie))goto 23 id=iep goto 25 23 continue 25 if(id.eq.0)call report('Did not find unique j') c move lj(id) to place ie do 26 ii=id,ie+1,-1 it=lj(ii) lj(ii)=lj(ii-1) lj(ii-1)=it 26 sf=-sf a(idiff)=li(ie) ap(idiff)=lj(ie) else c orbital lj(ic) exist in j same as li(ie): c move lj(ic) to place ie: do 17 ii=ic,ie+1,-1 it=lj(ii) lj(ii)=lj(ii-1) lj(ii-1)=it 17 sf=-sf endif 15 continue return end c ============================================================ function yyy(l1,m1,l2,m2,l3,m3) c implicit none real*8 yyy,pi,la,ma,lb,mb,lc,mc,rthrj,z0,one,two,four integer*4 l1,m1,l2,m2,l3,m3 data z0,one,two,pi,four/0.0d0,1.0d0,2.0d0, 13.14159265358979d0,4.0d0/ la=DBLE(l1) ma=DBLE(m1) lb=DBLE(l2) mb=DBLE(m2) lc=DBLE(l3) mc=DBLE(m3) c c int Y(*)l1m1 Yl2m2 Yl3m3 = dsqrt((2l1+1)(2l2+1)(2l3+1)/(4pi)) c c x (l1 l2 l3) (l1 l2 l3) x (-1)^m1 c (m1 m2 m3) (0 0 0 ) c yyy=dsqrt((two*la+one)*(two*lb+one)*(two*lc+one)/(four*pi)) 1*rthrj(la,lb,lc,ma,mb,mc)*rthrj(la,lb,lc,z0,z0,z0) if(mod(m1,2).ne.0)yyy=-yyy return end c ============================================================ subroutine inispec(i,f,s) implicit none integer*4 i character*(*)f,s open(i,file=f) write(i,*)' '//s write(i,3000) 3000 format(' n wavelength (nm) dipole strength (D^2)', 1' rotatory strength (cgs/10**-36)',/,80(1h-)) return end c ============================================================ subroutine re(e,n) integer*4 i,n real*8 e(*) open(9,file='energies') do 1 i=1,n 1 read(9,*)e(i) close(9) return end c ============================================================ subroutine readopt(lex,kev,nmax,pexc,lpolar,lopt) integer*4 kev,nmax logical lex,lpolar,lopt(*) real*8 pexc(4) character*5 s5 real*8 r2,r4,r6,t,plim common/rs/r2,r4,r6 c pexc - weighting factors of the f^n ->f^n-1 d pexc(1) c d^10 f^n ->d^9 f^n+1 pexc(2) c f^n ->f^n-1 g pexc(3) c excitations common/copt/t,plim c temperature in Kelvins: t=300.0d0 plim=0.01d0 do 2 i=1,10 2 lopt(i)=.false. lopt(1)=.true. lopt(2)=.true. lopt(3)=.true. c lopt(1) ... , integrals c lopt(2) ... magnetic moments - direct c lopt(3) ... electric moments - direct c lopt(4) ... magnetic moments - Racah r2=0.2648d0 r4=0.1678d0 r6=0.2123d0 pexc(1)=1.0d0 pexc(2)=0.0d0 pexc(3)=0.0d0 pexc(4)=0.0d0 kev =1 nmax=1 lpolar=.false. inquire(file='ANAL.OPT',exist=lex) if(lex)then open(9,FILE='ANAL.OPT') 1 read(9,90,end=99,err=99)s5 90 format(a5) if(s5(1:3).eq.'KEV')read(9,*)kev if(s5(1:4).eq.'NMAX')read(9,*)nmax if(s5(1:4).eq.'PEXC')read(9,*)i,pexc(i) if(s5(1:5).eq.'POLAR')read(9,*)lpolar if(s5(1:4).eq.'R246')read(9,*)r2,r4,r6 if(s5(1:4).eq.'L2S2')read(9,*)lopt(1) if(s5(1:4).eq.'MAGN')read(9,*)lopt(2) if(s5(1:4).eq.'ELEC')read(9,*)lopt(3) if(s5(1:4).eq.'RACA')read(9,*)lopt(4) if(s5(1:4).eq.'TEMP')read(9,*)t if(s5(1:4).eq.'PLIM')read(9,*)plim goto 1 99 close(9) write(6,*)'ANAL.OPT found' endif return end c ============================================================ function rt(N) integer*4 N real*8 rt real*8 r2,r4,r6 common/rs/r2,r4,r6 c rt=1.0d0 if(N.gt.0.and.N.le.3)rt=r2 if(N.gt.3.and.N.le.5)rt=r4 if(N.gt.5)rt=r6 return end c ============================================================ subroutine uvp(Y1z,Y1m,Y1p,vr,vi,l1,m1,l3,m3,almr,almi, 1polr,poli,c,nic,kic,mic,lpolar) c elements c elements c elements c elements , V is crystal field potential implicit none real*8 Y1z,Y1m,Y1p,yyy,fac,pi,pi4,vr,vi,almr(*),almi(*),rt, 1polr(*),poli(*),c(*) integer*4 l1,m1,l3,m3,ia,ial,iam,nic(*),kic(*),mic(*) logical lpolar Y1z=yyy(l1,m1,1, 0,l3,m3) Y1m=yyy(l1,m1,1,-1,l3,m3) Y1p=yyy(l1,m1,1, 1,l3,m3) pi=4.0d0*atan(1.0d0) pi4=4.0d0*pi ia=0 do 5 ial=0,7 c Clm = sqrt(4pi/(2l+1)) Ylm fac=dsqrt(pi4/dble(2*ial+1)) do 5 iam=-ial,ial ia=ia+1 c = vr=vr+yyy(l1,m1,ial,-iam,l3,m3)*fac*almr(ia)*rt(ial) 5 vi=vi-yyy(l1,m1,ial,-iam,l3,m3)*fac*almi(ia)*rt(ial) if(lpolar)then do 501 ia=1,120 c Clm = sqrt(4pi/(2l+1)) Ylm fac=dsqrt(pi4/dble(2*kic(ia)+1)) c = vr=vr+yyy(l1,m1,kic(ia),-mic(ia),l3,m3)*fac*polr(ia)* 1 c(ia)*rt(nic(ia)) 501 vi=vi-yyy(l1,m1,kic(ia),-mic(ia),l3,m3)*fac*poli(ia)* 1 c(ia)*rt(nic(ia)) endif return end c ============================================================ function timemy() integer*4 timemy,time8 timemy=0 timemy=time8() return end c ============================================================ function qpt(e,N,emin) c partition function real*8 e(*),t,qpt,s,f,tcm,emin,plim integer*4 i,N common/copt/t,plim emin=e(1) do 2 i=2,N 2 if(e(i).lt.emin) emin=e(i) s=0.0d0 c temperature in cm-1: tcm=t*0.695d0 do 1 i=1,N f=(e(i)-emin)/tcm 1 if(f.lt.40.0d0)s=s+exp(-f) qpt=s return end c ============================================================