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=50,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) character*2 atomname character*1 input data nl/mxtl*0/ open(unit=iscr1,status='scratch',form='unformatted') open(unit=iscr2,status='scratch',form='unformatted') now=-1 write(6,*) 'Enter atomic symbol' read(5,*) atomname open(unit=iout,file=TRIM(atomname)//'.energyresults',form='formatted') c 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,ecoul,etxc,etot) 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,ecoul,etxc,etot) goto 1111 endif c c list all orbitals and specify core or valence states c do io=1,norbit readloop: do write(6,*)'for the following state, enter "c" for core or "v" for valence' write(6,600) io,nlst(io),llst(io),occ(io),eig(io) read(5,*) input if (input=="c".or.input=="C") then iscore(io)=0 exit readloop elseif (input=="v".or.input=="V") then iscore(io)=1 exit readloop else write(6,*) 'Typo -- please try again' endif enddo readloop c(old) 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 If (now.eq.0) write(6,*)'Just completed frozen core calculation' If (now.eq.1) write(6,*)'Just completed relaxed core calculation' write(6,*)'enter 0 if you would like to run new frozen core calculation' write(6,*)'enter 1 if you would like to run new relaxed core calculation' read(5,*) now if (now.ne.0.and.now.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 (now.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,ecoul,etxc,etot) endif if (now.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,ecoul,etxc,etot) 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 (now.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 (now.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 c write(6,*) 'input rdenmax' c read(5,*) rdenmax c do i=1,n c if (r(i).lt.rdenmax) write(iden,'(1p2e15.7)') r(i),den(i) c enddo c write out all of the results write(6,*) write(6,*) ' Summary of all results ' write(6,*) ' energy results for # configurations =',m write(iout,*) write(iout,*) ' Summary of all results ' write(iout,*) ' energy results for # configurations =',m last=m next=1 do left=min(next+3,last) write(iout,800)0,(i,i=next,left) 800 format(/' Label ',7(6x,i1,7x)) c write(iout,801)tc0,(tc(i),i=1,m) c write(iout,802)tv0,(tv(i),i=1,m) c write(iout,803)ccoul0,(ccoul(i),i=1,m) c write(iout,804)vcoul0,(vcoul(i),i=1,m) c write(iout,805)cext0,(cext(i),i=1,m) c write(iout,806)vext0,(vext(i),i=1,m) c write(iout,807)cexc0,(cexc(i),i=1,m) c write(iout,808)vexc0,(vexc(i),i=1,m) c write(iout,809)cvcoul0,(cvcoul(i),i=1,m) c write(iout,810)texc0,(texc(i),i=1,m) c write(iout,811)ec0,(ec(i),i=1,m) c write(iout,812)ev0,(ev(i),i=1,m) c write(iout,813)ecv0,(ecv(i),i=1,m) write(iout,814)etotal0,(etotal(i),i=next,left) write(iout,815)evale0,(evale(i),i=next,left) 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(' Tot',1p5e14.7) 815 format(' Val',1p5e14.7) 914 format('Delta Tot',1p5e14.7) 915 format('Delta Val',1p5e14.7) c write(iout,816) c 816 format(//'differences -- ') c do i=1,m c tc(i)=tc(i)-tc0 c tv(i)=tv(i)-tv0 c cext(i)=cext(i)-cext0 c vext(i)=vext(i)-vext0 c ccoul(i)=ccoul(i)-ccoul0 c vcoul(i)=vcoul(i)-vcoul0 c cexc(i)=cexc(i)-cexc0 c vexc(i)=vexc(i)-vexc0 c cvcoul(i)=cvcoul(i)-cvcoul0 c texc(i)=texc(i)-texc0 c ec(i)=ec(i)-ec0 c ev(i)=ev(i)-ev0 c ecv(i)=ecv(i)-ecv0 c etotal(i)=etotal(i)-etotal0 c evale(i)=evale(i)-evale0 c enddo c write(iout,801)tc0,(tc(i),i=1,m) c write(iout,802)tv0,(tv(i),i=1,m) c write(iout,803)ccoul0,(ccoul(i),i=1,m) c write(iout,804)vcoul0,(vcoul(i),i=1,m) c write(iout,805)cext0,(cext(i),i=1,m) c write(iout,806)vext0,(vext(i),i=1,m) c write(iout,807)cexc0,(cexc(i),i=1,m) c write(iout,808)vexc0,(vexc(i),i=1,m) c write(iout,809)cvcoul0,(cvcoul(i),i=1,m) c write(iout,810)texc0,(texc(i),i=1,m) c write(iout,811)ec0,(ec(i),i=1,m) c write(iout,812)ev0,(ev(i),i=1,m) c write(iout,813)ecv0,(ecv(i),i=1,m) write(iout,914)etotal0,(etotal(i)-etotal0,i=next,left) write(iout,915)evale0,(evale(i)-evale0,i=next,left) next=left+1 if(next.gt.last) exit enddo stop end