!--------------------------------------------------------------------- subroutine exch_corr(nspin,rhog,rhor,exc) !--------------------------------------------------------------------- ! ! calculate exch-corr potential and energy ! use dft_mod use gvec use parm implicit none ! input integer nspin ! rhog contains the charge density in G space ! rhor contains the charge density in R space complex(kind=8) rhog(ng,nspin) ! output ! rhor contains the exchange-correlation potential real(kind=8) rhor(nnr,nspin), exc ! local real(kind=8), allocatable:: gradr(:,:,:) ! ! fillgrad calculates gradr=grad(rho) in real space using fft's ! if (dft.ne.lda) then allocate(gradr(nnr,3,nspin)) call fillgrad(nspin,rhog,gradr) end if ! exc=0.0 if (dft.eq.lda) then call expxc(nnr,nspin,rhor,exc) else if (dft.eq.pw91) then call ggapw(nspin,rhog,gradr,rhor,exc) !!! call ggapwold(nspin,rhog,gradr,rhor,exc) else if (dft.eq.blyp) then call ggablyp4(nspin,rhog,gradr,rhor,exc) else if (dft.eq.pbe) then call ggapbe(nspin,rhog,gradr,rhor,exc) else call error('exc-cor','no such exch-corr',dft) end if ! ! second part of the xc-potential for gradient-corrected dft: ! \sum_alpha (D / D r_alpha) ( D(rho*exc)/D(grad_alpha rho) ) ! if (dft.ne.lda) then call grad2(nspin,gradr,rhor) deallocate(gradr) end if ! exc=exc*omega/dfloat(nr1*nr2*nr3) ! return end ! !----------------------------------------------------------------------- subroutine force_cc(ei1,ei2,ei3,tsum) !----------------------------------------------------------------------- ! ! non linear core correction to fion (to be finished) ! use elct, only:ng0 use gvec use ions use parm use ncprm ! implicit none ! complex(kind=8) ei1(-nr1:nr1,nas,nsp), & & ei2(-nr2:nr2,nas,nsp), & & ei3(-nr3:nr3,nas,nsp) complex(kind=8) tsum(3,nax,nsx) ! integer is, ia, ix, ig, cc real(kind=8) wz complex(kind=8) eigrx, CSUM ! cc=0 do is=1,nsp cc=cc+ifpcor(is) end do if (cc.eq.0) return ! do is=1,nsp if(ifpcor(is).eq.1) then ! do ia=1,na(is) ! do ix=1,3 ! rhog(1)=0.0 ! do ig=ng0,ng ! eigrx = ei1(in1p(ig),ia,is) & ! & *ei2(in2p(ig),ia,is) & ! & *ei3(in3p(ig),ia,is) ! rhog(ig)= eigrx*v(np(ig))*rhoc(ig,is) & ! & *cmplx(0.d0,gx(ig,ix)) ! end do ! tsum(ix,ia,is) = tsum(ix,ia,is)+wz*CSUM(ng,rhog,1) ! end do ! end do call error('force_cc','not yet ready',is) endif end do return end ! !----------------------------------------------------------------------- subroutine force_ion(tau0,esr,fion) !----------------------------------------------------------------------- ! use cnst use ions use parm ! implicit none ! ! input ! real(kind=8) tau0(3,nax,nsx) ! ! output ! real(kind=8) fion(3,nax,nsx), esr, ERFC ! ! local variables ! integer i,j,k,l,m, lax, inf real(kind=8) rlm(3), rckj, rlmn, arg, addesr, addpre, repand, fxx ! ! esr=0.d0 ! desr=0.d0 ! do k=1,nsp do j=k,nsp rckj=sqrt(rcmax(k)**2+rcmax(j)**2)!ATTENZIONE: originale...... ! rckj=0.5*sqrt(rcmax(k)**2+rcmax(j)**2)!ATTENZIONE per i calcoli MgO e LiCl, solo 8 atomi per cella.... lax=na(k) if(k.eq.j) lax=lax-1 ! do l=1,lax inf=1 if(k.eq.j) inf=l+1 ! do m=inf,na(j) rlm(1)=tau0(1,l,k)-tau0(1,m,j) rlm(2)=tau0(2,l,k)-tau0(2,m,j) rlm(3)=tau0(3,l,k)-tau0(3,m,j) call pbc(rlm,a1,a2,a3,ainv,rlm) ! rlmn=sqrt(rlm(1)**2+rlm(2)**2+rlm(3)**2) ! arg=rlmn/rckj addesr=zv(k)*zv(j)*ERFC(arg)/rlmn esr=esr+addesr addpre=2.d0*zv(k)*zv(j)*exp(-arg*arg)/rckj/sqrt(pi) ! desr=desr-(addesr+addpre)/3.d0/omega repand=(addesr+addpre)/rlmn/rlmn ! do i=1,3 fxx=repand*rlm(i) fion(i,l,k)=fion(i,l,k)+fxx fion(i,m,j)=fion(i,m,j)-fxx end do end do end do end do end do ! return end ! !----------------------------------------------------------------------- subroutine force_ps(rhotemp,rhog,vtemp,ei1,ei2,ei3,tsum) !----------------------------------------------------------------------- ! use cnst use elct use gvec use gvecs use ions use parm use pseu ! implicit none complex(kind=8) rhotemp(ng), rhog(ng,nspin), vtemp(ng), & & ei1(-nr1:nr1,nas,nsp), & & ei2(-nr2:nr2,nas,nsp), & & ei3(-nr3:nr3,nas,nsp) complex(kind=8) tsum(3,nax,nsx) ! integer ig, is, isa, ism, ia, ix, iss, isup, isdw real(kind=8) wz complex(kind=8) eigrx, vcgs, cnvg, cvn, CSUM ! ! wz=2.0 do is=1,nsp isa=0 do ism=1,is-1 isa=isa+na(ism) end do do ia=1,na(is) isa=isa+1 do ix=1,3 if(nspin.eq.1)then iss=1 if (ng0.eq.2) vtemp(1)=0.0 do ig=ng0,ngs vcgs=conjg(rhotemp(ig))*fpi/(tpiba2*g(ig)) cnvg=rhops(ig,is)*vcgs cvn=vps(ig,is)*conjg(rhog(ig,iss)) eigrx=ei1(in1p(ig),ia,is)*ei2(in2p(ig),ia,is) & & *ei3(in3p(ig),ia,is) vtemp(ig)=eigrx*(cnvg+cvn)*cmplx(0.0,gx(ig,ix)) end do else isup=1 isdw=2 if (ng0.eq.2) vtemp(1)=0.0 do ig=ng0,ngs vcgs=conjg(rhotemp(ig))*fpi/(tpiba2*g(ig)) cnvg=rhops(ig,is)*vcgs cvn=vps(ig,is)*conjg(rhog(ig,isup) & & +rhog(ig,isdw)) eigrx=ei1(in1p(ig),ia,is)*ei2(in2p(ig),ia,is) & & *ei3(in3p(ig),ia,is) vtemp(ig)=eigrx*(cnvg+cvn)*cmplx(0.0,gx(ig,ix)) end do endif tsum(ix,ia,is) = tsum(ix,ia,is) + wz*CSUM(ngs,vtemp,1) end do end do end do ! return end ! !----------------------------------------------------------------------- subroutine formf !----------------------------------------------------------------------- !computes (a) the self-energy eself of the ionic pseudocharges; ! (b) the form factors of: (i) pseudopotential (vps), ! (ii) ionic pseudocharge (rhops) ! all quantities are returned in module pseu ! use gvec use gvecs use cnst use parm use ions use pseu use elct, only: ng0 use ener, only: eself use ncprm ! implicit none real(kind=8), allocatable:: f(:),vscr(:), figl(:) real(kind=8) vpsum, rhopsum, xg, r2new, fint, ERF integer is, ir, ig ! ! ================================================================== ! calculation of gaussian selfinteraction ! ================================================================== eself=0. do is=1,nsp eself=eself+float(na(is))*zv(is)*zv(is)/rcmax(is) end do eself=eself/sqrt(2.*pi) write(6,1200) eself 1200 format(2x,'formf: eself=',f10.5) ! allocate(figl(ngsl)) allocate(f(mmaxx)) allocate(vscr(mmaxx)) ! do is=1,nsp ! ================================================================== ! fourier transform of local pp ! ================================================================== if(ipp(is).ne.3) then ! ! local potential given numerically on logarithmic mesh ! ! vscr(ir) = r*vscr(r) with the longe-range tail -Z/r subracted ! ! do ir=1,mesh(is) if(r(ir,is).le.10.0)then ! ! for r > 10 the potential is assumed to be zero ! vscr(ir)=0.5*rucore(ir,1,is) + & & zv(is)*ERF(r(ir,is)/rcmax(is)) else vscr(ir)=0.0 end if end do ! g=0 ! if (ng0.eq.2) then do ir=1,mesh(is) f(ir)=vscr(ir)*r(ir,is) end do if (ipp(is).eq.0) then call herman_skillman_int(mesh(is),cmesh(is),f,fint) else call simpson(mesh(is),f,rab(1,is),fint) end if ! vps(1,is)=fpi*fint/omega end if ! ! g>0 ! do ig=ng0,ngsl xg=sqrt(gl(ig))*tpiba do ir=1,mesh(is) f(ir)=vscr(ir)*sin(r(ir,is)*xg) end do if (ipp(is).eq.0) then call herman_skillman_int & & (mesh(is),cmesh(is),f,figl(ig)) else call simpson(mesh(is),f,rab(1,is),figl(ig)) end if end do ! do ig=ng0,ngs xg=sqrt(g(ig))*tpiba vps(ig,is)=fpi*figl(igl(ig))/(omega*xg) end do ! else ! ! bhs pseudopotentials can be fourier transformed analytically ! call formf_bhs(is) endif ! ! ================================================================== ! fourier transform of gaussian nuclear charge ! ================================================================== ! if (ng0.eq.2) rhops(1,is)=-zv(is)/omega r2new=0.25*tpiba2*rcmax(is)**2 do ig=ng0,ngs rhops(ig,is)=-zv(is)*exp(-r2new*g(ig))/omega end do ! vpsum =0.0 rhopsum=0.0 do ig=1,ngs vpsum =vpsum +vps (ig,is) rhopsum=rhopsum+rhops(ig,is) end do ! #ifdef PARA call reduce(1,vpsum) call reduce(1,rhopsum) #endif write(6,1250) vps(1,is),rhops(1,is) write(6,1300) vpsum,rhopsum ! end do 1250 format(2x,'formf: vps(g=0)=',f12.7,' rhops(g=0)=',f12.7) 1300 format(2x,'formf: sum_g vps(g)=',f12.7,' sum_g rhops(g)=',f12.7) !ATTENZIONE i termini in G=0 non e' detto che siano sul processore=1 ! deallocate(vscr) deallocate(f) deallocate(figl) ! return end !----------------------------------------------------------------------- subroutine formf_bhs(is) !----------------------------------------------------------------------- ! computes analytic form factors for BHS pseudopotentials ! use bhs use gvec use gvecs use cnst use ions use pseu use elct, only: ng0 use parm, only: omega ! implicit none ! input integer is ! local real(kind=8) e1, e2, emax, el, ql, par, sp, gps, sfp, & & r2max, r21, r22, r2l integer ib, ig ! ! r2max=rcmax(is)**2 r21=rc1(is)**2 r22=rc2(is)**2 ! ! ------------------------------------------------------------------ ! g=0 ! ------------------------------------------------------------------ if (ng0.eq.2) then gps=-zv(is)*pi*(-wrc2(is)*r22-wrc1(is)*r21+r2max)/omega sfp=0. do ib=1,3 r2l=rcl(ib,is,lloc(is))**2 ql=0.25*r2l*g(1)*tpiba2 el=exp(-ql) par=al(ib,is,lloc(is))+bl(ib,is,lloc(is))*r2l*(1.5-ql) sp=(pi*r2l)**1.5*el/omega sfp=sp*par+sfp end do vps(1,is)=gps+sfp end if ! ! ------------------------------------------------------------------ ! g>0 ! ------------------------------------------------------------------ do ig=ng0,ngs emax=exp(-0.25*r2max*g(ig)*tpiba2) e1=exp(-0.25*r21*g(ig)*tpiba2) e2=exp(-0.25*r22*g(ig)*tpiba2) gps=-zv(is)*(wrc1(is)*e1-emax+wrc2(is)*e2)/omega gps=gps*fpi/(g(ig)*tpiba2) sfp=0. do ib=1,3 r2l=rcl(ib,is,lloc(is))**2 ql=0.25*r2l*g(ig)*tpiba2 par = al(ib,is,lloc(is))+bl(ib,is,lloc(is))*r2l & & *(1.5-ql) sp=(pi*r2l)**1.5*exp(-ql)/omega sfp=sp*par+sfp end do vps(ig,is)=sfp+gps end do ! return end ! !----------------------------------------------------------------------- subroutine init0(nbeg,ndr,ndw,nomore,iprint & & ,delt,emass,emaec & & ,frice,grease,twall & & ,tortho,eps,max,reps,ntop & & ,trane,ampre & & ,tfor,fricp,greasp & & ,trhor,trhow,tvlocw & & ,tnosep,qnp,tempw & & ,tnosee,qne,ekincw) !----------------------------------------------------------------------- ! this subroutine reads control variables from unit 5 on most machines ! excepted t3e and origin for which the input is read from unit 9 ! (this ensures that the input is read on all nodes) ! ------------------------------------------------------------------ implicit none real(kind=8) delt, emass, emaec, grease, eps, reps, ampre, & & frice, fricp, greasp, qnp, tempw, qne, ekincw integer nbeg, ndr,ndw,nomore, iprint, max, ntop integer unit logical tortho,trane,twall,tnosee logical tfor,tnosep,trhor,trhow,tvlocw ! #if defined(T3E) || defined(ORIGIN) || defined(OSF1) unit=9 #else unit=5 #endif read(unit,*) nbeg,ndr,ndw,nomore,iprint read(unit,*) delt,emass,emaec read(unit,*) frice,grease,twall read(unit,*) tortho,eps,max,reps,ntop read(unit,*) trane,ampre read(unit,*) tfor,fricp,greasp read(unit,*) trhor,trhow,tvlocw read(unit,*) tnosep,qnp,tempw read(unit,*) tnosee,qne,ekincw ! ! ------------------------------------------------------- ! correct sloppy inputs ! if (.not.tortho) tnosee=.false. if(.not.tfor) then tnosep=.false. tnosee=.false. end if if (nbeg.eq.-2) trane=.true. if (nbeg.eq.0.and..not.trhor) trane=.false. ! ! -------------------------------------------------------- ! print out heading ! write(6,400) write(6,410) write(6,420) write(6,410) write(6,400) write(6,500) nbeg,nomore,iprint write(6,505) delt write(6,510) emass,emaec ! write(6,511) eps,max ! if (tnosee) frice = 0. write(6,509) write(6,514) frice,grease ! if(trhor)then write(6,720) endif ! if(.not.trhor.and.trhow)then write(6,721) endif ! if(tvlocw)then write(6,722) endif ! if(trane) then write(6,515) ampre endif write(6,516) ! if(tfor) then if(tnosep) fricp=0. write(6,520) write(6,522) fricp,greasp else write(6,518) endif ! if(tfor) then if(tnosep) then write(6,562) tempw,qnp else write(6,550) end if if(tnosee) then write(6,566) ekincw,qne end if end if ! ! 400 format('************************************', & & '************************************') 410 format('**** ', & & ' ****') 420 format('**** ab-initio molecular dynamics: ', & & ' car-parrinello vanderbilt bhs ****') 500 format(// & & ' nbeg=',i3,' nomore=',i7,3x,' iprint=',i4,/ & & ' reads from unit 91 writes on unit 92') 505 format(' time step = ',f9.4/) 510 format(' parameters for electron dynamics:'/ & & ' emass= ',f10.2,2x,'emaec= ',f10.2,'ry') 511 format(' orthog. with lagrange multipliers: eps=',e10.2, & & ' max=',i3) 509 format(' verlet algorithm for electron dynamics') 514 format(' with friction frice = ',f7.4,' , grease = ',f7.4) 720 format(' charge density is read from unit 47',/) 721 format(' charge density is written in unit 47',/) 722 format(' local potential is written in unit 46',/) 515 format(' initial random displacement of el. coordinates with ', & & ' amplitude=',f10.6,/ & & ' trane not to be used with mass preconditioning') 516 format(/) 518 format(' ions are not allowed to move'/) 520 format(' ions are allowed to move') 522 format(' ion dynamics with fricp = ',f7.4,' and greasp = ',f7.4) 550 format(' ion dynamics: the temperature is not controlled'//) 562 format(' ion dynamics with nose` temp. control:'/ & & ' temperature required=',f10.5,'(kelvin)',' nose` mass = ',& & f10.3//) 566 format(' electronic dynamics with nose` temp. control:'/ & & ' elec. kin. en. required=',f10.5,'(hartree)', & & ' nose` mass = ',f10.3//) return end !----------------------------------------------------------------------- subroutine nlinit !----------------------------------------------------------------------- ! this routine initalizes arrays beta qrad qq rhoc ! ! beta(ig,l,is) = 4pi/sqrt(omega) y^r(l,q^) ! int_0^inf dr r^2 j_l(qr) beta(l,is,r) ! ! qrad(igl,l,k,is) = 4pi/omega int_0^r dr r^2 j_l(qr) q(r,l,k,is) ! ! beta(g)_lm,is = (-i)^l*beta(ig,l,is) ! ! qq_ij=int_0^r q_ij(r)=omega*qg(g=0) ! use gvec use cvan use core use cnst use ions use parm use elct use ncprm use qradb_mod use qgb_mod use gvecb use parmb ! implicit none integer lmax, is, il, l, lp, ig, ir, iv, jv, ijv, ind, lm real(kind=8) ylmr, xg, c, fac real(kind=8), allocatable:: fint(:), jl(:), betagl(:) complex(kind=8), allocatable:: qgbs(:) external ylmr ! ! read the local (pseudo)potential and wavefunctions ! call readpp ! ------------------------------------------------------------------ ! find number of beta functions per species, max dimensions, ! total number of beta functions (all and Vanderbilt only) ! ------------------------------------------------------------------ lmax=0 nhx=0 nhsa=0 nhsavb=0 do is=1,nsp ind=0 do iv=1,nbeta(is) lmax =max(lmax,lll(iv,is)) ind=ind+2*lll(iv,is)+1 end do nh(is)=ind nhx=max(nhx,nh(is)) ish(is)=nhsa nhsa=nhsa+na(is)*nh(is) if(ipp(is).le.1) nhsavb=nhsavb+na(is)*nh(is) end do lmax=lmax+1 if (lmax.gt.3) call error('nlinit ',' l > 3 ,l= ',lmax) if (nhsa.le.0) call error('nlinit ','not implemented ?',nhsa) ! ! initialize array ap ! call aainit(lmax,ap,lpx,lpl) ! allocate(beta(ngw,nhx,nsp)) allocate(qradb(nglb,nbrx,nbrx,lx,nsp)) allocate(qgb(ngb,nhx*(nhx+1)/2,nsp)) allocate(qq(nhx,nhx,nsp)) allocate(dvan(nhx,nhx,nsp)) allocate(rhoc(ngb,nsp)) allocate(nhtol(nhx,nsp)) allocate(indv (nhx,nsp)) allocate(indlm(nhx,nsp)) ! call zero(nglb*nbrx*nbrx*lx*nsp,qradb) call zero(nhx*nhx*nsp,qq) call zero(nhx*nhx*nsp,dvan) ! ! qgbs, betagl are used as work space ! allocate(qgbs(ngb)) allocate(betagl(ngl)) ! ! ------------------------------------------------------------------ ! definition of indices nhtol, indv, indlm ! ------------------------------------------------------------------ do is=1,nsp ind=0 do iv=1,nbeta(is) if(lll(iv,is).eq.0)then lm=0 else if (lll(iv,is).eq.1) then lm=1 else if (lll(iv,is).eq.2) then lm=4 endif do il=1,2*lll(iv,is)+1 lm=lm+1 ind=ind+1 indlm(ind,is)=lm nhtol(ind,is)=lll(iv,is) indv(ind,is)=iv end do end do end do ! do is=1,nsp ! write(6,*) ' nlinit: nh(is) = ',nh(is) ! do iv=1,nh(is) write(6,901) iv,indv(iv,is),nhtol(iv,is) end do 901 format(2x,i2,' indv= ',i2,' ang. mom= ',i2) ! allocate(fint(kkbeta(is))) allocate(jl(kkbeta(is))) ! if(ipp(is).le.1) then ! =============================================================== ! initialization for vanderbilt species ! =============================================================== ! ! qqq and beta are now indexed and taken in the same order ! as vanderbilts ppot-code prints them out ! ! --------------------------------------------------------------- ! calculation of array qradb(iglb,iv,jv,is) ! --------------------------------------------------------------- write(6,*) ' qradb ' c=fpi/omegab write(6,*) ' nlinit nh(is),nglb,is,kkbeta,nqlc = ', & & nh(is),nglb,is,kkbeta(is),nqlc(is) do l=1,nqlc(is) do il=1,nglb xg=sqrt(glb(il))*tpibab call bess(xg,l,kkbeta(is),r(1,is),jl) do iv= 1,nbeta(is) do jv=iv,nbeta(is) ! ! note qrl(r)=r^2*q(r) ! do ir=1,kkbeta(is) fint(ir)=qrl(ir,iv,jv,l,is)*jl(ir) end do if (ipp(is).eq.0) then call herman_skillman_int & & (kkbeta(is),cmesh(is),fint,qradb(il,iv,jv,l,is)) else call simpson & & (kkbeta(is),fint,rab(1,is),qradb(il,iv,jv,l,is)) end if qradb(il,iv,jv,l,is)=c*qradb(il,iv,jv,l,is) qradb(il,jv,iv,l,is)=qradb(il,iv,jv,l,is) end do end do end do end do ! ! --------------------------------------------------------------- ! calculation of array qq ! --------------------------------------------------------------- write(6,*) ' qq ' do iv= 1,nh(is) do jv=iv,nh(is) call qvan2b(1,iv,jv,is,qgbs) qq(iv,jv,is)=omegab*real(qgbs(1)) qq(jv,iv,is)=qq(iv,jv,is) end do end do ! ! --------------------------------------------------------------- ! stocking of qgb(igb,ijv,is) ! --------------------------------------------------------------- ijv=0 do iv= 1,nh(is) do jv=iv,nh(is) ! ! compact indices because qgb is symmetric: ! ivjv: 11 12 13 ... 22 23... ! ijv : 1 2 3 ... ! ijv=ijv+1 call qvan2b(ngb,iv,jv,is,qgbs) do ig=1,ngb qgb(ig,ijv,is)=qgbs(ig) end do end do end do ! write(6,*) write(6,*) write(6,'(20x,a)') ' qqq ' do iv=1,nbeta(is) write(6,'(8f9.4)') (qqq(iv,jv,is),jv=1,nbeta(is)) end do write(6,*) ! endif ! ! =============================================================== ! initialization that is common to all PPs ! =============================================================== ! fac converts ry to hartree ! if (ipp(is).eq.3) then fac=1.0 else fac=0.5 end if ! ! --------------------------------------------------------------- ! calculation of array beta(ig,iv,is) ! --------------------------------------------------------------- if(nh(is).gt.0) write(6,*) ' beta ' c=fpi/sqrt(omega) do iv=1,nh(is) l=nhtol(iv,is)+1 lp=indlm(iv,is) do ig=1,ngwl xg=sqrt(gl(ig))*tpiba call bess(xg,l,kkbeta(is),r(1,is),jl) ! ! betar(ir)=r*beta(r) ! do ir=1,kkbeta(is) fint(ir)=r(ir,is)*betar(ir,indv(iv,is),is)*jl(ir) end do if (ipp(is).eq.0) then call herman_skillman_int & & (kkbeta(is),cmesh(is),fint,betagl(ig)) else call simpson(kkbeta(is),fint,rab(1,is),betagl(ig)) end if end do do ig=1,ngw beta(ig,iv,is)=c*ylmr(lp,ig)*betagl(igl(ig)) end do end do ! ! --------------------------------------------------------------- ! calculate array dvan(iv,jv,is) ! --------------------------------------------------------------- do iv=1,nh(is) do jv=1,nh(is) if ( indlm(iv,is).eq.indlm(jv,is) ) then dvan(iv,jv,is)=fac*dion(indv(iv,is),indv(jv,is),is) endif end do end do ! write(6,*) write(6,'(20x,a)') ' dion ' do iv=1,nbeta(is) write(6,'(8f9.4)') (fac*dion(iv,jv,is),jv=1,nbeta(is)) end do ! --------------------------------------------------------------- ! non-linear core-correction ( rhoc(ig,is) ) ! ---------------------------------------------------------------- if(ifpcor(is).eq.1) then l=1 c=1.0/omegab do ig=1,nglb xg=sqrt(glb(ig))*tpibab call bess(xg,l,kkbeta(is),r(1,is),jl) do ir=1,kkbeta(is) fint(ir)=rscore(ir,is)*jl(ir) end do call simpson(kkbeta(is),fint,rab(1,is),qgbs(ig)) end do do ig=1,ngb rhoc(ig,is)=c*qgbs(iglb(ig)) end do write(6,'(a,f12.8)') & & ' integrated core charge= ',omegab*rhoc(1,is) endif ! deallocate(jl) deallocate(fint) end do ! deallocate(betagl) deallocate(qgbs) ! return end !------------------------------------------------------------------------- subroutine qvan2b(ngy,iv,jv,is,qg) !-------------------------------------------------------------------------- ! q(g,l,k) = sum_lm (-i)^l ap(lm,l,k) yr_lm(g^) qrad(g,l,l,k) ! use qradb_mod use cvan use gvecb ! implicit none integer ngy,iv,jv,is integer ivs, jvs, ivl, jvl, i, l, lp, ig complex(kind=8) qg(ngb), sig real(kind=8), allocatable:: ylm(:) ! ! iv = 1..8 s_1 p_x1 p_z1 p_y1 s_2 p_x2 p_z2 p_y2 ! ivs = 1..4 s_1 s_2 p_1 p_2 ! ivl = 1..4 s p_x p_z p_y ! ivs=indv(iv,is) jvs=indv(jv,is) ivl=indlm(iv,is) jvl=indlm(jv,is) ! allocate(ylm(ngb)) ! do i=1,ngb qg(i) =(0.,0.) end do ! ! lpx = max number of allowed y_lm ! lp = composite lm to identify them ! do i=1,lpx(ivl,jvl) lp=lpl(ivl,jvl,i) ! ! extraction of angular momentum l from lp: ! if (lp.eq.1) then l=1 else if ((lp.ge.2) .and. (lp.le.4)) then l=2 else if ((lp.ge.5) .and. (lp.le.9)) then l=3 else if ((lp.ge.10).and.(lp.le.16)) then l=4 else if ((lp.ge.17).and.(lp.le.25)) then l=5 else if (lp.ge.26) then call error(' qvanb ',' lp.ge.26 ',lp) endif ! ! sig= (-i)^l ! sig=(0.,-1.)**(l-1) call ylmr2b(lp,ngy,ngb,gxnb,ylm) sig=sig*ap(lp,ivl,jvl) do ig=1,ngy qg(ig)=qg(ig)+sig*ylm(ig)*qradb(iglb(ig),ivs,jvs,l,is) end do end do ! deallocate(ylm) ! return end !----------------------------------------------------------------------- subroutine readfile & & ( flag,nx,nax,nsp,ndr,nfi,ngw,n,nr,c0,cm,tau0,vol0,cdm0,cdmm,& & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) !----------------------------------------------------------------------- implicit none integer flag,nx,nax,nsp,ndr,nfi,ngw,n,nr complex(kind=8) c0(ngw,n), cm(ngw,n) real(kind=8) cdm0(3),cdmm(3),acc(10), a1(3),a2(3),a3(3), alat real(kind=8) lambda(nx,nx), lambdam(nx,nx) real(kind=8) xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,ekincm,volm,vol0 real(kind=8) taum(3,nax,nsp),tau0(3,nax,nsp),taui(3,nax,nsp) real(kind=8) vel(3,nax,nsp) ! integer i,ig,j,ngwr complex(kind=8) dummy ! rewind ndr if (flag.eq.0) then write(6,'((a,i3,a))') ' ## reading from file ',ndr,' only c0 ##' else write(6,'((a,i3))') ' ## reading from file ',ndr end if read(ndr) nfi,ngwr,nr if(flag.eq.0) then if(ngwr.gt.ngw) write(6,*) ' ## warning ## ngwr.gt.ngw ',ngwr,ngw else if(ngwr.gt.ngw) call error('readfile','ngwr.ne.ngw',ngwr) endif if(ngwr.le.ngw) then read(ndr) ((c0(ig,i),ig=1,ngwr),i=1,min(nr,n)) else read(ndr) ((c0(ig,i),ig=1,ngw),(dummy,ig=ngw+1,ngwr),i=1,min(nr,n)) endif if (flag.eq.0) return ! read(ndr) ((cm(ig,i),ig=1,ngwr),i=1,min(nr,n)) read(ndr) tau0,vol0,cdm0,cdmm read(ndr) taum,volm,acc,taui,a1,a2,a3,alat read(ndr) ((lambda(i,j),i=1,n),j=1,n) read(ndr) ((lambdam(i,j),i=1,n),j=1,n) read(ndr) xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm ! return end !----------------------------------------------------------------------- subroutine rhoofr (nfi,iprint,tbuff,trhor,trhow, & & c,irb,eigrb,bec,rhovan,rhor,rhog,rhos) !----------------------------------------------------------------------- ! the normalized electron density rhor in real space ! the kinetic energy ekin ! subroutine uses complex fft so it computes two ft's ! simultaneously ! ! rho_i,ij = sum_n < beta_i,i | psi_n >< psi_n | beta_i,j > ! < psi_n | beta_i,i > = c_n(0) beta_i,i(0) + ! 2 sum_g> re(c_n*(g) (-i)**l beta_i,i(g) e^-ig.r_i) ! ! e_v = sum_i,ij rho_i,ij d^ion_is,ji ! use ions use gvec use gvecs use gvecb, only: ngb use cvan use parm use parms use elct use cnst use pseu use ener use control use work1 ! implicit none ! real(kind=8) bec(nhsa,n), rhovan(nat,nhx,nhx,nspin) real(kind=8) rhor(nnr,nspin), rhos(nnrs,nspin) complex(kind=8) eigrb(ngb,nas,nsp), c(ngw,nx), rhog(ng,nspin) integer irb(3,nax,nsx), nfi, iprint logical tbuff,trhor,trhow ! local variables integer iss, isup, isdw, iss1, iss2, ios, i, ia, is, iv, jv, & & ir, ig, ism, isa, inl, jnl real(kind=8) sums(2), rsumr(2), rsumg(2), sum, sa1, sa2, SSUM real(kind=8) rnegsum, rmin, rmax, rsum real(kind=8) sk(nx) ! automatic array complex(kind=8), pointer:: psi(:), psis(:) external SSUM ! complex(kind=8) ci, fp, fm ! ! call tictac(4,0) psi => wrk1 psis=> wrk1 ci=(0.0,1.0) do iss=1,nspin call zero(nnr,rhor(1,iss)) call zero(nnrs,rhos(1,iss)) call zero(2*ng,rhog(1,iss)) end do ! ! ================================================================== ! calculation of kinetic energy ekin ! ================================================================== do i=1,n sk(i)=0.0 do ig=ng0,ngw sk(i)=sk(i)+real(conjg(c(ig,i))*c(ig,i))*g(ig) end do end do #ifdef PARA call reduce(n,sk) #endif ekin=0.0 do i=1,n ekin=ekin+f(i)*sk(i) end do ekin=ekin*tpiba2 ! ! ================================================================== ! calculation of non-local energy ! ================================================================== enl=0. do is=1,nsp do iv= 1,nh(is) do jv=iv,nh(is) isa=0 do ism=1,is-1 isa=isa+na(ism) end do do ia=1,na(is) inl=ish(is)+(iv-1)*na(is)+ia jnl=ish(is)+(jv-1)*na(is)+ia isa=isa+1 do iss=1,nspin sums(iss)=0. end do do i=1,n iss=ispin(i) sums(iss) = sums(iss) +f(i)*bec(inl,i)*bec(jnl,i) end do sum=0. do iss=1,nspin rhovan(isa,iv,jv,iss) = sums(iss) rhovan(isa,jv,iv,iss) = sums(iss) sum=sum+sums(iss) end do if(iv.ne.jv) sum=2.*sum enl=enl+sum*dvan(jv,iv,is) end do end do end do end do ! if(trhor) then ! ================================================================== ! charge density is read from unit 47 ! ================================================================== #ifdef PARA call read_rho(47,nspin,rhor) #else read(47) ((rhor(ir,iss),ir=1,nnr),iss=1,nspin) #endif rewind(47) if(nspin.eq.1) then iss=1 do ir=1,nnr psi(ir)=cmplx(rhor(ir,iss),0.) end do call fwfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ig=1,ng rhog(ig,iss)=psi(np(ig)) end do else isup=1 isdw=2 do ir=1,nnr psi(ir)=cmplx(rhor(ir,isup),rhor(ir,isdw)) end do call fwfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ig=1,ng fp = psi(np(ig)) + psi(nm(ig)) fm = psi(np(ig)) - psi(nm(ig)) rhog(ig,isup)=0.5*cmplx( real(fp),aimag(fm)) rhog(ig,isdw)=0.5*cmplx(aimag(fp),-real(fm)) end do endif ! else ! ================================================================== ! self-consistent charge ! ================================================================== ! ! important: if n is odd then c(*,n+1)=0. ! if (mod(n,2).ne.0) then do ig=1,ngw c(ig,n+1)=(0.,0.) end do endif ! call tictac(12,0) do i=1,n,2 call zero(2*nnrs,psis) do ig=1,ngw psis(nms(ig))=conjg(c(ig,i)) + ci*conjg(c(ig,i+1)) psis(nps(ig))=c(ig,i)+ci*c(ig,i+1) end do ! call ivfftw(psis,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx) ! ! wavefunctions in unit 21 ! !!!!#if defined(CRAYY) || defined(NEC) ! A.BONGIORNO #if defined(CRAYY) if(tbuff) buffer out(21,0) (psis(1),psis(nnrs)) #else if(tbuff) write(21,iostat=ios) psis #endif iss1=ispin(i) sa1=f(i)/omega if (i.ne.n) then iss2=ispin(i+1) sa2=f(i+1)/omega else iss2=iss1 sa2=0.0 end if do ir=1,nnrs rhos(ir,iss1)=rhos(ir,iss1) + sa1*( real(psis(ir)))**2 rhos(ir,iss2)=rhos(ir,iss2) + sa2*(aimag(psis(ir)))**2 end do ! ! buffer 21 ! if(tbuff) then !!!!#if defined(CRAYY) || defined(NEC) ! A.BONGIORNO #if defined(CRAYY) ios=unit(21) #endif if(ios.ne.0) call error & & (' rhoofr',' error in writing unit 21',ios) endif ! end do call tictac(12,1) ! if(tbuff) rewind 21 ! ! smooth charge in g-space is put into rhog(ig) ! if(nspin.eq.1)then iss=1 do ir=1,nnrs psis(ir)=cmplx(rhos(ir,iss),0.) end do call fwffts(psis,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx) do ig=1,ngs rhog(ig,iss)=psis(nps(ig)) end do else isup=1 isdw=2 do ir=1,nnrs psis(ir)=cmplx(rhos(ir,isup),rhos(ir,isdw)) end do call fwffts(psis,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx) do ig=1,ngs fp= psis(nps(ig)) + psis(nms(ig)) fm= psis(nps(ig)) - psis(nms(ig)) rhog(ig,isup)=0.5*cmplx( real(fp),aimag(fm)) rhog(ig,isdw)=0.5*cmplx(aimag(fp),-real(fm)) end do endif ! if(nspin.eq.1)then ! ================================================================== ! case nspin=1 ! ------------------------------------------------------------------ iss=1 call zero(2*nnr,psi) do ig=1,ngs psi(nm(ig))=conjg(rhog(ig,iss)) psi(np(ig))= rhog(ig,iss) end do call invfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ir=1,nnr rhor(ir,iss)=real(psi(ir)) end do else ! ================================================================== ! case nspin=2 ! ------------------------------------------------------------------ isup=1 isdw=2 call zero(2*nnr,psi) do ig=1,ngs psi(nm(ig))=conjg(rhog(ig,isup))+ ci*conjg(rhog(ig,isdw)) psi(np(ig))=rhog(ig,isup) + ci*rhog(ig,isdw) end do call invfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ir=1,nnr rhor(ir,isup)= real(psi(ir)) rhor(ir,isdw)=aimag(psi(ir)) end do endif ! if(iprsta.ge.3)then do iss=1,nspin rsumg(iss)=omega*real(rhog(1,iss)) rsumr(iss)=SSUM(nnr,rhor(1,iss),1)*omega/dfloat(nr1*nr2*nr3) end do #ifdef PARA if (ng0.ne.2) then ! in the parallel case, only one processor has G=0 ! do iss=1,nspin rsumg(iss)=0.0 end do end if call reduce(nspin,rsumg) call reduce(nspin,rsumr) #endif if (nspin.eq.1) then write(6,1) rsumg(1),rsumr(1) else write(6,2) (rsumg(iss),iss=1,nspin),(rsumr(iss),iss=1,nspin) endif endif ! ================================================================== ! ! vanderbilt contribution to the charge density ! call rhov(irb,eigrb,rhovan,rhog,rhor) ! #ifdef PARA if(trhow) call write_rho(47,nspin,rhor) #else if(trhow) write(47) ((rhor(ir,iss),ir=1,nnr),iss=1,nspin) #endif endif ! ======================================endif for trhor============= ! ! here to check the integral of the charge density ! if(iprsta.ge.2) then call checkrho(nnr,nspin,rhor,rmin,rmax,rsum,rnegsum) rnegsum=rnegsum*omega/dfloat(nr1*nr2*nr3) rsum=rsum*omega/dfloat(nr1*nr2*nr3) write(6,'(a,4(1x,f12.6))') & & ' rhoofr: rmin rmax rnegsum rsum ',rmin,rmax,rnegsum,rsum end if ! if(nfi.eq.0.or.mod(nfi-1,iprint).eq.0) then do iss=1,nspin rsumg(iss)=omega*real(rhog(1,iss)) rsumr(iss)=SSUM(nnr,rhor(1,iss),1)*omega/dfloat(nr1*nr2*nr3) end do #ifdef PARA if (ng0.ne.2) then ! in the parallel case, only one processor has G=0 ! do iss=1,nspin rsumg(iss)=0.0 end do end if call reduce(nspin,rsumg) call reduce(nspin,rsumr) #endif if (nspin.eq.1) then write(6,1) rsumg(1),rsumr(1) else if(iprsta.ge.3) & & write(6,2) rsumg(1),rsumg(2),rsumr(1),rsumr(2) write(6,1) rsumg(1)+rsumg(2),rsumr(1)+rsumr(2) endif endif ! call tictac(4,1) ! 2 format(//' subroutine rhoofr: total integrated electronic', & & ' density'/' in g-space =',f10.6,2x,f10.6,4x, & & ' in r-space =',f10.6,2x,f10.6) 1 format(//' subroutine rhoofr: total integrated electronic', & & ' density'/' in g-space =',f10.6,4x, & & ' in r-space =',f10.6) ! return end !------------------------------------------------------------------------- subroutine vofrho(nfi,iprint,rhor,rhog,rhos,tfor,tfirst,tlast, & & ei1,ei2,ei3,sfac,tau0,tvlocw,fion) !----------------------------------------------------------------------- ! computes: the one-particle potential v in real space, ! the total energy etot, ! the forces fion acting on the ions, ! rhor input : electronic charge on dense real space grid ! rhog input : electronic charge in g space till ng ! rhos input : electronic charge on smooth real space grid ! rhor output: total potential on dense real space grid ! rhos output: total potential on smooth real space grid ! use ions use gvec use gvecs use parm use parms use elct use cnst use ener use pseu use core use ncprm use gvecb use parmb use dft_mod use control use work1 use work_box ! implicit none ! logical tfor,tlast,tfirst,tvlocw integer nfi, iprint real(kind=8) rhor(nnr,nspin), rhos(nnrs,nspin), fion(3,nax,nsx) real(kind=8) tau0(3,nax,nsp) complex(kind=8) ei1(-nr1:nr1,nas,nsp), ei2(-nr2:nr2,nas,nsp), & & ei3(-nr3:nr3,nas,nsp), rhog(ng,nspin), sfac(ngs,nsp) integer iss, isup, isdw, ism, isa, ix, ig, ir, i,j,k, & & is, ia, l, m real(kind=8) vave, ebac, wz, SSUM complex(kind=8) tsum(3,nax,nsx), fp, fm, ci, eh, eps, CSUM complex(kind=8), pointer:: v(:), vs(:) complex(kind=8), allocatable:: rhotemp(:), vtemp(:) external SSUM ! call tictac(5,0) ci=(0.,1.) ! ! wz = factor for g.neq.0 because of c*(g)=c(-g) ! wz = 2.0 v => wrk1 vs=> wrk1 allocate(rhotemp(ng)) allocate(vtemp(ng)) ! ! first routine in which fion is calculated: annihilation ! do is=1,nsp do ia=1,nax do i=1,3 fion(i,ia,is)=0. tsum(i,ia,is)=(0.,0.) end do end do end do ! =================================================================== ! calculation forces, ionic term ! ------------------------------------------------------------------- if(tfor.or.tfirst) call force_ion(tau0,esr,fion) if(nspin.eq.1)then iss=1 do i=1,ng rhotemp(i)=rhog(i,iss) end do else isup=1 isdw=2 do i=1,ng rhotemp(i)=rhog(i,isup)+rhog(i,isdw) end do endif ! =================================================================== ! calculation local potential energy ! ------------------------------------------------------------------- do i=1,ng vtemp(i)=(0.,0.) end do do is=1,nsp do ig=1,ngs vtemp(ig)=vtemp(ig)+conjg(rhotemp(ig))*sfac(ig,is)*vps(ig,is) end do end do ! eps=wz*CSUM(ngs,vtemp,1) if (ng0.eq.2) eps=eps-vtemp(1) #ifdef PARA call reduce(1,eps) #endif ! =================================================================== ! calculation hartree energy ! ------------------------------------------------------------------- do is=1,nsp do ig=1,ngs rhotemp(ig)=rhotemp(ig)+sfac(ig,is)*rhops(ig,is) end do end do vtemp(1)=0.0 do ig=ng0,ng vtemp(ig)=conjg(rhotemp(ig))*rhotemp(ig)/g(ig) end do ! eh=CSUM(ng,vtemp,1)*wz*.5*fpi/tpiba2 #ifdef PARA call reduce(1,eh) #endif ! =================================================================== ! calculation forces, local pseudopotential contribution ! ------------------------------------------------------------------- if(tfor) call force_ps(rhotemp,rhog,vtemp,ei1,ei2,ei3,tsum) ! =================================================================== ! calculation hartree potential ! ------------------------------------------------------------------- ! vtemp = v_tot(g) - v_xc(g) ! if (ng0.eq.2) vtemp(1)=(0.,0.) do ig=ng0,ng vtemp(ig)=rhotemp(ig)*fpi/(tpiba2*g(ig)) end do ! =================================================================== ! calculation local pseudo potential ! ------------------------------------------------------------------- do is=1,nsp do ig=1,ngs vtemp(ig)=vtemp(ig)+sfac(ig,is)*vps(ig,is) end do end do ! =================================================================== ! xc energy and potential with nonlinear core-correction ! ------------------------------------------------------------------- ! (not yet implemented) ! =================================================================== ! calculation exchange and correlation energy and potential ! ------------------------------------------------------------------- ! call exch_corr(nspin,rhog,rhor,exc) ! ! rhor contains the xc potential in r-space ! ! =================================================================== ! fourier transform of xc potential to g-space (dense grid) ! ------------------------------------------------------------------- if(nspin.eq.1)then iss=1 do ig=1,ng rhog(ig,iss)=vtemp(ig) end do else isup=1 isdw=2 do ig=1,ng rhog(ig,isup)=vtemp(ig) rhog(ig,isdw)=vtemp(ig) end do endif ! if(nspin.eq.1)then iss=1 do ir=1,nnr v(ir)=cmplx(rhor(ir,iss),0.0) end do ! ! v_xc(r) --> v_xc(g) ! call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x) ! do ig=1,ng rhog(ig,iss)=rhog(ig,iss)+v(np(ig)) end do ! ! v_tot(g) = (v_tot(g) - v_xc(g)) +v_xc(g) ! rhog contains the total potential in g-space ! else isup=1 isdw=2 do ir=1,nnr v(ir)=cmplx(rhor(ir,isup),rhor(ir,isdw)) end do call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ig=1,ng fp=v(np(ig))+v(nm(ig)) fm=v(np(ig))-v(nm(ig)) rhog(ig,isup)=rhog(ig,isup)+0.5*cmplx( real(fp),aimag(fm)) rhog(ig,isdw)=rhog(ig,isdw)+0.5*cmplx(aimag(fp),-real(fm)) end do endif ! if(tfor) then call force_cc(ei1,ei2,ei3,tsum) #ifdef PARA call reduce(2*3*nax*nsp,tsum) #endif ! ! contribution to fion ! do k=1,3 do is=1,nsp do ia=1,na(is) fion(k,ia,is) = fion(k,ia,is) + & & real(tsum(k,ia,is))*omega*tpiba end do end do end do endif ! =================================================================== ! fourier transform of total potential to r-space (dense grid) ! ------------------------------------------------------------------- call zero(2*nnr,v) if(nspin.eq.1)then iss=1 do ig=1,ng v(nm(ig))=conjg(rhog(ig,iss)) v(np(ig))=rhog(ig,iss) end do ! ! v(g) --> v(r) ! call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x) ! do ir=1,nnr rhor(ir,iss)=real(v(ir)) end do ! ! calculation of average potential ! vave=SSUM(nnr,rhor(1,iss),1)/dfloat(nr1*nr2*nr3) else isup=1 isdw=2 do ig=1,ng v(nm(ig))=conjg(rhog(ig,isup)) & & + ci*conjg(rhog(ig,isdw)) v(np(ig))=rhog(ig,isup) + ci*rhog(ig,isdw) end do call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x) do ir=1,nnr rhor(ir,isup)= real(v(ir)) rhor(ir,isdw)=aimag(v(ir)) end do ! ! calculation of average potential ! vave=(SSUM(nnr,rhor(1,isup),1)+SSUM(nnr,rhor(1,isdw),1)) & & /2.0/dfloat(nr1*nr2*nr3) endif #ifdef PARA call reduce(1,vave) #endif ! =================================================================== ! fourier transform of total potential to r-space (smooth grid) ! ------------------------------------------------------------------- call zero(2*nnrs,vs) if(nspin.eq.1)then iss=1 do ig=1,ngs vs(nms(ig))=conjg(rhog(ig,iss)) vs(nps(ig))=rhog(ig,iss) end do ! call ivffts(vs,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx) ! do ir=1,nnrs rhos(ir,iss)=real(vs(ir)) end do else isup=1 isdw=2 do ig=1,ngs vs(nms(ig))=conjg(rhog(ig,isup)) & & + ci*conjg(rhog(ig,isdw)) vs(nps(ig))=rhog(ig,isup) + ci*rhog(ig,isdw) end do call ivffts(vs,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx) do ir=1,nnrs rhos(ir,isup)= real(vs(ir)) rhos(ir,isdw)=aimag(vs(ir)) end do endif ! ebac=vave*qbac ebac=0.0 ! eht=real(eh)*omega+esr-eself epseu=real(eps)*omega ! ! etot is the total energy ; enl was calculated in rhoofr ! etot=ekin+eht+epseu+enl+exc+ebac ! if(tvlocw.and.tlast)then #ifdef PARA call write_rho(46,nspin,rhor) #else write(46) ((rhor(ir,iss),ir=1,nnr),iss=1,nspin) #endif endif ! deallocate(vtemp) deallocate(rhotemp) ! if(nfi.eq.0.or.tfirst.or.tlast.or.mod(nfi-1,iprint).eq.0) & & write(6,1) etot,ekin,eht,esr,eself,epseu,enl,exc,vave call tictac(5,1) ! 1 format(//' total energy = ',f14.5,' a.u.'/ & & ' kinetic energy = ',f14.5,' a.u.'/ & & ' electrostatic energy = ',f14.5,' a.u.'/ & & ' esr = ',f14.5,' a.u.'/ & & ' eself = ',f14.5,' a.u.'/ & & ' pseudopotential energy = ',f14.5,' a.u.'/ & & ' n-l pseudopotential energy = ',f14.5,' a.u.'/ & & ' exchange-correlation energy = ',f14.5,' a.u.'/ & & ' average potential = ',f14.5,' a.u.'//) ! return end !----------------------------------------------------------------------- subroutine writefile & & ( nx,nax,nsp,ndw,nfi,ngw,n,c0,cm,tau0,vol0,cdm0,cdmm, & & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) !----------------------------------------------------------------------- ! implicit none integer nx,nax,nsp,ndw,nfi,ngw,n complex(kind=8) c0(ngw,n), cm(ngw,n) real(kind=8) cdm0(3),cdmm(3),acc(10), a1(3),a2(3),a3(3), alat real(kind=8) lambda(nx,nx), lambdam(nx,nx) real(kind=8) xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,ekincm,volm,vol0 real(kind=8) taum(3,nax,nsp),tau0(3,nax,nsp),taui(3,nax,nsp) real(kind=8) vel(3,nax,nsp) ! integer i,ig,j ! rewind ndw write(ndw) nfi,ngw,n write(ndw) ((c0(ig,i),ig=1,ngw),i=1,n) write(ndw) ((cm(ig,i),ig=1,ngw),i=1,n) write(ndw) tau0,vol0,cdm0,cdmm write(ndw) taum,volm,acc,taui,a1,a2,a3,alat write(ndw) ((lambda(i,j),i=1,n),j=1,n) write(ndw) ((lambdam(i,j),i=1,n),j=1,n) write(ndw) xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm ! return end ! !---------------------------------------------------------------------- subroutine print_times !---------------------------------------------------------------------- ! use timex_mod ! implicit none integer i dimension routine(maxclock) character*10 routine data routine / ' full ', & & ' phbox ', & & ' strucf ', & & ' rhoofr ', & & ' vofrho ', & & ' newd ', & & ' dforce ', & & ' nlfor ', & & ' calphi ', & & ' r_iter ', & & ' ortho ', & & ' wfs**2 ', & & ' rhov ', & & ' fft ', & & ' ffts+w ', & & ' fftb ', & & ' 2*nlsm1 ', & & ' v*wf ', & & ' abc ', & & ' rsg ', & & ' x_iter ', & & ' nlsm2 ', & & ' nl_dforce', & & ' becp-phi ', & & ' bhs ', & & ' prefor ', & & 'setfftpara', & & 'fftscatter', & & ' reduce ', & & ' '/ ! write(6,*) ' routine calls totcpu avgcpu tottime avgtime' write(6,*) do i=1, maxclock if (ntimes(i).gt.0) write(6,30) routine(i), ntimes(i), & & cputime(i), cputime(i)/ntimes(i), & & elapsed(i), elapsed(i)/ntimes(i) end do 30 format(a10,i7,f8.1,f8.3,f8.1,f8.3) ! return end