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,Rlim,Rlii, 1Xlim,Xlii real*8,allocatable::x(:),q(:),c(:),pol(:) integer*4,allocatable::ty(:), 1nic(:),nlc(:),kic(:),mic(:),mlc(:),klc(:) logical lpolar,les,lro,lex 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) ',/, 1 ' LAXYZ (optional) ',/, 1 ' density.cub (optional) ',/, 1 ' potential.cub (optional) ',/, 2 ' output: cfp.inp (V from charges)',/, 2 ' cfpp.inp (V from polarizabilities)',/, 2 ' cfq.inp (V from pot)',/, 2 ' cfq2.inp (V from pot2)',/, 2 ' cfro.inp (V from density)',/) c Read options if present call readopt(lpolar,lmax,kmax,n0,les,Rlim,Rlii,Xlim,Xlii,lro) 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.-1)call report(' More then one lanthanide found') if(ila.eq.0)then write(6,*)' No lanthanide found - virtual one put to origin' xa=0.0d0 ya=0.0d0 za=0.0d0 inquire(file='LAXYZ',exist=lex) if(lex)then write(6,*)'LAXYZ found, coordinates read from it.' open(9,file='LAXYZ') read(9,*)xa,xa,ya,za xa=xa/0.529177d0 ya=ya/0.529177d0 za=za/0.529177d0 close(9) endif else write(6,*)' Lanthanide: atom number ',ila,' type:', 1 ty(ila) xa=x(1+3*(ila-1)) ya=x(2+3*(ila-1)) za=x(3+3*(ila-1)) endif write(6,*) write(6,*)'Electrosctatic potential from charges' 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 601 format(i10,' atoms',/,i10,' 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,*) 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 604 format(i10,' atoms',/,i10,' center',/,' cfpp.inp written') endif if(lro)call ror(n,xa,ya,za,Xlim,Xlii,x,ty) if(les)call cles(n,xa,ya,za,Rlim,Rlii) if(les)call cles2(n,xa,ya,za,Rlim,Rlii) 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(*) real*8 qt character*50 FN logical lex NG=0 if(ic.eq.1)then do 7 i=1,n 7 q(i)=0.0d0 endif 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:37).EQ.'Mulliken charges and spin densities:'.or. 1 FN(2:25).EQ.'Mulliken atomic charges:'.or. 1 FN(2:18).EQ.'Mulliken charges:'). 1 and.ic.eq.1)then qt=0 read(2,*) do 6 i=1,n READ(2,2000)FN read(FN(11:21),*)q(i) 6 qt=qt+q(i) write(6,600)qt 600 format(' Mullican charges found, sum = ',f10.2) 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 8 l=1,n 8 READ(2,*)ty(l),(x(i+3*(l-1)),i=1,3) write(6,*)n,' atoms in FILE.X' do 9 i=1,n q(i)=0.0d0 if(ty(i).eq.1)q(i)= 0.39d0 9 if(ty(i).eq.8)q(i)=-0.78d0 2222 CLOSE(2) endif c ------------------------------------------------------------- return end c ============================================================ 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 c ============================================================ 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 c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ 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,les,Rlim,Rlii,Xlim, 1Xlii,lro) implicit none logical lpolar,lex,les,lro integer*4 lmax,n0,kmax real*8 Rlim,Rlii,Xlim,Xlii character*5 s5 lpolar=.false. les=.false. lro=.false. kmax=7 lmax=7 n0=1000 Rlii=0.1d0 Rlim=1.0d0 Xlii=1.8d0 Xlim=100.0d0 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.eq.'POTEN')read(9,*)les if(s5.eq.'DENSI')read(9,*)lro if(s5.eq.'RMAXI')read(9,*)Rlim if(s5.eq.'RMINI')read(9,*)Rlii if(s5.eq.'XMAXI')read(9,*)Xlim if(s5.eq.'XMINI')read(9,*)Xlii 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 c ============================================================ 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 c ============================================================ subroutine wmtr(A,n0,n,st,ic) 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(n0,n) character*(*) st write(6,*)st write(6,*) N1=1 1 N3=MIN(N1+4,N) WRITE(6,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(6,17)LN,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(6,*) return end c ============================================================ SUBROUTINE INV(a,ai,n) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(n,n),a(n,n) real*8,allocatable::e(:,:) TOL=1.0D-10 allocate(e(n,2*n)) C 10000 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 call wmtr(a,n,n,' original matrix:',1) call report('cannot invert matrix') c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END c ============================================================ subroutine EMA(A,n) c negative matrix A = -A implicit none integer*4 i,j,n real*8 A(n,n) do 1 i=1,n do 1 j=1,n 1 A(i,j)=-A(i,j) return end c ============================================================ subroutine EMV(A,n) c negative vector A = -A implicit none integer*4 i,n real*8 A(*) do 1 i=1,n 1 A(i)=-A(i) return end c ============================================================ subroutine MMA(A,B,C,n) c matrix multiplication A = B x C implicit none integer*4 i,j,k,n real*8 t,A(n,n),B(n,n),C(n,n) do 1 i=1,n do 1 j=1,n t=0.0d0 do 2 k=1,n 2 t=t+B(i,k)*C(k,j) 1 A(i,j)=t return end c ============================================================ subroutine MMV(A,B,C,n) c matrix vector multiplication A = B x C implicit none integer*4 i,k,n real*8 t,A(*),B(n,n),C(*) do 1 i=1,n t=0.0d0 do 2 k=1,n 2 t=t+B(i,k)*C(k) 1 A(i)=t return end c ============================================================ subroutine AMA(A,B,C,n) c matrix addition A = B + C implicit none integer*4 i,j,n real*8 A(n,n),B(n,n),C(n,n) do 1 i=1,n do 1 j=1,n 1 A(i,j)=B(i,j)+C(i,j) return end c ============================================================ subroutine AMV(A,B,C,n) c vector addition A = B + C implicit none integer*4 i,n real*8 A(n),B(n),C(n) do 1 i=1,n 1 A(i)=B(i)+C(i) return end c ============================================================ subroutine cles(n,xa,ya,za,Rlim,Rlii) implicit none integer*4 nat,i1,i2,i3,n,md,nc(3),l,m,lp,mp,i,j real*8 r0(3),rc(3,3),xa,ya,za,rect,imct,rectp,imctp, 1r,t,f,xc,yc,zc,xx,yy,zz,rl,rlp,Rlim,Rlimau,wai,war, 1Rlii,Rliiau,fi3 real*8,allocatable::fi(:),cr(:),ci(:),mi(:,:),mr(:,:), 1iri(:,:),t1(:),t2(:),ar(:),ai(:),rr(:,:),ri(:,:), 1b(:,:),ir(:,:) write(6,*) write(6,*)'ES parameters from potential on grid requested' c c setup matrix dimension md=0 do 1 l=0,7 do 1 m=-l,l 1 md=md+1 allocate(mr(md,md),cr(md),mi(md,md),ci(md), 1iri(md,md),t1(md),t2(md),ar(md),ai(md),rr(md,md),ri(md,md), 1b(md,md),ir(md,md)) do 2 i=1,md ci(i)=0.0d0 cr(i)=0.0d0 do 2 j=1,md mr(i,j)=0.0d0 2 mi(i,j)=0.0d0 Rlimau=Rlim/0.529177d0 Rliiau=Rlii/0.529177d0 c read gaussian cub file: write(6,*)'Reading potential.cub' open(9,file='potential.cub',status='old') read(9,*) read(9,*) c r0 and rc supposed to be in bohrs, no conversion now: read(9,*)nat,r0(1),r0(2),r0(3) if(nat.ne.n)call report('Number of atoms does not match.') read(9,*)nc(1),(rc(1,i),i=1,3) read(9,*)nc(2),(rc(2,i),i=1,3) read(9,*)nc(3),(rc(3,i),i=1,3) write(6,600)nc(1)*nc(2)*nc(3),Rlii,Rlim 600 format(i10,' grid points, R_min =',f6.2,' R_max =',f6.2) do 5 i=1,n 5 read(9,*) allocate(fi(nc(3))) do 6 i1=1,nc(1) do 6 i2=1,nc(2) read(9,4000)(fi(i3),i3=1,nc(3)) 4000 format(6E13.5) do 6 i3=1,nc(3) xc=r0(1)+rc(1,1)*dble(i1-1)+rc(2,1)*dble(i2-1)+rc(3,1)*dble(i3-1) yc=r0(2)+rc(1,2)*dble(i1-1)+rc(2,2)*dble(i2-1)+rc(3,2)*dble(i3-1) zc=r0(3)+rc(1,3)*dble(i1-1)+rc(2,3)*dble(i2-1)+rc(3,3)*dble(i3-1) xx=xc-xa yy=yc-ya zz=zc-za r=dsqrt(xx**2+yy**2+zz**2) if(r.le.Rlimau.and.r.ge.Rliiau)then t=acos(zz/r) f=atan2(yy,xx) fi3=fi(i3) i=0 do 61 l=0,7 rl=r**l*fi3 do 61 m=-l,l i=i+1 call ct(l,m,t,f,rect,imct) cr(i)=cr(i)+rl*rect ci(i)=ci(i)+rl*imct j=0 do 61 lp=0,7 rlp=r**(l+lp) do 61 mp=-lp,lp j=j+1 call ct(lp,mp,t,f,rectp,imctp) mr(i,j)=mr(i,j)+rlp*(rect*rectp-imct*imctp) 61 mi(i,j)=mi(i,j)+rlp*(rect*imctp+imct*rectp) endif 6 continue close(9) write(6,*)' .... done.' c M = Mr + i Mi c M^-1 = rr + i ri c M^-1 x M = 1 c ri=Mr^-1: call INV(mr,ri,md) c ir = Mi x Mr^-1 : call MMA(ir,mi,ri,md) c iri = ir x Mi = Mi x Mr^-1 x Mi: call MMA(iri,ir,mi,md) c B = Mr + Mi x Mr^-1 x Mi: call AMA(B,mr,iri,md) c rr = B^-1 = (Mr + Mi x Mr^-1 x Mi)^-1: call INV(B,rr,md) c ri = - rr x ir = - rr x Mi x Mr^-1 call MMA(ri,rr,ir,md) call EMA(ri,md) c A = M^-1 x c = Ar + i Ai c Ar = rr x cr - ri x ci c Ai = ri x cr + rr x ci call MMV(t1,rr,cr,md) call MMV(t2,ri,ci,md) call EMV(t2,md) call AMV(ar,t1,t2,md) call MMV(t1,ri,cr,md) call MMV(t2,rr,ci,md) call AMV(ai,t1,t2,md) open(9,file='cfq.inp') i=0 do 7 l=0,7 do 7 m=-l,l i=i+1 war=-ar(i)*219474.0d0 wai=-ai(i)*219474.0d0 7 write(9,900)l,m,war,wai 900 format(2i3,2f12.3) close(9) write(6,601) 601 format(' Grid potential transformed to cfq.inp.') return end c ============================================================ subroutine cles2(n,xa,ya,za,Rlim,Rlii) implicit none integer*4 nat,i1,i2,i3,n,md,nc(3),l,m,lp,mp,i,j real*8 r0(3),rc(3,3),xa,ya,za,rect,imct,rectp,imctp, 1r,t,f,xc,yc,zc,xx,yy,zz,rl,rlp,Rlim,Rlimau,wai,war, 1Rlii,Rliiau,fi3 real*8,allocatable::fi(:),cr(:),ci(:),mi(:,:),mr(:,:), 1iri(:,:),t1(:),t2(:),ar(:),ai(:),rr(:,:),ri(:,:), 1b(:,:),ir(:,:) write(6,*) write(6,*)'ES parameters from potential on grid requested' c c setup matrix dimension md=0 do 1 l=0,7 do 1 m=-l,l 1 md=md+1 allocate(mr(md,md),cr(md),mi(md,md),ci(md), 1iri(md,md),t1(md),t2(md),ar(md),ai(md),rr(md,md),ri(md,md), 1b(md,md),ir(md,md)) do 2 i=1,md ci(i)=0.0d0 cr(i)=0.0d0 do 2 j=1,md mr(i,j)=0.0d0 2 mi(i,j)=0.0d0 Rlimau=Rlim/0.529177d0 Rliiau=Rlii/0.529177d0 c read gaussian cub file: write(6,*)'Reading potential.cub' open(9,file='potential.cub',status='old') read(9,*) read(9,*) c r0 and rc supposed to be in bohrs, no conversion now: read(9,*)nat,r0(1),r0(2),r0(3) if(nat.ne.n)call report('Number of atoms does not match.') read(9,*)nc(1),(rc(1,i),i=1,3) read(9,*)nc(2),(rc(2,i),i=1,3) read(9,*)nc(3),(rc(3,i),i=1,3) write(6,600)nc(1)*nc(2)*nc(3),Rlii,Rlim 600 format(i10,' grid points, R_min =',f6.2,' R_max =',f6.2) do 5 i=1,n 5 read(9,*) allocate(fi(nc(3))) do 6 i1=1,nc(1) do 6 i2=1,nc(2) read(9,4000)(fi(i3),i3=1,nc(3)) 4000 format(6E13.5) do 6 i3=1,nc(3) xc=r0(1)+rc(1,1)*dble(i1-1)+rc(2,1)*dble(i2-1)+rc(3,1)*dble(i3-1) yc=r0(2)+rc(1,2)*dble(i1-1)+rc(2,2)*dble(i2-1)+rc(3,2)*dble(i3-1) zc=r0(3)+rc(1,3)*dble(i1-1)+rc(2,3)*dble(i2-1)+rc(3,3)*dble(i3-1) xx=xc-xa yy=yc-ya zz=zc-za r=dsqrt(xx**2+yy**2+zz**2) if(r.le.Rlimau.and.r.ge.Rliiau)then t=acos(zz/r) f=atan2(yy,xx) fi3=fi(i3) i=0 do 61 l=0,7 rl=r**l*fi3 do 61 m=-l,l i=i+1 call ct(l,m,t,f,rect,imct) cr(i)=cr(i)+rl*rect ci(i)=ci(i)-rl*imct j=0 do 61 lp=0,7 rlp=r**(l+lp) do 61 mp=-lp,lp j=j+1 call ct(lp,mp,t,f,rectp,imctp) mr(i,j)=mr(i,j)+rlp*( rect*rectp+imct*imctp) 61 mi(i,j)=mi(i,j)+rlp*( rect*imctp-imct*rectp) endif 6 continue close(9) write(6,*)' .... done.' c M = Mr + i Mi c M^-1 = rr + i ri c M^-1 x M = 1 c ri=Mr^-1: call INV(mr,ri,md) c ir = Mi x Mr^-1 : call MMA(ir,mi,ri,md) c iri = ir x Mi = Mi x Mr^-1 x Mi: call MMA(iri,ir,mi,md) c B = Mr + Mi x Mr^-1 x Mi: call AMA(B,mr,iri,md) c rr = B^-1 = (Mr + Mi x Mr^-1 x Mi)^-1: call INV(B,rr,md) c ri = - rr x ir = - rr x Mi x Mr^-1 call MMA(ri,rr,ir,md) call EMA(ri,md) c A = M^-1 x c = Ar + i Ai c Ar = rr x cr - ri x ci c Ai = ri x cr + rr x ci call MMV(t1,rr,cr,md) call MMV(t2,ri,ci,md) call EMV(t2,md) call AMV(ar,t1,t2,md) call MMV(t1,ri,cr,md) call MMV(t2,rr,ci,md) call AMV(ai,t1,t2,md) open(9,file='cfq2.inp') i=0 do 7 l=0,7 do 7 m=-l,l i=i+1 war=-ar(i)*219474.0d0 wai=-ai(i)*219474.0d0 7 write(9,900)l,m,war,wai 900 format(2i3,2f12.3) close(9) write(6,601) 601 format(' Grid potential transformed to cfq2.inp.') return end c ============================================================ subroutine ror(n,xa,ya,za,Rlim,Rlii,x,ty) implicit none integer*4 nat,i1,i2,i3,n,md,nc(3),l,m,i,ty(*),zsum real*8 r0(3),rc(3,3),xa,ya,za,rect,imct,rl, 1r,t,f,xc,yc,zc,xx,yy,zz,Rlim,Rlimau,wai,war, 1Rlii,Rliiau,ro3,dV,sgn,dnorm,x(*) real*8,allocatable::ro(:),ar(:),ai(:) write(6,*) write(6,*)'ES potential from density' c c setup matrix dimension md=0 do 1 l=0,7 do 1 m=-l,l 1 md=md+1 allocate(ar(md),ai(md)) do 2 i=1,md ai(i)=0.0d0 2 ar(i)=0.0d0 Rlimau=Rlim/0.529177d0 Rliiau=Rlii/0.529177d0 dnorm=0.0d0 zsum=0 c read gaussian cub file: write(6,*)'Reading density.cub' open(9,file='density.cub',status='old') read(9,*) read(9,*) c r0 and rc supposed to be in bohrs, no conversion now: read(9,*)nat,r0(1),r0(2),r0(3) if(nat.ne.n)call report('Number of atoms does not match.') read(9,*)nc(1),(rc(1,i),i=1,3) read(9,*)nc(2),(rc(2,i),i=1,3) read(9,*)nc(3),(rc(3,i),i=1,3) dV=rc(1,1)*rc(2,2)*rc(3,3) write(6,600)nc(1)*nc(2)*nc(3),Rlii,Rlim 600 format(i10,' grid points, X_min =',f6.2,' X_max =',f6.2) do 5 i=1,n 5 read(9,*) allocate(ro(nc(3))) do 6 i1=1,nc(1) do 6 i2=1,nc(2) read(9,4000)(ro(i3),i3=1,nc(3)) 4000 format(6E13.5) do 6 i3=1,nc(3) xc=r0(1)+rc(1,1)*dble(i1-1)+rc(2,1)*dble(i2-1)+rc(3,1)*dble(i3-1) yc=r0(2)+rc(1,2)*dble(i1-1)+rc(2,2)*dble(i2-1)+rc(3,2)*dble(i3-1) zc=r0(3)+rc(1,3)*dble(i1-1)+rc(2,3)*dble(i2-1)+rc(3,3)*dble(i3-1) xx=xc-xa yy=yc-ya zz=zc-za r=dsqrt(xx**2+yy**2+zz**2) if(r.le.Rlimau.and.r.gt.Rliiau)then t=acos(zz/r) f=atan2(yy,xx) ro3=-ro(i3)*dV dnorm=dnorm+ro3 i=0 do 61 l=0,7 rl=ro3/r**(l+1) do 61 m=-l,l i=i+1 call ct(l,-m,t,f,rect,imct) if(mod(m+1,2).eq.0)then sgn= rl else sgn=-rl endif ar(i)=ar(i)+rect*sgn 61 ai(i)=ai(i)+imct*sgn endif 6 continue close(9) do 8 i3=1,n xx=x(1+3*(i3-1))-xa yy=x(2+3*(i3-1))-ya zz=x(3+3*(i3-1))-za r=dsqrt(xx**2+yy**2+zz**2) if(r.le.Rlimau.and.r.ge.Rliiau)then t=acos(zz/r) f=atan2(yy,xx) c the sign is that e-e potential is positive, c electron-nuclear is negative c as put in the lanthanide.c program, c hopefully this is correct if(ty(i3).eq.1)then zsum=zsum+ty(i3) ro3=dble(ty(i3)) else c frozen core: zsum=zsum+ty(i3)-2 ro3=dble(ty(i3)-2) endif i=0 do 81 l=0,7 do 81 m=-l,l i=i+1 call ct(l,-m,t,f,rect,imct) if(mod(m+1,2).eq.0)then sgn= 1.0d0 else sgn=-1.0d0 endif ar(i)=ar(i)+ro3/r**(l+1)*rect*sgn 81 ai(i)=ai(i)+ro3/r**(l+1)*imct*sgn endif 8 continue write(6,*)' .... done.' write(6,*)'Frozen core density supposed' write(6,609)dnorm,zsum 609 format(' Dnorm = ',f9.3,' Zsum = ',i7) open(9,file='cfro.inp') i=0 do 7 l=0,7 do 7 m=-l,l i=i+1 war=ar(i)*219474.0d0 wai=ai(i)*219474.0d0 7 write(9,900)l,m,war,wai 900 format(2i3,2f12.3) close(9) write(6,601) 601 format(' Potential from density written to cfro.inp.') return end c ============================================================