TABLE OF CONTENTS
ABINIT/m_odamix [ Modules ]
NAME
m_odamix
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 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 51 implicit none 52 53 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
175 subroutine odamix(deltae,dtset,elast,energies,etotal,& 176 & gprimd,gsqcut,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 177 & nkxc,ntypat,nvresid,n3xccc,optres,paw_ij,& 178 & paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,& 179 & red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,usepaw,& 180 & usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,& 181 & taur,vxctau,add_tfw) ! optional arguments 182 183 !Arguments ------------------------------------ 184 !scalars 185 integer,intent(in) :: my_natom,n3xccc,nfft,nkxc,ntypat,optres 186 integer,intent(in) :: usepaw,usexcnhat 187 logical,intent(in),optional :: add_tfw 188 real(dp),intent(in) :: gsqcut,ucvol 189 real(dp),intent(inout) :: elast 190 real(dp),intent(out) :: deltae,etotal,vxcavg 191 type(MPI_type),intent(in) :: mpi_enreg 192 type(dataset_type),intent(in) :: dtset 193 type(energies_type),intent(inout) :: energies 194 type(pawang_type),intent(in) :: pawang 195 type(pseudopotential_type),intent(in) :: psps 196 !arrays 197 integer,intent(in) :: ngfft(18) 198 logical :: add_tfw_ 199 real(dp),intent(in) :: gprimd(3,3) 200 real(dp),intent(in) :: red_ptot(3),rprimd(3,3),vpsp(nfft),xred(3,dtset%natom) 201 real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden) 202 real(dp),intent(inout) :: kxc(nfft,nkxc),nhat(nfft,dtset%nspden*usepaw) 203 real(dp),intent(inout) :: nvresid(nfft,dtset%nspden),rhog(2,nfft) 204 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft) 205 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 206 real(dp),intent(inout) :: xccc3d(n3xccc) 207 real(dp),intent(out) :: strsxc(6) 208 real(dp),intent(inout),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 209 type(paw_an_type),intent(inout) :: paw_an(my_natom) 210 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 211 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 212 type(pawrad_type),intent(in) :: pawrad(ntypat) 213 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 214 type(pawtab_type),intent(in) :: pawtab(ntypat) 215 216 !Local variables------------------------------- 217 !scalars 218 integer :: cplex,iatom,ider,idir,ierr,ifft,ipert,irhoij,ispden,itypat,izero,iir,jjr,kkr 219 integer :: jrhoij,klmn,klmn1,kmix,nfftot,nhatgrdim,nzlmopt,nk3xc,option,optxc 220 logical :: nmxc,with_vxctau 221 real(dp) :: alphaopt,compch_fft,compch_sph,doti,e1t10,e_ksnm1,e_xcdc_vxctau 222 real(dp) :: eenth,fp0,gammp1,ro_dlt,ucvol_local 223 character(len=500) :: message 224 type(xcdata_type) :: xcdata 225 !arrays 226 real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3) 227 real(dp) :: gmet(3,3),gprimdlc(3,3),qpt(3),rmet(3,3),tsec(2) 228 real(dp),allocatable :: nhatgr(:,:,:),rhoijtmp(:,:) 229 230 ! ********************************************************************* 231 232 !DEBUG 233 !write(std_out,*)' odamix : enter' 234 !ENDDEBUG 235 236 call timab(80,1,tsec) 237 238 !Check that usekden is not 0 if want to use vxctau 239 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 240 241 !To be adjusted for the call to rhotoxc 242 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 243 nk3xc=1;nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 244 245 !faire un test sur optres=1, usewvl=0, nspden=1,nhatgrdim 246 if(optres/=1)then 247 write(message,'(a,i0,a)')' optres=',optres,', not allowed in oda => stop ' 248 ABI_ERROR(message) 249 end if 250 251 if(dtset%usewvl/=0)then 252 write(message,'(a,i0,a)')' usewvl=',dtset%usewvl,', not allowed in oda => stop ' 253 ABI_ERROR(message) 254 end if 255 256 if(dtset%nspden/=1)then 257 write(message,'(a,i0,a)')' nspden=',dtset%nspden,', not allowed in oda => stop ' 258 ABI_ERROR(message) 259 end if 260 261 if (my_natom>0) then 262 if(paw_ij(1)%has_dijhat==0)then 263 message = ' dijhat variable must be allocated in odamix ! ' 264 ABI_ERROR(message) 265 end if 266 if(paw_ij(1)%cplex_dij==2.or.paw_ij(1)%qphase==2)then 267 message = ' complex dij not allowed in odamix! ' 268 ABI_ERROR(message) 269 end if 270 end if 271 272 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 273 !!!!!!!!!! calculation of f'(0)= Eband_new-EH_old-E_xcdc_old-Ek_old-E_loc_old-E_nonloc_old 274 !!!!!!!!!! save previous energy E(rho_tild_n) 275 276 fp0=energies%e_eigenvalues-energies%h0-two*energies%e_hartree-energies%e_xcdc 277 if (usepaw==1) then 278 do iatom=1,my_natom 279 ABI_CHECK(pawrhoij(iatom)%qphase==1,'ODA mixing not allowed with a Q phase in PAW objects!') 280 itypat=pawrhoij(iatom)%itypat 281 do ispden=1,pawrhoij(iatom)%nspden 282 jrhoij=1 283 do irhoij=1,pawrhoij(iatom)%nrhoijsel 284 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 285 ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 286 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden)) 287 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 288 end do 289 klmn1=1 290 do klmn=1,pawrhoij(iatom)%lmn2_size 291 ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,ispden)*pawtab(itypat)%dltij(klmn) 292 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden)) 293 klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij 294 end do 295 end do 296 if (paw_ij(iatom)%ndij>=2.and.pawrhoij(iatom)%nspden==1) then 297 jrhoij=1 298 do irhoij=1,pawrhoij(iatom)%nrhoijsel 299 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 300 ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,1)*pawtab(itypat)%dltij(klmn) 301 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2)) 302 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 303 end do 304 klmn1=1 305 do klmn=1,pawrhoij(iatom)%lmn2_size 306 ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,1)*pawtab(itypat)%dltij(klmn) 307 e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2)) 308 klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij 309 end do 310 e1t10=half*e1t10 311 end if 312 end do 313 if (mpi_enreg%nproc_atom>1) then 314 call xmpi_sum(e1t10,mpi_enreg%comm_atom,ierr) 315 end if 316 fp0=fp0-e1t10 317 end if 318 e_ksnm1=etotal 319 320 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 321 !!!!! Calculation of quantities that do not depend on rho_n+1 322 323 !PAW: eventually recompute compensation density (and gradients) 324 nhatgrdim=0 325 if (usepaw==1) then 326 ider=-1;if (dtset%xclevel==2.or.usexcnhat==0) ider=0 327 if (dtset%xclevel==2.and.usexcnhat==1) ider=ider+2 328 if (ider>0) then 329 nhatgrdim=1 330 ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3)) 331 end if 332 if (ider>=0) then 333 ider=0;izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero 334 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 335 nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,& 336 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,& 337 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 338 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 339 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 340 end if 341 end if 342 343 !------Compute Hartree and xc potentials---------------------------------- 344 nfftot=PRODUCT(ngfft(1:3)) 345 346 call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,& 347 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 348 349 call xcdata_init(xcdata,dtset=dtset) 350 351 !Compute xc potential (separate up and down if spin-polarized) 352 optxc=1 353 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 354 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,strsxc,& 355 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 356 357 !------Compute parts of total energy depending on potentials-------- 358 359 ucvol_local=ucvol 360 361 !Compute Hartree energy energies%e_hartree 362 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 363 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 364 energies%e_hartree=half*energies%e_hartree 365 366 367 !Compute local psp energy energies%e_localpsp 368 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,& 369 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 370 371 !Compute double-counting XC energy energies%e_xcdc 372 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 373 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 374 if (with_vxctau) then 375 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,& 376 & vxctau(:,:,1),ucvol_local,mpi_comm_sphgrid=mpi_enreg%comm_fft) 377 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 378 end if 379 380 if (usepaw/=0) then 381 nzlmopt=dtset%pawnzlm; option=2 382 do iatom=1,my_natom 383 itypat=paw_ij(iatom)%itypat 384 ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size)) 385 paw_ij(iatom)%has_dijhartree=1 386 end do 387 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom,dtset%nspden,ntypat,& 388 & dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,& 389 & pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,& 390 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 391 do iatom=1,my_natom 392 ABI_FREE(paw_ij(iatom)%dijhartree) 393 paw_ij(iatom)%has_dijhartree=0 394 end do 395 end if 396 397 !When the finite-temperature VG broadening scheme is used, 398 !the total entropy contribution "tsmear*entropy" has a meaning, 399 !and gather the two last terms of Eq.8 of the VG paper 400 !Warning : might have to be changed for fixed moment calculations 401 if(dtset%occopt>=3 .and. dtset%occopt<=8) then 402 energies%e_entropy = - dtset%tsmear * energies%entropy 403 else 404 energies%e_entropy = zero 405 end if 406 !Turn it into an electric enthalpy,refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]] 407 ! the missing volume is added here 408 energies%e_elecfield = zero 409 if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG 410 411 energies%e_elecfield = -dot_product(dtset%red_efieldbar,red_ptot) 412 413 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 414 eenth = zero 415 do iir=1,3 416 do jjr=1,3 417 eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 418 end do 419 end do 420 eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth 421 energies%e_elecfield = energies%e_elecfield + eenth 422 423 end if 424 425 energies%e_magfield = zero 426 !if (dtset%berryopt == 5) then 427 !emag = dot_product(mag_cart,dtset%bfield) 428 !energies%e_magfield = emag 429 !end if 430 431 !HONG Turn it into an internal enthalpy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]], 432 !but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j 433 if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then 434 energies%e_elecfield=zero 435 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 436 do iir=1,3 437 do jjr=1,3 438 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 439 end do 440 end do 441 energies%e_elecfield = ucvol_local/(8.0d0*pi)*energies%e_elecfield 442 end if 443 444 !HONG calculate internal energy and electric enthalpy for mixed BC case. 445 if ( dtset%berryopt == 17 ) then 446 energies%e_elecfield = zero 447 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 448 A(:,:)=(4*pi/ucvol_local)*rmet(:,:) 449 A1(:,:)=A(:,:) 450 A_new(:,:)=A(:,:) 451 efield_new(:)=dtset%red_efield(:) 452 eenth = zero 453 454 do kkr=1,3 455 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 456 457 ! step 1 add -ebar*p 458 eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr) 459 460 ! step 2 chang to e_new (change e to ebar) 461 efield_new(kkr)=dtset%red_efieldbar(kkr) 462 463 ! step 3 chang matrix A to A1 464 465 do iir=1,3 466 do jjr=1,3 467 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 468 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 469 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 470 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 471 end do 472 end do 473 474 A(:,:)=A1(:,:) 475 A_new(:,:)=A1(:,:) 476 end if 477 478 end do ! end fo kkr 479 480 481 do iir=1,3 482 do jjr=1,3 483 eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 484 end do 485 end do 486 487 energies%e_elecfield=energies%e_elecfield+eenth 488 489 end if ! berryopt==17 490 491 etotal = energies%e_kinetic+ energies%e_hartree + energies%e_xc + & 492 & energies%e_localpsp + energies%e_nlpsp_vfock - energies%e_fock0 + energies%e_corepsp + & 493 & energies%e_entropy + energies%e_elecfield + energies%e_magfield + & 494 & energies%e_nucdip 495 !etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - & 496 !& energies%e_xcdc + energies%e_corepsp + & 497 !& energies%e_entropy + energies%e_elecfield 498 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 499 if (usepaw==1) then 500 etotal = etotal + energies%e_paw 501 end if 502 503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 504 !!!!!!!!!!!!!! now, compute mixed densities 505 506 gammp1=etotal-e_ksnm1-fp0 507 if (fp0>0.d0) then 508 write(std_out,*) "fp0 est positif" 509 ! stop 510 end if 511 write(std_out,*) "fp0 ",fp0 512 alphaopt=-fp0/two/gammp1 513 514 if (alphaopt>one.or.alphaopt<0.d0) alphaopt=one 515 if (abs(energies%h0)<=tol10) alphaopt=one 516 write(std_out,*) " alphaopt",alphaopt 517 518 energies%h0=(one-alphaopt)*energies%h0 + alphaopt*(energies%e_kinetic+energies%e_localpsp) 519 energies%h0=energies%h0 + alphaopt*energies%e_nlpsp_vfock 520 521 rhor= rhor+(alphaopt-one)*nvresid 522 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,0) 523 524 if (usepaw==1) then 525 if (my_natom>0) then 526 ABI_MALLOC(rhoijtmp,(pawrhoij(1)%cplex_rhoij*pawrhoij(1)%lmn2_size,pawrhoij(1)%nspden)) 527 end if 528 do iatom=1,my_natom 529 rhoijtmp=zero 530 if (pawrhoij(iatom)%cplex_rhoij==1) then 531 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 532 do ispden=1,pawrhoij(iatom)%nspden 533 do irhoij=1,pawrhoij(iatom)%nrhoijsel 534 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 535 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 536 end do 537 end do 538 end if 539 do ispden=1,pawrhoij(iatom)%nspden 540 do kmix=1,pawrhoij(iatom)%lmnmix_sz 541 klmn=pawrhoij(iatom)%kpawmix(kmix) 542 rhoijtmp(klmn,ispden)=rhoijtmp(klmn,ispden)+(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn,ispden) 543 end do 544 end do 545 else 546 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 547 jrhoij=1 548 do ispden=1,pawrhoij(iatom)%nspden 549 do irhoij=1,pawrhoij(iatom)%nrhoijsel 550 klmn=2*pawrhoij(iatom)%rhoijselect(irhoij)-1 551 rhoijtmp(klmn:klmn+1,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden) 552 jrhoij=jrhoij+2 553 end do 554 end do 555 end if 556 do ispden=1,pawrhoij(iatom)%nspden 557 do kmix=1,pawrhoij(iatom)%lmnmix_sz 558 klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 559 rhoijtmp(klmn:klmn+1,ispden)=rhoijtmp(klmn:klmn+1,ispden) & 560 & +(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 561 end do 562 end do 563 end if 564 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,& 565 & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,pawrhoij(iatom)%lmn2_size,& 566 & pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 567 end do ! iatom 568 if (allocated(rhoijtmp)) then 569 ABI_FREE(rhoijtmp) 570 end if 571 end if ! usepaw 572 573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 574 !!!!! Calcul des quantites qui dependent de rho_tilde_n+1 (rho apres mixing) 575 576 if (usepaw==1) then 577 if (ider>=0) then 578 izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero 579 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 580 & nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,nhatgr,& 581 & nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,& 582 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 583 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 584 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 585 end if 586 end if 587 588 !------Compute Hartree and xc potentials---------------------------------- 589 590 call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,& 591 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 592 593 !Compute xc potential (separate up and down if spin-polarized) 594 optxc=1;if (nkxc>0) optxc=2 595 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 596 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,strsxc,& 597 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 598 599 if (nhatgrdim>0) then 600 ABI_FREE(nhatgr) 601 end if 602 603 !------Compute parts of total energy depending on potentials-------- 604 605 ucvol_local = ucvol 606 607 !Compute Hartree energy energies%e_hartree 608 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 609 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 610 energies%e_hartree=half*energies%e_hartree 611 612 !Compute double-counting XC energy energies%e_xcdc 613 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 614 & mpi_comm_sphgrid=mpi_enreg%comm_fft) 615 616 etotal=energies%h0+energies%e_hartree+energies%e_xc+energies%e_corepsp + & 617 & energies%e_entropy + energies%e_elecfield + energies%e_magfield 618 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 619 if (usepaw==1) then 620 do iatom=1,my_natom 621 itypat=paw_ij(iatom)%itypat 622 ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size)) 623 paw_ij(iatom)%has_dijhartree=1 624 end do 625 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom, & 626 & dtset%nspden,ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang, & 627 & dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,& 628 & dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,& 629 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 630 do iatom=1,my_natom 631 ABI_FREE(paw_ij(iatom)%dijhartree) 632 paw_ij(iatom)%has_dijhartree=0 633 end do 634 etotal=etotal+energies%e_paw 635 end if 636 !Compute energy residual 637 deltae=etotal-elast 638 elast=etotal 639 640 do ispden=1,min(dtset%nspden,2) 641 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vpsp,vxc) 642 do ifft=1,nfft 643 vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden) 644 end do 645 end do 646 if(dtset%nspden==4) vtrial(:,3:4)=vxc(:,3:4) 647 648 call timab(80,2,tsec) 649 650 !DEBUG 651 !write(std_out,*) 'eeig-ehart+enxc-enxcdc+eew+eii+eent+enefield+epawdc',eeig,ehart,enxc,enxcdc,eew,eii,eent,enefield,epawdc 652 !ENDEBUG 653 654 end subroutine odamix