TABLE OF CONTENTS
- ABINIT/e_eigen
- ABINIT/vtorho
- ABINIT/wvl_comm_eigen
- ABINIT/wvl_nscf_loop
- ABINIT/wvl_nscf_loop_bigdft
- ABINIT/wvl_occ
- ABINIT/wvl_occ_bigdft
ABINIT/e_eigen [ Functions ]
NAME
e_eigen
FUNCTION
Computes eigenvalues energy from eigen, occ, kpt, wtk
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
eigen(nkpt*nsppol)=eigenvalues mband= maximum number of bands nband(nkpt*nsppol)= number of bands for each k-point and spin nkpt= number of k-points nsppol= number of spin polarization occ(mband*nkpt*nsppol)=occupations wtk(nkpt)= k-point weights
OUTPUT
e_eigenvalues= eigenvalues energy
SIDE EFFECTS
NOTES
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
2053 subroutine e_eigen(eigen,e_eigenvalues,mband,nband,nkpt,nsppol,occ,wtk) 2054 2055 use defs_basis 2056 use m_errors 2057 2058 !This section has been created automatically by the script Abilint (TD). 2059 !Do not modify the following lines by hand. 2060 #undef ABI_FUNC 2061 #define ABI_FUNC 'e_eigen' 2062 !End of the abilint section 2063 2064 implicit none 2065 2066 !Arguments ------------------------------------ 2067 integer , intent(in) :: mband,nkpt,nsppol 2068 integer , intent(in) :: nband(nkpt*nsppol) 2069 real(dp) , intent(in) :: eigen(mband*nkpt*nsppol) 2070 real(dp) , intent(in) :: occ(mband*nkpt*nsppol) 2071 real(dp) , intent(in) :: wtk(nkpt) 2072 real(dp) , intent(out) :: e_eigenvalues 2073 2074 !Local variables------------------------------- 2075 integer :: ib,iband,ii,ikpt,isppol,nband_k 2076 real(dp) :: wtk_k 2077 2078 ! ************************************************************************* 2079 2080 DBG_ENTER("COLL") 2081 ii=0;ib=0 2082 do isppol=1,nsppol 2083 do ikpt=1,nkpt 2084 ii=ii+1 2085 nband_k=nband(ii) ; wtk_k=wtk(ii) 2086 do iband=1,nband_k 2087 ib=ib+1 2088 if(abs(occ(ib)) > tol8) then 2089 e_eigenvalues = e_eigenvalues + wtk_k*occ(ib)*eigen(ib) 2090 end if 2091 end do 2092 end do 2093 end do 2094 2095 DBG_EXIT("COLL") 2096 2097 end subroutine e_eigen
ABINIT/vtorho [ Functions ]
NAME
vtorho
FUNCTION
This routine compute the new density from a fixed potential (vtrial) but might also simply compute eigenvectors and eigenvalues. The main part of it is a wf update over all k points.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MF, AR, MM, MT, FJ, MB, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
afford=used to dimension susmat atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cpus= cpu time limit in seconds dbl_nnsclo=if 1, will double the value of dtset%nnsclo dielop= if positive, the dielectric matrix must be computed. dielstrt=number of the step at which the dielectric preconditioning begins. dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only) dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem =number of k points treated by this node. | mpw=maximum dimensioned size of npw | nfft=(effective) number of FFT grid points (for this processor) | nkpt=number of k points. | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetry elements in space group | typat= array of types of the natoms electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation etotal=total energy (Ha) - only needed for tddft fock <type(fock_type)>= quantities to calculate Fock exact exchange gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations (3x3 tensor) and grads wrt atomic coordinates (3*natom) gsqcut=cutoff on (k+G)^2 (bohr^-2) hdr <type(hdr_type)>=the header of wf, den and pot files indsym(4,nsym,natom)=indirect indexing array for atom labels irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for diel matrix nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise istep=index of the number of steps in the routine scfcv istep_mix=index of the number of steps for the SCF mixing (can be <istep) kg(3,mpw*mkmem)=reduced planewave coordinates. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. kxc(nfftf,nkxc)=exchange-correlation kernel, needed only if nkxc/=0 . lmax_diel=1+max. value of l angular momentum used for dielectric matrix mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix mpi_enreg=informations about MPI parallelization my_natom=number of atoms treated by current processor natom=number of atoms in cell. nattyp(ntypat)= # atoms of each type. nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs) nfftdiel=number of fft grid points for the computation of the diel matrix ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix, see ~abinit/doc/variables/vargs.htm#ngfft nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwarr(nkpt)=number of planewaves in basis at this k point npwdiel=size of the susmat array. ntypat=number of types of atoms in unit cell. optforces=option for the computation of forces (0: no force;1: forces) optres=0: the new value of the density is computed in place of the input value 1: only the density residual is computed ; the input density is kept paw_dmft <type(paw_dmft_type)>= paw+dmft related data paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang)>=paw angular mesh and related data pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise phnonsdiel(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases, for diel matr nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information for the dielectric matrix psps <type(pseudopotential_type)>=variables related to pseudopotentials pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations (see initberry.f) rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive vectors symrec(3,3,nsym)=symmetry operations in reciprocal space ucvol=unit cell volume in bohr**3. usecprj=1 if cprj datastructure is stored in memory wffnew,unit numbers for wf disk files. vtrial(nfftf,nspden)=INPUT potential Vtrial(r). vxctau=(only for meta-GGA): derivative of XC energy density with respect to kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also the gradient of vxctau (gvxctau) in vxctau(:,:,2:4) xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point for the dielectric matrix
OUTPUT
compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid dphase(3) : dphase(idir) = accumulated change in the string-averaged Zak phase along the idir-th direction caused by the update of all the occupied Bloch states at all the k-points (only if finite electric field) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) resid(mband*nkpt*nsppol)=residuals for each band over all k points. residm=maximum value from resid array (except for nbdbuf highest bands) susmat(2,npwdiel*afford,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space === if optforces>0 === grnl(3*natom)=stores grads of nonlocal energy wrt length scales ==== if optres==1 nres2=square of the norm of the residual nvresid(nfftf,nspden)=density residual ==== if psps%usepaw==1 cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors: cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector. nhat(nfftf,nspden*psps%usepaw)=compensation charge density on rectangular grid in real space
SIDE EFFECTS
cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions. At output contains updated wavefunctions coefficients; if nkpt>1, these are kept in a disk file. energies <type(energies_type)>=storage for energies computed here : | e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree) | e_kinetic=kinetic energy part of total energy | e_nonlocalpsp=nonlocal pseudopotential part of total energy | e_fermie=fermi energy (Hartree) occ(mband*nkpt*nsppol)=occupation number for each band for each k. (input if insulator - occopt<3 - ; output if metallic) pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data rhog(2,nfftf)=Fourier transform of total electron density rhor(nfftf,nspden)=total electron density (el/bohr**3) taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density wvl <type(wvl_data)>=wavelets structures in case of wavelets basis. ==== if (usepaw==1) ==== cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors: cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
PARENTS
scfcv
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
NOTES
Be careful to the meaning of nfft (size of FFT grids): - In case of norm-conserving calculations the FFT grid is the usual FFT grid. - In case of PAW calculations: Two FFT grids are used; one with nfft points (coarse grid) for the computation of wave functions ; one with nfftf points (fine grid) for the computation of total density. The total electronic density (rhor,rhog) is divided into two terms: - The density related to WFs =Sum[Psi**2] - The compensation density (nhat) - only in PAW The parallelisation needed for the electric field should be made an independent subroutine, so that this routine could be put back in the 95_drive directory.
SOURCE
180 #if defined HAVE_CONFIG_H 181 #include "config.h" 182 #endif 183 184 #include "abi_common.h" 185 186 subroutine vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,& 187 & dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,& 188 & eigen,electronpositron,energies,etotal,gbound_diel,& 189 & gmet,gprimd,grnl,gsqcut,hdr,indsym,irrzon,irrzondiel,& 190 & istep,istep_mix,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,& 191 & my_natom,natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,& 192 & npwarr,npwdiel,nres2,ntypat,nvresid,occ,optforces,& 193 & optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,& 194 & phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,& 195 & pwind,pwind_alloc,pwnsfac,resid,residm,rhog,rhor,& 196 & rmet,rprimd,susmat,symrec,taug,taur,& 197 & ucvol,usecprj,wffnew,vtrial,vxctau,wvl,xred,ylm,ylmgr,ylmdiel) 198 199 use defs_basis 200 use defs_datatypes 201 use defs_abitypes 202 use defs_wvltypes 203 use m_profiling_abi 204 use m_xmpi 205 use m_ab7_mixing 206 use m_errors 207 use m_wffile 208 use m_efield 209 use m_cgtools 210 use m_gemm_nonlop 211 212 use m_geometry, only : xred2xcart 213 use m_occ, only : newocc 214 use m_dtset, only : testsusmat 215 use m_pawang, only : pawang_type 216 use m_pawtab, only : pawtab_type 217 use m_paw_ij, only : paw_ij_type 218 use m_pawfgrtab, only : pawfgrtab_type 219 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_get_nspden 220 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim 221 use m_pawfgr, only : pawfgr_type 222 use m_energies, only : energies_type 223 use m_hamiltonian, only : init_hamiltonian,destroy_hamiltonian,& 224 & load_spin_hamiltonian,load_k_hamiltonian,gs_hamiltonian_type 225 use m_bandfft_kpt, only : bandfft_kpt,bandfft_kpt_type,bandfft_kpt_set_ikpt,& 226 & bandfft_kpt_savetabs,bandfft_kpt_restoretabs 227 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 228 use m_paw_dmft, only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft,saveocc_dmft 229 use m_crystal, only : crystal_init, crystal_free, crystal_t 230 use m_oper, only : oper_type,init_oper,destroy_oper 231 use m_io_tools, only : flush_unit 232 use m_abi2big, only : wvl_occ_abi2big,wvl_rho_abi2big,wvl_occopt_abi2big 233 use m_fock, only : fock_type,fock_ACE_type,fock_updateikpt,fock_calc_ene 234 use m_invovl, only : make_invovl 235 use m_tddft, only : tddft 236 use m_kg, only : mkkin 237 use m_suscep_stat, only : suscep_stat 238 239 #if defined HAVE_BIGDFT 240 use BigDFT_API, only : last_orthon,evaltoocc,write_energies 241 #endif 242 243 !This section has been created automatically by the script Abilint (TD). 244 !Do not modify the following lines by hand. 245 #undef ABI_FUNC 246 #define ABI_FUNC 'vtorho' 247 use interfaces_14_hidewrite 248 use interfaces_18_timing 249 use interfaces_32_util 250 use interfaces_53_ffts 251 use interfaces_56_recipspace 252 use interfaces_62_wvl_wfs 253 use interfaces_65_paw 254 use interfaces_66_nonlocal 255 use interfaces_66_wfs 256 use interfaces_67_common 257 use interfaces_68_dmft 258 use interfaces_79_seqpar_mpi, except_this_one => vtorho 259 !End of the abilint section 260 261 implicit none 262 263 !Arguments ------------------------------- 264 integer, intent(in) :: afford,dbl_nnsclo,dielop,dielstrt,istep,istep_mix,lmax_diel,mcg,mcprj,mgfftdiel 265 integer, intent(in) :: my_natom,natom,nfftf,nfftdiel,nkxc,npwdiel 266 integer, intent(in) :: ntypat,optforces,optres,pwind_alloc,usecprj 267 real(dp), intent(in) :: cpus,etotal,gsqcut,ucvol 268 real(dp), intent(out) :: compch_fft,nres2,residm 269 type(MPI_type), intent(inout) :: mpi_enreg 270 type(datafiles_type), intent(in) :: dtfil 271 type(dataset_type), intent(inout) :: dtset 272 type(efield_type), intent(inout) :: dtefield 273 type(electronpositron_type),pointer :: electronpositron 274 type(energies_type), intent(inout) :: energies 275 type(hdr_type), intent(inout) :: hdr 276 type(paw_dmft_type), intent(inout) :: paw_dmft 277 type(pawang_type), intent(in) :: pawang 278 type(pawfgr_type), intent(in) :: pawfgr 279 type(pseudopotential_type), intent(in) :: psps 280 type(fock_type),pointer, intent(inout) :: fock 281 type(wffile_type), intent(inout) :: wffnew 282 type(wvl_data), intent(inout) :: wvl 283 integer, intent(in) :: atindx(natom),atindx1(natom),gbound_diel(2*mgfftdiel+8,2) 284 integer, intent(in) :: indsym(4,dtset%nsym,natom) 285 integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 286 integer, intent(in) :: irrzondiel(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 287 integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),kg_diel(3,npwdiel),nattyp(ntypat),ngfftdiel(18),npwarr(dtset%nkpt) 288 integer, intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym) 289 real(dp), intent(in) :: dmatpawu(:,:,:,:),gmet(3,3),gprimd(3,3),ph1d(2,3*(2*dtset%mgfft+1)*natom) 290 real(dp), intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*psps%usepaw) 291 real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 292 real(dp), intent(in) :: phnonsdiel(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 293 real(dp), intent(in) :: pwnsfac(2,pwind_alloc),rmet(3,3),rprimd(3,3) 294 real(dp), intent(inout) :: vtrial(nfftf,dtset%nspden) 295 real(dp), intent(inout) :: xred(3,natom) 296 real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 297 real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 298 real(dp), intent(in) :: ylmdiel(npwdiel,lmax_diel**2) 299 real(dp), intent(out) :: dphase(3),grnl(3*natom) 300 real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 301 real(dp), intent(out) :: nhat(nfftf,dtset%nspden*psps%usepaw) 302 real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 303 real(dp), intent(out) :: nvresid(nfftf,dtset%nspden) 304 real(dp), intent(out) :: susmat(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden) 305 real(dp), intent(inout) :: cg(2,mcg) 306 real(dp), intent(inout) :: kxc(nfftf,nkxc),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 307 real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden) 308 real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden) 309 real(dp), intent(inout),optional :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 310 type(pawcprj_type),allocatable,intent(inout) :: cprj(:,:) 311 type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw) 312 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 313 type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw) 314 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 315 316 !Local variables------------------------------- 317 !scalars 318 integer,parameter :: level=111,tim_mkrho=2 319 integer,save :: nwarning=0 320 integer :: bdtot_index,counter,cplex,dimffnl,enunit,iband,iband1,ibdkpt 321 integer :: ibg,icg,ider,idir,ierr,ifft,ifor,ifor1,ii,ikg,ikpt 322 integer :: ikpt_loc,ikpt1,my_ikpt,ikxc,ilm,imagn,index1,iorder_cprj,ipert,iplex 323 integer :: iscf,ispden,isppol,istwf_k,mband_cprj,mbdkpsp,mb2dkpsp 324 integer :: mcgq,mcprj_local,mcprj_tmp,me_distrb,mkgq,mpi_comm_sphgrid 325 integer :: mwarning,my_nspinor,n1,n2,n3,n4,n5,n6,nband_eff 326 integer :: nband_k,nband_cprj_k,nbuf,neglect_pawhat,nfftot,nkpg,nkpt1,nnn,nnsclo_now 327 integer :: nproc_distrb,npw_k,nspden_rhoij,option,prtvol 328 integer :: spaceComm_distrb,usecprj_local,usefock_ACE,usetimerev 329 logical :: berryflag,computesusmat,fixed_occ 330 logical :: locc_test,paral_atom,remove_inv,usefock,with_vxctau 331 logical :: do_last_ortho,wvlbigdft=.false. 332 real(dp) :: dmft_ldaocc 333 real(dp) :: edmft,ebandlda,ebanddmft,ebandldatot,ekindmft,ekindmft2,ekinlda 334 real(dp) :: min_occ,vxcavg_dum,strsxc(6) 335 character(len=500) :: message 336 type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null() 337 type(gs_hamiltonian_type) :: gs_hamk 338 !arrays 339 integer,allocatable :: kg_k(:,:) 340 real(dp) :: dielar(7),dphase_k(3),kpoint(3),qpt(3),rhodum(1),tsec(2),ylmgr_dum(0,0,0) 341 real(dp),allocatable :: EigMin(:,:),buffer1(:),buffer2(:),cgq(:,:) 342 real(dp),allocatable :: cgrkxc(:,:),cgrvtrial(:,:),doccde(:) 343 real(dp),allocatable :: dphasek(:,:),eig_k(:),ek_k(:),ek_k_nd(:,:,:),eknk(:),eknk_nd(:,:,:,:,:) 344 real(dp),allocatable :: enl_k(:),enlnk(:),focknk(:),fockfornk(:,:,:),ffnl(:,:,:,:),grnl_k(:,:), xcart(:,:) 345 real(dp),allocatable :: grnlnk(:,:),kinpw(:),kpg_k(:,:),occ_k(:),ph3d(:,:,:) 346 real(dp),allocatable :: pwnsfacq(:,:),resid_k(:),rhoaug(:,:,:,:),rhowfg(:,:),rhowfr(:,:) 347 real(dp),allocatable :: vlocal(:,:,:,:),vlocal_tmp(:,:,:),vxctaulocal(:,:,:,:,:),ylm_k(:,:),zshift(:) 348 complex(dpc),allocatable :: nucdipmom_k(:) 349 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 350 type(oper_type) :: lda_occup 351 type(pawrhoij_type),pointer :: pawrhoij_unsym(:) 352 type(crystal_t) :: cryst_struc 353 integer :: idum1(0),idum3(0,0,0) 354 real(dp) :: rdum2(0,0),rdum4(0,0,0,0) 355 356 !Variables for BigDFT 357 #if defined HAVE_BIGDFT 358 integer :: occopt_bigdft 359 #endif 360 361 ! ********************************************************************* 362 363 DBG_ENTER("COLL") 364 365 !Keep track of total time spent in vtorho 366 call timab(980,1,tsec) 367 call timab(981,1,tsec) 368 369 !Structured debugging if prtvol==-level 370 prtvol=dtset%prtvol 371 372 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 373 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 374 375 !Several inits 376 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 377 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 378 usecprj_local=0;if (psps%usepaw==1) usecprj_local=1 379 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 380 paral_atom=(my_natom/=natom) 381 compch_fft=-1.d5 382 383 !Check that usekden is not 0 if want to use vxctau 384 with_vxctau = (dtset%usekden/=0) 385 386 !Check that fock is present if want to use fock option 387 usefock = (dtset%usefock==1 .and. associated(fock)) 388 usefock_ACE=0 389 if (usefock) usefock_ACE=fock%fock_common%use_ACE 390 391 392 !Init MPI 393 spaceComm_distrb=mpi_enreg%comm_cell 394 if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt 395 if (mpi_enreg%paral_hf ==1) spaceComm_distrb=mpi_enreg%comm_kpt 396 nproc_distrb=xmpi_comm_size(spaceComm_distrb) 397 me_distrb=xmpi_comm_rank(spaceComm_distrb) 398 mpi_comm_sphgrid=mpi_enreg%comm_fft 399 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 400 if (mpi_enreg%me_img/=0) nwarning=nwarning+1 401 402 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 403 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 404 MSG_BUG('wrong values for nfft, nfftf!') 405 end if 406 407 !Test optforces (to prevent memory overflow) 408 if (optforces/=0.and.optforces/=1) then 409 write(message,'(a,i0)')' wrong value for optforces = ',optforces 410 MSG_BUG(message) 411 end if 412 413 iscf=dtset%iscf 414 fixed_occ=(dtset%occopt<3.or.electronpositron_calctype(electronpositron)==1) 415 if(.not. wvlbigdft) then 416 energies%e_eigenvalues = zero 417 energies%e_kinetic = zero 418 energies%e_nonlocalpsp = zero 419 if (usefock) then 420 energies%e_fock=zero 421 energies%e_fockdc=zero 422 end if 423 grnl(:)=zero 424 resid(:) = zero ! JWZ 13 May 2010. resid and eigen need to be fully zeroed each time before use 425 eigen(:) = zero 426 bdtot_index=0 427 ibg=0;icg=0 428 mbdkpsp=dtset%mband*dtset%nkpt*dtset%nsppol 429 if(paw_dmft%use_dmft==1) mb2dkpsp=2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol 430 end if 431 432 if(dtset%usewvl==0) then 433 ABI_ALLOCATE(eknk,(mbdkpsp)) 434 ABI_ALLOCATE(enlnk,(mbdkpsp)) 435 ABI_ALLOCATE(eknk_nd,(dtset%nsppol,dtset%nkpt,2,dtset%mband,dtset%mband*paw_dmft%use_dmft)) 436 ABI_ALLOCATE(EigMin,(2,dtset%mband)) 437 ABI_ALLOCATE(grnlnk,(3*natom,mbdkpsp*optforces)) 438 if (usefock) then 439 ABI_ALLOCATE(focknk,(mbdkpsp)) 440 focknk=zero 441 if (optforces>0)then 442 ABI_ALLOCATE(fockfornk,(3,natom,mbdkpsp)) 443 fockfornk=zero 444 end if 445 end if 446 eknk(:)=zero;enlnk(:)=zero 447 if (optforces>0) grnlnk(:,:)=zero 448 if(paw_dmft%use_dmft==1) eknk_nd=zero 449 end if !usewvl==0 450 451 !Initialize rhor if needed; store old rhor 452 if(iscf>=0 .or. iscf==-3) then 453 if (optres==1) nvresid=rhor 454 ! NC and plane waves 455 if (psps%usepaw==0 .and. dtset%usewvl==0) then 456 rhor=zero 457 ! PAW 458 elseif(psps%usepaw==1) then 459 ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden)) 460 ABI_ALLOCATE(rhowfg,(2,dtset%nfft)) 461 rhowfr(:,:)=zero 462 end if 463 end if 464 465 !Set max number of non-self-consistent loops nnsclo_now for use in vtowfk 466 if(iscf<0)then 467 ! ===== Non self-consistent ===== 468 nnsclo_now=dtset%nstep 469 else 470 ! ===== Self-consistent ===== 471 if(dtset%nnsclo>0) then 472 ! ===== Self-consistent + imposed ===== 473 nnsclo_now=dtset%nnsclo 474 else 475 ! ===== Self-consistent + default ===== 476 nnsclo_now=1 477 if (dtset%usewvl==0) then 478 ! ----- Plane waves ----- 479 if (istep<=2.and.iscf/=0) nnsclo_now=2 480 else 481 ! ----- Wavelets ----- 482 if (iscf==0) then 483 nnsclo_now=0 484 else if (istep<=2) then 485 nnsclo_now=3 486 else if (istep<=4) then 487 nnsclo_now=2 488 end if 489 end if 490 end if 491 ! ===== Double is required ===== 492 if(dbl_nnsclo==1)then 493 nnsclo_now=nnsclo_now*2 494 end if 495 end if 496 if(dtset%wfoptalg==2)nnsclo_now=40 ! UNDER DEVELOPMENT 497 498 write(message, '(a,i0,a,3(i0,1x))' ) ' vtorho : nnsclo_now=',nnsclo_now,& 499 & ', note that nnsclo,dbl_nnsclo,istep=',dtset%nnsclo,dbl_nnsclo,istep 500 call wrtout(std_out,message,'COLL') 501 502 !==== Initialize most of the Hamiltonian ==== 503 !Allocate all arrays and initialize quantities that do not depend on k and spin. 504 call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,& 505 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 506 & paw_ij=paw_ij,ph1d=ph1d,usecprj=usecprj_local,electronpositron=electronpositron,fock=fock,& 507 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 508 & nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda) 509 510 !Initializations for PAW (projected wave functions) 511 mcprj_local=0 ; mband_cprj=0 512 if (psps%usepaw==1) then 513 mband_cprj=dtset%mband 514 if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band 515 iorder_cprj=0 ; mcprj_local=mcprj 516 if (usecprj==0) then 517 mcprj_local=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol 518 !This is a check but should always be true since scfcv allocated cprj anyway 519 if (allocated(cprj)) then 520 !Was allocated in scfcv so we just destroy and reconstruct it as desired 521 call pawcprj_free(cprj) 522 ABI_DATATYPE_DEALLOCATE(cprj) 523 end if 524 ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj_local)) 525 call pawcprj_alloc(cprj,0,gs_hamk%dimcprj) 526 end if 527 end if 528 529 call timab(981,2,tsec) 530 531 !=================================================================== 532 ! WAVELETS - Branching with a separate VTORHO procedure 533 !=================================================================== 534 535 if (dtset%usewvl == 1) then 536 #ifndef HAVE_BIGDFT 537 BIGDFT_NOTENABLED_ERROR() 538 #else 539 540 ! do_last_ortho in case of diagonalization scheme 541 if ( wvlbigdft) do_last_ortho=(dtset%iscf/=0) 542 if (.not.wvlbigdft) do_last_ortho=(.true.) 543 544 ABI_ALLOCATE(xcart,(3, dtset%natom)) 545 call xred2xcart(dtset%natom, rprimd, xcart, xred) 546 547 if(wvlbigdft) then 548 ! NSCF loop for wvlbigdt: 549 call wvl_nscf_loop_bigdft() 550 else 551 ! NSCF loop for WVL: (not wvlbigdft) 552 call wvl_nscf_loop() 553 end if 554 555 ! Eventually orthogonalize WFs now 556 if (do_last_ortho) then 557 call write_energies(ii,0,wvl%e%energs,0.d0,0.d0,"final") 558 call last_orthon(me_distrb, nproc_distrb, ii, wvl%wfs%ks, wvl%e%energs%evsum, .true.) 559 if(wvlbigdft) energies%e_xcdc = wvl%e%energs%evxc 560 ! If occupation numbers are not changed... 561 if (fixed_occ .or. (iscf<0 .and. iscf/=-3)) then 562 call wvl_comm_eigen() 563 end if 564 ! ... or update occupations: 565 if( ( .not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then 566 if(wvlbigdft) then 567 call wvl_occ_bigdft() 568 else 569 ! Communicate eigenvalues: 570 call wvl_comm_eigen() 571 ! Update occ and Fermi level 572 call wvl_occ() 573 end if 574 end if 575 ! This might accelerate convergence: 576 wvl%wfs%ks%diis%energy_min=one 577 wvl%wfs%ks%diis%alpha=two 578 end if !do_last_ortho 579 580 ! Compute eigenvalues energy 581 if(.not. wvlbigdft .and. nnsclo_now>0) then 582 call e_eigen(eigen,energies%e_eigenvalues,dtset%mband,dtset%nband,dtset%nkpt,& 583 & dtset%nsppol,occ,dtset%wtk) 584 else 585 energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp & 586 & + energies%e_xcdc + two*energies%e_hartree +energies%e_nonlocalpsp 587 end if 588 589 if (optforces == 1) then ! not compatible with iscf=0 and wvlbigdftcomp=1 + PAW 590 call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart) 591 end if 592 593 ! For iscf<0 we do not update the density 594 if (dtset%iscf>=0) then !(dtset%iscf>=0 .and. .not. wvlbigdft ) then 595 call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den) 596 end if 597 ABI_DEALLOCATE(xcart) 598 599 ! Note in WVL+NC: the rest will be skipped. 600 ! For PAW: we will compute Rho_ij at the end. 601 !if(wvlbigdft) return 602 #endif 603 else 604 605 !=================================================================== 606 ! PLANE WAVES - Standard VTORHO procedure 607 !=================================================================== 608 609 ! Electric fields: set flag to turn on various behaviors 610 berryflag = (dtset%berryopt == 4 .or. dtset%berryopt == 14 .or. & 611 & dtset%berryopt == 6 .or. dtset%berryopt == 16 .or. & 612 & dtset%berryopt == 7 .or. dtset%berryopt == 17) 613 614 ! Electric field: allocate dphasek 615 nkpt1 = dtset%nkpt 616 if ( berryflag ) then 617 ABI_ALLOCATE(dphasek,(3,dtset%nkpt*dtset%nsppol)) 618 dphasek(:,:) = zero 619 nkpt1 = dtefield%mkmem_max 620 end if 621 622 ABI_ALLOCATE(rhoaug,(n4,n5,n6,gs_hamk%nvloc)) 623 ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamk%nvloc)) 624 if(with_vxctau) then 625 ABI_ALLOCATE(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4)) 626 end if 627 628 ! LOOP OVER SPINS 629 do isppol=1,dtset%nsppol 630 call timab(982,1,tsec) 631 632 ikpt_loc = 0 633 ikg=0 634 635 ! Set up local potential vlocal with proper dimensioning, from vtrial 636 ! Also take into account the spin. 637 if(dtset%nspden/=4)then 638 if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then 639 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal,2) 640 if(with_vxctau) then 641 do ii=1,4 642 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,& 643 & vxctau(:,:,ii),vxctaulocal(:,:,:,:,ii),2) 644 end do 645 end if 646 else 647 ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden)) 648 call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 649 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal,2) 650 ABI_DEALLOCATE(cgrvtrial) 651 end if 652 else 653 ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6)) 654 if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then 655 do ispden=1,dtset%nspden 656 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal_tmp,2) 657 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 658 end do 659 else 660 ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden)) 661 call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 662 do ispden=1,dtset%nspden 663 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal_tmp,2) 664 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 665 end do 666 ABI_DEALLOCATE(cgrvtrial) 667 end if 668 ABI_DEALLOCATE(vlocal_tmp) 669 end if ! nspden 670 671 rhoaug(:,:,:,:)=zero 672 673 ! Continue to initialize the Hamiltonian 674 call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.) 675 if (with_vxctau) then 676 call load_spin_hamiltonian(gs_hamk,isppol,vxctaulocal=vxctaulocal) 677 end if 678 679 call timab(982,2,tsec) 680 681 ! BIG FAT k POINT LOOP 682 ! MVeithen: I had to modify the structure of this loop in order to implement MPI // of the electric field 683 ! Note that the loop here differs from the similar one in berryphase_new.F90. 684 ! here, ikpt_loc numbers the kpts treated by the current processor. 685 ! in berryphase_new.F90, ikpt_loc ALSO includes info about value of isppol. 686 687 ikpt = 0 688 do while (ikpt_loc < nkpt1) 689 690 call timab(997,1,tsec) 691 692 if ( .not.berryflag ) then 693 ikpt_loc = ikpt_loc + 1 694 ikpt = ikpt_loc 695 my_ikpt = mpi_enreg%my_kpttab(ikpt) 696 else 697 if (ikpt_loc < dtset%mkmem) ikpt = ikpt + 1 698 if ((ikpt > dtset%nkpt).and.(ikpt_loc < dtset%mkmem)) exit 699 my_ikpt=ikpt 700 end if 701 702 dphase_k(:) = zero 703 counter=100*ikpt+isppol 704 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 705 nband_cprj_k=nband_k/mpi_enreg%nproc_band 706 istwf_k=dtset%istwfk(ikpt) 707 npw_k=npwarr(ikpt) 708 709 mcgq=1 ; mkgq=1 710 if (.not.berryflag) then 711 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 712 eigen(1+bdtot_index : nband_k+bdtot_index) = zero 713 resid(1+bdtot_index : nband_k+bdtot_index) = zero 714 bdtot_index=bdtot_index+nband_k 715 cycle 716 end if 717 else 718 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) & 719 & .and.(ikpt_loc <= dtset%mkmem)) then 720 eigen(1+bdtot_index : nband_k+bdtot_index) = zero 721 resid(1+bdtot_index : nband_k+bdtot_index) = zero 722 bdtot_index = bdtot_index + nband_k 723 cycle 724 end if 725 ikpt_loc = ikpt_loc + 1 726 mcgq = dtset%mpw*my_nspinor*nband_k*dtefield%nneigh(ikpt) 727 ikg = dtefield%kgindex(ikpt) 728 mkgq = 6*dtset%mpw 729 end if ! berryflag 730 731 call timab(997,2,tsec) 732 733 ! In case of MPI // of a finite field calculation 734 ! build the cgq array that stores the wavefunctions for the 735 ! neighbours of ikpt, and the pwnsfacq array that stores the 736 ! corresponding phase factors (in case of tnons) 737 ABI_ALLOCATE(cgq,(2,mcgq)) 738 ABI_ALLOCATE(pwnsfacq,(2,mkgq)) 739 if ( berryflag ) then 740 call cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,& 741 & me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,& 742 & npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb) 743 if (ikpt_loc > dtset%mkmem) then 744 ABI_DEALLOCATE(cgq) 745 ABI_DEALLOCATE(pwnsfacq) 746 cycle 747 end if 748 end if !berryopt 749 750 call timab(984,1,tsec) 751 752 if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt) 753 call bandfft_kpt_set_ikpt(ikpt,mpi_enreg) 754 ! my_ikpt = ikpt 755 ! if (mpi_enreg%paral_kgb==1) then 756 ! my_ikpt = mpi_enreg%my_kpttab(ikpt) 757 ! my_bandfft_kpt => bandfft_kpt(my_ikpt) 758 ! call bandfft_kpt_set_ikpt(ikpt,mpi_enreg) 759 ! end if 760 761 ABI_ALLOCATE(eig_k,(nband_k)) 762 ABI_ALLOCATE(ek_k,(nband_k)) 763 ABI_ALLOCATE(enl_k,(nband_k)) 764 ABI_ALLOCATE(ek_k_nd,(2,nband_k,nband_k*paw_dmft%use_dmft)) 765 ABI_ALLOCATE(occ_k,(nband_k)) 766 ABI_ALLOCATE(resid_k,(nband_k)) 767 ABI_ALLOCATE(zshift,(nband_k)) 768 ABI_ALLOCATE(grnl_k,(3*natom,nband_k*optforces)) 769 770 eig_k(:)=zero 771 ek_k(:)=zero 772 enl_k(:)=zero 773 if(paw_dmft%use_dmft==1) ek_k_nd(:,:,:)=zero 774 if (optforces>0) grnl_k(:,:)=zero 775 kpoint(:)=dtset%kptns(:,ikpt) 776 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 777 resid_k(:)=zero 778 zshift(:)=dtset%eshift 779 780 ABI_ALLOCATE(kg_k,(3,npw_k)) 781 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 782 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 783 if (psps%useylm==1) then 784 do ilm=1,psps%mpsang*psps%mpsang 785 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 786 end do 787 end if 788 789 ! Set up remaining of the Hamiltonian 790 791 ! Compute (1/2) (2 Pi)**2 (k+G)**2: 792 ABI_ALLOCATE(kinpw,(npw_k)) 793 ! call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k) 794 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0) 795 796 ! Compute (k+G) vectors (only if useylm=1) 797 nkpg=3*optforces*dtset%nloalg(3) 798 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 799 if ((mpi_enreg%paral_kgb/=1.or.istep<=1).and.nkpg>0) then 800 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 801 end if 802 803 ! Compute nonlocal form factors ffnl at all (k+G): 804 ider=0;idir=0;dimffnl=1 805 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 806 if (mpi_enreg%paral_kgb/=1.or.istep<=1) then 807 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,& 808 & gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 809 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 810 & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,& 811 & psps%usepaw,psps%useylm,ylm_k,ylmgr) 812 end if 813 814 ! compute and load nuclear dipole Hamiltonian at current k point 815 if(any(abs(gs_hamk%nucdipmom)>0.0)) then 816 if(allocated(nucdipmom_k)) then 817 ABI_DEALLOCATE(nucdipmom_k) 818 end if 819 ABI_ALLOCATE(nucdipmom_k,(npw_k*(npw_k+1)/2)) 820 call mknucdipmom_k(gmet,kg_k,kpoint,natom,gs_hamk%nucdipmom,nucdipmom_k,npw_k,rprimd,ucvol,xred) 821 call load_k_hamiltonian(gs_hamk,nucdipmom_k=nucdipmom_k) 822 end if 823 824 825 ! Load k-dependent part in the Hamiltonian datastructure 826 ! - Compute 3D phase factors 827 ! - Prepare various tabs in case of band-FFT parallelism 828 ! - Load k-dependent quantities in the Hamiltonian 829 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk)) 830 831 if(usefock_ACE/=0) then 832 call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,& 833 & kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,fockACE_k=fock%fockACE(ikpt,isppol),ph3d_k=ph3d,& 834 & compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),& 835 & compute_gbound=(mpi_enreg%paral_kgb/=1)) 836 else 837 call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,& 838 & kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,& 839 & compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),& 840 & compute_gbound=(mpi_enreg%paral_kgb/=1)) 841 end if 842 843 ! Load band-FFT tabs (transposed k-dependent arrays) 844 if (mpi_enreg%paral_kgb==1) then 845 if (istep<=1) then 846 call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg) 847 end if 848 call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, & 849 & gbound_k =my_bandfft_kpt%gbound, & 850 & kinpw_k =my_bandfft_kpt%kinpw_gather, & 851 & kg_k =my_bandfft_kpt%kg_k_gather, & 852 & kpg_k =my_bandfft_kpt%kpg_k_gather, & 853 ffnl_k =my_bandfft_kpt%ffnl_gather, & 854 ph3d_k =my_bandfft_kpt%ph3d_gather) 855 end if 856 857 ! Build inverse of overlap matrix for chebfi 858 if(psps%usepaw == 1 .and. dtset%wfoptalg == 1 .and. istep <= 1) then 859 call make_invovl(gs_hamk, dimffnl, ffnl, ph3d, mpi_enreg) 860 end if 861 862 ! Setup gemm_nonlop 863 if (gemm_nonlop_use_gemm) then 864 !set the global variable indicating to gemm_nonlop where to get its data from 865 gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt 866 if (istep <= 1) then 867 !Init the arrays 868 call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax, & 869 & gs_hamk%ntypat, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, gs_hamk%ffnl_k,& 870 & gs_hamk%ph3d_k) 871 end if 872 end if 873 874 #if defined HAVE_GPU_CUDA 875 if (dtset%use_gpu_cuda==1) then 876 call gpu_update_ffnl_ph3d(ph3d,size(ph3d),ffnl,size(ffnl)) 877 end if 878 #endif 879 880 call timab(984,2,tsec) 881 882 ! Update the value of ikpt,isppol in fock_exchange and allocate the memory space to perform HF calculation. 883 if (usefock) then 884 call fock_updateikpt(fock%fock_common,ikpt,isppol) 885 end if 886 if ((psps%usepaw==1).and.(usefock)) then 887 if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then 888 fock%fock_common%forces_ikpt=zero 889 end if 890 end if 891 892 ! Compute the eigenvalues, wavefunction, residuals, 893 ! contributions to kinetic energy, nonlocal energy, forces, 894 ! and update of rhor to this k-point and this spin polarization. 895 call vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,& 896 & dtset,eig_k,ek_k,ek_k_nd,enl_k,fixed_occ,grnl_k,gs_hamk,& 897 & ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,mband_cprj,mcg,mcgq,mcprj_local,mkgq,& 898 & mpi_enreg,dtset%mpw,natom,nband_k,dtset%nkpt,nnsclo_now,npw_k,npwarr,& 899 & occ_k,optforces,prtvol,pwind,pwind_alloc,pwnsfac,pwnsfacq,resid_k,& 900 & rhoaug,paw_dmft,dtset%wtk(ikpt),zshift) 901 902 call timab(985,1,tsec) 903 904 #if defined HAVE_GPU_CUDA 905 if(dtset%use_gpu_cuda==1) then 906 call gpu_finalize_ffnl_ph3d() 907 end if 908 #endif 909 ABI_DEALLOCATE(ffnl) 910 ABI_DEALLOCATE(kinpw) 911 ABI_DEALLOCATE(kg_k) 912 ABI_DEALLOCATE(kpg_k) 913 if(allocated(nucdipmom_k)) then 914 ABI_DEALLOCATE(nucdipmom_k) 915 end if 916 ABI_DEALLOCATE(ylm_k) 917 ABI_DEALLOCATE(ph3d) 918 ABI_DEALLOCATE(cgq) 919 ABI_DEALLOCATE(pwnsfacq) 920 921 ! electric field 922 if (berryflag) then 923 dphasek(:,ikpt + (isppol - 1)*dtset%nkpt) = dphase_k(:) 924 925 ! The overlap matrices for all first neighbours of ikpt are no more up to date 926 do idir = 1, 3 927 do ifor = 1, 2 928 ikpt1 = dtefield%ikpt_dk(dtefield%i2fbz(ikpt),ifor,idir) 929 ikpt1 = dtefield%indkk_f2ibz(ikpt1,1) 930 ifor1 = -1*ifor + 3 ! ifor = 1 -> ifor1 = 2 & ifor = 2 -> ifor1 = 1 931 dtefield%sflag(:,ikpt1+(isppol-1)*dtset%nkpt,ifor1,idir) = 0 932 end do 933 end do 934 end if ! berryflag 935 936 ! Save eigenvalues (hartree), residuals (hartree**2) 937 eigen(1+bdtot_index : nband_k+bdtot_index) = eig_k(:) 938 eknk (1+bdtot_index : nband_k+bdtot_index) = ek_k (:) 939 if(usefock) then 940 focknk (1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%eigen_ikpt (:) 941 if (optforces>0) fockfornk(:,:,1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%forces_ikpt(:,:,:) 942 end if 943 if(paw_dmft%use_dmft==1) eknk_nd(isppol,ikpt,:,:,:) = ek_k_nd(:,:,:) 944 resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:) 945 if (optforces>0) grnlnk(:,1+bdtot_index : nband_k+bdtot_index) = grnl_k(:,:) 946 enlnk(1+bdtot_index : nband_k+bdtot_index) = enl_k(:) 947 948 if(iscf>0 .or. iscf==-3)then 949 ! Accumulate sum over k points for band, nonlocal and kinetic energies, 950 ! also accumulate gradients of Enonlocal: 951 do iband=1,nband_k 952 if (abs(occ_k(iband))>tol8) then 953 energies%e_kinetic = energies%e_kinetic + dtset%wtk(ikpt)*occ_k(iband)*ek_k(iband) 954 energies%e_eigenvalues = energies%e_eigenvalues + dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband) 955 energies%e_nonlocalpsp = energies%e_nonlocalpsp + dtset%wtk(ikpt)*occ_k(iband)*enl_k(iband) 956 if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ_k(iband)*grnl_k(:,iband) 957 if (usefock) energies%e_fock=energies%e_fock + half*fock%fock_common%eigen_ikpt(iband)*occ_k(iband)*dtset%wtk(ikpt) 958 end if 959 end do 960 961 ! Calculate Fock contribution to the total energy if required 962 if ((psps%usepaw==1).and.(usefock)) then 963 if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then 964 call fock_calc_ene(dtset,fock%fock_common,energies%e_exactX,ikpt,nband_k,occ_k) 965 end if 966 end if 967 end if 968 969 call timab(985,2,tsec) 970 971 ABI_DEALLOCATE(eig_k) 972 ABI_DEALLOCATE(ek_k) 973 ABI_DEALLOCATE(ek_k_nd) 974 ABI_DEALLOCATE(grnl_k) 975 ABI_DEALLOCATE(occ_k) 976 ABI_DEALLOCATE(resid_k) 977 ABI_DEALLOCATE(zshift) 978 ABI_DEALLOCATE(enl_k) 979 980 ! Keep track of total number of bands (all k points so far, even for k points not treated by me) 981 bdtot_index=bdtot_index+nband_k 982 983 ! Also shift array memory if dtset%mkmem/=0 984 if (dtset%mkmem/=0) then 985 ibg=ibg+my_nspinor*nband_cprj_k 986 icg=icg+npw_k*my_nspinor*nband_k 987 ikg=ikg+npw_k 988 end if 989 990 end do ! End big k point loop 991 992 call timab(986,1,tsec) 993 994 if (fixed_occ .and. mpi_enreg%paral_kgb==1) then 995 call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) !Sum the contributions over bands/FFT/spinors 996 end if 997 998 ! Transfer density on augmented fft grid to normal fft grid in real space 999 ! Also take into account the spin. 1000 if(iscf>0.or.iscf==-3)then 1001 if (psps%usepaw==0) then 1002 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,1),1) 1003 if(dtset%nspden==4)then 1004 do imagn=2,4 1005 call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,imagn),1) 1006 end do 1007 end if 1008 else 1009 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,1),1) 1010 if(dtset%nspden==4)then 1011 do imagn=2,4 1012 call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,imagn),1) 1013 end do 1014 end if 1015 end if 1016 end if 1017 1018 call timab(986,2,tsec) 1019 1020 end do ! End loop over spins 1021 1022 call timab(988,1,tsec) 1023 if (usefock) then 1024 if(fock%fock_common%optfor) then 1025 call xmpi_sum(fock%fock_common%forces,mpi_enreg%comm_kpt,ierr) 1026 end if 1027 end if 1028 ! Electric field: compute string-averaged change in Zak phase 1029 ! along each direction, store it in dphase(idir) 1030 ! ji: it is not convenient to do this anymore. Remove. Set dphase(idir)=0.0_dp. 1031 ! eventually, dphase(idir) will have to go... 1032 if (berryflag) then 1033 dphase(:) = zero 1034 ! In case of MPI // of a finite field calculation, send dphasek to all cpus 1035 call xmpi_sum(dphasek,spaceComm_distrb,ierr) 1036 ABI_DEALLOCATE(dphasek) 1037 end if ! berryflag 1038 ABI_DEALLOCATE(rhoaug) 1039 ABI_DEALLOCATE(vlocal) 1040 if(with_vxctau) then 1041 ABI_DEALLOCATE(vxctaulocal) 1042 end if 1043 1044 call timab(988,2,tsec) 1045 1046 ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1047 doccde(:)=zero 1048 1049 ! Treat now varying occupation numbers, in the self-consistent case 1050 if((.not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then 1051 1052 ! Parallel case 1053 if (mpi_enreg%nproc_kpt>1) then 1054 1055 call timab(989,1,tsec) 1056 1057 ! If needed, exchange the values of eigen,resid,eknk,enlnk,grnlnk 1058 ABI_ALLOCATE(buffer1,((4+3*natom*optforces+dtset%usefock+3*natom*dtset%usefock*optforces)*mbdkpsp)) 1059 if(paw_dmft%use_dmft==1) then 1060 ABI_ALLOCATE(buffer2,(mb2dkpsp*paw_dmft%use_dmft)) 1061 end if 1062 ! Pack eigen,resid,eknk,enlnk,grnlnk in buffer1 1063 buffer1(1 : mbdkpsp)=eigen(:) 1064 buffer1(1+ mbdkpsp:2*mbdkpsp)=resid(:) 1065 buffer1(1+2*mbdkpsp:3*mbdkpsp)=eknk(:) 1066 buffer1(1+3*mbdkpsp:4*mbdkpsp)=enlnk(:) 1067 index1=4*mbdkpsp 1068 if (optforces>0) then 1069 buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(grnlnk,(/(3*natom)*mbdkpsp/) ) 1070 index1=index1+3*natom*mbdkpsp 1071 end if 1072 if (usefock) then 1073 buffer1(1+index1:index1+mbdkpsp)=focknk(:) 1074 if (optforces>0) then 1075 index1=index1+mbdkpsp 1076 buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(fockfornk,(/(3*natom)*mbdkpsp/) ) 1077 end if 1078 end if 1079 if(paw_dmft%use_dmft==1) then 1080 nnn=0 1081 do ikpt=1,dtset%nkpt 1082 do isppol=1,dtset%nsppol 1083 do iband=1,dtset%mband 1084 do iband1=1,dtset%mband 1085 do iplex=1,2 1086 nnn=nnn+1 1087 buffer2(nnn)=eknk_nd(isppol,ikpt,iplex,iband,iband1) 1088 end do 1089 end do 1090 end do 1091 end do 1092 end do 1093 if(nnn.ne.mb2dkpsp) then 1094 write(message,*)' BUG in vtorho2, buffer2',nnn,mb2dkpsp 1095 MSG_BUG(message) 1096 end if 1097 end if 1098 ! Build sum of everything 1099 call timab(48,1,tsec) 1100 call xmpi_sum(buffer1,mpi_enreg%comm_kpt,ierr) 1101 if(mpi_enreg%paral_kgb/=1.and.paw_dmft%use_dmft==1) then 1102 call xmpi_sum(buffer2,mpi_enreg%comm_kpt,ierr) 1103 end if 1104 call timab(48,2,tsec) 1105 1106 ! Unpack eigen,resid,eknk,enlnk,grnlnk in buffer1 1107 eigen(:) =buffer1(1 : mbdkpsp) 1108 resid(:) =buffer1(1+ mbdkpsp:2*mbdkpsp) 1109 eknk(:) =buffer1(1+2*mbdkpsp:3*mbdkpsp) 1110 enlnk(:) =buffer1(1+3*mbdkpsp:4*mbdkpsp) 1111 index1=4*mbdkpsp 1112 if (optforces>0) then 1113 grnlnk(:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3*natom,mbdkpsp/) ) 1114 end if 1115 if (usefock) then 1116 focknk(:)=buffer1(1+index1:index1+mbdkpsp) 1117 if (optforces>0) then 1118 index1=index1+mbdkpsp 1119 fockfornk(:,:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3,natom,mbdkpsp/) ) 1120 end if 1121 end if 1122 if(paw_dmft%use_dmft==1) then 1123 nnn=0 1124 do ikpt=1,dtset%nkpt 1125 do isppol=1,dtset%nsppol 1126 do iband=1,dtset%mband 1127 do iband1=1,dtset%mband 1128 do iplex=1,2 1129 nnn=nnn+1 1130 eknk_nd(isppol,ikpt,iplex,iband,iband1)=buffer2(nnn) 1131 end do 1132 end do 1133 end do 1134 end do 1135 end do 1136 end if 1137 if(allocated(buffer2)) then 1138 ABI_DEALLOCATE(buffer2) 1139 end if 1140 ABI_DEALLOCATE(buffer1) 1141 call timab(989,2,tsec) 1142 1143 end if ! nproc_kpt>1 1144 1145 ! Compute the new occupation numbers from eigen 1146 call timab(990,1,tsec) 1147 call newocc(doccde,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,& 1148 & dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,& 1149 & dtset%nsppol,occ,dtset%occopt,prtvol,dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk) 1150 call timab(990,2,tsec) 1151 1152 ! !========= DMFT call begin ============================================ 1153 dmft_ldaocc=0 1154 if(paw_dmft%use_dmft==1.and.psps%usepaw==1.and.dtset%nbandkss==0) then 1155 call timab(991,1,tsec) 1156 1157 if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(pawtab(:)%upawu)>=tol8.or. & 1158 & sum(pawtab(:)%jpawu)>tol8).and.dtset%dmft_entropy==0) energies%entropy=zero 1159 1160 ! == 0 to a dmft calculation and do not use lda occupations 1161 ! == 1 to a lda calculation with the dmft loop 1162 if(dtset%dmftcheck==-1) dmft_ldaocc=1 1163 1164 ! == initialise occnd 1165 paw_dmft%occnd=zero 1166 1167 if(dmft_ldaocc==0) then 1168 if(dtset%occopt/=3) then 1169 message = ' occopt should be equal to 3 in dmft' 1170 MSG_ERROR(message) 1171 end if 1172 ! == initialise edmft 1173 if(paw_dmft%use_dmft>=1) edmft = zero 1174 1175 ! Compute residm to check the value 1176 ibdkpt=1 1177 residm=zero 1178 do isppol=1,dtset%nsppol 1179 do ikpt=1,dtset%nkpt 1180 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1181 nband_eff=max(1,nband_k-dtset%nbdbuf) 1182 residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1))) 1183 ibdkpt=ibdkpt+nband_k 1184 end do 1185 end do 1186 1187 ! Test residm 1188 if (paw_dmft%use_dmft>0 .and. residm>tol4 .and. dtset%dmftcheck>=0) then 1189 if(dtset%dmft_entropy>0) then 1190 write(message,'(a,e12.3)')& 1191 & ' WARNING: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm 1192 call wrtout(std_out,message,'COLL') 1193 else 1194 write(message,'(a,e12.3)')& 1195 & ' ERROR: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm 1196 call wrtout(std_out,message,'COLL') 1197 write(message,'(a,i0)')' Action: increase nline and nnsclo',dtset%nstep 1198 MSG_ERROR(message) 1199 end if 1200 1201 else if (paw_dmft%use_dmft>0 .and. residm>tol10.and. dtset%dmftcheck>=0) then 1202 write(message,'(3a)')ch10,& 1203 & ' Wavefunctions not converged : DFT+DMFT calculation might not be carried out safely ',ch10 1204 MSG_WARNING(message) 1205 end if 1206 1207 ! == allocate paw_dmft%psichi and paw_dmft%eigen_lda 1208 call init_dmft(dmatpawu,dtset,energies%e_fermie,dtfil%fnameabo_app,dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat) 1209 call print_dmft(paw_dmft,dtset%pawprtvol) 1210 1211 ! == gather crystal structure date into data "cryst_struc" 1212 remove_inv=.false. 1213 if(dtset%nspden==4) remove_inv=.true. 1214 call crystal_init(dtset%amu_orig(:,1),cryst_struc,dtset%spgroup,natom,dtset%npsp,ntypat, & 1215 & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,& 1216 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,& 1217 & dtset%symrel,dtset%tnons,dtset%symafm) 1218 1219 ! == compute psichi 1220 call xmpi_barrier(spaceComm_distrb) 1221 call init_oper(paw_dmft,lda_occup) 1222 call flush_unit(std_out) 1223 call timab(620,1,tsec) 1224 call datafordmft(cryst_struc,cprj,gs_hamk%dimcprj,dtset,eigen,energies%e_fermie,& 1225 & lda_occup,dtset%mband,mband_cprj,dtset%mkmem,mpi_enreg,& 1226 & dtset%nkpt,my_nspinor,dtset%nsppol,occ,& 1227 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj_local,dtfil%unpaw) 1228 call timab(620,2,tsec) 1229 call flush_unit(std_out) 1230 1231 ! == solve dmft loop 1232 call xmpi_barrier(spaceComm_distrb) 1233 1234 call dmft_solve(cryst_struc,istep,lda_occup,paw_dmft,pawang,pawtab,dtset%pawprtvol) 1235 edmft=paw_dmft%edmft 1236 energies%e_paw=energies%e_paw+edmft 1237 energies%e_pawdc=energies%e_pawdc+edmft 1238 call flush_unit(std_out) 1239 ! paw_dmft%occnd(:,:,:,:,:)=0.5_dp 1240 1241 ! call print_dmft(paw_dmft,dtset%pawprtvol) 1242 if(dtset%paral_kgb==1) then 1243 write(message,'(5a)')ch10,& 1244 & ' Parallelization over bands is not yet compatible with self-consistency in DMFT ',ch10,& 1245 & ' Calculation of density does not taken into account non diagonal occupations',ch10 1246 call wrtout(std_out,message,'COLL') 1247 call wrtout(ab_out,message,'COLL') 1248 ! MSG_ERROR(message) 1249 if(dtset%nstep>1) then 1250 write(message,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep 1251 MSG_ERROR(message) 1252 end if 1253 residm=zero 1254 end if 1255 if(dtset%nspinor==2) then 1256 ! call flush_unit(ab_out) 1257 ! write(message,'(3a)')& 1258 ! & ' Self consistent DFT+DMFT with nspinor==2 is not possible yet ',ch10,& 1259 ! & ' Calculation are restricted to nstep =1' 1260 ! ! MSG_ERROR(message) 1261 ! if(dtset%nstep>1) then 1262 ! write(message,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep 1263 ! ! MSG_ERROR(message) 1264 ! endif 1265 end if 1266 1267 if(me_distrb==0) then 1268 call saveocc_dmft(paw_dmft) 1269 end if 1270 call destroy_dmft(paw_dmft) 1271 1272 ! == destroy crystal_t cryst_struc 1273 call crystal_free(cryst_struc) 1274 call destroy_oper(lda_occup) 1275 end if ! dmft_ldaocc 1276 call timab(991,2,tsec) 1277 1278 end if ! usedmft 1279 1280 if(dtset%nbandkss/=0) then 1281 write(message,'(a,i3,2a,i3,4a)') & 1282 & " dtset%nbandkss ==",dtset%nbandkss,ch10,& 1283 & " and dtset%usedmft ==",dtset%usedmft,ch10,& 1284 & " a DFT loop is carried out without DMFT.",ch10,& 1285 & " Only psichi's will be written at convergence of the DFT loop." 1286 call wrtout(std_out,message,'COLL') 1287 end if 1288 ! !========= DMFT call end ============================================ 1289 1290 call timab(992,1,tsec) 1291 1292 ! Compute eeig, ek,enl and grnl from the new occ, and the shared eknk,enlnk,grnlnk 1293 energies%e_eigenvalues = zero 1294 energies%e_kinetic = zero 1295 energies%e_nonlocalpsp = zero 1296 if (usefock) then 1297 energies%e_fock = zero 1298 if (optforces>0) fock%fock_common%forces=zero 1299 end if 1300 if (optforces>0) grnl(:)=zero 1301 if(paw_dmft%use_dmft>=1) then 1302 ebandlda = zero 1303 ebanddmft = zero 1304 ebandldatot = zero 1305 ekindmft = zero 1306 ekindmft2 = zero 1307 ekinlda = zero 1308 end if 1309 1310 ! Compute new energy terms due to non diagonal occupations and DMFT. 1311 bdtot_index=1 1312 do isppol=1,dtset%nsppol 1313 do ikpt=1,dtset%nkpt 1314 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1315 do iband=1,nband_k 1316 1317 locc_test = abs(occ(bdtot_index))>tol8 1318 ! dmft 1319 if(paw_dmft%use_dmft>=1.and.dtset%nbandkss==0) then 1320 if(paw_dmft%band_in(iband)) then 1321 if( paw_dmft%use_dmft == 1 .and. dmft_ldaocc == 1 ) then ! test of the code 1322 paw_dmft%occnd(1,iband,iband,ikpt,isppol)=occ(bdtot_index) 1323 end if 1324 locc_test = abs(paw_dmft%occnd(1,iband,iband,ikpt,isppol))+& 1325 & abs(paw_dmft%occnd(2,iband,iband,ikpt,isppol))>tol8 1326 end if 1327 end if 1328 1329 if (locc_test) then 1330 ! dmft 1331 if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then 1332 ebandldatot=ebandldatot+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index) 1333 if(paw_dmft%band_in(iband)) then 1334 ebandlda=ebandlda+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index) 1335 ekinlda=ekinlda+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index) 1336 occ(bdtot_index)=paw_dmft%occnd(1,iband,iband,ikpt,isppol) 1337 ebanddmft=ebanddmft+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index) 1338 ekindmft=ekindmft+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index) 1339 end if 1340 end if 1341 1342 energies%e_eigenvalues = energies%e_eigenvalues + & 1343 & dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index) 1344 energies%e_kinetic = energies%e_kinetic + & 1345 & dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index) 1346 energies%e_nonlocalpsp = energies%e_nonlocalpsp + & 1347 & dtset%wtk(ikpt)*occ(bdtot_index)*enlnk(bdtot_index) 1348 if (usefock) then 1349 energies%e_fock=energies%e_fock + half*focknk(bdtot_index)*occ(bdtot_index)*dtset%wtk(ikpt) 1350 if (optforces>0) fock%fock_common%forces(:,:)=fock%fock_common%forces(:,:)+& 1351 & dtset%wtk(ikpt)*occ(bdtot_index)*fockfornk(:,:,bdtot_index) 1352 end if 1353 if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ(bdtot_index)*grnlnk(:,bdtot_index) 1354 end if 1355 bdtot_index=bdtot_index+1 1356 if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then 1357 do iband1=1,nband_k 1358 if(paw_dmft%band_in(iband).and.paw_dmft%band_in(iband1)) then 1359 ! write(std_out,*) "II+", isppol,ikpt,iband,iband1 1360 ekindmft2=ekindmft2 + dtset%wtk(ikpt)*paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*& 1361 & eknk_nd(isppol,ikpt,1,iband,iband1) 1362 if(dtset%nspinor==2) then 1363 ekindmft2=ekindmft2 - dtset%wtk(ikpt)*paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*& 1364 & eknk_nd(isppol,ikpt,2,iband,iband1) 1365 end if 1366 ! write(std_out,*) "II", occnd(1,iband,iband1,ikpt,isppol),eknk_nd(isppol,ikpt,iband,iband1) 1367 end if 1368 end do 1369 end if 1370 end do 1371 end do 1372 end do 1373 1374 if(paw_dmft%use_dmft==1) then 1375 energies%e_kinetic = energies%e_kinetic -ekindmft+ekindmft2 1376 if(abs(dtset%pawprtvol)>=2) then 1377 write(message,'(4a,7(2x,a,2x,e14.7,a),a)') & 1378 & "-----------------------------------------------",ch10,& 1379 & "--- Energy for DMFT and tests (in Ha) ",ch10,& 1380 & "--- Ebandldatot (Ha.) = ",ebandldatot,ch10,& 1381 & "--- Ebandlda (Ha.) = ",ebandlda,ch10,& 1382 & "--- Ebanddmft (Ha.) = ",ebanddmft,ch10,& 1383 & "--- Ekinlda (Ha.) = ",ekinlda,ch10, & 1384 & "--- Ekindmftdiag (Ha.) = ",ekindmft,ch10,& 1385 & "--- Ekindmftnondiag(Ha.) = ",ekindmft2,ch10,& 1386 & "--- Edmft= (Ha.) = ",edmft,ch10,& 1387 & "-----------------------------------------------" 1388 call wrtout(std_out,message,'COLL') 1389 end if 1390 if(paw_dmft%use_dmft==1.and.mpi_enreg%paral_kgb==1) paw_dmft%use_dmft=0 1391 end if 1392 1393 if (psps%usepaw==0) then 1394 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1395 & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 1396 else 1397 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1398 & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 1399 end if 1400 call timab(992,2,tsec) 1401 1402 ! Treat fixed occupation numbers or non-self-consistent case 1403 else 1404 1405 if (mpi_enreg%nproc_kpt>1) then 1406 1407 call timab(989,1,tsec) 1408 1409 nbuf=2*mbdkpsp+dtset%nfft*dtset%nspden+3+3*natom*optforces 1410 ! * If Hartree-Fock calculation, the exact exchange energy is k-dependent. 1411 if(dtset%usefock==1) then 1412 nbuf=nbuf+1 1413 if (optforces>0) nbuf=nbuf+3*natom 1414 end if 1415 if(iscf==-1 .or. iscf==-2)nbuf=2*mbdkpsp 1416 ABI_ALLOCATE(buffer1,(nbuf)) 1417 ! Pack eigen,resid,rho[wf]r,grnl,enl,ek 1418 buffer1(1:mbdkpsp)=eigen(:) 1419 buffer1(1+mbdkpsp:2*mbdkpsp)=resid(:) 1420 index1=2*mbdkpsp 1421 if(iscf/=-1 .and. iscf/=-2)then 1422 if (psps%usepaw==0) then 1423 buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhor,(/dtset%nfft*dtset%nspden/)) 1424 else 1425 buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhowfr,(/dtset%nfft*dtset%nspden/)) 1426 end if 1427 index1=index1+dtset%nfft*dtset%nspden 1428 buffer1(index1+1) = energies%e_kinetic 1429 buffer1(index1+2) = energies%e_eigenvalues 1430 buffer1(index1+3) = energies%e_nonlocalpsp 1431 index1=index1+3 1432 ! * If Hartree-Fock calculation, save e_fock in buffer1 1433 if (dtset%usefock==1) then 1434 buffer1(index1+1) = energies%e_fock 1435 index1=index1+1 1436 if (optforces>0)then 1437 buffer1(index1+1:index1+3*natom)=reshape(fock%fock_common%forces,(/3*natom/)) 1438 index1=index1+3*natom 1439 end if 1440 end if 1441 if (optforces>0) buffer1(index1+1:index1+3*natom)=grnl(1:3*natom) 1442 end if 1443 1444 ! Build sum of everything 1445 call timab(48,1,tsec) 1446 call xmpi_sum(buffer1,nbuf,mpi_enreg%comm_kpt ,ierr) 1447 call timab(48,2,tsec) 1448 1449 ! Unpack the final result 1450 eigen(:)=buffer1(1:mbdkpsp) 1451 resid(:)=buffer1(1+mbdkpsp:2*mbdkpsp) 1452 index1=2*mbdkpsp 1453 if(iscf/=-1 .and. iscf/=-2)then 1454 if (psps%usepaw==0) then 1455 ii=1 1456 do ispden=1,dtset%nspden 1457 do ifft=1,dtset%nfft 1458 rhor(ifft,ispden)=buffer1(index1+ii) 1459 ii=ii+1 1460 end do 1461 end do 1462 else 1463 ii=1 1464 do ispden=1,dtset%nspden 1465 do ifft=1,dtset%nfft 1466 rhowfr(ifft,ispden)=buffer1(index1+ii) 1467 ii=ii+1 1468 end do 1469 end do 1470 end if 1471 index1=index1+dtset%nfft*dtset%nspden 1472 energies%e_kinetic = buffer1(index1+1) 1473 energies%e_eigenvalues = buffer1(index1+2) 1474 energies%e_nonlocalpsp = buffer1(index1+3) 1475 index1=index1+3 1476 ! * If Hartree-Fock calculation, save e_fock in buffer1 1477 if (dtset%usefock==1) then 1478 energies%e_fock = buffer1(index1+1) 1479 index1=index1+1 1480 if (optforces>0) then 1481 fock%fock_common%forces(:,:)=reshape(buffer1(index1+1:index1+3*natom),(/3,natom/)) 1482 index1=index1+3*natom 1483 end if 1484 end if 1485 if (optforces>0) grnl(1:3*natom)=buffer1(index1+1:index1+3*natom) 1486 end if 1487 ABI_DEALLOCATE(buffer1) 1488 call timab(989,2,tsec) 1489 1490 end if ! nproc_kpt>1 1491 1492 1493 ! Compute the highest occupied eigenenergy 1494 if(iscf/=-1 .and. iscf/=-2)then 1495 call timab(993,1,tsec) 1496 energies%e_fermie = -huge(one) 1497 bdtot_index=1 1498 do isppol=1,dtset%nsppol 1499 do ikpt=1,dtset%nkpt 1500 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1501 do iband=1,nband_k 1502 if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermie+tol10) then 1503 energies%e_fermie=eigen(bdtot_index) 1504 end if 1505 bdtot_index=bdtot_index+1 1506 end do 1507 end do 1508 end do 1509 call xmpi_max(energies%e_fermie,spaceComm_distrb,ierr) 1510 call timab(993,2,tsec) 1511 end if 1512 1513 ! If needed, compute rhog, and symmetrizes the density 1514 if (iscf > 0 .or. iscf==-3 ) then 1515 ! energies%e_fermie=zero ! Actually, should determine the maximum of the valence band XG20020802 1516 nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3) 1517 1518 call timab(994,1,tsec) 1519 if (psps%usepaw==0) then 1520 call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,& 1521 & dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhog ,rhor ,rprimd,dtset%symafm,dtset%symrel) 1522 else 1523 call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,& 1524 & dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhowfg,rhowfr,rprimd,dtset%symafm,dtset%symrel) 1525 end if 1526 call timab(994,2,tsec) 1527 ! We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2 1528 ! we also have the spin-up density, symmetrized, in rhor(:,2). 1529 end if 1530 1531 end if ! End of test on varying or fixed occupation numbers 1532 1533 call timab(994,1,tsec) 1534 1535 1536 ! Compute the kinetic energy density 1537 if(dtset%usekden==1 .and. (iscf > 0 .or. iscf==-3 ) )then 1538 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1539 & taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1540 end if 1541 1542 ABI_DEALLOCATE(eknk) 1543 if (usefock) then 1544 ABI_DEALLOCATE(focknk) 1545 if (optforces>0)then 1546 ABI_DEALLOCATE(fockfornk) 1547 end if 1548 end if 1549 ABI_DEALLOCATE(eknk_nd) 1550 ABI_DEALLOCATE(grnlnk) 1551 ABI_DEALLOCATE(enlnk) 1552 1553 ! In the non-self-consistent case, print eigenvalues and residuals 1554 if(iscf<=0)then 1555 option=2 ; enunit=1 ; vxcavg_dum=zero 1556 call prteigrs(eigen,enunit,energies%e_fermie,dtfil%fnameabo_app_eig,& 1557 & ab_out,iscf,dtset%kptns,dtset%kptopt,dtset%mband,dtset%nband,& 1558 & dtset%nkpt,nnsclo_now,dtset%nsppol,occ,dtset%occopt,option,& 1559 & dtset%prteig,prtvol,resid,dtset%tolwfr,vxcavg_dum,dtset%wtk) 1560 end if 1561 1562 ! Find largest residual over bands, k points, and spins, except for nbdbuf highest bands 1563 ibdkpt=1 1564 residm=zero 1565 do isppol=1,dtset%nsppol 1566 do ikpt=1,dtset%nkpt 1567 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1568 nband_eff=max(1,nband_k-dtset%nbdbuf) 1569 residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1))) 1570 ibdkpt=ibdkpt+nband_k 1571 end do 1572 end do 1573 1574 end if !usewvl==0 1575 1576 !=================================================================== 1577 ! End of PLANE WAVES section 1578 !=================================================================== 1579 1580 !In the self-consistent case, diagnose lack of unoccupied state (for each spin and k-point). 1581 !Print a warning if the number of such messages already written does not exceed mwarning. 1582 mwarning=5 1583 if(nwarning<mwarning .and. iscf>=0)then 1584 nwarning=nwarning+1 1585 bdtot_index=1 1586 do isppol=1,dtset%nsppol 1587 do ikpt=1,dtset%nkpt 1588 min_occ=two 1589 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1590 do iband=1,nband_k 1591 if(occ(bdtot_index)<min_occ)min_occ=occ(bdtot_index) 1592 bdtot_index=bdtot_index+1 1593 end do 1594 if(min_occ>0.01_dp)then 1595 if(dtset%nsppol==1)then 1596 write(message, '(a,i0,3a,f7.3,5a)' )& 1597 & 'For k-point number ',ikpt,',',ch10,& 1598 & 'The minimal occupation factor is',min_occ,'.',ch10,& 1599 & 'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,& 1600 & 'Action: increase slightly the number of bands.' 1601 else 1602 write(message, '(a,i0,3a,i0,a,f7.3,5a)' )& 1603 & 'For k-point number ',ikpt,', and',ch10,& 1604 & 'for spin polarization',isppol,'the minimal occupation factor is',min_occ,'.',ch10,& 1605 & 'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,& 1606 & 'Action: increase slightly the number of bands.' 1607 end if 1608 MSG_WARNING(message) 1609 exit ! It is enough if one lack of adequate occupation is identified, so exit. 1610 end if 1611 end do 1612 end do 1613 end if 1614 1615 if (iscf>0.or.iscf==-3 .or. (dtset%usewvl==1 .and. iscf==0)) then 1616 1617 ! PAW: Build new rhoij quantities from new occ then symetrize them 1618 ! Compute and add the compensation density to rhowfr to get the total density 1619 if (psps%usepaw==1) then 1620 call timab(555,1,tsec) 1621 if (paral_atom) then 1622 ABI_DATATYPE_ALLOCATE(pawrhoij_unsym,(natom)) 1623 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 1624 call pawrhoij_alloc(pawrhoij_unsym,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,& 1625 & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0) 1626 else 1627 pawrhoij_unsym => pawrhoij 1628 end if 1629 if (usecprj_local==1) then 1630 call pawmkrhoij(atindx,atindx1,cprj,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,& 1631 & mcprj_local,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,& 1632 & occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,dtfil%unpaw,& 1633 & dtset%usewvl,dtset%wtk) 1634 else 1635 mcprj_tmp=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol 1636 ABI_DATATYPE_ALLOCATE(cprj_tmp,(natom,mcprj_tmp)) 1637 call pawcprj_alloc(cprj_tmp,0,gs_hamk%dimcprj) 1638 call ctocprj(atindx,cg,1,cprj_tmp,gmet,gprimd,0,0,0,dtset%istwfk,kg,dtset%kptns,& 1639 & mcg,mcprj_tmp,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,& 1640 & dtset%natom,nattyp,dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg,& 1641 & npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,& 1642 & ucvol,dtfil%unpaw,xred,ylm,ylmgr_dum) 1643 call pawmkrhoij(atindx,atindx1,cprj_tmp,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,& 1644 & dtset%mband,mband_cprj,mcprj_tmp,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,& 1645 & dtset%nspinor,dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,& 1646 & dtfil%unpaw,dtset%usewvl,dtset%wtk) 1647 call pawcprj_free(cprj_tmp) 1648 ABI_DATATYPE_DEALLOCATE(cprj_tmp) 1649 end if 1650 call timab(555,2,tsec) 1651 ! Build symetrized packed rhoij and compensated pseudo density 1652 cplex=1;ipert=0;idir=0;qpt(:)=zero 1653 if(dtset%usewvl==0) then 1654 call pawmkrho(compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,& 1655 & my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,& 1656 & dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,& 1657 & symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog) 1658 else 1659 ! here do not pass rhog, we do not use it 1660 call pawmkrho(compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,& 1661 & my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,& 1662 & dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,& 1663 & symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat) 1664 ! In WVL: copy density to BigDFT object: 1665 call wvl_rho_abi2big(1,rhor,wvl%den) 1666 end if 1667 if (paral_atom) then 1668 call pawrhoij_free(pawrhoij_unsym) 1669 ABI_DATATYPE_DEALLOCATE(pawrhoij_unsym) 1670 end if 1671 end if 1672 1673 if(paw_dmft%use_dmft==1) then 1674 ! == check noccmmp 1675 call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,dtset%nsppol,0,ntypat,& 1676 & paw_ij,pawang,dtset%pawprtvol,pawrhoij,pawtab,rdum2,idum1,dtset%typat,0,dtset%usepawu,& 1677 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1678 end if 1679 1680 ! Find and print minimum and maximum total electron density and locations 1681 ! Compute density residual (if required) and its squared norm 1682 if (iscf>=0) then 1683 if (psps%usepaw==0) then 1684 call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,ucvol=ucvol) 1685 else 1686 call prtrhomxmn(std_out,mpi_enreg,nfftf,pawfgr%ngfft,dtset%nspden,1,rhor,ucvol=ucvol) 1687 end if 1688 if (optres==1) then 1689 nvresid=rhor-nvresid 1690 call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid,mpi_comm_sphgrid=mpi_comm_sphgrid) 1691 end if 1692 end if 1693 1694 end if ! iscf>0 or iscf=-3 1695 1696 if(psps%usepaw==1.and.(iscf>=0.or.iscf==-3)) then 1697 ABI_DEALLOCATE(rhowfr) 1698 ABI_DEALLOCATE(rhowfg) 1699 end if 1700 1701 call timab(994,2,tsec) 1702 1703 if (iscf==-1) then 1704 ! Eventually compute the excited states within tddft 1705 call timab(995,1,tsec) 1706 if (psps%usepaw==1) then 1707 ! In case of PAW calculation, have to transfer kxc from the fine to the coarse grid: 1708 ABI_ALLOCATE(cgrkxc,(dtset%nfft,nkxc)) 1709 do ikxc=1,nkxc 1710 call transgrid(1,mpi_enreg,1,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrkxc(:,ikxc),kxc(:,ikxc)) 1711 end do 1712 call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,& 1713 & kg,cgrkxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,& 1714 & ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew) 1715 ABI_DEALLOCATE(cgrkxc) 1716 else 1717 call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,& 1718 & kg,kxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,& 1719 & ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew) 1720 end if 1721 call timab(995,2,tsec) 1722 1723 else 1724 ! Eventually compute the susceptibility matrix and the 1725 ! dielectric matrix when istep_mix is equal to 1 or dielstrt 1726 call timab(996,1,tsec) 1727 call testsusmat(computesusmat,dielop,dielstrt,dtset,istep_mix) !test if the matrix is to be computed 1728 if(computesusmat) then 1729 dielar(1)=dtset%diecut;dielar(2)=dtset%dielng 1730 dielar(3)=dtset%diemac;dielar(4)=dtset%diemix 1731 dielar(5)=dtset%diegap;dielar(6)=dtset%dielam 1732 dielar(7)=dtset%diemix;if (iscf>=10) dielar(7)=dtset%diemixmag 1733 usetimerev=1 1734 if (psps%usepaw==1.and.dtset%pawspnorb>0.and.dtset%kptopt/=1.and.dtset%kptopt/=2) usetimerev=0 1735 neglect_pawhat=1-dtset%pawsushat 1736 call suscep_stat(atindx,atindx1,cg,cprj,dielar,& 1737 & gs_hamk%dimcprj,doccde,eigen,gbound_diel,gprimd,& 1738 & irrzondiel,dtset%istwfk,kg,kg_diel,lmax_diel,& 1739 & dtset%mband,mcg,mcprj_local,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,natom,dtset%nband,& 1740 & neglect_pawhat,nfftdiel,ngfftdiel,& 1741 & dtset%nkpt,npwarr,npwdiel,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,& 1742 & occ,dtset%occopt,pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,& 1743 & susmat,dtset%symafm,dtset%symrel,dtset%tnons,dtset%typat,ucvol,& 1744 & dtfil%unpaw,usecprj_local,psps%usepaw,usetimerev,dtset%wtk,ylmdiel) 1745 end if 1746 call timab(996,2,tsec) 1747 1748 end if ! end condition on iscf 1749 1750 call destroy_hamiltonian(gs_hamk) 1751 1752 if(dtset%usewvl==0) then 1753 ABI_DEALLOCATE(EigMin) 1754 ABI_DEALLOCATE(doccde) 1755 #if defined HAVE_GPU_CUDA 1756 if(dtset%use_gpu_cuda==1) then 1757 call gpu_finalize_ham_data() 1758 end if 1759 #endif 1760 end if 1761 1762 call timab(980,2,tsec) 1763 1764 DBG_EXIT("COLL") 1765 1766 contains
ABINIT/wvl_comm_eigen [ Functions ]
NAME
wvl_comm_eigen
FUNCTION
Computes occupations for the wavelet case Using BigDFT routines
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (D.Caliste, T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
2269 subroutine wvl_comm_eigen() 2270 2271 use defs_basis 2272 use m_errors 2273 2274 use m_abi2big, only : wvl_eigen_abi2big 2275 #if defined HAVE_BIGDFT 2276 use BigDFT_API, only: eigensystem_info 2277 #endif 2278 2279 !This section has been created automatically by the script Abilint (TD). 2280 !Do not modify the following lines by hand. 2281 #undef ABI_FUNC 2282 #define ABI_FUNC 'wvl_comm_eigen' 2283 !End of the abilint section 2284 2285 implicit none 2286 2287 !Arguments ------------------------------------ 2288 2289 !Local variables------------------------------- 2290 #if defined HAVE_BIGDFT 2291 integer:: ikpt,norb,shift 2292 #endif 2293 2294 ! ************************************************************************* 2295 2296 DBG_ENTER("COLL") 2297 2298 #if defined HAVE_BIGDFT 2299 if(wvlbigdft) then 2300 ! Communicates eigenvalues to all procs. 2301 ! This will print out the eigenvalues and Fermi level. 2302 call eigensystem_info(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl,0.d0,& 2303 & wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,& 2304 & wvl%wfs%ks%orbs,wvl%wfs%ks%psi) 2305 else 2306 ! Send all eigenvalues to all procs. 2307 ! I simply communicate eigenvalues: I do not print them into screen, nor calculate Fermi-level. 2308 norb=wvl%wfs%ks%orbs%norb 2309 if (mpi_enreg%nproc_wvl > 1) then 2310 shift=1 2311 do ikpt = 1, wvl%wfs%ks%orbs%nkpts 2312 call xmpi_bcast(wvl%wfs%ks%orbs%eval(shift:shift+norb-1),wvl%wfs%ks%orbs%ikptproc(ikpt),mpi_enreg%comm_wvl,ierr) 2313 shift=shift+norb 2314 end do 2315 end if 2316 end if 2317 2318 !Copy eigenvalues from BigDFT object to "eigen" 2319 call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs) 2320 2321 #else 2322 BIGDFT_NOTENABLED_ERROR() 2323 #endif 2324 2325 DBG_EXIT("COLL") 2326 2327 end subroutine wvl_comm_eigen 2328 2329 end subroutine vtorho
ABINIT/wvl_nscf_loop [ Functions ]
NAME
wvl_nscf_loop
FUNCTION
Non-self-consistent field cycle in Wavelets See also "wvl_nscf_loop_bigdft"
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
nnsclo= number of non-self consistent field iterations
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
1800 subroutine wvl_nscf_loop() 1801 1802 use defs_basis 1803 use defs_abitypes 1804 use defs_wvltypes 1805 1806 use m_errors 1807 use m_energies, only : energies_type 1808 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 1809 1810 !This section has been created automatically by the script Abilint (TD). 1811 !Do not modify the following lines by hand. 1812 #undef ABI_FUNC 1813 #define ABI_FUNC 'wvl_nscf_loop' 1814 use interfaces_62_wvl_wfs 1815 !End of the abilint section 1816 1817 implicit none 1818 1819 !Arguments ------------------------------------ 1820 ! integer, intent(in) :: istep,mcprj,nfft,nnsclo 1821 ! real(dp), intent(inout) :: residm 1822 ! type(dataset_type), intent(in) :: dtset 1823 ! type(MPI_type), intent(in) :: mpi_enreg 1824 ! type(energies_type), intent(inout) :: energies 1825 ! type(wvl_data), intent(inout) :: wvl 1826 ! !arrays 1827 ! real(dp), intent(inout) :: xcart(3, dtset%natom) 1828 ! real(dp), dimension(6), intent(out) :: strsxc 1829 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj 1830 1831 !Local variables------------------------------- 1832 integer :: inonsc,ii 1833 integer,parameter :: iscf_=-1 !do not do a SCF cycle 1834 logical,parameter :: do_scf=.false. !do not do a SCF cycle 1835 logical,parameter :: wvlbigdft=.false. 1836 real(dp) :: dum,eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc 1837 1838 ! ************************************************************************* 1839 1840 DBG_ENTER("COLL") 1841 1842 if(nnsclo_now>0) then 1843 do inonsc=1,nnsclo_now 1844 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,& 1845 & istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,& 1846 & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,& 1847 & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc) 1848 call wvl_hpsitopsi(cprj,dtset,energies,inonsc,mcprj_local,mpi_enreg, & 1849 & residm,wvl,xcart) 1850 if(residm<dtset%tolwfr) exit !Exit loop if converged 1851 end do 1852 1853 else 1854 do ii=1, dtset%nline 1855 ! Direct minimization technique: no diagonalization 1856 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,& 1857 & istep,ii,iscf_,mpi_enreg%me_wvl,dtset%natom,& 1858 & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,& 1859 & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc) 1860 call wvl_hpsitopsi(cprj,dtset,energies,ii,mcprj_local,mpi_enreg, & 1861 & residm,wvl,xcart) 1862 if(residm<dtset%tolwfr) exit !Exit loop if converged 1863 end do 1864 end if 1865 1866 ! Update energies depending on new WF 1867 energies%e_kinetic=ekin 1868 energies%e_nonlocalpsp=enl 1869 energies%e_exactX=eexctx 1870 energies%e_sicdc=esicdc 1871 1872 ! Eventually update energies depending on density 1873 if (dtset%iscf<10) then 1874 energies%e_localpsp=eloc 1875 energies%e_hartree=eh 1876 energies%e_xc=exc ; energies%e_xcdc=evxc 1877 else if (nnsclo_now==0) then 1878 energies%e_localpsp=eloc 1879 end if 1880 1881 ! End of nscf iterations 1882 if (do_last_ortho) then 1883 ! !Don't update energies (nscf cycle has been done); just recompute potential 1884 inonsc=nnsclo_now;if (nnsclo_now==0) inonsc=dtset%nline 1885 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,& 1886 & istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,& 1887 & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,& 1888 & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc) 1889 end if 1890 1891 DBG_EXIT("COLL") 1892 1893 end subroutine wvl_nscf_loop
ABINIT/wvl_nscf_loop_bigdft [ Functions ]
NAME
wvl_nscf_loop_bigdft
FUNCTION
Non-self-consistent field cycle in Wavelets It follows the BigDFT scheme. See also "wvl_nscf_loop"
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (T. Rangel, D. Caliste) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
nnsclo= number of non-self consistent field iterations
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
1929 #if defined HAVE_CONFIG_H 1930 #include "config.h" 1931 #endif 1932 1933 #include "abi_common.h" 1934 1935 subroutine wvl_nscf_loop_bigdft() 1936 1937 use defs_basis 1938 use defs_abitypes 1939 use defs_wvltypes 1940 use m_errors 1941 1942 use m_energies, only : energies_type 1943 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 1944 1945 !This section has been created automatically by the script Abilint (TD). 1946 !Do not modify the following lines by hand. 1947 #undef ABI_FUNC 1948 #define ABI_FUNC 'wvl_nscf_loop_bigdft' 1949 use interfaces_62_wvl_wfs 1950 !End of the abilint section 1951 1952 implicit none 1953 1954 !Arguments ------------------------------------ 1955 ! integer, intent(in) :: istep,mcprj,nfft,nnsclo 1956 ! real(dp), intent(inout) :: residm 1957 ! real(dp), intent(out) :: nres2 1958 ! type(dataset_type), intent(in) :: dtset 1959 ! type(MPI_type), intent(in) :: mpi_enreg 1960 ! type(energies_type), intent(inout) :: energies 1961 ! type(wvl_data), intent(inout) :: wvl 1962 !arrays 1963 ! real(dp), intent(inout) :: xcart(3, dtset%natom) 1964 ! real(dp), dimension(6), intent(out) :: strsxc 1965 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj 1966 1967 !Local variables------------------------------- 1968 integer :: inonsc 1969 integer,parameter :: iscf_=-1 !do not do a SCF cycle 1970 logical,parameter :: do_scf=.false. !do not do a SCF cycle 1971 logical,parameter :: wvlbigdft=.true. 1972 real(dp) :: eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc 1973 1974 ! ************************************************************************* 1975 1976 DBG_ENTER("COLL") 1977 1978 call wvl_hpsitopsi(cprj,dtset, energies, istep, mcprj_local,mpi_enreg, & 1979 & residm, wvl,xcart) 1980 1981 if (nnsclo_now>2) then 1982 do inonsc = 2, nnsclo_now-1 1983 call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, & 1984 & energies%e_hartree, energies%e_kinetic, energies%e_localpsp, & 1985 & energies%e_nonlocalpsp, energies%e_sicdc, istep, inonsc, iscf_, & 1986 & mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,& 1987 & dtset%nspden, nres2, do_scf,energies%e_xcdc, & 1988 & wvl, wvlbigdft, xcart, strsxc) 1989 call wvl_hpsitopsi(cprj,dtset, energies, inonsc, mcprj_local,mpi_enreg, & 1990 & residm, wvl,xcart) 1991 end do 1992 end if 1993 1994 ! End of nscf iterations 1995 if (do_last_ortho.and.nnsclo_now<=1) then 1996 ! !Don't update energies (nscf cycle has been done); just recompute potential 1997 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc, & 1998 & istep, 1, iscf_, mpi_enreg%me_wvl, dtset%natom, nfftf, & 1999 & mpi_enreg%nproc_wvl,dtset%nspden, nres2, do_scf,evxc, & 2000 & wvl, wvlbigdft, xcart, strsxc) 2001 else if (do_last_ortho.and.nnsclo_now>1) then 2002 ! !Update energies and potential (nscf cycles are not finished) 2003 call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, & 2004 & energies%e_hartree,energies%e_kinetic, energies%e_localpsp, & 2005 & energies%e_nonlocalpsp, energies%e_sicdc, istep, nnsclo_now, iscf_, & 2006 & mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,& 2007 & dtset%nspden, nres2, do_scf,energies%e_xcdc, & 2008 & wvl, wvlbigdft, xcart, strsxc) 2009 end if 2010 2011 DBG_EXIT("COLL") 2012 2013 end subroutine wvl_nscf_loop_bigdft
ABINIT/wvl_occ [ Functions ]
NAME
wvl_occ
FUNCTION
Computes occupations for the wavelet case
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
2130 subroutine wvl_occ() 2131 2132 use defs_basis 2133 use m_errors 2134 use defs_wvltypes 2135 2136 !This section has been created automatically by the script Abilint (TD). 2137 !Do not modify the following lines by hand. 2138 #undef ABI_FUNC 2139 #define ABI_FUNC 'wvl_occ' 2140 !End of the abilint section 2141 2142 implicit none 2143 2144 !Local variables------------------------------- 2145 real(dp):: doccde_(dtset%mband*dtset%nkpt*dtset%nsppol) 2146 ! ************************************************************************* 2147 2148 DBG_ENTER("COLL") 2149 2150 ! Compute the new occupation numbers from eigen 2151 call newocc(doccde_,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,& 2152 & dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,& 2153 & dtset%nsppol,occ,dtset%occopt,prtvol,& 2154 & dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk) 2155 2156 ! Copy occupations and efermi to BigDFT variables 2157 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs) 2158 2159 #if defined HAVE_BIGDFT 2160 ! Copy Fermi level to BigDFT variable: 2161 wvl%wfs%ks%orbs%efermi=energies%e_fermie 2162 #endif 2163 2164 DBG_EXIT("COLL") 2165 2166 end subroutine wvl_occ
ABINIT/wvl_occ_bigdft [ Functions ]
NAME
wvl_occ_bigdft
FUNCTION
Computes occupations for the wavelet case Using BigDFT routines
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (D.Caliste, T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
PARENTS
vtorho
CHILDREN
eigensystem_info,wvl_eigen_abi2big,xmpi_bcast
SOURCE
2200 subroutine wvl_occ_bigdft() 2201 2202 use defs_basis 2203 use m_errors 2204 2205 !This section has been created automatically by the script Abilint (TD). 2206 !Do not modify the following lines by hand. 2207 #undef ABI_FUNC 2208 #define ABI_FUNC 'wvl_occ_bigdft' 2209 !End of the abilint section 2210 2211 implicit none 2212 2213 ! ************************************************************************* 2214 2215 DBG_ENTER("COLL") 2216 2217 ! Transfer occopt from ABINIT to BigDFT 2218 #if defined HAVE_BIGDFT 2219 occopt_bigdft=dtset%occopt 2220 call wvl_occopt_abi2big(occopt_bigdft,occopt_bigdft,1) 2221 2222 !Calculate occupations using BigDFT routine 2223 call evaltoocc(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, .false., & 2224 & dtset%tsmear, wvl%wfs%ks%orbs, occopt_bigdft) 2225 2226 !Pass e_fermi from BigDFT object to ABINIT variable: 2227 energies%e_fermie = wvl%wfs%ks%orbs%efermi 2228 2229 !Copy occupations from BigDFT to ABINIT variables 2230 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs) 2231 #endif 2232 2233 DBG_EXIT("COLL") 2234 2235 end subroutine wvl_occ_bigdft