From 5df79c53745fde5d6c3340a2979b1429cd5892c1 Mon Sep 17 00:00:00 2001 From: Henrik Rydberg Date: Sat, 8 Oct 2011 20:30:28 +0200 Subject: Initial import of htcd system 1.0 Signed-off-by: Henrik Rydberg --- src/labat/labat.f | 587 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 587 insertions(+) create mode 100644 src/labat/labat.f (limited to 'src/labat/labat.f') diff --git a/src/labat/labat.f b/src/labat/labat.f new file mode 100644 index 0000000..da1f7fe --- /dev/null +++ b/src/labat/labat.f @@ -0,0 +1,587 @@ + program atom + implicit double precision (a-h,o-z) + dimension dval(2,561),nval(2,7,10),ddval(2,561),ndval(2,7,10) + dimension dcore(2,561),vhv(561),fpu(561) +c +c free atom program by m.p. summer 1985 (valence energy and core +c density are calculated) +c +c GGA implemented by T. Holmquist and U. Yxklinten. +c +c files used: +c 6: terminal output +c 22: density and potential output +c 33: density and potential input +c 44: potential and density double prec. output +c 46: potential and density double prec. input +c 55: control input +c 66: printer output +c 77: density (total and core) output (J.H.) +c + dimension r(561),veff(2,561),vc(561),vcold(561),vv(561) + dimension gr(561),lm(2),rnocc(2,7,10),ux1(561),ux2(561) + dimension u(561),roo(561),spl(561),nb(7) + dimension db(2,561),v(561),vold(2,561) + dimension snlo(561),wj(301),ekin(2) + dimension ebound(2,7,10),nbound(2,7),de(2,7,10) + dimension apu(561),bpu(561),cpu(561),dpu(561),nlp(2,7,10) + dimension epu(561) +c=================================================================== +c gga + dimension excgga(561),uxcup(561),uxcdn(561) + logical doned,donela,donegr,donez,donec +c gga +c=================================================================== + common/sc/gr,r,snlo,nbl + common/mess/wj,dx,nblock,jblock + common/pot/v + frs(x)=(3.d0/(4.d0*pi*x))**(1.d0/3.d0) + pi=4.d0*datan(1.d0) + pi2=pi/2.d0 + a=(4.d0/9.d0/pi)**(1.d0/3.d0) +c + open(22,file='adensout.dat',status='unknown') + open(33,file='adensin.dat',status='unknown') + open(44,file='apotout.dat',status='unknown') + open(46,file='apotin.dat',status='unknown') + open(55,file='atomctrl.dat',status='old') + open(66,file='aprtout.dat',status='unknown') + open(77,file='atomdens.dat',status='unknown') + write(6,*) '=== Its ===== Energy (Ha) ======' + call flush(6) + read(55,*)z,zion +c z : atomic number +c zion : ionicity +c + write(66,2) z,zion +2 format(/' atom number:',f5.0,' ionicity:',f5.0) + sf=4.d0*pi + read(55,*)jblock,nbl,c +c parameters of the hermann-skillmann mesh + read(55,*)fback,thresh,qsc +C fback : feedback +c thresh : error allowed for the bound state eigenenergy +c e.g. 0.00001 +c qsc : parameter for the screened green's function (0.005) + read(55,*)itmax,inopt,ipr,iout +c itmax : max. no. of iterations +c iopt : =0 starting potential generated; =1 starting potential +c read in; =2 starting potential read in from unit 46. +c (double precision format) +c ipr, iout : print out parameters +c +c set up the herman-skillman mesh +c + nblock=jblock + mesh=nblock*nbl+1 + n=mesh + mest=mesh + dx=c*0.0025d0 + en0=-0.01d0 + i=1 + r(i)=0.d0 + deltax=dx + do 251 j=1,nblock + do 241 jk=1,nbl + i=i+1 + 241 r(i)=r(i-1)+deltax + deltax=2.d0*deltax + 251 continue +c + read(55,*)ne,iys,ixc +c ne: no. bound states +c iys: =1 for spin compensated; =2 spin polarized +c ixc: = 1 b-h xc; =0 c-a xc + if(ixc.eq.1)write(66,2149) + if(ixc.eq.0)write(66,2148) + if(ixc.eq.2)write(66,2147) +2149 format(/' von barth - hedin xc') +2148 format(/' ceperley - alder xc') +2147 format(/' perdew - wang xc') + dxx=r(n)-r(n-1) + write(66,5)jblock,nbl,c + write(66,5849)dx,dxx,r(n) + 5 format(/' hermann-skillmann mesh',/,' blocks, nbl, c:',2(i5,','), + j f10.7) +5849 format(' first interval, last interval, final r:',2(f10.6,','), + j f10.6) + write(66,2312) thresh,qsc,fback +2312 format(/' energy eigenvalue threshold, screening parameter, ', + j 'feedback:',2(f12.6,','),f6.3) + write(66,9) + do 5273 i=1,7 +5273 nb(i)=0 + lmax=0 +9 format(//' initial bound state configuration',/, + j ' nlm energy occup. val.el. d-band el.') + zv=0.d0 + do 5595 ispin=1,iys + do 5595 ii=1,ne + read(55,*)nnlz,eb0,rnoc,nv,ndv +c +c nnlz : e.g. 100 +c eb0 : energy eigenvalue guess +c rnoc : occupation number +c nv : =1 when orbital is a valence orbital; =0 for core orbitals +c ndv : =1 when orbital is a d-valence orbital; =0 for other orbitals +c + nnn=nnlz/100 + lll=(nnlz-nnn*100)/10 + nnn=nnn-lll + if(nnn.gt.nb(lll+1))nb(lll+1)=nnn + if(lll.gt.lmax)lmax=lll + rnocc(ispin,lll+1,nnn)=rnoc + nval(ispin,lll+1,nnn)=nv + ndval(ispin,lll+1,nnn)=ndv + if(iys.eq.1)nval(2,lll+1,nnn)=nv + if(iys.eq.1)ndval(2,lll+1,nnn)=ndv + if(iys.eq.1)rnocc(2,lll+1,nnn)=rnoc + if(iys.eq.1)rnoc=rnoc*2 + write(66,19)nnlz,eb0,rnoc,nv,ndv + zv=zv+rnoc*(nv+ndv) +5595 ebound(ispin,lll+1,nnn)=eb0 + write(66,5120)zv +5120 format(/' number of valence electrons: ',f5.1) + 19 format(i6,e12.4,f7.2,i7,i10) + nb1=nbl+1 + iter=0 + itr=itmax-iout +c + 16 format(7(f12.5,1x)) + 17 format(/' charge and spin density profiles '/) + 18 format(1h0/' potentials vs. distance after ',i3,' iterations'/) +c +c tabulate screened green's function + do 10 i=1,n + 10 gr(i)=exp(-qsc*r(i)) +c +c integration weigths + do 121 i=1,nb1 + 121 wj(i)=(3.d0+(-1.d0)**i)/3.d0 + wj(1)=1.d0/3.d0 + wj(nb1)=1.d0/3.d0 +c +c initial values of the potentials + if(inopt.eq.1) goto 302 + if(inopt.eq.2) goto 306 +c + itt=0 +c thomas-fermi potential + do 30 i=2,n + x=r(i)/0.88534135d0*z**(1.d0/3.d0) + xx=sqrt(x) + vc(i)=-z/r(i)/(1.+0.02747d0*xx+1.243d0*x-0.1486d0*x*xx + j+0.2302d0*x*x+0.007298d0*x*x*xx+0.006944d0*x*x*x) + do 30 ispin=1,2 + 30 veff(ispin,i)=vc(i) + goto 311 +c read in an old potential for the input +302 read(33,2645)itt +2645 format(i3) +1645 format(i3,' z,ion:',2f4.0,' jblock,nbl,c:',2i5,f10.7) + do 303 i=1,n + 303 read(33,907) vc(i),veff(1,i),veff(2,i),du1,du2 + goto 311 + 306 read(46,2645) itt + do 304 i = 1, n + read(46,*) du1, vc(i), vve, du2 + veff(1,i) = vve + veff(2,i) = vve + 304 continue + 311 do 312 i=1,n + vcold(i)=vc(i) + do 312 ispin=1,2 + 312 vold(ispin,i)=veff(ispin,i) + 5555 format(e15.8) + 7771 format(/' initial potentials'/) + 7773 format(5f13.5) + if(iout.lt.-1)go to 1437 + write(66,7771) + do 7772 i=2,n,ipr + 7772 write(66,7773) r(i),vc(i),veff(1,i),veff(2,i) +1437 eold=0.d0 + etot=-1.d0 + itec=0 +c +c iteration +c + 100 iter=iter+1 +c write(6,*) 'NEW ITERATION' +c call flush(6) +c convergency check + if(abs(etot-eold).lt.5.d-5)itec=itec+1 + if(abs(etot-eold).lt.5.d-5.and.itec.eq.1)itmax=iter + eold=etot + if(iter-itmax)999,999,400 +999 itt=itt+1 + write(66,6147)itt + do 117 i=1,n + dval(2,i)=0.d0 + ddval(2,i)=0.d0 + dcore(1,i)=0.d0 + dcore(2,i)=0.d0 + dval(1,i)=0.d0 + 117 ddval(1,i)=0.d0 + 6147 format(///' ********************* iteration:', + 1 i4,' ***********************'/) + do 917 ispin=1,iys + do 116 i=1,n + vv(i)=veff(ispin,i) + v(i)=2*vv(i) + 116 db(ispin,i)=0.d0 + lm(ispin)=-1 + do 807 i=0,4 + ik=i + if(i)817,817,818 +817 do 819 jj=1,4 +819 u(jj)=r(jj)-z*r(jj)**2 + go to 822 +818 do 821 jj=1,4 +821 u(jj)=r(jj)**(i+1) +822 call schrhs(vv,0.d0,ik,u) + merkki=1 +c determine the number of bound states + ncross=0 + do 826 ii=2,n + if(merkki)823,823,824 +823 if(u(ii))826,826,825 +824 if(u(ii))825,826,826 +825 merkki=-merkki + ncross=ncross+1 +826 continue + dlo=(u(n)-u(n-1))*u(n) + if(dlo.ge.0)nbound(ispin,i+1)=ncross + if(dlo.lt.0)nbound(ispin,i+1)=ncross+1 + if(nbound(ispin,i+1).eq.0)go to 827 + lm(ispin)=i + write(66,678)i,nbound(ispin,i+1) + nbound(ispin,i+1)=min0(nbound(ispin,i+1),nb(i+1)) +678 format(' l = ',i4,',',i4,' bound states ') +807 continue +827 if(lm(ispin).lt.0.and.ispin.eq.2)go to 674 + if(lm(ispin).lt.0.and.ispin.eq.1)go to 917 + lm(ispin)=min0(lm(ispin),lmax) + lmm=lm(ispin) + do 918 i=0,lmm + nii=min0(nbound(ispin,i+1),8) + ik=i + do 918 ii=1,nii + if(rnocc(ispin,i+1,ii).lt.0.1d0)go to 918 + nn=i+ii + en=2.d0*ebound(ispin,i+1,ii) + if(en.ge.0.d0)en=en0 + en1=en + dde=de(ispin,i+1,ii) + if(iter.le.2)dde=0.d0 +c determine the bound state energy and eigenfunction + call scheq(z,en,ik,nn,mest,mesh,c,thresh,iflag,npr,dde) +c scheq operates in rydberg units + if(iter.ge.2)de(ispin,i+1,ii)=abs(en-en1) + if(iflag.eq.1) goto 400 + 163 ebound(ispin,i+1,ii)=en/2.d0 + nlp(ispin,i+1,ii)=npr + do 9188 ji=1,n + dval(ispin,ji)=dval(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/sf + j*nval(ispin,i+1,ii) + ddval(ispin,ji)=ddval(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/ + j sf*ndval(ispin,i+1,ii) + dcore(ispin,ji)=dcore(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/ + j sf*(1-ndval(ispin,i+1,ii))*(1-nval(ispin,i+1,ii)) + 9188 db(ispin,ji)=db(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/sf + 918 continue + do 919 ji=2,n + dcore(ispin,ji)=dcore(ispin,ji)/r(ji)**2 + 919 db(ispin,ji)=db(ispin,ji)/r(ji)**2 + db(ispin,1)=db(ispin,2) + dcore(ispin,1)=dcore(ispin,2) + do 5739 ji=2,n + ddval(ispin,ji)=ddval(ispin,ji)/r(ji)**2 +5739 dval(ispin,ji)=dval(ispin,ji)/r(ji)**2 + dval(ispin,1)=dval(ispin,2) + ddval(ispin,1)=ddval(ispin,2) + 917 continue + if(iys.eq.2)go to 923 + do 921 i=1,n + dval(2,i)=dval(1,i) + ddval(2,i)=ddval(1,i) + dcore(2,i)=dcore(1,i) +921 db(2,i)=db(1,i) + lmm=lm(1) + lm(2)=lm(1) + do 922 i=0,lmm + nii=nbound(1,i+1) + nbound(2,i+1)=nii + do 922 ii=1,nii +922 ebound(2,i+1,ii)=ebound(1,i+1,ii) +923 continue + 8 format(' ',i6,' energy: ',f12.5,' iteration loops: ',i4) +c +c total bound state energy +c +c write(6,*) 'Total energy starts' +c call flush(6) + 661 eb=0.d0 + ebv=0.d0 + write(66,7496) +7496 format(//' bound states'/) + do 673 ispin=1,2 + lmm=lm(ispin) + do 673 i=0,lmm + nii=nbound(ispin,i+1) + do 673 j=1,nii + nn=i+j + nnlz=100*nn+10*i + eb=eb+ebound(ispin,i+1,j)*rnocc(ispin,i+1,j) + ebv=ebv+ebound(ispin,i+1,j)*rnocc(ispin,i+1,j)* + j (nval(ispin,i+1,j)+ndval(ispin,i+1,j)) + if(rnocc(ispin,i+1,j).lt.0.1d0) go to 673 + write(66,8) nnlz,ebound(ispin,i+1,j),nlp(ispin,i+1,j) + 673 continue + go to 679 +674 write(66,676) +676 format('0no bound states') +679 continue +c +c total charge and spin densities +c computing energy integrals +c +506 continue +c========================================= +c gga +c set the "done" controle variables to false in the beginning of +c each iteration + + if (ixc.eq.2) then + doned = .false. + donela = .false. + donegr = .false. + donez = .false. + donec = .false. + endif +c gga +c======================================== + do 660 i=2,n + roo(i)=db(1,i)+db(2,i) + rr=roo(i) +c When db = (0,0) then the spin-polarization, spl = 0, and not 0/0. + if (rr.gt.0.0000000001) then + spl(i)=(db(1,i)-db(2,i))/rr + else + spl(i)=0.0 + endif + x1=r(i)**2 + apu(i)=x1*veff(1,i)*db(1,i)*sf + bpu(i)=x1*veff(2,i)*db(2,i)*sf + if (ixc.eq.2) then +c================================================================ +c gga +c calculate the gga exchange-correlation energy + + call ggaexc(r,db,i,1,n,doned,excgga(i)) + cpu(i) = x1*rr*excgga(i)*sf +c gga +c================================================================ + else + cpu(i)=x1*rr*exc(rr,spl(i),ixc)*sf + endif + dpu(i)=x1*spl(i)*rr*sf + 660 continue + 26 format(' induced moment: ',f10.5) + call simpsh(apu,v1) + call simpsh(bpu,v2) + call simpsh(cpu,eexc) + call simpsh(dpu,smom) +c +c compute coulomb energy +c + cz=z + do 683 i=1,n + bpu(i)=roo(i)*r(i)**2*sf + dpu(i)=(dval(1,i)+dval(2,i))*r(i)**2*sf + apu(i)=(ddval(1,i)+ddval(2,i))*r(i)**2*sf + cpu(i)=roo(i)*(r(i)**2*vc(i)-z*r(i))*sf/2.d0 +683 continue + call simpsh(dpu,sum1) + call simpsh(apu,sum3) + call simpsh(bpu,sum2) + call simpsh(cpu,ec) + 24 format(' kin: ',2(e14.6,1x),' coul: ',e14.6,' exc: ',e14.6) + 25 format(/' total energy: ',e15.7) + write(66,192) sum2,sum1,sum3 + write(66,26) smom + ekin(1)=-v1 + ekin(2)=-v2 + etot=ekin(1)+ekin(2)+eb+ec+eexc + write(66,7453) +7453 format(/' energy terms') + write(66,622)eb + 622 format(' energy eigenvalue sum: ',e14.6) + write(66,24) ekin(1),ekin(2),ec,eexc + write(66,25) etot +c write(6,2573) itt,sum2,etot +c call flush(6) +c2573 format(' iter:',i4,' total ch:',f6.2,' total en:',f15.7) + write(6,2573) itt,etot + call flush(6) +2573 format(' ',i4,' ',f15.5) + + if(iter.ne.1 .and. iter.lt.itr) goto 167 + if(iout.lt.0)go to 167 + write(66,17) + 192 format(/' integr charges: total ',f10.7,' valence: ',2f10.7) + do 166 i=1,n,ipr + r22=r(i)*r(i) + db1=db(1,i) + db2=db(2,i) + dr=db1+db2 + dcc=dcore(1,i)+dcore(2,i) + dvv=dval(1,i)+dval(2,i)+ddval(1,i)+ddval(2,i) + ddb=db1-db2 +c166 write(66,1624) r(i),db1,db2,dval(i),ddval(i),dr,ddb + 166 write(66,1624) r(i),db1,db2,dr,dcc,dvv +c +c +c compute coulomb potential +c + 167 do 55 i=2,n + 55 vc(i)=r(i)*vc(i)+zion +c write(6,*) 'Coulomb potential' +c call flush(6) + zi=z-zion + vc(1)=-zi + call scrhs(roo,vc,zi,qsc,n) + do 50 i=2,n + 50 vc(i)=(vc(i)-zion)/r(i) +c +c set up the total potential +c + do 60 i=2,n + dr=roo(i) + vc(i)=(1.d0-fback)*vcold(i)+fback*vc(i) + vcold(i)=vc(i) + if (ixc.eq.2) then +c================================================================ +c gga +c calculate the gga exchange-correlation potential + + call ggauxc(r,db,spl,i,1,n,donela,donegr,donez, + j uxcup(i),uxcdn(i)) + ux1(i) = uxcup(i) + ux2(i) = uxcdn(i) +c gga +c================================================================ + else + ux1(i)=uxc(dr,spl(i),1,ixc) + ux2(i)=uxc(dr,spl(i),2,ixc) + endif + 6848 veff(1,i)=vc(i)+ux1(i) + veff(2,i)=vc(i)+ux2(i) + do 60 ispin=1,2 + veff(ispin,i)=(1.d0-fback)*vold(ispin,i)+fback*veff(ispin,i) + 60 vold(ispin,i)=veff(ispin,i) +c calculate the valence energy + do 5100 i=2,mesh + do 5110 j=1,mesh + apu(j)=0.d0 + if(i.le.j)apu(j)=(dval(1,j)+dval(2,j)+ddval(1,j)+ddval(2,j))* + j (r(j)**2/r(i)-r(j)) +5110 continue + call simpsh(apu,vp) +5100 vhv(i)=-vp*4.d0*pi+zv/r(i) + do 5130 i=1,mesh + x1=r(i)**2 + rr=db(1,i)+db(2,i) + spt=0.d0 + if(rr.gt.1.d-10)spt=(db(1,i)-db(2,i))/rr + rrv=dval(1,i)+dval(2,i)+ddval(1,i)+ddval(2,i) + spv=0.d0 + if(rrv.gt.1.d-10)spv=(dval(1,i)-dval(2,i)+ddval(1,i)-ddval(2,i)) + j /rrv + rrc=dcore(1,i)+dcore(2,i) + spc=0.d0 + if(rrc.gt.1.d-10)spc=(dcore(1,i)-dcore(2,i))/rrc + apu(i)=-0.5d0*x1*sf*rrv*vhv(i) + if (ixc.eq.2) then +c============================================================= +c gga +c calculate the gga exchange-correlation energy for core electrons + + if (rrc.lt.1.d-100) then + excgac = 0.d0 + else + call ggaexc(r,dcore,i,1,n,donec,excgac) + endif + bpu(i)=x1*rr*excgga(i)*sf + cpu(i)=-x1*rrc*excgac*sf + fpu(i)=(dcore(1,i)+dcore(2,i))*excgga(i)*x1*sf + dpu(i)=(-(dval(1,i)+ddval(1,i))*uxcup(i) + j -(dval(2,i)+ddval(2,i))*uxcdn(i))*sf*x1 +c gga +c================================================================== + else + bpu(i)=x1*rr*exc(rr,spt,ixc)*sf + cpu(i)=-x1*rrc*exc(rrc,spc,ixc)*sf + fpu(i)=(dcore(1,i)+dcore(2,i))*exc(rr,spt,ixc)*x1*sf + dpu(i)=(-(dval(1,i)+ddval(1,i))*uxc(rr,spt,1,ixc) + j -(dval(2,i)+ddval(2,i))*uxc(rr,spt,2,ixc))*sf*x1 + endif + epu(i)=apu(i)+bpu(i)+cpu(i)+dpu(i)+fpu(i) +5130 continue +c call simpsh(apu,ec) +c call simpsh(bpu,eext) +c call simpsh(fpu,eexx) +c call simpsh(cpu,eexc) +c call simpsh(dpu,eexv) +c etv=ebv+ec+eext+eexc+eexv + call simpsh(epu,etv) + etv = etv + ebv + write(66,5140)etv +5140 format(' valence energy:',e15.7) + if(iter.ne.1 .and. iter.lt.itr) go to 100 + if(iout.lt.0)go to 100 + write(66,18) itt + do 75 i=2,n,ipr + dr=roo(i) + 75 write(66,1624) r(i),vc(i),ux1(i),ux2(i),veff(1,i),veff(2,i) +1624 format(1x,f12.5,6e14.5) + goto 100 + 400 continue +c +c +c + do 1140 i=2,mesh + do 1130 j=1,mesh + apu(j)=0.d0 + if(i.le.j)apu(j)=(db(1,j)+db(2,j))*(r(j)**2/r(i)-r(j)) +1130 continue + call simpsh(apu,vp) +1140 vc(i)=-vp*sf + write(22,1645)itt,z,zion,jblock,nbl,c + write(44,1645)itt,z,zion,jblock,nbl,c + do 900 i=1,n + dd=db(1,i)+db(2,i) + dc=dcore(1,i)+dcore(2,i) + dddc = dd + dc + dv1=dval(1,i)+dval(2,i) + ddv1=ddval(1,i)+ddval(2,i) + vve = 0.5d00*(veff(1,i) + veff(2,i)) + write(44,*) r(i), vc(i), vve, dddc +c The line below writes total and core densities to atomdens.dat (J.H.) + write(77,906) r(i), 4*pi*dd*r(i)*r(i), 4*pi*dc*r(i)*r(i) + write(22,907)vc(i),veff(1,i),veff(2,i),dddc,r(i) + 900 continue + 907 format(6e15.7) +c900 write(22,906)vc(i),veff(1,i),veff(2,i),dd,dc +c900 write(22,906)vc(i),veff(1,i),veff(2,i),dd,dv1,ddv1 + 906 format(6e20.12) + write(6,*) '================================' + call flush(6) +c ieeer=ieee_flags('clear','exception','all',ieeeout) + close(6) + close(22) + close(33) + close(55) + close(66) + close(77) + stop + end -- cgit v1.2.3