program frozencore c program to calculate the self-consistent density functional c atom ground state for atom with atomic number nz c then -- states are partitioned in terms of core and valence c contributions c additional configurations are calculated within the frozen core c approximation c only states that have been calculated during initial scf calculation c can be used in new configurations implicit real*8 (a-h,o-z) parameter (mxnl=10,mxtl=mxnl*mxnl,nmx=40000,mxo=30,mxloop=200) parameter (iscr1=30,iscr2=31,iout=32,ncs=6,iden=7) dimension occ(mxo),wfn(nmx,mxo),eig(mxo),w(nmx,mxo),e(mxo) dimension r(nmx),den(nmx),rv(nmx),denout(nmx),rwfn(nmx,mxo) dimension coreden(nmx),valeden(nmx),totden(nmx) dimension tc(ncs),tv(ncs),ccoul(ncs),vcoul(ncs) dimension cext(ncs),vext(ncs),cexc(ncs),vexc(ncs) dimension cvcoul(ncs),texc(ncs),ec(ncs),ev(ncs),ecv(ncs) dimension etotal(ncs),evale(ncs) dimension nl(mxnl,mxnl),nlst(mxo),llst(mxo),no(mxnl) dimension iscore(mxo) data nl/mxtl*0/ open(unit=iscr1,status='scratch',form='unformatted') open(unit=iscr2,status='scratch',form='unformatted') open(unit=iout,file='energyresults',form='formatted') open(unit=iden,file='density',form='formatted') c call dfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,norbit, x nlst,llst,no,qf,nmx,r,wfn,eig,den,denout,rv,itype,mxloop,iscr1, x iscr2,mxo,w,e,small) c c optionally change base configuration (useful for poorly c convergent cases like Cu c 1111 write(6,*) 'enter 0 if wish to recalculate atomic configuration' read(5,*) itis if (itis.eq.0) then write(6,*) 'for each of the following states enter new occ' do io=1,norbit write(6,600) io,nlst(io),llst(io),occ(io),eig(io) 600 format(3i5,1p2e15.7) read(5,*) occ(io) enddo call rcdfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,norbit, x nlst,llst,no,qf,nmx,r,wfn,eig,den,denout,rv,itype,mxloop,iscr1, x iscr2,mxo,w,e,small) goto 1111 endif c c list all orbitals and specify core or valence states c do io=1,norbit write(6,*) 'for the following state - enter 0 for core or 1 for valence' write(6,600) io,nlst(io),llst(io),occ(io),eig(io) read(5,*) iscore(io) enddo call energyanalysis(nz,n,h,mxo,occ,norbit,qcore,qvale, x nlst,llst,iscore,nmx,r,wfn,itype,coreden,valeden,den,tc0,tv0, x ccoul0,vcoul0,cext0,vext0,cexc0,vexc0,cvcoul0,texc0,ec0,ev0, x ecv0,etotal0,evale0,mxo,w) c write(iout,*) ' Calculation for nz = ',nz write(iout,*) ' with rmax =',r(n),' , n = ',n write(iout,*) ' Convergence parameter = ', small write(iout,*) ' Reference configuration:' write(iout,*) ' Frozen core:' do io=1,norbit if (iscore(io).eq.0) then write(iout,600) io,nlst(io),llst(io),occ(io),eig(io) endif enddo write(iout,*) ' Reference valence:' do io=1,norbit if (iscore(io).eq.1) then write(iout,600) io,nlst(io),llst(io),occ(io),eig(io) endif enddo write(iout,*) ' Etotal = ',etotal0 write(iout,*) ' Evale = ',evale0 m=0 1001 continue if (m.ge.ncs) go to 4001 write(6,*) 'enter 0 if you would like to run a frozen core calculation' write(6,*) 'enter 1 if you would like to run a relaxed core calculation' read(5,*) ii if (ii.ne.0.and.ii.ne.1) then write(6,*) 'terminating frozencore calculation' go to 4001 endif m=m+1 write(6,*) 'for each old n l occ eig listed, enter new occ' do io=1,norbit if (iscore(io).eq.1) then write(6,600) io,nlst(io),llst(io),occ(io),eig(io) read(5,*) occ(io) endif enddo c write(6,*) ' new orbital listings completed' write(iout,*) ' Modified valence configuration #', m if (ii.eq.0) then write(iout,*)' -- frozen core ' call fcdfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ, x norbit,nlst,llst,no,qf,nmx,r,wfn,eig,valeden,denout,rv,itype,mxloop, x qcore,coreden,totden,iscore,iscr1,iscr2,mxo,w,e,small) endif if (ii.eq.1) then write(iout,*)' -- relaxed core ' call rcdfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,norbit, x nlst,llst,no,qf,nmx,r,rwfn,eig,den,denout,rv,itype,mxloop,iscr1, x iscr2,mxo,w,e,small) endif do io=1,norbit if (iscore(io).eq.1) then write(iout,600) io,nlst(io),llst(io),occ(io),eig(io) endif enddo if (ii.eq.0) then call energyanalysis(nz,n,h,mxo,occ,norbit,qcore,qvale, x nlst,llst,iscore,nmx,r,wfn,itype,coreden,valeden,den,tc(m),tv(m), x ccoul(m),vcoul(m),cext(m),vext(m),cexc(m),vexc(m),cvcoul(m), x texc(m),ec(m),ev(m),ecv(m),etotal(m),evale(m),mxo,w) write(iout,*) ' Etotal = ',etotal(m) write(iout,*) ' Evale = ',evale(m) endif if (ii.eq.1) then call energyanalysis(nz,n,h,mxo,occ,norbit,qcore,qvale, x nlst,llst,iscore,nmx,r,rwfn,itype,coreden,valeden,den,tc(m),tv(m), x ccoul(m),vcoul(m),cext(m),vext(m),cexc(m),vexc(m),cvcoul(m), x texc(m),ec(m),ev(m),ecv(m),etotal(m),evale(m),mxo,w) write(iout,*) ' Etotal = ',etotal(m) write(iout,*) ' Evale = ',evale(m) endif go to 1001 4001 continue c write out total density write(6,*) 'input rdenmax' read(5,*) rdenmax do i=1,n if (r(i).lt.rdenmax) write(iden,'(1p2e15.7)') r(i),den(i) enddo c write out all of the results write(6,*) ' energy results for # configurations =',m write(iout,800)0,(i,i=1,m) 800 format(/'type ',7(6x,i1,7x)) write(iout,801)tc0,(tc(i),i=1,m) write(iout,802)tv0,(tv(i),i=1,m) write(iout,803)ccoul0,(ccoul(i),i=1,m) write(iout,804)vcoul0,(vcoul(i),i=1,m) write(iout,805)cext0,(cext(i),i=1,m) write(iout,806)vext0,(vext(i),i=1,m) write(iout,807)cexc0,(cexc(i),i=1,m) write(iout,808)vexc0,(vexc(i),i=1,m) write(iout,809)cvcoul0,(cvcoul(i),i=1,m) write(iout,810)texc0,(texc(i),i=1,m) write(iout,811)ec0,(ec(i),i=1,m) write(iout,812)ev0,(ev(i),i=1,m) write(iout,813)ecv0,(ecv(i),i=1,m) write(iout,814)etotal0,(etotal(i),i=1,m) write(iout,815)evale0,(evale(i),i=1,m) 801 format('tc ',1p5e14.7) 802 format('tv ',1p5e14.7) 803 format('ccoul ',1p5e14.7) 804 format('vcoul ',1p5e14.7) 805 format('cext ',1p5e14.7) 806 format('vext ',1p5e14.7) 807 format('cexc ',1p5e14.7) 808 format('vexc ',1p5e14.7) 809 format('cvcoul ',1p5e14.7) 810 format('texc ',1p5e14.7) 811 format('ec ',1p5e14.7) 812 format('ev ',1p5e14.7) 813 format('ecv ',1p5e14.7) 814 format('etotal ',1p5e14.7) 815 format('evale ',1p5e14.7) write(iout,816) 816 format(//'differences -- ') do i=1,m tc(i)=tc(i)-tc0 tv(i)=tv(i)-tv0 cext(i)=cext(i)-cext0 vext(i)=vext(i)-vext0 ccoul(i)=ccoul(i)-ccoul0 vcoul(i)=vcoul(i)-vcoul0 cexc(i)=cexc(i)-cexc0 vexc(i)=vexc(i)-vexc0 cvcoul(i)=cvcoul(i)-cvcoul0 texc(i)=texc(i)-texc0 ec(i)=ec(i)-ec0 ev(i)=ev(i)-ev0 ecv(i)=ecv(i)-ecv0 etotal(i)=etotal(i)-etotal0 evale(i)=evale(i)-evale0 enddo write(iout,801)tc0,(tc(i),i=1,m) write(iout,802)tv0,(tv(i),i=1,m) write(iout,803)ccoul0,(ccoul(i),i=1,m) write(iout,804)vcoul0,(vcoul(i),i=1,m) write(iout,805)cext0,(cext(i),i=1,m) write(iout,806)vext0,(vext(i),i=1,m) write(iout,807)cexc0,(cexc(i),i=1,m) write(iout,808)vexc0,(vexc(i),i=1,m) write(iout,809)cvcoul0,(cvcoul(i),i=1,m) write(iout,810)texc0,(texc(i),i=1,m) write(iout,811)ec0,(ec(i),i=1,m) write(iout,812)ev0,(ev(i),i=1,m) write(iout,813)ecv0,(ecv(i),i=1,m) write(iout,814)etotal0,(etotal(i),i=1,m) write(iout,815)evale0,(evale(i),i=1,m) stop end