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