1 files changed, 587 insertions, 0 deletions

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