program cff implicit none integer*4 l,m,i,ila,ild,n,lmax,kmax,nc,n0 real*8 t,f,xa,ya,za,xx,yy,zz,r,ar,ai,sgn,rect,imct real*8,allocatable::x(:),q(:),c(:),pol(:) integer*4,allocatable::ty(:), 1nic(:),nlc(:),kic(:),mic(:),mlc(:),klc(:) logical lpolar c write(6,600) 600 format(/,' Crystal field perturbation parameters',/,/, 1 ' (for the lanthanide program)',/,/, 1 ' Input: G.OUT ',/, 1 ' CFF.OPT (optional) ',/, 1 ' POL.LST (optional) ',/, 2 ' output: cfp.inp (charges)',/, 2 ' cfpp.inp (polarizabilities)',/) c Read options if present call readopt(lpolar,lmax,kmax,n0) c initialize single and double factorials: call InitFac c read number of atoms: call rd(0,x,q,n,ty) allocate(q(n),x(3*n),ty(n)) c read geometry and charges: call rd(1,x,q,n,ty) c transform to atomic units: call au(x,n) c identify Lanthanide: ila=ild(ty,n) if(ila.eq.0)call report(' No lanthanide found') if(ila.eq.-1)call report(' More then one lanthanide found') xa=x(1+3*(ila-1)) ya=x(2+3*(ila-1)) za=x(3+3*(ila-1)) open(9,file='cfp.inp') do 1 l=0,7 do 1 m=-l,l ar=0.0d0 ai=0.0d0 do 2 i=1,n if(i.ne.ila)then xx=x(1+3*(i-1))-xa yy=x(2+3*(i-1))-ya zz=x(3+3*(i-1))-za r=dsqrt(xx**2+yy**2+zz**2) t=acos(zz/r) f=atan2(yy,xx) call ct(l,-m,t,f,rect,imct) ar=ar+q(i)/r**(l+1)*rect ai=ai+q(i)/r**(l+1)*imct endif 2 continue if(mod(m+1,2).eq.0)then sgn= 1.0d0 else sgn=-1.0d0 endif ar=ar*sgn*219474.0d0 ai=ai*sgn*219474.0d0 1 write(9,900)l,m,ar,ai 900 format(2i3,2f12.3) close(9) write(6,601)n,ila,ty(ila) 601 format(i10,' atoms',/,i10,i4,' center',/,' cfp.inp written') if(lpolar)then allocate(nic(n0),nlc(n0),kic(n0),mic(n0),mlc(n0),klc(n0),c(n0), 1 pol(n)) write(6,*)'Polarizability parameters requested' c prepare r^-4 polynomial expansion call r4m(n0,nc,c ,nic,nlc,kic,mic,klc,mlc,lmax,kmax) c prepare polarizabilities call getpol(n,pol,ty) open(9,file='cfpp.inp') do 3 l=1,nc ar=0.0d0 ai=0.0d0 do 4 i=1,n if(i.ne.ila)then xx=x(1+3*(i-1))-xa yy=x(2+3*(i-1))-ya zz=x(3+3*(i-1))-za r=dsqrt(xx**2+yy**2+zz**2) t=acos(zz/r) f=atan2(yy,xx) call ct(klc(l),mlc(l),t,f,rect,imct) ar=ar-pol(i)/r**nlc(l)*rect ai=ai-pol(i)/r**nlc(l)*imct endif 4 continue ar=ar*219474.0d0 ai=ai*219474.0d0 3 write(9,944)l,nic(l),nlc(l),kic(l),mic(l),klc(l),mlc(l),c(l), 1 ar,ai 944 format(i5,3(1x,2i3),1x,3f15.6) close(9) write(6,604)n,ila,ty(ila) 604 format(i10,' atoms',/,i10,i4,' center',/,' cfpp.inp written') endif end function ild(ty,n) implicit none integer*4 ild,n,ty(*),i,ii,ic ii=0 do 1 i=1,n if(ty(i).ge.57.and.ty(i).le.70)then ii=ii+1 ic=i endif 1 continue if(ii.gt.1)then ild=-1 else if(ii.eq.0)then ild=0 else ild=ic endif endif return end subroutine au(x,n) implicit none real*8 x(*) integer*4 n,i do 1 i=1,3*n 1 x(i)=x(i)/0.529177d0 return end subroutine rd(ic,x,q,n,ty) implicit none integer*4 ic,n,ty(*),ig98,I,l,KA,NG real*8 x(*),q(*) character*50 FN logical lex NG=0 inquire(file='G.OUT',exist=lex) c ------------------------------------------------------------- if(lex)then open(2,file='G.OUT',status='old') 1 READ(2,2000,END=1000,ERR=1000)FN 2000 FORMAT(A50) IF(FN(19:39).EQ.'Z-Matrix orientation:'.OR. 1 FN(26:46).EQ.'Z-Matrix orientation:'.OR. 1 FN(20:37).EQ.'Input orientation:'.OR. 1 FN(27:44).EQ.'Input orientation:'.OR. 2 FN(20:40).EQ.'Standard orientation:'.OR. 2 FN(26:46).EQ.'Standard orientation:')THEN NG=NG+1 ig98=0 if(FN(26:46).EQ.'Z-Matrix orientation:'.OR. 1 FN(27:44).EQ.'Input orientation:'.OR. 1 FN(26:46).EQ.'Standard orientation:')ig98=1 DO 4 I=1,4 4 READ(2,*) l=0 5 READ(2,2000)FN IF(FN(2:4).NE.'---')THEN l=l+1 if(ic.eq.1)then BACKSPACE 2 if(ig98.eq.0)then READ(2,*)KA,KA,(x(i+3*(l-1)),i=1,3) else READ(2,*)KA,KA,x(1+3*(l-1)),(x(i+3*(l-1)),i=1,3) endif ty(l)=KA endif GOTO 5 ENDIF n=l if(ic.eq.0)goto 1000 ENDIF IF(FN(2:17).EQ.'Mulliken charges'.and.ic.eq.1)then read(2,*) do 6 i=1,n READ(2,2000)FN 6 read(FN(11:21),*)q(i) write(6,*)'Mullican charges found' goto 1000 ENDIF goto 1 1000 CLOSE(2) write(6,*)NG,' geometries in G.OUT' c ------------------------------------------------------------- else open(2,file='FILE.X',status='old') read(2,*) read(2,*)n if(ic.eq.0)goto 2222 do 7 l=1,n 7 READ(2,*)ty(l),(x(i+3*(l-1)),i=1,3) write(6,*)n,' atoms in FILE.X' do 8 i=1,n q(i)=0.0d0 if(ty(i).eq.1)q(i)= 0.39d0 8 if(ty(i).eq.8)q(i)=-0.78d0 2222 CLOSE(2) endif c ------------------------------------------------------------- return end subroutine sph(l,m,t,f,re,im) c spherical function Yl,m(t,f) ... r, im, real and imaginaty parts implicit none integer*4 l,m real*8 t,f,re,im,n,norm,pol,p norm=n(l,m) pol=p(m,l,t) re=norm*pol*cos(dble(m)*f) im=norm*pol*sin(dble(m)*f) return end subroutine ct(l,m,t,f,rect,imct) c tensor operators implicit none integer*4 l,m real*8 t,f,rect,imct,re,im,pi,fac pi=4.0d0*datan(1.0d0) call sph(l,m,t,f,re,im) fac=dsqrt(4.0d0*pi/dble(2*l+1)) rect=fac*re imct=fac*im return end subroutine report(s) character*(*) s write(6,*)s stop end function p(m,l,t) c associated Legendre polynomial implicit none integer*4 nfac0 parameter (nfac0=100) integer*4 l,m,ll,i,mp real*8 t,z,p,sgn,pmm,pmmp,po real*8 dm,d0,df common/dpar/dm,d0,df(nfac0) real*8 f0,ff common/npar/f0,ff(nfac0) real*8,allocatable::pp(:) if(l.eq.0.and.m.eq.0)then p=1.0d0 return endif z=cos(t) c pp is p for positive m: mp=abs(m) c Pm,m pmm=df(2*mp-1)*(1-z**2)**(dble(mp)/2.0d0) if(mp.eq.l)then po=pmm else c Pm,m+1 pmmp=z*dble(2*mp+1)*pmm if(mp+1.eq.l)then po=pmmp else allocate(pp(l-mp+1)) pp(1)=pmm pp(2)=pmmp do 1 ll=mp+2,l i=ll-mp+1 1 pp(i)= 1 (z*dble(2*ll-1)*pp(i-1)-dble(ll+mp-1)*pp(i-2))/dble(ll-mp) po=pp(l-mp+1) endif endif if(m.gt.0)then p=po else if(mod(mp,2).eq.0)then sgn= 1.0d0 else sgn=-1.0d0 endif p=sgn*ff(l-mp)/ff(l+mp)*po endif return end function n(l,m) implicit none real*8 n,fp integer*4 l,m integer*4 nfac0 parameter (nfac0=100) real*8 f0,factorial common/npar/f0,factorial(nfac0) fp=4.0d0*4.0d0*datan(1.0d0) c n=dsqrt(dble(2*l+1)*factorial(l-abs(m))/fp/factorial(l+abs(m))) n=dsqrt(dble(2*l+1)*factorial(l-m)/fp/factorial(l+m)) 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 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 subroutine readopt(lpolar,lmax,kmax,n0) implicit none logical lpolar,lex integer*4 lmax,n0,kmax character*5 s5 lpolar=.false. kmax=7 lmax=7 n0=1000 inquire(file='CFF.OPT',exist=lex) if(lex)then write(6,*)'CFF.OPT read' open(9,file='CFF.OPT') 1 read(9,90,end=99,err=99)s5 write(6,90)s5 90 format(a5) if(s5.eq.'POLAR')read(9,*)lpolar if(s5(1:4).eq.'KMAX')read(9,*)kmax if(s5(1:4).eq.'LMAX')read(9,*)lmax if(s5(1:2).eq.'N0')read(9,*)n0 goto 1 99 close(9) endif return end c ============================================================ subroutine r4m(n0,nc,c ,nic,nlc,kic,mic,klc,mlc,lmax,kmax) implicit none integer*4 i,n0,kmax,lmax,k,m,n ,np, 1nc,nic(*),nlc(*),kic(*),mic(*),mlc(*),klc(*) real*8 pi,c(*) real*8,allocatable:: p(:),pp(:) integer*4,allocatable:: 1ni( :),nl( :),ki( :),mi( :),ml( :),kl( :), 1nip(:),nlp(:),kip(:),mip(:),mlp(:),klp(:) pi=4.0d0*atan(1.0d0) c sum(k=0..7,m=-k,k)4pi(2k+1) ri^ni/rl(^nl+1)Yki,miYklml c Ykm - spherical functions allocate(ni( n0),nl( n0),ki( n0),mi( n0),ml( n0),kl( n0),p (n0), 1 nip(n0),nlp(n0),kip(n0),mip(n0),mlp(n0),klp(n0),pp(n0)) write(6,*)'preparing 1/r^4 expansion' i=0 do 1 k=0,kmax do 1 m=-k,k i=i+1 if(i.gt.n0)call report('too many terms') p(i)=4.0d0*pi/dble(2*k+1) if(mod(m,2).ne.0)p(i)=-p(i) ni(i)=k nl(i)=k+1 ki(i)=k mi(i)=m kl(i)=k 1 ml(i)=-m n=i write(6,*)n,' terms' call pq( n0,n ,p ,ni ,nl, ki, mi, kl, ml, 1 n ,p ,ni ,nl, ki, mi, kl, ml, 1 np,pp,nip,nlp,kip,mip,klp,mlp,lmax,kmax) write(6,*)'squared -',np,' terms' call pq( n0,np,pp,nip,nlp,kip,mip,klp,mlp, 1 np,pp,nip,nlp,kip,mip,klp,mlp, 1 nc,c ,nic,nlc,kic,mic,klc,mlc,lmax,kmax) write(6,*)'fourth power -',nc,' terms' return end c ============================================================ subroutine pq(n0,np,p ,nip,nlp,kip,mip,klp,mlp, 1 nq,q ,niq,nlq,kiq,miq,klq,mlq, 1 k,s ,nis,nls,kis,mis,kls,mls,lmax,kmax) c s = p x q ;special polynoms: c p = sum (j) ri^nij/rL^nLj Ykjmj(i)YkLjmL(j(L) integer*4 n0,lmax,i,j,lp,mp,lq,mq,kk,kmax, 1np,nip(*),nlp(*),kip(*),mip(*),klp(*),mlp(*), 1nq,niq(*),nlq(*),kiq(*),miq(*),klq(*),mlq(*), 1 k,nis(*),nls(*),kis(*),mis(*),kls(*),mls(*) real*8 p(*),q(*),s(*),c,yyy,cpq,tol c do 1 i=1,n0 nis(i)=0 nls(i)=0 kis(i)=0 mis(i)=0 kls(i)=0 mls(i)=0 1 s(i)=0.0d0 c generate all possible Yi x YL terms and digest k=0 do 7 i=1,np do 7 j=1,nq c Y(i)p x Y(i)q = sum YYY Ylm(i), etc cpq=p(i)*q(j) if(cpq.ne.0.0d0.and.nip(i)+niq(j).le.lmax)then do 71 lp=0,kmax do 71 mp=-lp,lp do 71 lq=0,kmax do 71 mq=-lq,lq c=yyy(lp,mp,kip(i),mip(i),kiq(j),miq(j)) 1 * yyy(lq,mq,klp(i),mlp(i),klq(j),mlq(j)) if(c.ne.0.0d0)then c try to find this term: do 8 kk=1,k if(nis(kk).eq.nip(i)+niq(j).and. 1 nls(kk).eq.nlp(i)+nlq(j).and. 1 kis(kk).eq.lp.and. 1 mis(kk).eq.mp.and. 1 kls(kk).eq.lq.and. 1 mls(kk).eq.mq)then c this term exists, add it and go to next: s(kk)=s(kk)+cpq*c goto 71 endif 8 continue c this term was not found, add it: k=k+1 if(k.gt.n0)call report('too many terms - Y stage') nis(k)=nip(i)+niq(j) nls(k)=nlp(i)+nlq(j) kis(k)=lp mis(k)=mp kls(k)=lq mls(k)=mq s(kk)=cpq*c endif 71 continue endif 7 continue write(6,*)k,' final unique terms' c delete small ones tol=1.0d-5 99 do 9 i=1,k if(dabs(s(i)).lt.tol)then do 91 kk=i,k-1 s(kk)=s(kk+1) nis(kk)=nis(kk+1) nls(kk)=nls(kk+1) kis(kk)=kis(kk+1) mis(kk)=mis(kk+1) kls(kk)=kls(kk+1) 91 mls(kk)=mls(kk+1) k=k-1 goto 99 endif 9 continue write(6,*)k,' large terms left' 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 ============================================================ 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 getpol(n,pol,ty) implicit none integer*4 n,ty(*),i logical lex real*8 pol(*) inquire(file='POL.LST',exist=lex) if(lex)then write(6,*)' Polarizabilities read from POL.LST.' open(9,file='POL.LST') do 1 i=1,n 1 read(9,*)pol(i),pol(i) close(9) else write(6,*)' Polarizabilities roughly estimated' do 2 i=1,n if(ty(i).eq.1)then pol(i)=1.0d0 else pol(i)=8.0d0 endif 2 continue endif return end