1!! Copyright (C) 2016-2019 N. Tancogne-Dejean, X. Andrade 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19 20subroutine X(lda_u_apply)(this, d, mesh, psib, hpsib) 21 type(lda_u_t), intent(in) :: this 22 type(mesh_t), intent(in) :: mesh 23 type(wfs_elec_t), intent(in) :: psib 24 type(wfs_elec_t), intent(inout) :: hpsib 25 type(states_elec_dim_t), intent(in) :: d 26 27 integer :: ibatch, ios, imp, im, ispin, bind1, bind2, inn, ios2 28 R_TYPE, allocatable :: dot(:,:,:,:), reduced(:,:,:), psi(:,:) 29 type(orbitalset_t), pointer :: os 30 type(profile_t), save :: prof 31 integer :: el_per_state 32 R_TYPE :: weight 33 34 call profiling_in(prof, "DFTU_APPLY") 35 36 PUSH_SUB(lda_u_apply) 37 38 SAFE_ALLOCATE(reduced(1:this%maxnorbs,1:psib%nst_linear, 1:this%norbsets)) 39 SAFE_ALLOCATE(dot(1:d%dim,1:this%maxnorbs, 1:this%norbsets, 1:psib%nst)) 40 SAFE_ALLOCATE(psi(1:mesh%np, 1:d%dim)) 41 42 ispin = states_elec_dim_get_spin_index(d, psib%ik) 43 if(d%ispin == UNPOLARIZED) then 44 el_per_state = 2 45 else 46 el_per_state = 1 47 end if 48 49 50 ! We have to compute 51 ! hpsi> += sum_m |phi m> sum_mp Vmmp <phi mp | psi > 52 ! 53 ! We first compute <phi m | psi> for all orbitals of the atom 54 ! 55 do ibatch = 1, psib%nst 56 call batch_get_state(psib, ibatch, mesh%np, psi) 57 do ios = 1, this%norbsets 58 os => this%orbsets(ios) 59 call X(orbitalset_get_coefficients)(os, d%dim, psi, psib%ik, psib%has_phase, this%basisfromstates, & 60 dot(1:d%dim,1:os%norbs,ios,ibatch)) 61 end do 62 end do 63 64 reduced(:,:,:) = R_TOTYPE(M_ZERO) 65 66 do ios = 1, this%norbsets 67 ! 68 os => this%orbsets(ios) 69 do ibatch = 1, psib%nst 70 bind1 = psib%ist_idim_to_linear((/ibatch, 1/)) 71 bind2 = psib%ist_idim_to_linear((/ibatch, 2/)) 72 do im = 1,os%norbs 73 ! sum_mp Vmmp <phi mp | psi > 74 do imp = 1, os%norbs 75 !Note here that V_{mmp} =U/2(delta_{mmp}-2n_{mpm}) 76 if(d%ispin /= SPINORS) then 77 reduced(im, ibatch, ios) = reduced(im, ibatch, ios) + this%X(V)(im, imp, ispin, ios)*dot(1,imp, ios, ibatch) 78 else 79 reduced(im, bind1, ios) = reduced(im, bind1, ios) + this%X(V)(im, imp, 1, ios)*dot(1, imp, ios, ibatch) 80 reduced(im, bind1, ios) = reduced(im, bind1, ios) + this%X(V)(im, imp, 3, ios)*dot(2, imp, ios, ibatch) 81 reduced(im, bind2, ios) = reduced(im, bind2, ios) + this%X(V)(im, imp, 4, ios)*dot(1, imp, ios, ibatch) 82 reduced(im, bind2, ios) = reduced(im, bind2, ios) + this%X(V)(im, imp, 2, ios)*dot(2, imp, ios, ibatch) 83 end if 84 end do 85 end do 86 end do !ibatch 87 88 !We add the intersite interaction 89 if(this%intersite) then 90 !Loop over nearest neighbors 91 do inn = 1, os%nneighbors 92 ios2 = os%map_os(inn) 93 94 if(psib%has_phase) then 95#ifdef R_TCOMPLEX 96 weight = os%phase_shift(inn, psib%ik)*os%V_ij(inn, 0)/el_per_state 97#endif 98 else 99 weight = os%V_ij(inn, 0)/el_per_state 100 end if 101 102 do ibatch = 1, psib%nst 103 bind1 = psib%ist_idim_to_linear((/ibatch, 1/)) 104 bind2 = psib%ist_idim_to_linear((/ibatch, 2/)) 105 do im = 1, os%norbs 106 do imp = 1, this%orbsets(ios2)%norbs 107 if(d%ispin /= SPINORS) then 108 reduced(im, ibatch, ios) = reduced(im, ibatch, ios) - dot(1, imp, ios2, ibatch) & 109 *R_CONJ(this%X(n_ij)(im, imp, ispin, ios, inn))*weight 110 else ! Spinors 111 reduced(im, bind1, ios) = reduced(im, bind1, ios) - dot(1, imp, ios2, ibatch) & 112 *R_CONJ(this%X(n_ij)(im, imp, 1, ios, inn))*weight 113 reduced(im, bind1, ios) = reduced(im, bind1, ios) - dot(2, imp, ios2, ibatch) & 114 *R_CONJ(this%X(n_ij)(im, imp, 4, ios, inn))*weight 115 reduced(im, bind2, ios) = reduced(im, bind2, ios) - dot(1, imp, ios2, ibatch) & 116 *R_CONJ(this%X(n_ij)(im, imp, 3, ios, inn))*weight 117 reduced(im, bind2, ios) = reduced(im, bind2, ios) - dot(2, imp, ios2, ibatch) & 118 *R_CONJ(this%X(n_ij)(im, imp, 2, ios, inn))*weight 119 end if 120 end do !imp 121 end do !im 122 end do !ibatch 123 end do !inn 124 125 end if 126 end do 127 128 !We add the orbitals properly weighted to hpsi 129 do ios = 1, this%norbsets 130 os => this%orbsets(ios) 131 call X(orbitalset_add_to_batch)(os, d%dim, hpsib, this%basisfromstates, reduced(:,:,ios)) 132 end do 133 134 SAFE_DEALLOCATE_A(psi) 135 SAFE_DEALLOCATE_A(dot) 136 SAFE_DEALLOCATE_A(reduced) 137 138 POP_SUB(lda_u_apply) 139 call profiling_out(prof) 140 141end subroutine X(lda_u_apply) 142 143! --------------------------------------------------------- 144!> This routine computes the values of the occupation matrices 145! --------------------------------------------------------- 146subroutine X(update_occ_matrices)(this, namespace, mesh, st, lda_u_energy, phase) 147 type(lda_u_t), intent(inout) :: this 148 type(namespace_t), intent(in) :: namespace 149 type(mesh_t), intent(in) :: mesh 150 type(states_elec_t), intent(in) :: st 151 FLOAT, intent(inout) :: lda_u_energy 152 CMPLX, optional, pointer :: phase(:,:) 153 154 integer :: ios, im, ik, ist, ispin, norbs, idim, inn, im2, ios2 155 R_TYPE, allocatable :: psi(:,:) 156 R_TYPE, allocatable :: dot(:,:,:) 157 FLOAT :: weight 158 R_TYPE :: renorm_weight 159 type(orbitalset_t), pointer :: os, os2 160 type(profile_t), save :: prof 161 integer :: spec_ind 162 163 call profiling_in(prof, "DFTU_OCC_MATRICES") 164 165 PUSH_SUB(update_occ_matrices) 166 167 this%X(n)(1:this%maxnorbs, 1:this%maxnorbs, 1:st%d%nspin, 1:this%norbsets) = R_TOTYPE(M_ZERO) 168 SAFE_ALLOCATE(psi(1:mesh%np, 1:st%d%dim)) 169 SAFE_ALLOCATE(dot(1:st%d%dim, 1:this%maxnorbs, 1:this%norbsets)) 170 if(this%level == DFT_U_ACBN0) then 171 this%X(n_alt)(1:this%maxnorbs, 1:this%maxnorbs, 1:st%d%nspin, 1:this%norbsets) = R_TOTYPE(M_ZERO) 172 if(this%intersite) then 173 this%X(n_ij)(1:this%maxnorbs, 1:this%maxnorbs, 1:st%d%nspin, & 174 1:this%norbsets, 1:this%maxneighbors) = R_TOTYPE(M_ZERO) 175 this%X(n_alt_ij)(1:this%maxnorbs, 1:this%maxnorbs, 1:st%d%nspin, & 176 1:this%norbsets, 1:this%maxneighbors) = R_TOTYPE(M_ZERO) 177 end if 178 end if 179 180 !TODO: use symmetries of the occupation matrices 181 do ik = st%d%kpt%start, st%d%kpt%end 182 ispin = states_elec_dim_get_spin_index(st%d,ik) 183 184 do ist = st%st_start, st%st_end 185 186 weight = st%d%kweights(ik)*st%occ(ist, ik) 187 if(weight < CNST(1.0e-10)) cycle 188 189 call states_elec_get_state(st, mesh, ist, ik, psi ) 190 191 if(present(phase)) then 192#ifdef R_TCOMPLEX 193 call states_elec_set_phase(st%d, psi, phase(:,ik), mesh%np, .false.) 194#endif 195 end if 196 197 do ios = 1, this%norbsets 198 os => this%orbsets(ios) 199 !We first compute the matrix elemets <orb_m |\psi> 200 !taking into account phase correction if needed 201 call X(orbitalset_get_coefficients)(os, st%d%dim, psi, ik, present(phase), & 202 this%basisfromstates, dot(1:st%d%dim,1:os%norbs,ios)) 203 end do !ios 204 205 !We compute the on-site occupation of the site, if needed 206 if(this%level == DFT_U_ACBN0) then 207 this%X(renorm_occ)(:,:,:,ist,ik) = R_TOTYPE(M_ONE)*(M_ONE-this%acbn0_screening) 208 do ios = 1, this%norbsets 209 os => this%orbsets(ios) 210 if(this%basisfromstates) then 211 spec_ind = 1 212 else 213 spec_ind = species_index(os%spec) 214 end if 215 norbs = os%norbs 216 do im = 1, norbs 217 do idim = 1, st%d%dim 218 this%X(renorm_occ)(spec_ind,os%nn,os%ll,ist,ik) = & 219 this%X(renorm_occ)(spec_ind,os%nn,os%ll,ist,ik) & 220 + abs(dot(idim,im,ios))**2*this%acbn0_screening 221 end do 222 end do 223 end do 224 end if 225 226 227 if(st%d%ispin /= SPINORS) then !Collinear case 228 229 !We can compute the (renormalized) occupation matrices 230 do ios = 1, this%norbsets 231 os => this%orbsets(ios) 232 if(this%basisfromstates) then 233 spec_ind = 1 234 else 235 spec_ind = species_index(os%spec) 236 end if 237 norbs = os%norbs 238 do im = 1, norbs 239 this%X(n)(1:norbs, im, ispin, ios) = this%X(n)(1:norbs, im, ispin, ios) & 240 + weight*dot(1, 1:norbs, ios)*R_CONJ(dot(1, im, ios)) 241 !We compute the renomalized occupation matrices 242 if(this%level == DFT_U_ACBN0) then 243 renorm_weight = this%X(renorm_occ)(spec_ind,os%nn,os%ll,ist,ik)*weight 244 this%X(n_alt)(1:norbs,im,ispin,ios) = this%X(n_alt)(1:norbs,im,ispin,ios) & 245 + renorm_weight*dot(1,1:norbs,ios)*R_CONJ(dot(1,im,ios)) 246 247 !Generalized occupation matrices 248 if(this%intersite) then 249 do inn = 1, os%nneighbors 250 ios2 = os%map_os(inn) 251 os2 => this%orbsets(ios2) 252 253 renorm_weight = sqrt(this%X(renorm_occ)(species_index(os%spec), os%nn, os%ll, ist, ik) & 254 *this%X(renorm_occ)(species_index(os2%spec), os2%nn, os2%ll, ist, ik))*weight 255 do im2 = 1, os2%norbs 256 if(present(phase)) then 257 this%X(n_ij)(im, im2, ispin, ios, inn) = this%X(n_ij)(im, im2, ispin, ios, inn) & 258 + weight*dot(1, im, ios)*R_CONJ(dot(1, im2, ios2)*os%phase_shift(inn, ik)) 259 this%X(n_alt_ij)(im, im2, ispin, ios, inn) = this%X(n_alt_ij)(im, im2, ispin, ios, inn) & 260 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(1, im2,ios2)*os%phase_shift(inn,ik)) 261 else 262 this%X(n_ij)(im, im2, ispin, ios, inn) = this%X(n_ij)(im, im2, ispin, ios, inn) & 263 + weight*dot(1, im, ios)*R_CONJ(dot(1, im2, ios2)) 264 this%X(n_alt_ij)(im, im2, ispin, ios, inn) = this%X(n_alt_ij)(im, im2, ispin, ios, inn) & 265 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(1, im2,ios2)) 266 end if 267 end do !im2 268 end do !inn 269 end if 270 end if 271 end do !im 272 end do !ios 273 274 else !Noncollinear case 275 276 !We can compute the (renormalized) occupation matrices 277 do ios = 1, this%norbsets 278 os => this%orbsets(ios) 279 if(this%basisfromstates) then 280 spec_ind = 1 281 else 282 spec_ind = species_index(os%spec) 283 end if 284 285 norbs = os%norbs 286 do im = 1, norbs 287 this%X(n)(1:norbs, im, 1, ios) = this%X(n)(1:norbs, im, 1, ios) & 288 + weight*dot(1, 1:norbs, ios)*R_CONJ(dot(1, im, ios)) 289 this%X(n)(1:norbs, im, 2, ios) = this%X(n)(1:norbs, im, 2, ios) & 290 + weight*dot(2, 1:norbs, ios)*R_CONJ(dot(2, im, ios)) 291 this%X(n)(1:norbs, im, 3, ios) = this%X(n)(1:norbs, im, 3, ios) & 292 + weight*dot(1, 1:norbs, ios)*R_CONJ(dot(2, im, ios)) 293 this%X(n)(1:norbs, im, 4, ios) = this%X(n)(1:norbs, im, 4, ios) & 294 + weight*dot(2, 1:norbs, ios)*R_CONJ(dot(1, im, ios)) 295 !We compute the renomalized occupation matrices 296 if(this%level == DFT_U_ACBN0) then 297 renorm_weight = this%X(renorm_occ)(spec_ind,os%nn,os%ll,ist,ik)*weight 298 this%X(n_alt)(1:norbs,im,1,ios) = this%X(n_alt)(1:norbs,im,1,ios) & 299 + renorm_weight*dot(1,1:norbs,ios)*R_CONJ(dot(1,im,ios)) 300 this%X(n_alt)(1:norbs,im,2,ios) = this%X(n_alt)(1:norbs,im,2,ios) & 301 + renorm_weight*dot(2,1:norbs,ios)*R_CONJ(dot(2,im,ios)) 302 this%X(n_alt)(1:norbs,im,3,ios) = this%X(n_alt)(1:norbs,im,3,ios) & 303 + renorm_weight*dot(1,1:norbs,ios)*R_CONJ(dot(2,im,ios)) 304 this%X(n_alt)(1:norbs,im,4,ios) = this%X(n_alt)(1:norbs,im,4,ios) & 305 + renorm_weight*dot(2,1:norbs,ios)*R_CONJ(dot(1,im,ios)) 306 307 if(this%intersite) then 308 do inn = 1, os%nneighbors 309 ios2 = os%map_os(inn) 310 os2 => this%orbsets(ios2) 311 312 renorm_weight = sqrt(this%X(renorm_occ)(species_index(os%spec), os%nn, os%ll, ist, ik) & 313 *this%X(renorm_occ)(species_index(os2%spec), os2%nn, os2%ll, ist, ik))*weight 314 do im2 = 1, os2%norbs 315 if(present(phase)) then 316 this%X(n_ij)(im, im2, 1, ios, inn) = this%X(n_ij)(im, im2, 1, ios, inn) & 317 + weight*dot(1, im, ios)*R_CONJ(dot(1, im2, ios2)*os%phase_shift(inn, ik)) 318 this%X(n_ij)(im, im2, 2, ios, inn) = this%X(n_ij)(im, im2, 2, ios, inn) & 319 + weight*dot(2, im, ios)*R_CONJ(dot(2, im2, ios2)*os%phase_shift(inn, ik)) 320 this%X(n_ij)(im, im2, 3, ios, inn) = this%X(n_ij)(im, im2, 3, ios, inn) & 321 + weight*dot(1, im, ios)*R_CONJ(dot(2, im2, ios2)*os%phase_shift(inn, ik)) 322 this%X(n_ij)(im, im2, 4, ios, inn) = this%X(n_ij)(im, im2, 4, ios, inn) & 323 + weight*dot(2, im, ios)*R_CONJ(dot(1, im2, ios2)*os%phase_shift(inn, ik)) 324 325 326 this%X(n_alt_ij)(im, im2, 1, ios, inn) = this%X(n_alt_ij)(im, im2, 1, ios, inn) & 327 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(1, im2,ios2)*os%phase_shift(inn,ik)) 328 this%X(n_alt_ij)(im, im2, 2, ios, inn) = this%X(n_alt_ij)(im, im2, 2, ios, inn) & 329 + renorm_weight*dot(2, im, ios)*R_CONJ(dot(2, im2,ios2)*os%phase_shift(inn,ik)) 330 this%X(n_alt_ij)(im, im2, 3, ios, inn) = this%X(n_alt_ij)(im, im2, 3, ios, inn) & 331 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(2, im2,ios2)*os%phase_shift(inn,ik)) 332 this%X(n_alt_ij)(im, im2, 4, ios, inn) = this%X(n_alt_ij)(im, im2, 4, ios, inn) & 333 + renorm_weight*dot(2, im, ios)*R_CONJ(dot(1, im2,ios2)*os%phase_shift(inn,ik)) 334 else 335 this%X(n_ij)(im, im2, 1, ios, inn) = this%X(n_ij)(im, im2, 1, ios, inn) & 336 + weight*dot(1, im, ios)*R_CONJ(dot(1, im2, ios2)) 337 this%X(n_ij)(im, im2, 2, ios, inn) = this%X(n_ij)(im, im2, 2, ios, inn) & 338 + weight*dot(2, im, ios)*R_CONJ(dot(2, im2, ios2)) 339 this%X(n_ij)(im, im2, 3, ios, inn) = this%X(n_ij)(im, im2, 3, ios, inn) & 340 + weight*dot(1, im, ios)*R_CONJ(dot(2, im2, ios2)) 341 this%X(n_ij)(im, im2, 4, ios, inn) = this%X(n_ij)(im, im2, 4, ios, inn) & 342 + weight*dot(2, im, ios)*R_CONJ(dot(1, im2, ios2)) 343 344 345 this%X(n_alt_ij)(im, im2, 1, ios, inn) = this%X(n_alt_ij)(im, im2, 1, ios, inn) & 346 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(1, im2,ios2)) 347 this%X(n_alt_ij)(im, im2, 2, ios, inn) = this%X(n_alt_ij)(im, im2, 2, ios, inn) & 348 + renorm_weight*dot(2, im, ios)*R_CONJ(dot(2, im2,ios2)) 349 this%X(n_alt_ij)(im, im2, 3, ios, inn) = this%X(n_alt_ij)(im, im2, 3, ios, inn) & 350 + renorm_weight*dot(1, im, ios)*R_CONJ(dot(2, im2,ios2)) 351 this%X(n_alt_ij)(im, im2, 4, ios, inn) = this%X(n_alt_ij)(im, im2, 4, ios, inn) & 352 + renorm_weight*dot(2, im, ios)*R_CONJ(dot(1, im2,ios2)) 353 end if 354 end do !im2 355 end do !inn 356 end if 357 end if 358 end do !im 359 end do !ios 360 end if !Spinors 361 end do 362 end do 363 364 SAFE_DEALLOCATE_A(dot) 365 SAFE_DEALLOCATE_A(psi) 366 367#if defined(HAVE_MPI) 368 if(st%parallel_in_states .or. st%d%kpt%parallel) then 369 call comm_allreduce(st%st_kpt_mpi_grp%comm, this%X(n)) 370 if(this%level == DFT_U_ACBN0) then 371 call comm_allreduce(st%st_kpt_mpi_grp%comm, this%X(n_alt)) 372 if(this%intersite) then 373 call comm_allreduce(st%st_kpt_mpi_grp%comm, this%X(n_ij)) 374 call comm_allreduce(st%st_kpt_mpi_grp%comm, this%X(n_alt_ij)) 375 end if 376 end if 377 end if 378#endif 379 380 if(this%level == DFT_U_ACBN0 .and. .not.this%freeze_u) then 381 if(this%nspins > 1 ) then 382 do ios = 1, this%norbsets 383 if(this%orbsets(ios)%ndim == 1) then 384 call X(compute_ACBNO_U)(this, ios, namespace) 385 if(this%intersite) call X(compute_ACBNO_V)(this, ios) 386 else 387 call compute_ACBNO_U_noncollinear(this, ios, namespace) 388 if(this%intersite) then 389 call messages_not_implemented("Intersite interaction with spinors orbitals.", namespace=namespace) 390 end if 391 end if 392 end do 393 else 394 call X(compute_ACBNO_U_restricted)(this) 395 call X(compute_ACBNO_V_restricted)(this) 396 end if 397 end if 398 399 400 call X(compute_dftu_energy)(this, lda_u_energy, st) 401 call X(lda_u_update_potential)(this,st) 402 403 POP_SUB(update_occ_matrices) 404 call profiling_out(prof) 405end subroutine X(update_occ_matrices) 406 407! --------------------------------------------------------- 408!> This routine computes the value of the double counting term in the LDA+U energy 409! --------------------------------------------------------- 410subroutine X(compute_dftu_energy)(this, energy, st) 411 type(lda_u_t), intent(inout) :: this 412 FLOAT, intent(inout) :: energy 413 type(states_elec_t), intent(in) :: st 414 415 integer :: ios, imp, im, ispin, inn 416 type(orbitalset_t), pointer :: os 417 FLOAT :: nsigma 418 419 PUSH_SUB(compute_dftu_energy) 420 421 energy = M_ZERO 422 423 do ios = 1, this%norbsets 424 do ispin = 1, this%nspins 425 !TODO: These are matrix operations, that could be optimized 426 do im = 1, this%orbsets(ios)%norbs 427 do imp = 1, this%orbsets(ios)%norbs 428 energy = energy - M_HALF*this%orbsets(ios)%Ueff & 429 *abs(this%X(n)(im, imp, ispin, ios))**2/st%smear%el_per_state 430 end do 431 if(ispin <= this%spin_channels) & 432 energy = energy + M_HALF*this%orbsets(ios)%Ueff*TOFLOAT(this%X(n)(im, im, ispin, ios)) 433 end do 434 end do 435 end do 436 437 if(this%intersite) then 438 do ios = 1, this%norbsets 439 os => this%orbsets(ios) 440 do inn = 1, os%nneighbors 441 do ispin = 1, this%nspins 442 do im = 1, os%norbs 443 do imp = 1, this%orbsets(os%map_os(inn))%norbs 444 energy = energy - M_HALF*os%V_ij(inn,0)/st%smear%el_per_state & 445 *abs(this%X(n_ij)(im, imp, ispin, ios, inn))**2 446 end do 447 end do 448 end do 449 end do 450 end do 451 end if 452 453 !See for instance PRB 67, 153106 (2003), Eq.(4) 454 if(this%double_couting == DFT_U_AMF) then 455 ASSERT(st%d%ispin /= SPINORS) 456 do ios = 1, this%norbsets 457 do ispin = 1, this%nspins 458 nsigma = M_ZERO 459 do im = 1, this%orbsets(ios)%norbs 460 nsigma = nsigma + R_REAL(this%X(n)(im,im,ispin,ios))/st%smear%el_per_state 461 end do 462 463 do im = 1, this%orbsets(ios)%norbs 464 energy = energy + M_HALF*this%orbsets(ios)%Ueff*nsigma*(M_ONE-nsigma/this%orbsets(ios)%norbs) 465 energy = energy - M_HALF*this%orbsets(ios)%Ueff*R_REAL(this%X(n)(im, im, ispin, ios)) 466 end do 467 end do 468 end do 469 end if 470 471 POP_SUB(compute_dftu_energy) 472end subroutine X(compute_dftu_energy) 473 474 475! --------------------------------------------------------- 476!> This routine computes the potential that, once multiplied 477!> by the projector Pmm' and summed over m and m' for all the atoms 478!> gives the full Hubbard potential 479! --------------------------------------------------------- 480subroutine X(lda_u_update_potential)(this, st) 481 type(lda_u_t), intent(inout) :: this 482 type(states_elec_t), intent(in) :: st 483 484 integer :: ios, im, ispin, norbs 485 type(profile_t), save :: prof 486 FLOAT :: nsigma 487 488 call profiling_in(prof, "DFTU_POTENTIAL") 489 490 PUSH_SUB(lda_u_update_potential) 491 492 this%X(V) = M_ZERO 493 494 do ios = this%orbs_dist%start, this%orbs_dist%end 495 norbs = this%orbsets(ios)%norbs 496 do ispin = 1, this%nspins 497 do im = 1, norbs 498 this%X(V)(1:norbs,im,ispin,ios) = - this%orbsets(ios)%Ueff*this%X(n)(1:norbs,im,ispin,ios)/st%smear%el_per_state 499 ! Only the diagonal part in spin space (for spinors) 500 if(ispin <= this%spin_channels) & 501 this%X(V)(im,im,ispin,ios) = this%X(V)(im,im,ispin,ios) + M_HALF*this%orbsets(ios)%Ueff 502 503 if(this%orbsets(ios)%alpha > CNST(1.0e-6)) then 504 this%X(V)(im,im,ispin,ios) = this%X(V)(im,im,ispin,ios) + this%orbsets(ios)%alpha 505 end if 506 end do 507 end do 508 509 !See for instance PRB 67, 153106 (2003), Eq.(4) 510 if(this%double_couting == DFT_U_AMF) then 511 ASSERT(st%d%ispin /= SPINORS) 512 do ispin = 1, this%nspins 513 nsigma = M_ZERO 514 do im = 1, norbs 515 nsigma = nsigma + R_REAL(this%X(n)(im,im,ispin,ios))/st%smear%el_per_state 516 end do 517 do im = 1, norbs 518 this%X(V)(im,im,ispin,ios) = this%X(V)(im,im,ispin,ios) + this%orbsets(ios)%Ueff & 519 *(nsigma/norbs - M_HALF) 520 end do 521 end do 522 end if 523 end do 524 525 526 if(this%orbs_dist%parallel) then 527 call comm_allreduce(this%orbs_dist%mpi_grp%comm, this%X(V)) 528 end if 529 530 POP_SUB(lda_u_update_potential) 531 call profiling_out(prof) 532end subroutine X(lda_u_update_potential) 533 534! --------------------------------------------------------- 535!> This routine computes the effective U following the expression 536!> given in Agapito et al., Phys. Rev. X 5, 011006 (2015) 537! --------------------------------------------------------- 538subroutine X(compute_ACBNO_U)(this, ios, namespace) 539 type(lda_u_t), intent(inout) :: this 540 integer, intent(in) :: ios 541 type(namespace_t), intent(in) :: namespace 542 543 integer :: im, imp, impp, imppp, ispin1, ispin2, norbs 544 FLOAT :: numU, numJ, denomU, denomJ, tmpU, tmpJ 545 546 PUSH_SUB(compute_ACBNO_U) 547 548 ASSERT(this%orbsets(ios)%ndim == 1) 549 550 norbs = this%orbsets(ios)%norbs 551 numU = M_ZERO 552 numJ = M_ZERO 553 denomU = M_ZERO 554 denomJ = M_ZERO 555 556 if(norbs > 1) then 557 558 do im = 1, norbs 559 do imp = 1,norbs 560 do impp = 1, norbs 561 do imppp = 1, norbs 562 ! We first compute the terms 563 ! sum_{alpha,beta} P^alpha_{mmp}P^beta_{mpp,mppp} 564 ! sum_{alpha} P^alpha_{mmp}P^alpha_{mpp,mppp} 565 tmpU = M_ZERO 566 tmpJ = M_ZERO 567 568 do ispin1 = 1, this%spin_channels 569 do ispin2 = 1, this%spin_channels 570 tmpU = tmpU + R_REAL(this%X(n_alt)(im,imp,ispin1,ios)*this%X(n_alt)(impp,imppp,ispin2,ios)) 571 end do 572 tmpJ = tmpJ + R_REAL(this%X(n_alt)(im,imp,ispin1,ios)*this%X(n_alt)(impp,imppp,ispin1,ios)) 573 end do 574 if(this%nspins>this%spin_channels) then !Spinors 575 tmpJ = tmpJ + R_REAL(this%X(n_alt)(im,imp,3,ios)*this%X(n_alt)(impp,imppp,4,ios)) & 576 +R_REAL(this%X(n_alt)(im,imp,4,ios)*this%X(n_alt)(impp,imppp,3,ios)) 577 end if 578 ! These are the numerator of the ACBN0 U and J 579 numU = numU + tmpU*this%coulomb(im,imp,impp,imppp,ios) 580 numJ = numJ + tmpJ*this%coulomb(im,imppp,impp,imp,ios) 581 end do 582 end do 583 584 ! We compute the term 585 ! sum_{alpha} sum_{m,mp/=m} N^alpha_{m}N^alpha_{mp} 586 tmpJ = M_ZERO 587 tmpU = M_ZERO 588 if(imp/=im) then 589 do ispin1 = 1, this%spin_channels 590 tmpJ = tmpJ + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin1,ios)) 591 tmpU = tmpU + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin1,ios)) 592 end do 593 if(this%nspins>this%spin_channels) then !Spinors 594 tmpJ = tmpJ + R_REAL(this%X(n)(im,im,3,ios)*this%X(n)(imp,imp,4,ios)) & 595 +R_REAL(this%X(n)(im,im,4,ios)*this%X(n)(imp,imp,3,ios)) 596 end if 597 end if 598 599 ! We compute the term 600 ! sum_{alpha,beta} sum_{m,mp} N^alpha_{m}N^beta_{mp} 601 do ispin1 = 1, this%spin_channels 602 do ispin2 = 1, this%spin_channels 603 if(ispin1 /= ispin2) then 604 tmpU = tmpU + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin2,ios)) 605 end if 606 end do 607 end do 608 609 if(this%nspins>this%spin_channels) then !Spinors 610 if(im == imp) then 611 tmpU = tmpU -(R_REAL(this%X(n)(im,im,3,ios)*this%X(n)(im,im,4,ios)) & 612 +R_REAL(this%X(n)(im,im,4,ios)*this%X(n)(im,im,3,ios))) 613 end if 614 end if 615 616 if(this%rot_inv) then 617 !Rotationally invariance term 618 !sum_{alpha} sum_{m,mp/=m} n^alpha_{mmp}n^alpha_{mpm} 619 if(imp/=im) then 620 do ispin1 = 1, this%spin_channels 621 tmpJ = tmpJ + R_REAL(this%X(n)(im,imp,ispin1,ios)*this%X(n)(imp,im,ispin1,ios)) 622 tmpU = tmpU + R_REAL(this%X(n)(im,imp,ispin1,ios)*this%X(n)(imp,im,ispin1,ios)) 623 end do 624 ASSERT(this%nspins==this%spin_channels)!Spinors not yet implemented 625 end if 626 end if 627 628 denomJ = denomJ + tmpJ 629 denomU = denomU + tmpU 630 631 end do 632 end do 633 634 this%orbsets(ios)%Ueff = numU/denomU - numJ/denomJ 635 this%orbsets(ios)%Ubar = numU/denomU 636 this%orbsets(ios)%Jbar = numJ/denomJ 637 638 else !In the case of s orbitals, the expression is different 639 ! sum_{alpha/=beta} P^alpha_{mmp}P^beta_{mpp,mppp} 640 ! sum_{alpha,beta} sum_{m,mp} N^alpha_{m}N^beta_{mp} 641 numU = M_ZERO 642 denomU = M_ZERO 643 do ispin1 = 1, this%spin_channels 644 do ispin2 = 1, this%spin_channels 645 if(ispin1 /= ispin2) then 646 numU = numU + R_REAL(this%X(n_alt)(1,1,ispin1,ios))*R_REAL(this%X(n_alt)(1,1,ispin2,ios)) 647 denomU = denomU + R_REAL(this%X(n)(1,1,ispin1,ios))*R_REAL(this%X(n)(1,1,ispin2,ios)) 648 end if 649 end do 650 end do 651 652 if(this%nspins>this%spin_channels) then !Spinors 653 denomU = denomU - (R_REAL(this%X(n)(1,1,3,ios)*this%X(n)(1,1,4,ios)) & 654 +R_REAL(this%X(n)(1,1,4,ios)*this%X(n)(1,1,3,ios))) 655 end if 656 657 ! We have to be careful in the case of hydrogen atom for instance 658 if(abs(denomU)> CNST(1.0e-3)) then 659 this%orbsets(ios)%Ubar = (numU/denomU)*R_REAL(this%coulomb(1,1,1,1,ios)) 660 else 661 if( abs(numU-denomU) < CNST(1.0e-3)) then 662 this%orbsets(ios)%Ubar = R_REAL(this%coulomb(1,1,1,1,ios)) 663 else 664 this%orbsets(ios)%Ubar = (numU/denomU) 665 write(message(1),'(a,a)')' Small denominator value for the s orbital ', this%orbsets(ios)%Ubar 666 write(message(2),'(a,a)')' to be multiplied by ', this%coulomb(1,1,1,1,ios) 667 call messages_warning(2, namespace=namespace) 668 this%orbsets(ios)%Ubar = this%orbsets(ios)%Ubar*this%coulomb(1,1,1,1,ios) 669 end if 670 end if 671 672 this%orbsets(ios)%Jbar = 0 673 this%orbsets(ios)%Ueff = this%orbsets(ios)%Ubar 674 end if 675 676 POP_SUB(compute_ACBNO_U) 677end subroutine X(compute_ACBNO_U) 678 679 680! --------------------------------------------------------- 681!> This routine computes the effective Uin the spin-unpolarised case 682! --------------------------------------------------------- 683subroutine X(compute_ACBNO_U_restricted)(this) 684 type(lda_u_t), intent(inout) :: this 685 686 integer :: ios, im, imp, impp, imppp, norbs 687 FLOAT :: numU, numJ, denomU, denomJ 688 689 PUSH_SUB(compute_ACBNO_U_restricted) 690 691 do ios = 1, this%norbsets 692 norbs = this%orbsets(ios)%norbs 693 numU = M_ZERO 694 numJ = M_ZERO 695 denomU = M_ZERO 696 denomJ = M_ZERO 697 698 if(norbs > 1) then 699 700 do im = 1, norbs 701 do imp = 1,norbs 702 do impp = 1, norbs 703 do imppp = 1, norbs 704 numU = numU + R_REAL(this%X(n_alt)(im,imp,1,ios)*this%X(n_alt)(impp,imppp,1,ios)) & 705 *this%coulomb(im,imp,impp,imppp,ios) 706 numJ = numJ + R_REAL(this%X(n_alt)(im,imp,1,ios)*this%X(n_alt)(impp,imppp,1,ios)) & 707 *this%coulomb(im,imppp,impp,imp,ios) 708 end do 709 end do 710 ! We compute the term 711 ! sum_{m,mp/=m} N_{m}N_{mp} 712 if(imp/=im) then 713 denomJ = denomJ + R_REAL(this%X(n)(im,im,1,ios))*R_REAL(this%X(n)(imp,imp,1,ios)) 714 denomU = denomU + R_REAL(this%X(n)(im,im,1,ios))*R_REAL(this%X(n)(imp,imp,1,ios)) 715 end if 716 ! We compute the term 717 ! sum_{m,mp} N_{m}N_{mp} 718 denomU = denomU + R_REAL(this%X(n)(im,im,1,ios))*R_REAL(this%X(n)(imp,imp,1,ios)) 719 720 if(this%rot_inv) then 721 !Rotationally invariance term 722 !sum_{m,mp/=m} n^alpha_{mmp}n^alpha_{mpm} 723 if(imp/=im) then 724 denomJ = denomJ + R_REAL(this%X(n)(im,imp,1,ios)*this%X(n)(imp,im,1,ios)) 725 denomU = denomU + R_REAL(this%X(n)(im,imp,1,ios)*this%X(n)(imp,im,1,ios)) 726 end if 727 end if 728 end do 729 end do 730 731 this%orbsets(ios)%Ueff = M_TWO*numU/denomU - numJ/denomJ 732 this%orbsets(ios)%Ubar = M_TWO*numU/denomU 733 this%orbsets(ios)%Jbar = numJ/denomJ 734 735 else !In the case of s orbitals, the expression is different 736 ! P_{mmp}P_{mpp,mppp}(m,mp|mpp,mppp) 737 numU = R_REAL(this%X(n_alt)(1,1,1,ios))**2*this%coulomb(1,1,1,1,ios) 738 739 ! We compute the term 740 ! sum_{alpha,beta} sum_{m,mp} N^alpha_{m}N^beta_{mp} 741 denomU = R_REAL(this%X(n)(1,1,1,ios))**2 742 743 this%orbsets(ios)%Ueff = numU/denomU 744 this%orbsets(ios)%Ubar = numU/denomU 745 this%orbsets(ios)%Jbar = 0 746 end if 747 end do 748 749 POP_SUB(compute_ACBNO_U_restricted) 750end subroutine X(compute_ACBNO_U_restricted) 751 752! --------------------------------------------------------- 753!> This routine computes the effective V in the spin-polarized case 754! --------------------------------------------------------- 755subroutine X(compute_ACBNO_V)(this, ios) 756 type(lda_u_t), intent(inout) :: this 757 integer, intent(in) :: ios 758 759 integer :: ios2, im, imp 760 integer :: inn, norbs, norbs2 761 integer :: ispin1, ispin2 762 FLOAT :: numV, denomV 763 764 if(.not.this%intersite) return 765 766 PUSH_SUB(compute_ACBNO_V) 767 768 norbs = this%orbsets(ios)%norbs 769 do inn = 1, this%orbsets(ios)%nneighbors 770 numV = M_ZERO 771 denomV = M_ZERO 772 773 ios2 = this%orbsets(ios)%map_os(inn) 774 norbs2 = this%orbsets(ios2)%norbs 775 776 do im = 1, norbs 777 do imp= 1, norbs2 778 do ispin1 = 1, this%spin_channels 779 do ispin2 = 1, this%spin_channels 780 numV = numV + R_REAL(this%X(n_alt)(im,im,ispin1,ios))*R_REAL(this%X(n_alt)(imp,imp,ispin2,ios2)) & 781 *this%orbsets(ios)%coulomb_IIJJ(im,im,imp,imp,inn) 782 if(ispin1 == ispin2) then 783 numV = numV - R_REAL(this%X(n_alt_ij)(im,imp,ispin1,ios,inn)*R_CONJ(this%X(n_alt_ij)(im,imp,ispin1,ios,inn)))& 784 *this%orbsets(ios)%coulomb_IIJJ(im,im,imp,imp,inn) 785 end if 786 end do 787 end do 788 if(this%nspins>this%spin_channels) then !Spinors 789 numV = numV -(R_REAL(this%X(n_alt_ij)(im,imp,3,ios,inn)*R_CONJ(this%X(n_alt_ij)(im,imp,3,ios,inn))) & 790 +R_REAL(this%X(n_alt_ij)(im,imp,4,ios,inn)*R_CONJ(this%X(n_alt_ij)(im,imp,4,ios,inn)))) & 791 *this%orbsets(ios)%coulomb_IIJJ(im,im,imp,imp,inn) 792 end if 793 end do 794 end do 795 796 do im = 1, norbs 797 do imp = 1,norbs2 798 ! We compute the term 799 ! sum_{m,mp} ( 2*N_{m}N_{mp} - n_{mmp}n_{mpm}) 800 do ispin1 = 1, this%spin_channels 801 do ispin2 = 1, this%spin_channels 802 denomV = denomV + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin2,ios2)) 803 if(ispin1 == ispin2) denomV = denomV - abs(this%X(n_ij)(im,imp,ispin1,ios,inn))**2 804 end do 805 end do 806 if(this%nspins>this%spin_channels) then !Spinors 807 denomV = denomV - abs(this%X(n_ij)(im,imp,3,ios,inn))**2 - abs(this%X(n_ij)(im,imp,4,ios,inn))**2 808 end if 809 end do 810 end do 811 812 this%orbsets(ios)%V_ij(inn,0) = numV/denomV*M_HALF 813 end do 814 815 POP_SUB(compute_ACBNO_V) 816end subroutine X(compute_ACBNO_V) 817 818 819 820! --------------------------------------------------------- 821!> This routine computes the effective V in the spin-unpolarised case 822! --------------------------------------------------------- 823subroutine X(compute_ACBNO_V_restricted)(this) 824 type(lda_u_t), intent(inout) :: this 825 826 integer :: ios, ios2, im, imp 827 integer :: inn, norbs, norbs2 828 FLOAT :: numV, denomV 829 830 if(.not.this%intersite) return 831 832 PUSH_SUB(compute_ACBNO_V_restricted) 833 834 do ios = 1, this%norbsets 835 norbs = this%orbsets(ios)%norbs 836 do inn = 1, this%orbsets(ios)%nneighbors 837 numV = M_ZERO 838 denomV = M_ZERO 839 840 ios2 = this%orbsets(ios)%map_os(inn) 841 norbs2 = this%orbsets(ios2)%norbs 842 843 do im = 1, norbs 844 do imp = 1, norbs2 845 numV = numV + (M_TWO*R_REAL(this%X(n_alt)(im,im,1,ios)*this%X(n_alt)(imp,imp,1,ios2)) & 846 - R_REAL(this%X(n_alt_ij)(im,imp,1,ios,inn)*R_CONJ(this%X(n_alt_ij)(im,imp,1,ios,inn))))& 847 *this%orbsets(ios)%coulomb_IIJJ(im,im,imp,imp,inn) 848 end do 849 850 do imp = 1, norbs2 851 ! We compute the term 852 ! sum_{m,mp} ( 2*N_{m}N_{mp} - n_{mmp}n_{mpm}) 853 denomV = denomV + M_TWO*R_REAL(this%X(n)(im,im,1,ios))*R_REAL(this%X(n)(imp,imp,1,ios2)) & 854 - abs(this%X(n_ij)(im,imp,1,ios,inn))**2 855 end do 856 end do 857 858 this%orbsets(ios)%V_ij(inn,0) = numV/denomV*M_HALF 859 end do !inn 860 end do !ios 861 862 POP_SUB(compute_ACBNO_V_restricted) 863end subroutine X(compute_ACBNO_V_restricted) 864 865! --------------------------------------------------------- 866!> This routine computes the Kanamori U, Up, and J 867! --------------------------------------------------------- 868subroutine X(compute_ACBNO_U_kanamori)(this, kanamori) 869 type(lda_u_t), intent(in) :: this 870 FLOAT, intent(out) :: kanamori(:,:) 871 872 integer :: ios, im, imp, impp, imppp, norbs 873 integer :: ispin1, ispin2 874 FLOAT :: numU, numUp, numJ, denomU, denomUp, denomJ 875 FLOAT :: tmpU, tmpUp, tmpJ 876 877 PUSH_SUB(compute_ACBNO_U_kanamori) 878 879 ASSERT(this%nspins == this%spin_channels) 880 881 do ios = 1, this%norbsets 882 norbs = this%orbsets(ios)%norbs 883 numU = M_ZERO 884 denomU = M_ZERO 885 numUp = M_ZERO 886 denomUp = M_ZERO 887 numJ = M_ZERO 888 denomJ = M_ZERO 889 890 if(norbs > 1) then 891 892 do im = 1, norbs 893 do imp = 1,norbs 894 do impp = 1, norbs 895 do imppp = 1, norbs 896 tmpU = M_ZERO 897 tmpUp = M_ZERO 898 tmpJ = M_ZERO 899 900 do ispin1 = 1, this%spin_channels 901 do ispin2 = 1, this%spin_channels 902 tmpUp = tmpUp + R_REAL(this%X(n_alt)(im,imp,ispin1,ios)*this%X(n_alt)(impp,imppp,ispin2,ios)) 903 if(ispin1 /= ispin2) then 904 tmpU = tmpU + R_REAL(this%X(n_alt)(im,imp,ispin1,ios)*this%X(n_alt)(impp,imppp,ispin2,ios)) 905 end if 906 end do 907 tmpJ = tmpJ + R_REAL(this%X(n_alt)(im,imp,ispin1,ios)*this%X(n_alt)(impp,imppp,ispin1,ios)) 908 end do 909 910 if(im == imp .and. impp == imppp .and. im == impp) then 911 numU = numU + tmpU*this%coulomb(im,imp,impp,imppp,ios) 912 else 913 numUp = numUp + tmpUp*this%coulomb(im,imp,impp,imppp,ios) 914 numJ = numJ + tmpJ*this%coulomb(im,imppp,impp,imp,ios) 915 end if 916 end do 917 end do 918 919 tmpU = M_ZERO 920 tmpUp = M_ZERO 921 tmpJ = M_ZERO 922 if(imp/=im) then 923 do ispin1 = 1, this%spin_channels 924 tmpUp = tmpUp + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin1,ios)) 925 tmpJ = tmpJ + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin1,ios)) 926 end do 927 end if 928 929 do ispin1 = 1, this%spin_channels 930 do ispin2 = 1, this%spin_channels 931 if(ispin1 /= ispin2) then 932 if(im /= imp) then 933 tmpUp = tmpUp + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin2,ios)) 934 else 935 tmpU = tmpU + R_REAL(this%X(n)(im,im,ispin1,ios))*R_REAL(this%X(n)(imp,imp,ispin2,ios)) 936 end if 937 end if 938 end do 939 end do 940 941 denomU = denomU + tmpU 942 denomUp = denomUp + tmpUp 943 denomJ = denomJ + tmpJ 944 945 end do 946 end do 947 948 kanamori(1,ios) = numU/denomU 949 kanamori(2,ios) = numUp/denomUp 950 kanamori(3,ios) = numJ/denomJ 951 952 else !In the case of s orbitals, the expression is different 953 kanamori(1,ios) = this%orbsets(ios)%Ubar 954 kanamori(2,ios) = M_ZERO 955 kanamori(3,ios) = M_ZERO 956 end if 957 958 959 end do 960 961 POP_SUB(compute_ACBNO_U_kanamori) 962end subroutine X(compute_ACBNO_U_kanamori) 963 964! --------------------------------------------------------- 965!> This routine computes the Kanamori U, Up, and J 966! --------------------------------------------------------- 967subroutine X(compute_ACBNO_U_kanamori_restricted)(this, kanamori) 968 type(lda_u_t), intent(in) :: this 969 FLOAT, intent(out) :: kanamori(3) 970 971 integer :: ios, im, imp, impp, imppp, norbs 972 FLOAT :: numU, numUp, numJ, denomU, denomUp, denomJ 973 FLOAT :: tmpU, tmpUp, tmpJ 974 975 PUSH_SUB(compute_ACBNO_U_kanamori_restricted) 976 977 ASSERT(this%nspins == 1) 978 979 do ios = 1, this%norbsets 980 norbs = this%orbsets(ios)%norbs 981 numU = M_ZERO 982 denomU = M_ZERO 983 numUp = M_ZERO 984 denomUp = M_ZERO 985 numJ = M_ZERO 986 denomJ = M_ZERO 987 988 if(norbs > 1) then 989 990 do im = 1, norbs 991 do imp = 1,norbs 992 do impp = 1, norbs 993 do imppp = 1, norbs 994 tmpU = M_ZERO 995 tmpUp = M_ZERO 996 tmpJ = M_ZERO 997 998 tmpUp = tmpUp + M_TWO*R_REAL(this%X(n_alt)(im,imp,1,ios)*this%X(n_alt)(impp,imppp,1,ios)) 999 tmpU = tmpU + R_REAL(this%X(n_alt)(im,imp,1,ios)*this%X(n_alt)(impp,imppp,1,ios)) 1000 tmpJ = tmpJ + R_REAL(this%X(n_alt)(im,imp,1,ios)*this%X(n_alt)(impp,imppp,1,ios)) 1001 1002 ! These are the numerator of the ACBN0 U and J 1003 if(im == imp .and. impp == imppp .and. im == impp) then 1004 numU = numU + tmpU*this%coulomb(im,imp,impp,imppp,ios) 1005 else 1006 numUp = numUp + tmpUp*this%coulomb(im,imp,impp,imppp,ios) 1007 numJ = numJ + tmpU*this%coulomb(im,imppp,impp,imp,ios) 1008 end if 1009 end do 1010 end do 1011 1012 if(im /= imp) then 1013 denomUp = denomUp + M_TWO*R_REAL(this%X(n)(im,im,1,ios)*this%X(n)(imp,imp,1,ios)) 1014 denomJ = denomJ + R_REAL(this%X(n)(im,im,1,ios))*R_REAL(this%X(n)(imp,imp,1,ios)) 1015 else 1016 denomU = denomU + R_REAL(this%X(n)(im,im,1,ios)*this%X(n)(imp,imp,1,ios)) 1017 end if 1018 1019 end do 1020 end do 1021 1022 kanamori(1) = numU/denomU 1023 kanamori(2) = numUp/denomUp 1024 kanamori(3) = numJ/denomJ 1025 1026 1027 else !In the case of s orbitals, the expression is different 1028 kanamori(1) = this%orbsets(ios)%Ubar 1029 kanamori(2) = M_ZERO 1030 kanamori(3) = M_ZERO 1031 end if 1032 1033 1034 end do 1035 1036 POP_SUB(compute_ACBNO_U_kanamori_restricted) 1037end subroutine X(compute_ACBNO_U_kanamori_restricted) 1038 1039 1040! --------------------------------------------------------- 1041! TODO: Merge this with the two_body routine in system/output_me_inc.F90 1042subroutine X(compute_coulomb_integrals) (this, namespace, mesh, der, psolver) 1043 type(lda_u_t), intent(inout) :: this 1044 type(namespace_t), intent(in) :: namespace 1045 type(mesh_t), intent(in) :: mesh 1046 type(derivatives_t), intent(in) :: der 1047 type(poisson_t), intent(in) :: psolver 1048 1049 integer :: ist, jst, kst, lst, ijst, klst 1050 integer :: norbs, np_sphere, ios, ip 1051 integer :: idone, ntodo 1052 FLOAT, allocatable :: tmp(:), vv(:), nn(:) 1053 type(orbitalset_t), pointer :: os 1054 type(profile_t), save :: prof 1055 1056 call profiling_in(prof, "DFTU_COULOMB_INTEGRALS") 1057 1058 PUSH_SUB(X(compute_coulomb_integrals)) 1059 1060 SAFE_ALLOCATE(nn(1:this%max_np)) 1061 SAFE_ALLOCATE(vv(1:this%max_np)) 1062 SAFE_ALLOCATE(tmp(1:this%max_np)) 1063 1064 SAFE_ALLOCATE(this%coulomb(1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs, 1:this%norbsets)) 1065 this%coulomb(1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%norbsets) = M_ZERO 1066 1067 !Lets counts the number of orbital to treat, to display a progress bar 1068 ntodo = 0 1069 do ios = this%orbs_dist%start, this%orbs_dist%end 1070 norbs = this%orbsets(ios)%norbs 1071 ntodo= ntodo + ((norbs+1)*norbs/2)*((norbs+1)*norbs/2+1)/2 1072 end do 1073 idone = 0 1074 if(mpi_world%rank == 0) call loct_progress_bar(-1, ntodo) 1075 1076 1077 do ios = this%orbs_dist%start, this%orbs_dist%end 1078 os => this%orbsets(ios) 1079 norbs = os%norbs 1080 np_sphere = os%sphere%np 1081 1082 call submesh_build_global(os%sphere) 1083 1084 select case (this%sm_poisson) 1085 case(DFT_U_POISSON_DIRECT) 1086 call poisson_init_sm(os%poisson, namespace, psolver, der, os%sphere, method = POISSON_DIRECT_SUM) 1087 case(DFT_U_POISSON_ISF) 1088 call poisson_init_sm(os%poisson, namespace, psolver, der, os%sphere, method = POISSON_ISF) 1089 end select 1090 1091 ijst=0 1092 do ist = 1, norbs 1093 1094 do jst = 1, norbs 1095 if(jst > ist) cycle 1096 ijst=ijst+1 1097 1098 !$omp parallel do 1099 do ip = 1,np_sphere 1100 nn(ip) = TOFLOAT(os%X(orb)(ip,1,ist))*TOFLOAT(os%X(orb)(ip,1,jst)) 1101 end do 1102 !$omp end parallel do 1103 1104 !Here it is important to use a non-periodic poisson solver, e.g. the direct solver 1105 call dpoisson_solve_sm(os%poisson, os%sphere, vv(1:np_sphere), nn(1:np_sphere)) 1106 1107 klst=0 1108 do kst = 1, norbs 1109 do lst = 1, norbs 1110 if(lst > kst) cycle 1111 klst=klst+1 1112 if(klst > ijst) cycle 1113 1114 !$omp parallel do 1115 do ip = 1,np_sphere 1116 tmp(ip) = vv(ip)*TOFLOAT(os%X(orb)(ip,1,kst))*TOFLOAT(os%X(orb)(ip,1,lst)) 1117 end do 1118 !$omp end parallel do 1119 1120 this%coulomb(ist,jst,kst,lst,ios) = dsm_integrate(mesh, os%sphere, tmp(1:np_sphere), reduce = .false.) 1121 1122 if(abs(this%coulomb(ist,jst,kst,lst,ios))<CNST(1.0e-12)) then 1123 this%coulomb(ist,jst,kst,lst,ios) = M_ZERO 1124 else 1125 this%coulomb(kst,lst,ist,jst,ios) = this%coulomb(ist,jst,kst,lst,ios) 1126 this%coulomb(jst,ist,lst,kst,ios) = this%coulomb(ist,jst,kst,lst,ios) 1127 this%coulomb(lst,kst,jst,ist,ios) = this%coulomb(ist,jst,kst,lst,ios) 1128 this%coulomb(jst,ist,kst,lst,ios) = this%coulomb(ist,jst,kst,lst,ios) 1129 this%coulomb(lst,kst,ist,jst,ios) = this%coulomb(ist,jst,kst,lst,ios) 1130 this%coulomb(ist,jst,lst,kst,ios) = this%coulomb(ist,jst,kst,lst,ios) 1131 this%coulomb(kst,lst,jst,ist,ios) = this%coulomb(ist,jst,kst,lst,ios) 1132 end if 1133 1134 !Update the progress bar 1135 idone = idone + 1 1136 if(mpi_world%rank == 0) call loct_progress_bar(idone, ntodo) 1137 end do !lst 1138 end do !kst 1139 end do !jst 1140 end do !ist 1141 call poisson_end(os%poisson) 1142 call submesh_end_cube_map(os%sphere) 1143 call submesh_end_global(os%sphere) 1144 end do !iorb 1145 1146 if(mesh%parallel_in_domains) then 1147 call comm_allreduce(mesh%mpi_grp%comm, this%coulomb) 1148 end if 1149 1150 if(this%orbs_dist%parallel) then 1151 call comm_allreduce(this%orbs_dist%mpi_grp%comm, this%coulomb) 1152 end if 1153 1154 SAFE_DEALLOCATE_A(nn) 1155 SAFE_DEALLOCATE_A(vv) 1156 SAFE_DEALLOCATE_A(tmp) 1157 1158 POP_SUB(X(compute_coulomb_integrals)) 1159 call profiling_out(prof) 1160end subroutine X(compute_coulomb_integrals) 1161 1162subroutine X(compute_periodic_coulomb_integrals)(this, namespace, der, mc) 1163 type(lda_u_t), intent(inout) :: this 1164 type(namespace_t), intent(in) :: namespace 1165 type(derivatives_t), intent(in) :: der 1166 type(multicomm_t), intent(in) :: mc 1167 1168 integer :: ist, jst, kst, lst, ijst, klst 1169 integer :: norbs, np, ip 1170 integer :: idone, ntodo 1171 FLOAT, allocatable :: tmp(:), vv(:), nn(:) 1172 type(orbitalset_t), pointer :: os 1173 type(profile_t), save :: prof 1174 1175 call profiling_in(prof, "DFTU_PER_COULOMB") 1176 1177 !At the moment the basis is not spin polarized 1178 ASSERT(this%nspins == 1) 1179 1180 PUSH_SUB(X(compute_periodic_coulomb_integrals)) 1181 1182 SAFE_ALLOCATE(nn(1:der%mesh%np)) 1183 SAFE_ALLOCATE(vv(1:der%mesh%np)) 1184 SAFE_ALLOCATE(tmp(1:der%mesh%np)) 1185 1186 SAFE_ALLOCATE(this%coulomb(1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs, 1:this%norbsets)) 1187 this%coulomb(1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%norbsets) = M_ZERO 1188 1189 !Lets counts the number of orbital to treat, to display a progress bar 1190 ntodo = 0 1191 norbs = this%orbsets(1)%norbs 1192 ntodo= ntodo + ((norbs+1)*norbs/2)*((norbs+1)*norbs/2+1)/2 1193 idone = 0 1194 if(mpi_world%rank == 0) call loct_progress_bar(-1, ntodo) 1195 1196 os => this%orbsets(1) 1197 norbs = os%norbs 1198 np = der%mesh%np 1199 1200 call poisson_init(os%poisson, namespace, der, mc, M_ZERO, solver=POISSON_DIRECT_SUM) !POISSON_ISF) 1201 1202 ijst=0 1203 do ist = 1, norbs 1204 do jst = 1, norbs 1205 if(jst > ist) cycle 1206 ijst=ijst+1 1207 1208 !$omp parallel do 1209 do ip = 1,np 1210 nn(ip) = TOFLOAT(os%X(orb)(ip,1,ist))*TOFLOAT(os%X(orb)(ip,1,jst)) 1211 end do 1212 !$omp end parallel do 1213 1214 !Here it is important to use a non-periodic poisson solver, e.g. the direct solver 1215 call dpoisson_solve(os%poisson, vv, nn, all_nodes=.true.) 1216 1217 klst=0 1218 do kst = 1, norbs 1219 do lst = 1, norbs 1220 if(lst > kst) cycle 1221 klst=klst+1 1222 if(klst > ijst) cycle 1223 1224 !$omp parallel do 1225 do ip = 1,np 1226 tmp(ip) = vv(ip)*TOFLOAT(os%X(orb)(ip,1,lst))*TOFLOAT(os%X(orb)(ip,1,kst)) 1227 end do 1228 !$omp end parallel do 1229 1230 this%coulomb(ist,jst,kst,lst,1) = dmf_integrate(der%mesh, tmp) 1231 1232 if(abs(this%coulomb(ist,jst,kst,lst,1))<CNST(1.0e-12)) then 1233 this%coulomb(ist,jst,kst,lst,1) = M_ZERO 1234 else 1235 this%coulomb(kst,lst,ist,jst,1) = this%coulomb(ist,jst,kst,lst,1) 1236 this%coulomb(jst,ist,lst,kst,1) = this%coulomb(ist,jst,kst,lst,1) 1237 this%coulomb(lst,kst,jst,ist,1) = this%coulomb(ist,jst,kst,lst,1) 1238 this%coulomb(jst,ist,kst,lst,1) = this%coulomb(ist,jst,kst,lst,1) 1239 this%coulomb(lst,kst,ist,jst,1) = this%coulomb(ist,jst,kst,lst,1) 1240 this%coulomb(ist,jst,lst,kst,1) = this%coulomb(ist,jst,kst,lst,1) 1241 this%coulomb(kst,lst,jst,ist,1) = this%coulomb(ist,jst,kst,lst,1) 1242 end if 1243 1244 !Update the progress bar 1245 idone = idone + 1 1246 if(mpi_world%rank == 0) call loct_progress_bar(idone, ntodo) 1247 end do !lst 1248 end do !kst 1249 end do !jst 1250 end do !ist 1251 1252 call poisson_end(os%poisson) 1253 1254 SAFE_DEALLOCATE_A(nn) 1255 SAFE_DEALLOCATE_A(vv) 1256 SAFE_DEALLOCATE_A(tmp) 1257 1258 POP_SUB(X(compute_periodic_coulomb_integrals)) 1259 call profiling_out(prof) 1260end subroutine X(compute_periodic_coulomb_integrals) 1261 1262 ! --------------------------------------------------------- 1263 !> This routine computes [r,V_lda+u]. 1264 ! --------------------------------------------------------- 1265 subroutine X(lda_u_commute_r)(this, mesh, d, namespace, ik, psi, gpsi, has_phase) 1266 type(lda_u_t), intent(in) :: this 1267 type(mesh_t), intent(in) :: mesh 1268 type(states_elec_dim_t), intent(in) :: d 1269 type(namespace_t), intent(in) :: namespace 1270 R_TYPE, intent(in) :: psi(:,:) 1271 integer, intent(in) :: ik 1272 R_TYPE, intent(inout) :: gpsi(:, :, :) 1273 logical, intent(in) :: has_phase !True if the wavefunction has an associated phase 1274 1275 integer :: ios, idim, idir, im, imp, is, ispin 1276 integer :: idim_orb, inn, ios2, el_per_state 1277 R_TYPE, allocatable :: dot(:,:,:), epsi(:,:), reduced(:,:,:) 1278 type(orbitalset_t), pointer :: os 1279 type(profile_t), save :: prof 1280 1281 call profiling_in(prof, "DFTU_COMMUTE_R") 1282 1283 PUSH_SUB(lda_u_commute_r) 1284 1285 if(this%double_couting /= DFT_U_FLL) then 1286 call messages_not_implemented("AMF double couting and commutator [r,V_u]", namespace=namespace) 1287 end if 1288 1289 if(this%intersite .and. d%ispin == SPINORS) then 1290 call messages_not_implemented("Intersite interaction, spinors, and commutator [r,V_u]", namespace=namespace) 1291 end if 1292 1293 if((simul_box_is_periodic(mesh%sb) .and. .not. this%basis%submeshforperiodic) & 1294 .or. this%basisfromstates) then 1295 SAFE_ALLOCATE(epsi(1:mesh%np,1:d%dim)) 1296 else 1297 SAFE_ALLOCATE(epsi(1:this%max_np,1:d%dim)) 1298 end if 1299 SAFE_ALLOCATE(dot(1:d%dim,1:this%maxnorbs, 1:this%norbsets)) 1300 SAFE_ALLOCATE(reduced(1:d%dim,1:this%maxnorbs, 1:this%norbsets)) 1301 1302 ispin = states_elec_dim_get_spin_index(d, ik) 1303 if(d%ispin == UNPOLARIZED) then 1304 el_per_state = 2 1305 else 1306 el_per_state = 1 1307 end if 1308 1309 do ios = 1, this%norbsets 1310 ! We have to compute 1311 ! hpsi> += r sum_m |phi m> sum_mp Vmmp <phi mp | psi > 1312 ! 1313 ! We first compute <phi m | psi> for all orbitals of the atom 1314 ! 1315 os => this%orbsets(ios) 1316 ! 1317 call X(orbitalset_get_coefficients)(os, d%dim, psi, ik, has_phase, this%basisfromstates, dot(:,:,ios)) 1318 end do 1319 ! 1320 reduced(:,:,:) = M_ZERO 1321 ! 1322 do ios = 1, this%norbsets 1323 os => this%orbsets(ios) 1324 do im = 1, os%norbs 1325 ! sum_mp Vmmp <phi mp | psi > 1326 do imp = 1, os%norbs 1327 if(d%ispin /= SPINORS) then 1328 reduced(1,im, ios) = reduced(1,im, ios) + this%X(V)(im,imp,ispin,ios)*dot(1,imp, ios) 1329 else 1330 reduced(1,im, ios) = reduced(1,im, ios) + this%X(V)(im,imp,1,ios)*dot(1,imp, ios) 1331 reduced(1,im, ios) = reduced(1,im, ios) + this%X(V)(im,imp,3,ios)*dot(2,imp, ios) 1332 reduced(2,im, ios) = reduced(2,im, ios) + this%X(V)(im,imp,4,ios)*dot(1,imp, ios) 1333 reduced(2,im, ios) = reduced(2,im, ios) + this%X(V)(im,imp,2,ios)*dot(2,imp, ios) 1334 end if 1335 end do 1336 end do 1337 1338 !We add the intersite interaction 1339 if(this%intersite) then 1340 ASSERT(d%ispin /= SPINORS) 1341 1342 !Loop over nearest neighbors 1343 do inn = 1, os%nneighbors 1344 ios2 = os%map_os(inn) 1345 1346 do im = 1,os%norbs 1347 do imp = 1, this%orbsets(ios2)%norbs 1348 if(has_phase) then 1349 reduced(1,im,ios) = reduced(1,im,ios) - dot(1,imp,ios2)*os%phase_shift(inn, ik) & 1350 *this%X(n_ij)(im,imp,ispin,ios,inn)*M_HALF*os%V_ij(inn,0)/el_per_state 1351 reduced(1,imp,ios2) = reduced(1,imp,ios2) - dot(1, im, ios)*R_CONJ(os%phase_shift(inn, ik)) & 1352 *R_CONJ(this%X(n_ij)(im,imp,ispin,ios,inn))*M_HALF*os%V_ij(inn,0)/el_per_state 1353 1354 else 1355 reduced(1,im,ios) = reduced(1,im,ios) - dot(1,imp,ios2) & 1356 *this%X(n_ij)(im,imp,ispin,ios,inn)*M_HALF*os%V_ij(inn,0)/el_per_state 1357 reduced(1,imp,ios2) = reduced(1,imp,ios2) - dot(1,im,ios) & 1358 *R_CONJ(this%X(n_ij)(im,imp,ispin,ios,inn))*M_HALF*os%V_ij(inn,0)/el_per_state 1359 end if 1360 end do !imp 1361 end do !im 1362 end do !inn 1363 end if 1364 end do !ios 1365 1366 1367 do ios = 1, this%norbsets 1368 os => this%orbsets(ios) 1369 do idir = 1, mesh%sb%dim 1370 do idim = 1, d%dim 1371 idim_orb = min(idim,os%ndim) 1372 do im = 1, os%norbs 1373 !In case of phase, we have to apply the conjugate of the phase here 1374 if(has_phase) then 1375#ifdef R_TCOMPLEX 1376 if(simul_box_is_periodic(mesh%sb) .and. .not. this%basis%submeshforperiodic) then 1377 !If we orthogonalize, the orbital is not anymore zorb*phase 1378 if(.not.this%basis%orthogonalization) then 1379 epsi(:,idim) = R_TOTYPE(M_ZERO) 1380 do is = 1, os%sphere%np 1381 epsi(os%sphere%map(is),idim) = epsi(os%sphere%map(is),idim) & 1382 + os%sphere%x(is,idir)*os%zorb(is,idim_orb,im)*os%phase(is,ik) 1383 end do 1384 call lalg_axpy(mesh%np, reduced(idim,im,ios), epsi(1:mesh%np,idim), & 1385 gpsi(1:mesh%np,idir,idim)) 1386 else 1387 !$omp parallel do 1388 do is = 1, mesh%np 1389 epsi(is,idim) = os%sphere%x(is,idir)*os%eorb_mesh(is,im,idim_orb, ik) 1390 end do 1391 call lalg_axpy(mesh%np, reduced(idim,im,ios), epsi(1:mesh%np,idim), & 1392 gpsi(1:mesh%np,idir,idim)) 1393 end if 1394 else 1395 !$omp parallel do 1396 do is = 1, os%sphere%np 1397 epsi(is,idim) = os%sphere%x(is,idir)*os%eorb_submesh(is,idim_orb,im,ik) 1398 end do 1399 !$omp end parallel do 1400 call submesh_add_to_mesh(os%sphere, epsi(1:os%sphere%np,idim), & 1401 gpsi(1:mesh%np,idir,idim), reduced(idim,im,ios)) 1402 end if 1403#endif 1404 else 1405 if(this%basisfromstates) then 1406 !$omp parallel do 1407 do is = 1, mesh%np 1408 epsi(is,idim) = os%sphere%x(is,idir)*os%X(orb)(is,idim_orb,im) 1409 end do 1410 call lalg_axpy(mesh%np, reduced(idim,im,ios), epsi(1:mesh%np,idim), & 1411 gpsi(1:mesh%np,idir,idim)) 1412 else 1413 !$omp parallel do 1414 do is = 1, os%sphere%np 1415 epsi(is,idim) = os%sphere%x(is,idir)*os%X(orb)(is,idim_orb,im) 1416 end do 1417 call submesh_add_to_mesh(os%sphere, epsi(1:os%sphere%np,idim), & 1418 gpsi(1:mesh%np,idir,idim), reduced(idim,im,ios)) 1419 end if 1420 end if 1421 end do !im 1422 end do !idim 1423 end do !idir 1424 end do !ios 1425 1426 do idir = 1, mesh%sb%dim 1427 reduced = M_ZERO 1428 1429 ! We have to compute 1430 ! hpsi> -= sum_m |phi m> sum_mp Vmmp <phi mp| r | psi > 1431 ! 1432 ! We first compute <phi m| r | psi> for all orbitals of the atom 1433 ! 1434 ! 1435 if(simul_box_is_periodic(mesh%sb) .and. .not. this%basis%submeshforperiodic) then 1436 do is = 1, mesh%np 1437 epsi(is,1) = mesh%x(is,idir)*psi(is,1) 1438 end do 1439 else 1440 !$omp parallel do 1441 do is = 1, os%sphere%np 1442 epsi(is,1) = os%sphere%x(is,idir)*psi(os%sphere%map(is),1) 1443 end do 1444 !$omp end parallel do 1445 end if 1446 1447 do ios = 1, this%norbsets 1448 os => this%orbsets(ios) 1449 1450 if(has_phase) then 1451#ifdef R_TCOMPLEX 1452 if(simul_box_is_periodic(mesh%sb) .and. .not. this%basis%submeshforperiodic) then 1453 do im = 1, os%norbs 1454 do idim = 1, d%dim 1455 idim_orb = min(idim,os%ndim) 1456 dot(idim,im,ios) = X(mf_dotp)(mesh, os%eorb_mesh(1:mesh%np,im, idim_orb,ik),& 1457 epsi(1:mesh%np,idim), reduce = .false.) 1458 end do 1459 end do 1460 else 1461 do im = 1, os%norbs 1462 do idim = 1, d%dim 1463 idim_orb = min(idim,os%ndim) 1464 dot(idim,im,ios) = X(mf_dotp)(mesh, os%eorb_submesh(1:os%sphere%np,idim_orb,im,ik),& 1465 epsi(1:os%sphere%np,idim), reduce = .false., np = os%sphere%np) 1466 end do 1467 end do 1468 end if 1469#endif 1470 else 1471 if(this%basisfromstates) then 1472 do idim = 1, d%dim 1473 idim_orb = min(idim,os%ndim) 1474 dot(idim,im,ios) = X(mf_dotp)(mesh, os%X(orb)(1:mesh%np,idim_orb,im),& 1475 epsi(1:mesh%np,idim), reduce=.false.) 1476 end do 1477 else 1478 do im = 1, os%norbs 1479 do idim = 1, d%dim 1480 idim_orb = min(idim,os%ndim) 1481 dot(idim,im,ios) = X(mf_dotp)(mesh, os%X(orb)(1:os%sphere%np,idim_orb,im),& 1482 epsi(1:os%sphere%np,idim), reduce = .false., np = os%sphere%np) 1483 end do 1484 end do 1485 end if 1486 end if 1487 end do !ios 1488 1489 if(mesh%parallel_in_domains) call comm_allreduce(mesh%mpi_grp%comm, dot) 1490 1491 do ios = 1, this%norbsets 1492 os => this%orbsets(ios) 1493 do im = 1, os%norbs 1494 ! sum_mp Vmmp <phi mp|r| psi > 1495 do imp = 1, os%norbs 1496 if(d%ispin /= SPINORS) then 1497 reduced(1,im,ios) = reduced(1,im,ios) - this%X(V)(im,imp,ispin,ios)*dot(1,imp,ios) 1498 else 1499 reduced(1,im,ios) = reduced(1,im,ios) - this%X(V)(im,imp,1,ios)*dot(1,imp,ios) 1500 reduced(1,im,ios) = reduced(1,im,ios) - this%X(V)(im,imp,3,ios)*dot(2,imp,ios) 1501 reduced(2,im,ios) = reduced(2,im,ios) - this%X(V)(im,imp,4,ios)*dot(1,imp,ios) 1502 reduced(2,im,ios) = reduced(2,im,ios) - this%X(V)(im,imp,2,ios)*dot(2,imp,ios) 1503 end if 1504 end do 1505 end do 1506 1507 !We add the intersite interaction 1508 if(this%intersite) then 1509 ASSERT(d%ispin /= SPINORS) 1510 1511 !Loop over nearest neighbors 1512 do inn = 1, os%nneighbors 1513 1514 ios2 = os%map_os(inn) 1515 do im = 1,os%norbs 1516 do imp = 1, this%orbsets(ios2)%norbs 1517 if(has_phase) then 1518 reduced(1,im, ios) = reduced(1,im,ios) - dot(1,imp,ios2)*os%phase_shift(inn, ik) & 1519 *this%X(n_ij)(im,imp,ispin,ios,inn)*os%V_ij(inn,0)/el_per_state 1520 1521 else 1522 reduced(1,im,ios) = reduced(1,im,ios) - dot(1,imp,ios2) & 1523 *this%X(n_ij)(im,imp,ispin,ios,inn)*os%V_ij(inn,0)/el_per_state 1524 end if 1525 end do !imp 1526 end do !im 1527 end do !inn 1528 end if 1529 end do 1530 1531 do ios = 1, this%norbsets 1532 os => this%orbsets(ios) 1533 call X(orbitalset_add_to_psi)(os, d%dim, gpsi(1:mesh%np,idir,1:d%dim), ik, has_phase, & 1534 this%basisfromstates, reduced(1:d%dim,1:os%norbs, ios)) 1535 end do 1536 end do !idir 1537 1538 SAFE_DEALLOCATE_A(epsi) 1539 SAFE_DEALLOCATE_A(dot) 1540 1541 POP_SUB(lda_u_commute_r) 1542 call profiling_out(prof) 1543 end subroutine X(lda_u_commute_r) 1544 1545 subroutine X(lda_u_force)(this, namespace, mesh, st, iq, ndim, psib, grad_psib, force, phase) 1546 type(lda_u_t), intent(in) :: this 1547 type(namespace_t), intent(in) :: namespace 1548 type(mesh_t), intent(in) :: mesh 1549 type(states_elec_t), intent(in) :: st 1550 integer, intent(in) :: iq, ndim 1551 type(wfs_elec_t), intent(in) :: psib 1552 type(wfs_elec_t), intent(in) :: grad_psib(:) 1553 FLOAT, intent(inout) :: force(:, :) 1554 logical, intent(in) :: phase 1555 1556 integer :: ios, iatom, ibatch, ist, im, imp, ispin, idir 1557 type(orbitalset_t), pointer :: os 1558 R_TYPE :: ff(1:ndim) 1559 R_TYPE, allocatable :: psi(:,:), gpsi(:,:) 1560 R_TYPE, allocatable :: dot(:,:), gdot(:,:,:), gradn(:,:,:,:) 1561 FLOAT :: weight 1562 type(profile_t), save :: prof 1563 1564 if(this%level == DFT_U_NONE) return 1565 !In this casem there is no contribution to the force 1566 if(this%basisfromstates) return 1567 1568 if(this%double_couting /= DFT_U_FLL) then 1569 call messages_not_implemented("AMF double couting with forces", namespace=namespace) 1570 end if 1571 1572 PUSH_SUB(X(lda_u_force)) 1573 1574 call profiling_in(prof, "FORCES_DFTU") 1575 1576 !TODO: Implement 1577 if(this%intersite) then 1578 message(1) = "Intersite V forces are not implemented." 1579 call messages_warning(1, namespace=namespace) 1580 end if 1581 1582 SAFE_ALLOCATE(psi(1:mesh%np, 1:st%d%dim)) 1583 SAFE_ALLOCATE(gpsi(1:mesh%np, 1:st%d%dim)) 1584 SAFE_ALLOCATE(dot(1:st%d%dim, 1:this%maxnorbs)) 1585 SAFE_ALLOCATE(gdot(1:st%d%dim, 1:this%maxnorbs,1:ndim)) 1586 SAFE_ALLOCATE(gradn(1:this%maxnorbs,1:this%maxnorbs,1:this%nspins,1:ndim)) 1587 1588 ispin = states_elec_dim_get_spin_index(st%d, iq) 1589 1590 do ios = 1, this%norbsets 1591 os => this%orbsets(ios) 1592 iatom = os%iatom 1593 1594 gradn(1:os%norbs,1:os%norbs,1:this%nspins,1:ndim) = R_TOTYPE(M_ZERO) 1595 1596 do ibatch = 1, psib%nst 1597 ist = psib%ist(ibatch) 1598 weight = st%d%kweights(iq)*st%occ(ist, iq) 1599 if(weight < CNST(1.0e-10)) cycle 1600 1601 call batch_get_state(psib, ibatch, mesh%np, psi) 1602 !No phase here, this is already added 1603 1604 !We first compute the matrix elemets <\psi | orb_m> 1605 !taking into account phase correction if needed 1606 ! 1607 call X(orbitalset_get_coefficients)(os, st%d%dim, psi, iq, phase, this%basisfromstates, dot) 1608 1609 do idir = 1, ndim 1610 call batch_get_state(grad_psib(idir), ibatch, mesh%np, gpsi) 1611 !We first compute the matrix elemets <\psi | orb_m> 1612 !taking into account phase correction if needed 1613 ! 1614 !No phase here, this is already added 1615 1616 call X(orbitalset_get_coefficients)(os, st%d%dim, gpsi, iq, phase, & 1617 this%basisfromstates, gdot(1:st%d%dim,1:os%norbs,idir)) 1618 1619 if(st%d%ispin /= SPINORS) then 1620 do im = 1, os%norbs 1621 gradn(1:os%norbs,im,ispin,idir) = gradn(1:os%norbs,im,ispin,idir) & 1622 + weight*(gdot(1,1:os%norbs,idir)*R_CONJ(dot(1,im)) & 1623 +R_CONJ(gdot(1,im,idir))*dot(1,1:os%norbs)) 1624 end do 1625 else 1626 do im = 1, os%norbs 1627 do ispin = 1, this%spin_channels 1628 gradn(1:os%norbs,im,ispin,idir) = gradn(1:os%norbs,im,ispin,idir) & 1629 + weight*(gdot(ispin,1:os%norbs,idir)*R_CONJ(dot(ispin,im)) & 1630 +R_CONJ(gdot(ispin,im,idir))*dot(ispin,1:os%norbs)) 1631 end do 1632 gradn(1:os%norbs,im,3,idir) = gradn(1:os%norbs,im,3,idir) & 1633 + weight*(gdot(1,1:os%norbs,idir)*R_CONJ(dot(2,im)) & 1634 +R_CONJ(gdot(2,im,idir))*dot(1,1:os%norbs)) 1635 gradn(1:os%norbs,im,4,idir) = gradn(1:os%norbs,im,4,idir) & 1636 + weight*(gdot(2,1:os%norbs,idir)*R_CONJ(dot(1,im)) & 1637 +R_CONJ(gdot(1,im,idir))*dot(2,1:os%norbs)) 1638 end do 1639 1640 end if 1641 end do !idir 1642 1643 end do !ibatch 1644 1645 if(st%d%ispin /= SPINORS) then 1646 ff(1:ndim) = M_ZERO 1647 do im = 1, os%norbs 1648 do imp = 1, os%norbs 1649 ff(1:ndim) = ff(1:ndim) - this%X(n)(imp,im,ispin,ios)/st%smear%el_per_state*gradn(im,imp,ispin,1:ndim) 1650 end do !imp 1651 ff(1:ndim) = ff(1:ndim) + CNST(0.5)*gradn(im, im, ispin,1:ndim) 1652 end do !im 1653 else 1654 ff(1:ndim) = M_ZERO 1655 do ispin = 1, st%d%nspin 1656 do im = 1, os%norbs 1657 do imp = 1, os%norbs 1658 !We use R_CONJ to get n(imp,im, sigmap, sigma) from n(im,imp, sigma,sigmap) 1659 ff(1:ndim) = ff(1:ndim) - R_CONJ(this%X(n)(im,imp,ispin,ios))/st%smear%el_per_state & 1660 *gradn(im,imp,ispin,1:ndim) 1661 end do !imp 1662 if(ispin <= this%spin_channels) & 1663 ff(1:ndim) = ff(1:ndim) + CNST(0.5)*gradn(im, im, ispin,1:ndim) 1664 end do !im 1665 end do !ispin 1666 end if 1667 1668 ! We convert the force to Cartesian coordinates 1669 ! Grad_xyw = Bt Grad_uvw, see Chelikowsky after Eq. 10 1670 if (simul_box_is_periodic(mesh%sb) .and. mesh%sb%nonorthogonal ) then 1671 ff(1:ndim) = matmul(mesh%sb%klattice_primitive(1:ndim, 1:ndim), ff(1:ndim)) 1672 end if 1673 1674 force(1:ndim, iatom) = force(1:ndim, iatom) - os%Ueff*TOFLOAT(ff(1:ndim)) 1675 end do !ios 1676 1677 SAFE_DEALLOCATE_A(psi) 1678 SAFE_DEALLOCATE_A(gpsi) 1679 SAFE_DEALLOCATE_A(dot) 1680 SAFE_DEALLOCATE_A(gdot) 1681 SAFE_DEALLOCATE_A(gradn) 1682 1683 call profiling_out(prof) 1684 1685 POP_SUB(X(lda_u_force)) 1686 end subroutine X(lda_u_force) 1687 1688 ! --------------------------------------------------------- 1689 subroutine X(lda_u_set_occupations)(this, occ) 1690 type(lda_u_t), intent(inout) :: this 1691 R_TYPE, intent(in) :: occ(:) 1692 1693 integer :: ios, ispin, im, imp, ind, norbs, inn, ios2 1694 1695 PUSH_SUB(X(lda_u_set_occupations)) 1696 1697 ind = 0 1698 do ios = 1, this%norbsets 1699 norbs = this%orbsets(ios)%norbs 1700 do ispin = 1, this%nspins 1701 do im = 1, norbs 1702 do imp = 1,norbs 1703 ind = ind + 1 1704 this%X(n)(im,imp,ispin,ios) = occ(ind) 1705 end do 1706 end do 1707 end do 1708 end do 1709 1710 if(this%level == DFT_U_ACBN0) then 1711 do ios = 1, this%norbsets 1712 norbs = this%orbsets(ios)%norbs 1713 do ispin = 1, this%nspins 1714 do im = 1, norbs 1715 do imp = 1,norbs 1716 ind = ind + 1 1717 this%X(n_alt)(im,imp,ispin,ios) = occ(ind) 1718 end do 1719 end do 1720 end do 1721 end do 1722 1723 if(this%intersite) then 1724 do ios = 1, this%norbsets 1725 do inn = 1, this%orbsets(ios)%nneighbors 1726 ios2 = this%orbsets(ios)%map_os(inn) 1727 do ispin = 1, this%nspins 1728 do im = 1, this%orbsets(ios)%norbs 1729 do imp = 1,this%orbsets(ios2)%norbs 1730 ind = ind + 1 1731 this%X(n_ij)(im,imp,ispin,ios,inn) = occ(ind) 1732 ind = ind + 1 1733 this%X(n_alt_ij)(im,imp,ispin,ios,inn) = occ(ind) 1734 end do 1735 end do 1736 end do 1737 end do 1738 end do 1739 end if 1740 end if 1741 1742 POP_SUB(X(lda_u_set_occupations)) 1743 end subroutine X(lda_u_set_occupations) 1744 1745 ! --------------------------------------------------------- 1746 subroutine X(lda_u_get_occupations)(this, occ) 1747 type(lda_u_t), intent(in) :: this 1748 R_TYPE, intent(inout) :: occ(:) 1749 1750 integer :: ios, ispin, im, imp, ind, norbs, inn, ios2 1751 1752 PUSH_SUB(X(lda_u_get_occupations)) 1753 1754 ind = 0 1755 do ios = 1, this%norbsets 1756 norbs = this%orbsets(ios)%norbs 1757 do ispin = 1, this%nspins 1758 do im = 1, norbs 1759 do imp = 1,norbs 1760 ind = ind + 1 1761 occ(ind) = this%X(n)(im,imp,ispin,ios) 1762 end do 1763 end do 1764 end do 1765 end do 1766 1767 if(this%level == DFT_U_ACBN0) then 1768 do ios = 1, this%norbsets 1769 norbs = this%orbsets(ios)%norbs 1770 do ispin = 1, this%nspins 1771 do im = 1, norbs 1772 do imp = 1,norbs 1773 ind = ind + 1 1774 occ(ind) = this%X(n_alt)(im,imp,ispin,ios) 1775 end do 1776 end do 1777 end do 1778 end do 1779 1780 if(this%intersite) then 1781 do ios = 1, this%norbsets 1782 do inn = 1, this%orbsets(ios)%nneighbors 1783 ios2 = this%orbsets(ios)%map_os(inn) 1784 do ispin = 1, this%nspins 1785 do im = 1, this%orbsets(ios)%norbs 1786 do imp = 1,this%orbsets(ios2)%norbs 1787 ind = ind + 1 1788 occ(ind) = this%X(n_ij)(im,imp,ispin,ios,inn) 1789 ind = ind + 1 1790 occ(ind) = this%X(n_alt_ij)(im,imp,ispin,ios,inn) 1791 end do 1792 end do 1793 end do 1794 end do 1795 end do 1796 end if 1797 end if 1798 1799 POP_SUB(X(lda_u_get_occupations)) 1800 end subroutine X(lda_u_get_occupations) 1801 1802 ! --------------------------------------------------------- 1803 subroutine X(lda_u_allocate)(this, st) 1804 type(lda_u_t), intent(inout) :: this 1805 type(states_elec_t), intent(in):: st 1806 1807 integer :: maxorbs, nspin 1808 1809 PUSH_SUB(X(lda_u_allocate)) 1810 1811 maxorbs = this%maxnorbs 1812 nspin = this%nspins 1813 1814 SAFE_ALLOCATE(this%X(n)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets)) 1815 this%X(n)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets) = R_TOTYPE(M_ZERO) 1816 SAFE_ALLOCATE(this%X(V)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets)) 1817 this%X(V)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets) = R_TOTYPE(M_ZERO) 1818 1819 !In case we use the ab-initio scheme, we need to allocate extra resources 1820 if(this%level == DFT_U_ACBN0) then 1821 SAFE_ALLOCATE(this%X(n_alt)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets)) 1822 this%X(n_alt)(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets) = R_TOTYPE(M_ZERO) 1823 SAFE_ALLOCATE(this%X(renorm_occ)(this%nspecies,0:5,0:(MAX_L-1),st%st_start:st%st_end,st%d%kpt%start:st%d%kpt%end)) 1824 this%X(renorm_occ)(this%nspecies,0:5,0:(MAX_L-1),st%st_start:st%st_end,st%d%kpt%start:st%d%kpt%end) = R_TOTYPE(M_ZERO) 1825 end if 1826 1827 POP_SUB(X(lda_u_allocate)) 1828 end subroutine X(lda_u_allocate) 1829