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