TABLE OF CONTENTS
ABINIT/m_odamix [ Modules ]
NAME
m_odamix
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (FJ, 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_odamix 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_xcdata 29 use m_dtset 30 31 32 use defs_datatypes, only : pseudopotential_type 33 use defs_abitypes, only : MPI_type 34 use m_time, only : timab 35 use m_geometry, only : metric 36 use m_cgtools, only : dotprod_vn 37 use m_pawang, only : pawang_type 38 use m_pawrad, only : pawrad_type 39 use m_pawtab, only : pawtab_type 40 use m_paw_an, only : paw_an_type 41 use m_paw_ij, only : paw_ij_type 42 use m_pawfgrtab, only : pawfgrtab_type 43 use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter 44 use m_paw_nhat, only : pawmknhat 45 use m_paw_denpot, only : pawdenpot 46 use m_energies, only : energies_type 47 use m_spacepar, only : hartre 48 use m_rhotoxc, only : rhotoxc 49 use m_fft, only : fourdp 50 use m_xc_tb09, only : xc_tb09_update_c 51 52 implicit none 53 54 private
ABINIT/odamix [ Functions ]
NAME
odamix
FUNCTION
This routine is called to compute the total energy and various parts of it. The routine computes -if requested- the forces.
INPUTS
[add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy dtset <type(dataset_type)>=all input variables in this dataset berryopt = 4/14: electric field is on -> add the contribution of the -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy = 6/16, or 7/17: electric displacement field is on -> add the contribution of the Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy | berryopt = 5: magnetic field is on -> add the contribution of the | - \Omega B.M term to the total energy | /= 5: magnetic field is off | bfield = cartesian coordinates of the magnetic field in atomic units | dfield = cartesian coordinates of the electric displacement field in atomic units (berryopt==6/7) | efield = cartesian coordinates of the electric field in atomic units (berryopt==4) | red_dfield = reduced the electric displacement field (berryopt==16/17) | red_efieldbar = reduced the electric field (ebar) (berryopt==14) | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen | ionmov=governs the movement of atoms (see help file) | natom=number of atoms in cell. | nconeq=number of atomic constraint equations | nspden=number of spin-density components | nsym=number of symmetry elements in space group | occopt=option for occupancies | prtvol=integer controlling volume of printed output | tsmear=smearing energy or temperature (if metal) | wtatcon(3,natom,nconeq)=weights for atomic constraints | xclevel= XC functional level gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkxc=second dimension of the array kxc, see rhotoxc.f for a description ntypat=number of types of atoms in unit cell. nvresid(nfft,nspden)=potential or density residual n3xccc=dimension of the xccc3d array (0 or nfft). optres=0 if residual array (nvresid) contains the potential residual =1 if residual array (nvresid) contains the density residual paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh pawang <type(pawang_type)>=paw angular mesh and related data pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrad pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=array for Fourier transform of electron density rhor(nfft,nspden)=array for electron density in electrons/bohr**3 rprimd(3,3)=dimensional primitive translations in real space (bohr) [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density ucvol = unit cell volume (Bohr**3) usepaw= 0 for non paw calculation; =1 for paw calculation vhartr(nfft)=array for holding Hartree potential vpsp(nfft)=array for holding local psp vxc(nfft,nspden)=array for holding XC potential xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
deltae=change in total energy between the previous and present SCF cycle etotal=total energy (hartree)
SIDE EFFECTS
Input/Output: elast=previous value of the energy, needed to compute deltae, then updated. energies <type(energies_type)>=all part of total energy. | entropy(IN)=entropy due to the occupation number smearing (if metal) | e_localpsp(IN)=local psp energy (hartree) | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree) | e_chempot(IN)=energy from spatially varying chemical potential (hartree) | e_ewald(IN)=Ewald energy (hartree) | e_vdw_dftd(IN)=VdW DFT-D energy | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_kinetic(IN)=kinetic energy part of total energy. | e_nucdip(IN)=energy of nuclear dipole array | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy. | e_xc(IN)=exchange-correlation energy (hartree) | e_xcdc(IN)=exchange-correlation double-counting energy (hartree) | e_paw(IN)=PAW spherical part energy | e_pawdc(IN)=PAW spherical part double-counting energy | e_elecfield(OUT)=the term of the energy functional that depends explicitely | on the electric field: enefield = -ucvol*E*P | e_magfield(OUT)=the term of the energy functional that depends explicitely | on the magnetic field: enmagfield = -ucvol*B*M | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal) | this value is %entropy * dtset%tsmear (hartree). kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to kinetic energy density (metaGGA cases) (optional output) ===== if psps%usepaw==1 pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (gradients of rhoij for each atom with respect to atomic positions are computed here)
NOTES
In case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
! Developpers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In case of norm-conserving calculations the FFT grid is the usual FFT grid.
SOURCE
176 subroutine odamix(deltae,dtset,elast,energies,etotal,& 177 & gprimd,gsqcut,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 178 & nkxc,ntypat,nvresid,n3xccc,optres,paw_ij,& 179 & paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,& 180 & red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,usepaw,& 181 & usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,& 182 & taur,vxctau,add_tfw) ! optional arguments 183 184 !Arguments ------------------------------------ 185 !scalars 186 integer,intent(in) :: my_natom,n3xccc,nfft,nkxc,ntypat,optres 187 integer,intent(in) :: usepaw,usexcnhat 188 logical,intent(in),optional :: add_tfw 189 real(dp),intent(in) :: gsqcut,ucvol 190 real(dp),intent(inout) :: elast 191 real(dp),intent(out) :: deltae,etotal,vxcavg 192 type(MPI_type),intent(in) :: mpi_enreg 193 type(dataset_type),intent(in) :: dtset 194 type(energies_type),intent(inout) :: energies 195 type(pawang_type),intent(in) :: pawang 196 type(pseudopotential_type),intent(in) :: psps 197 !arrays 198 integer,intent(in) :: ngfft(18) 199 logical :: add_tfw_ 200 real(dp),intent(in) :: gprimd(3,3) 201 real(dp),intent(in) :: red_ptot(3),rprimd(3,3),vpsp(nfft),xred(3,dtset%natom) 202 real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden) 203 real(dp),intent(inout) :: kxc(nfft,nkxc),nhat(nfft,dtset%nspden*usepaw) 204 real(dp),intent(inout) :: nvresid(nfft,dtset%nspden),rhog(2,nfft) 205 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft) 206 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 207 real(dp),intent(inout) :: xccc3d(n3xccc) 208 real(dp),intent(out) :: strsxc(6) 209 real(dp),intent(inout),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 210 type(paw_an_type),intent(inout) :: paw_an(my_natom) 211 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 212 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 213 type(pawrad_type),intent(in) :: pawrad(ntypat) 214 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 215 type(pawtab_type),intent(in) :: pawtab(ntypat) 216 217 !Local variables------------------------------- 218 !scalars 219 integer :: cplex,iatom,ider,idir,ierr,ifft,ipert,irhoij,ispden,itypat,izero,iir,jjr,kkr 220 integer :: jrhoij,klmn,klmn1,kmix,nfftot,nhatgrdim,nzlmopt,nk3xc,option,optxc 221 logical :: nmxc,with_vxctau 222 real(dp) :: alphaopt,compch_fft,compch_sph,doti,e1t10,e_ksnm1,e_xcdc_vxctau 223 real(dp) :: eenth,fp0,gammp1,ro_dlt,ucvol_local 224 character(len=500) :: message 225 type(xcdata_type) :: xcdata 226 !arrays 227 real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3) 228 real(dp) :: gmet(3,3),gprimdlc(3,3),qpt(3),rmet(3,3),tsec(2) 229 real(dp),allocatable :: nhatgr(:,:,:),rhoijtmp(:,:) 230 231 ! ********************************************************************* 232 233 !DEBUG 234 !write(std_out,*)' odamix : enter' 235 !ENDDEBUG 236 237 call timab(80,1,tsec) 238 239 !Check that usekden is not 0 if want to use vxctau 240 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 241 242 !To be adjusted for the call to rhotoxc 243 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 244 nk3xc=1;nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 245 246 !faire un test sur optres=1, usewvl=0, nspden=1,nhatgrdim 247 if(optres/=1)then 248 write(message,'(a,i0,a)')' optres=',optres,', not allowed in oda => stop ' 249 ABI_ERROR(message) 250 end if 251 252 if(dtset%usewvl/=0)then 253 write(message,'(a,i0,a)')' usewvl=',dtset%usewvl,', not allowed in oda => stop ' 254 ABI_ERROR(message) 255 end if 256 257 if(dtset%nspden/=1)then 258 write(message,'(a,i0,a)')' nspden=',dtset%nspden,', not allowed in oda => stop ' 259 ABI_ERROR(message) 260 end if 261 262 if (my_natom>0) then 263 if(paw_ij(1)%has_dijhat==0)then 264 message = ' dijhat variable must be allocated in odamix ! ' 265 ABI_ERROR(message) 266 end if 267 if(paw_ij(1)%cplex_dij==2.or.paw_ij(1)%qphase==2)then 268 message = ' complex dij not allowed in odamix! ' 269 ABI_ERROR(message) 270 end if 271 end if 272 273 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 274 !!!!!!!!!! calculation of f'(0)= Eband_new-EH_old-E_xcdc_old-Ek_old-E_loc_old-E_nonloc_old 275 !!!!!!!!!! save previous energy E(rho_tild_n) 276 277 fp0=energies%e_eigenvalues-energies%h0-two*energies%e_hartree-energies%e_xcdc 278 if (usepaw==1) then 279 do iatom=1,my_natom 280 ABI_CHECK(pawrhoij(iatom)%qphase==1,'ODA mixing not allowed with a Q phase in PAW objects!') 281 itypat=pawrhoij(iatom)%itypat 282 do ispden=1,pawrhoij(iatom)%nspden 283 jrhoij=1 284 do irhoij=1,pawrhoij(iatom)%nrhoijsel 285 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 286 ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 287 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden)) 288 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 289 end do 290 klmn1=1 291 do klmn=1,pawrhoij(iatom)%lmn2_size 292 ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,ispden)*pawtab(itypat)%dltij(klmn) 293 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden)) 294 klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij 295 end do 296 end do 297 if (paw_ij(iatom)%ndij>=2.and.pawrhoij(iatom)%nspden==1) then 298 jrhoij=1 299 do irhoij=1,pawrhoij(iatom)%nrhoijsel 300 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 301 ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,1)*pawtab(itypat)%dltij(klmn) 302 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2)) 303 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 304 end do 305 klmn1=1 306 do klmn=1,pawrhoij(iatom)%lmn2_size 307 ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,1)*pawtab(itypat)%dltij(klmn) 308 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2)) 309 klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij 310 end do 311 e1t10=half*e1t10 312 end if 313 end do 314 if (mpi_enreg%nproc_atom>1) then 315 call xmpi_sum(e1t10,mpi_enreg%comm_atom,ierr) 316 end if 317 fp0=fp0-e1t10 318 end if 319 e_ksnm1=etotal 320 321 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 322 !!!!! Calculation of quantities that do not depend on rho_n+1 323 324 !PAW: eventually recompute compensation density (and gradients) 325 nhatgrdim=0 326 if (usepaw==1) then 327 ider=-1;if (dtset%xclevel==2.or.usexcnhat==0) ider=0 328 if (dtset%xclevel==2.and.usexcnhat==1) ider=ider+2 329 if (ider>0) then 330 nhatgrdim=1 331 ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3)) 332 end if 333 if (ider>=0) then 334 ider=0;izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero 335 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 336 nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,& 337 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,& 338 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 339 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 340 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 341 end if 342 end if 343 344 !------Compute Hartree and xc potentials---------------------------------- 345 nfftot=PRODUCT(ngfft(1:3)) 346 347 call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,& 348 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 349 350 call xcdata_init(xcdata,dtset=dtset) 351 352 !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value 353 if (dtset%xc_tb09_c>99._dp) then 354 call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, & 355 & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, & 356 & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, & 357 & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 358 & computation_type='all') 359 end if 360 361 !Compute xc potential (separate up and down if spin-polarized) 362 optxc=1 363 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 364 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,strsxc,& 365 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 366 367 !------Compute parts of total energy depending on potentials-------- 368 369 ucvol_local=ucvol 370 371 !Compute Hartree energy energies%e_hartree 372 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 373 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 374 energies%e_hartree=half*energies%e_hartree 375 376 377 !Compute local psp energy energies%e_localpsp 378 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,& 379 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 380 381 !Compute double-counting XC energy energies%e_xcdc 382 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 383 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 384 if (with_vxctau) then 385 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,& 386 & vxctau(:,:,1),ucvol_local,mpi_comm_sphgrid=mpi_enreg%comm_fft) 387 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 388 end if 389 390 if (usepaw/=0) then 391 nzlmopt=dtset%pawnzlm; option=2 392 do iatom=1,my_natom 393 itypat=paw_ij(iatom)%itypat 394 ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size)) 395 paw_ij(iatom)%has_dijhartree=1 396 end do 397 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom,dtset%nspden,ntypat,& 398 & dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,& 399 & pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp,& 400 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 401 do iatom=1,my_natom 402 ABI_FREE(paw_ij(iatom)%dijhartree) 403 paw_ij(iatom)%has_dijhartree=0 404 end do 405 end if 406 407 !When the finite-temperature VG broadening scheme is used, 408 !the total entropy contribution "tsmear*entropy" has a meaning, 409 !and gather the two last terms of Eq.8 of the VG paper 410 !Warning : might have to be changed for fixed moment calculations 411 if(dtset%occopt>=3 .and. dtset%occopt<=8) then 412 energies%e_entropy = - dtset%tsmear * energies%entropy 413 else 414 energies%e_entropy = zero 415 end if 416 !Turn it into an electric enthalpy,refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]] 417 ! the missing volume is added here 418 energies%e_elecfield = zero 419 if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG 420 421 energies%e_elecfield = -dot_product(dtset%red_efieldbar,red_ptot) 422 423 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 424 eenth = zero 425 do iir=1,3 426 do jjr=1,3 427 eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 428 end do 429 end do 430 eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth 431 energies%e_elecfield = energies%e_elecfield + eenth 432 433 end if 434 435 energies%e_magfield = zero 436 !if (dtset%berryopt == 5) then 437 !emag = dot_product(mag_cart,dtset%bfield) 438 !energies%e_magfield = emag 439 !end if 440 441 !HONG Turn it into an internal enthalpy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]], 442 !but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j 443 if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then 444 energies%e_elecfield=zero 445 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 446 do iir=1,3 447 do jjr=1,3 448 energies%e_elecfield = energies%e_elecfield + gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 449 end do 450 end do 451 energies%e_elecfield = ucvol_local/(8.0d0*pi)*energies%e_elecfield 452 end if 453 454 !HONG calculate internal energy and electric enthalpy for mixed BC case. 455 if ( dtset%berryopt == 17 ) then 456 energies%e_elecfield = zero 457 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 458 A(:,:)=(4*pi/ucvol_local)*rmet(:,:) 459 A1(:,:)=A(:,:) 460 A_new(:,:)=A(:,:) 461 efield_new(:)=dtset%red_efield(:) 462 eenth = zero 463 464 do kkr=1,3 465 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 466 467 ! step 1 add -ebar*p 468 eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr) 469 470 ! step 2 chang to e_new (change e to ebar) 471 efield_new(kkr)=dtset%red_efieldbar(kkr) 472 473 ! step 3 chang matrix A to A1 474 475 do iir=1,3 476 do jjr=1,3 477 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 478 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 479 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 480 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 481 end do 482 end do 483 484 A(:,:)=A1(:,:) 485 A_new(:,:)=A1(:,:) 486 end if 487 488 end do ! end fo kkr 489 490 491 do iir=1,3 492 do jjr=1,3 493 eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 494 end do 495 end do 496 497 energies%e_elecfield=energies%e_elecfield+eenth 498 499 end if ! berryopt==17 500 501 etotal = energies%e_kinetic+ energies%e_hartree + energies%e_xc + & 502 & energies%e_localpsp + energies%e_nlpsp_vfock - energies%e_fock0 + energies%e_corepsp + & 503 & energies%e_entropy + energies%e_elecfield + energies%e_magfield + & 504 & energies%e_nucdip 505 !etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - & 506 !& energies%e_xcdc + energies%e_corepsp + & 507 !& energies%e_entropy + energies%e_elecfield 508 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 509 if (usepaw==1) then 510 etotal = etotal + energies%e_paw 511 end if 512 513 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 514 !!!!!!!!!!!!!! now, compute mixed densities 515 516 gammp1=etotal-e_ksnm1-fp0 517 if (fp0>0.d0) then 518 write(std_out,*) "fp0 est positif" 519 ! stop 520 end if 521 write(std_out,*) "fp0 ",fp0 522 alphaopt=-fp0/two/gammp1 523 524 if (alphaopt>one.or.alphaopt<0.d0) alphaopt=one 525 if (abs(energies%h0)<=tol10) alphaopt=one 526 write(std_out,*) " alphaopt",alphaopt 527 528 energies%h0=(one-alphaopt)*energies%h0 + alphaopt*(energies%e_kinetic+energies%e_localpsp) 529 energies%h0=energies%h0 + alphaopt*energies%e_nlpsp_vfock 530 531 rhor= rhor+(alphaopt-one)*nvresid 532 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,0) 533 534 if (usepaw==1) then 535 if (my_natom>0) then 536 ABI_MALLOC(rhoijtmp,(pawrhoij(1)%cplex_rhoij*pawrhoij(1)%lmn2_size,pawrhoij(1)%nspden)) 537 end if 538 do iatom=1,my_natom 539 rhoijtmp=zero 540 if (pawrhoij(iatom)%cplex_rhoij==1) then 541 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 542 do ispden=1,pawrhoij(iatom)%nspden 543 do irhoij=1,pawrhoij(iatom)%nrhoijsel 544 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 545 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 546 end do 547 end do 548 end if 549 do ispden=1,pawrhoij(iatom)%nspden 550 do kmix=1,pawrhoij(iatom)%lmnmix_sz 551 klmn=pawrhoij(iatom)%kpawmix(kmix) 552 rhoijtmp(klmn,ispden)=rhoijtmp(klmn,ispden)+(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn,ispden) 553 end do 554 end do 555 else 556 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 557 jrhoij=1 558 do ispden=1,pawrhoij(iatom)%nspden 559 do irhoij=1,pawrhoij(iatom)%nrhoijsel 560 klmn=2*pawrhoij(iatom)%rhoijselect(irhoij)-1 561 rhoijtmp(klmn:klmn+1,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden) 562 jrhoij=jrhoij+2 563 end do 564 end do 565 end if 566 do ispden=1,pawrhoij(iatom)%nspden 567 do kmix=1,pawrhoij(iatom)%lmnmix_sz 568 klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 569 rhoijtmp(klmn:klmn+1,ispden)=rhoijtmp(klmn:klmn+1,ispden) & 570 & +(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 571 end do 572 end do 573 end if 574 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,& 575 & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,pawrhoij(iatom)%lmn2_size,& 576 & pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 577 end do ! iatom 578 if (allocated(rhoijtmp)) then 579 ABI_FREE(rhoijtmp) 580 end if 581 end if ! usepaw 582 583 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 584 !!!!! Calcul des quantites qui dependent de rho_tilde_n+1 (rho apres mixing) 585 586 if (usepaw==1) then 587 if (ider>=0) then 588 izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero 589 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 590 & nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,nhatgr,& 591 & nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,& 592 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 593 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 594 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 595 end if 596 end if 597 598 !------Compute Hartree and xc potentials---------------------------------- 599 600 call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,& 601 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 602 603 !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value 604 if (dtset%xc_tb09_c>99._dp) then 605 call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, & 606 & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, & 607 & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, & 608 & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 609 & computation_type='all') 610 end if 611 612 !Compute xc potential (separate up and down if spin-polarized) 613 optxc=1;if (nkxc>0) optxc=2 614 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 615 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,strsxc,& 616 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 617 618 if (nhatgrdim>0) then 619 ABI_FREE(nhatgr) 620 end if 621 622 !------Compute parts of total energy depending on potentials-------- 623 624 ucvol_local = ucvol 625 626 !Compute Hartree energy energies%e_hartree 627 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 628 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 629 energies%e_hartree=half*energies%e_hartree 630 631 !Compute double-counting XC energy energies%e_xcdc 632 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 633 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 634 635 etotal=energies%h0+energies%e_hartree+energies%e_xc+energies%e_corepsp + & 636 & energies%e_entropy + energies%e_elecfield + energies%e_magfield 637 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 638 if (usepaw==1) then 639 do iatom=1,my_natom 640 itypat=paw_ij(iatom)%itypat 641 ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size)) 642 paw_ij(iatom)%has_dijhartree=1 643 end do 644 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom, & 645 & dtset%nspden,ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang, & 646 & dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,& 647 & dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp,& 648 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 649 do iatom=1,my_natom 650 ABI_FREE(paw_ij(iatom)%dijhartree) 651 paw_ij(iatom)%has_dijhartree=0 652 end do 653 etotal=etotal+energies%e_paw 654 end if 655 !Compute energy residual 656 deltae=etotal-elast 657 elast=etotal 658 659 do ispden=1,min(dtset%nspden,2) 660 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vpsp,vxc) 661 do ifft=1,nfft 662 vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden) 663 end do 664 end do 665 if(dtset%nspden==4) vtrial(:,3:4)=vxc(:,3:4) 666 667 call timab(80,2,tsec) 668 669 !DEBUG 670 !write(std_out,*) 'eeig-ehart+enxc-enxcdc+eew+eii+eent+enefield+epawdc',eeig,ehart,enxc,enxcdc,eew,eii,eent,enefield,epawdc 671 !ENDEBUG 672 673 end subroutine odamix