1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8! 9 module m_dftu 10 11 use precision, only: dp 12 13 implicit none 14 15 private 16 17 public :: hubbard_term 18 19 ! Last change in population of DFT+U occupations 20 real(dp), save, public :: DFTU_dPop = huge(1._dp) 21 22 ! These are private variables that is used below 23 integer, save :: DFTU_pop_iter = 0 24 25 ! Maximum number of DFT+U projectors 26 integer, save :: maxdftu = 0 27 real(dp), dimension(:,:,:,:), pointer, save :: occu ! Array used to 28 ! store the 29 ! occupations of 30 ! the DFT+U proj. 31 real(dp), dimension(:,:,:,:), pointer, save :: occu_old ! Same as occu 32 ! but in the 33 ! previous step 34 35 CONTAINS 36 37 subroutine hubbard_term( scell, nua, na, isa, xa, indxua, 38 . maxnh, maxnd, lasto, iphorb, no_u, no_l, 39 . numd, listdptr, listd, numh, 40 . listhptr, listh, nspin, Dscf, 41 . Edftu, DEdftu, Hdftu, 42 . fa, stress, H, iscf, 43 . matrix_elements_only ) 44C ********************************************************************* 45C Calculates Hubbard-like U term contribution to total 46C energy, atomic forces, stress and hamiltonian matrix elements. 47C Energies in Ry. Lengths in Bohr. 48C Merged into the trunk by J. Junquera, March 2016. 49C Written by D. Sanchez-Portal, October 2008, 50C after subroutine nlefsm by J.Soler and P.Ordejon (June 1997). 51C Based on a previous version by S. Riikonen and D. Sanchez-Portal (2005) 52C **************************** INPUT ********************************** 53C real*8 scell(3,3) : Supercell vectors SCELL(IXYZ,IVECT) 54C integer nua : Number of atoms in unit cell 55C integer na : Number of atoms in supercell 56C integer isa(na) : Species index of each atom 57C real*8 xa(3,na) : Atomic positions in cartesian coordinates 58C integer indxua(na) : Index of equivalent atom in unit cell 59C integer maxnh : First dimension of H and listh 60C integer maxnd : Maximum number of elements of the 61C density matrix 62C integer lasto(0:na) : Position of last orbital of each atom 63C integer iphorb(no) : Orbital index of each orbital in its atom, 64C where no=lasto(na) 65C integer no_u : Number of orbitals in unit cell 66C integer no_l : Number of orbitals (local) 67C integer numd(nuo) : Number of nonzero elements of each row of the 68C density matrix 69C integer listdptr(nuo) : Pointer to the start of each row (-1) of the 70C density matrix 71C integer listd(maxnd) : Nonzero hamiltonian-matrix element column 72C indexes for each matrix row 73C integer numh(nuo) : Number of nonzero elements of each row of the 74C hamiltonian matrix 75C integer listhptr(nuo) : Pointer to the start of each row (-1) of the 76C hamiltonian matrix 77C integer listh(maxnh) : Nonzero hamiltonian-matrix element column 78C indexes for each matrix row 79C integer nspin : Number of spin components 80C integer iscf : Counter of the cycles of SCF iterations 81C real*8 Dscf(maxnd,nspin): Density matrix 82C logical matrix_elements_only: Determine whether only the matrix elements 83C of the Hamiltonian are computed, or also the 84C forces and stresses 85C ******************* INPUT and OUTPUT ********************************* 86C real*8 fa(3,na) : NL forces (added to input fa) 87C real*8 stress(3,3) : NL stress (added to input stress) 88C real*8 H(maxnh,nspin) : NL Hamiltonian (added to input H) 89C **************************** OUTPUT ********************************* 90C real*8 Edftu : U-Hubbard energy 1 91C real*8 DEdftu : U-hubbard energy 2 92C real*8 Hdftu(maxnh,nspin): Hamiltonian elements from DFT+U 93C********************************************************************* 94C 95C Modules 96C 97#ifdef MPI 98 use m_mpi_utils, only: globalize_sum 99#endif 100 101 use parallel, only : Node, Nodes 102 use parallelsubs, only : GetNodeOrbs, LocalToGlobalOrb 103 use parallelsubs, only : GlobalToLocalOrb 104 use atmfuncs, only : rcut, orb_gindex, dftu_gindex 105 use atmfuncs, only : nofis 106 use neighbour , only : iana=>jan, r2ki=>r2ij, xki=>xij 107 use neighbour , only : mneighb 108 use alloc, only : re_alloc, de_alloc 109 use m_new_matel, only : new_matel 110 use radial, only : rad_func 111 use atm_types, only : nspecies 112 use atm_types, only : species_info ! Derived type with all the info 113 ! about the radial functions 114 ! (PAOs, KB projectors, 115 ! DFT+U proj, 116 ! VNA potentials, etc) 117 ! for a given atomic specie 118 use atm_types, only : species ! Actual array where the 119 ! previous information is 120 ! stored 121 use dftu_specs, only : dDmax_threshold ! Parameter that defines the 122 ! criterium required to start 123 ! or update the calculation of 124 ! the populations of 125 ! the DFT+U projections 126 use dftu_specs, only : dtol_dftupop ! Parameter that defines the 127 ! convergence criterium 128 ! for the DFT+U local 129 ! population 130 use dftu_specs, only : dftu_init ! Flag that determines whether 131 ! the local populations are 132 ! calculated on the 133 ! first iteration 134 use dftu_specs, only : dftu_shift ! Flag that determines whether 135 ! the U parameter 136 ! is interpreted 137 ! as a local potential shift 138 use m_compute_max_diff, only: dDmax_current 139 140 integer, intent(in) :: 141 . maxnh, na, maxnd, nspin, nua, iscf, no_u, no_l 142 143 integer, intent(in) :: 144 . indxua(na), iphorb(*), isa(na), 145 . lasto(0:na), listd(maxnd), listh(maxnh), 146 . numd(no_l), numh(no_l), listdptr(no_l), listhptr(no_l) 147 real(dp), intent(inout) :: Hdftu(maxnh,nspin) 148 149 real(dp), intent(in) :: scell(3,3), Dscf(maxnd,nspin), 150 . xa(3,na) 151 real(dp), intent(inout) :: fa(3,nua), stress(3,3) 152 real(dp), intent(inout) :: H(maxnh,nspin) 153 real(dp), intent(out) :: Edftu, DEdftu 154 logical, intent(in) :: matrix_elements_only 155 156 real(dp) :: volcel 157 external :: timer, volcel 158 159C Internal variables ................................................ 160C maxno = maximum number of basis orbitals overlapping a KB projector 161 integer, save :: maxno = 500 162 163 integer 164 . ia, ishdftu, ina, ind, ino, 165 . io, iio, ioa, is, ispin, ix, ikb, 166 . j, jno, jo, jx, ka, ko, ks, kua, dftuidx, 167 . ndftuproj, nna, nno, no, nuotot, ja, 168 . lko, lkb, nprin_ko, nprin_kb, kg, ig 169 170 integer, dimension(:), pointer :: iano, iono 171 172 real(dp) 173 . Cijk, fik, rki, rmax, rmaxdftu, rmaxo, 174 . Sik, Sjk, volume, oc(nspin), Ueff, dn, 175 . Dij, rci 176 177 real(dp), dimension(:,:), pointer :: Vi, Di 178 real(dp), dimension(:,:), pointer :: Ski, xno 179 180 real(dp), dimension(:,:,:), pointer :: grSki 181 182#ifdef MPI 183 ! Reduction operations 184 real(dp), dimension(:,:), pointer :: buffer1 => null() 185#endif 186 187 logical :: first ! For a given set of 188 ! atomic positions 189 ! determine whether this 190 ! is the first step in the 191 ! SCF iterations 192 logical :: within, pop_conv 193 logical, dimension(:), pointer :: listed, listedall 194 logical, save :: firstime = .true. ! First time that this 195 ! subroutine is called? 196 type(species_info), pointer :: spp 197 type(rad_func), pointer :: pp 198C ...................... 199 200! Start time counter 201 call timer( 'hubbard_term', 1 ) 202 203 ! Nullify pointers 204 nullify( grSki, Ski, xno, iono, iano ) 205 nullify( listedall, listed, Vi, Di ) 206 207 ! Make sure the energies are zero 208 Edftu = 0.0_dp 209 DEdftu = 0.0_dp 210 211! Determine whether this is the first SCF step for a 212! given atomic configuration 213 first = (iscf == 1) 214 ! Reset population iteration 215 if ( first ) DFTU_pop_iter = 0 216 217! Initialization and allocation of matrices 218 if( firstime ) then 219! Find maximum number of DFT+U projectors on a given atom 220 maxdftu = 0 221 do ka = 1, na 222 is = isa(ka) 223 spp => species(is) 224 maxdftu = max(maxdftu,spp%nprojsdftu) 225 enddo 226 227! Allocate local array to store the occupations of the 228! DFT+U projectors 229 nullify( occu, occu_old ) 230 allocate( occu(maxdftu,maxdftu,nua,nspin) ) 231 call memory( 'A', 'D', size(occu), 'hubbard_term' ) 232 allocate( occu_old(maxdftu,maxdftu,nua,nspin) ) 233 call memory( 'A', 'D', size(occu_old), 'hubbard_term' ) 234 occu_old = 0.0_dp 235 236 firstime = .false. 237 endif 238 239! End initialization 240 241! Here we determine whether the occupations are computed on the first 242! SCF step or not. 243 if( first .and. (.not. dftu_init) ) then 244 if ( Node == 0 ) then 245 write(6,'(2a)') 'hubbard_term: ', 246 & 'not computing occupations in the first SCF step' 247 end if 248 call timer( 'hubbard_term', 2 ) 249 return 250 endif 251 252 253 occupations: if( first .or. 254 & (dDmax_current .lt. dDmax_threshold) .or. 255 & dftu_shift ) then 256 257 258! Find unit cell volume 259 volume = volcel( scell ) * nua / na 260 261! Find maximum range of the atomic orbitals (rmaxo) 262! and of the DFT+U projectors (rmaxdftu) 263 rmaxo = 0.0_dp 264 rmaxdftu = 0.0_dp 265 do is = 1, nspecies 266 267 ! Species orbital range 268 do io = 1, nofis(is) 269 rmaxo = max(rmaxo, rcut(is,io)) 270 enddo 271 272 ! Species DFTU range 273 spp => species(is) 274 do io = 1, spp%n_pjdftunl 275 pp => spp%pjdftu(io) 276 rmaxdftu = max(rmaxdftu, pp%cutoff) 277 end do 278 enddo 279 280 ! Total range of the projector is Oo 281 rmax = rmaxo + rmaxdftu 282 283! Initialize arrays Di and Vi only once 284 no = lasto(na) 285 nuotot = lasto(nua) 286 287! Allocate local memory 288 call re_alloc( Di, 1, no, 1, nspin, 289 & name='Di', routine='hubbard_term' ) 290 Di = 0.0_dp 291 292 call re_alloc( Vi, 1, no, 1, nspin, 293 & name='Vi', routine='hubbard_term' ) 294 Vi = 0.0_dp 295 296 call re_alloc( listed, 1, no, 297 & name='listed', routine='hubbard_term') 298 listed(1:no) = .false. 299 300 call re_alloc( listedall, 1, no, 301 & name='listedall', routine='hubbard_term' ) 302 listedall(1:no) = .false. 303 304! Make list of all orbitals needed for this node 305 do io = 1, no_l 306 call LocalToGlobalOrb(io,Node,Nodes,iio) 307 listedall(iio) = .true. 308 do j = 1, numh(io) 309 jo = listh(listhptr(io)+j) 310 listedall(jo) = .true. 311 enddo 312 enddo 313 314! Allocate local arrays that depend on saved parameters 315 call re_alloc( iano, 1, maxno, 316 & name='iano', routine='hubbard_term' ) 317 call re_alloc( iono, 1, maxno, 318 & name='iono', routine='hubbard_term' ) 319 call re_alloc( xno, 1, 3, 1, maxno, 320 & name='xno', routine='hubbard_term' ) 321 call re_alloc( Ski, 1, maxdftu, 1, maxno, 322 & name='Ski', routine='hubbard_term' ) 323 call re_alloc( grSki, 1, 3, 1, maxdftu, 1, maxno, 324 & name='grSki', routine='hubbard_term' ) 325 326! Counter for the SCF loops to converge the population of the 327! DFT+U projectors 328 DFTU_pop_iter = DFTU_pop_iter + 1 329 if( Node == 0 ) write(6,'(a,i4)') 330 & 'hubbard_term: recalculating local occupations ', DFTU_pop_iter 331 332! Initialize occupations 333 occu = 0.0_dp 334 335! Initialize neighb subroutine 336 call mneighb( scell, rmax, na, xa, 0, 0, nna ) 337 338! Loop on atoms with DFT+U projectors 339 do ka = 1, na 340 kua = indxua(ka) 341 ks = isa(ka) 342 spp => species(ks) 343 ndftuproj = spp%nprojsdftu 344 if( ndftuproj == 0 ) cycle 345 346! Find neighbour atoms 347 call mneighb( scell, rmax, na, xa, ka, 0, nna ) 348 349! Find neighbour orbitals 350 nno = 0 351 do ina = 1, nna 352 ia = iana(ina) 353 is = isa(ia) 354 rki = sqrt(r2ki(ina)) 355 do io = lasto(ia-1)+1, lasto(ia) 356 357! Only calculate if needed locally 358 if (listedall(io)) then 359 ioa = iphorb(io) 360 ig = orb_gindex(is,ioa) 361 rci = rcut(is,ioa) 362 363! Find if orbital is within range 364 within = .false. 365 do ko = 1, ndftuproj 366 dftuidx = spp%pjdftu_index(ko) 367 pp => spp%pjdftu(dftuidx) 368 if ( rci + pp%cutoff > rki ) 369 & within = .true. 370 end do 371 372! Find overlap between neighbour orbitals and DFT+U projectors 373 if (within) then 374! Check maxno - if too small then increase array sizes 375 if (nno.eq.maxno) then 376 maxno = maxno + 10 377 call re_alloc( iano, 1, maxno, name='iano', 378 & copy=.true., routine='hubbard_term' ) 379 call re_alloc( iono, 1, maxno, name='iono', 380 & copy=.true., routine='hubbard_term' ) 381 call re_alloc( xno, 1, 3, 1, maxno, name='xno', 382 & copy=.true., routine='hubbard_term' ) 383 call re_alloc( Ski,1, maxdftu, 1, maxno, name='Ski', 384 & copy=.true., routine='hubbard_term' ) 385 call re_alloc( grSki, 1, 3, 1, maxdftu, 1, maxno, 386 & name='grSki', routine='hubbard_term', 387 & copy=.true. ) 388 endif 389 nno = nno + 1 390! The nno-eme neighbour orbital of the DFT+U projector is io, 391! where io runs between 1 and the total number of orbitals 392! in the supercell 393 iono(nno) = io 394 395! The nno-eme neighbour orbital of the DFT+U projector belongs 396! to atom ia, 397! where ia runs between 1 and the total number of atoms 398! in the supercell 399 iano(nno) = ia 400 401! The relative position between the center of the DFT+U proj. 402! and the center of the nno-eme neighbour orbital 403! is xno 404 do ix = 1,3 405 xno(ix,nno) = xki(ix,ina) 406 enddo 407 408! Here we compute the overlap between a DFT+U projector 409! and a neighbour atomic orbital 410 do ko = 1, ndftuproj 411 kg = dftu_gindex(ks,ko) 412 call new_matel( 'S', kg, ig, xki(1:3,ina), 413 & Ski(ko,nno), grSki(1:3,ko,nno) ) 414 enddo 415 416 endif ! If on orbitals within range 417 418 endif ! If on orbitals within the local node 419 420 enddo ! End loop on neighbour orbitals 421 422 enddo ! End loop on neighbour atoms 423 424! Loop on neighbour orbitals 425 do ino = 1,nno 426 io = iono(ino) 427 ia = iano(ino) 428 429 ! Only atoms in the unit cell 430 if ( ia > nua ) cycle 431 432 ! Only local orbitals 433 call GlobalToLocalOrb(io,Node,Nodes,iio) 434 if ( iio < 1 ) cycle 435 436! Scatter filter and density matrix row of orbital io 437 do j = 1, numd(iio) 438 ind = listdptr(iio) + j 439 jo = listd(ind) 440 listed(jo) = .true. 441 do ispin = 1, nspin 442 Di(jo,ispin) = Di(jo,ispin) + Dscf(ind,ispin) 443 enddo 444 enddo 445 446! Find matrix elements with other neighbour orbitals 447 do jno = 1,nno 448 jo = iono(jno) 449 ja = iano(jno) 450 if ( listed(jo) ) then 451 452! Loop on DFT+U projectors 453 do ko = 1, ndftuproj 454 Sik = Ski(ko,ino) 455 do ikb = 1, ndftuproj 456 Sjk = Ski(ikb,jno) 457 do ispin = 1, nspin 458 Dij = Di(jo,ispin) 459 occu(ko,ikb,kua,ispin) = 460 & occu(ko,ikb,kua,ispin) + 461 & Dij * Sik * Sjk/(3.0_dp-dble(nspin)) 462 enddo 463 enddo 464 enddo 465 endif 466 enddo 467 468! Restore Di and listed 469 do j = 1, numh(iio) 470 jo = listh(listhptr(iio)+j) 471 listed(jo) = .false. 472 do ispin=1,nspin 473 Di(jo,ispin) = 0.0_dp 474 enddo 475 enddo 476 477 enddo ! End loop on neighbour orbitals (ino) 478 479 enddo ! End of the loop on atoms with DFT+U projectors 480 481#ifdef MPI 482! Global reduction of occupation 483 call re_alloc(buffer1, 1, maxdftu, 1, maxdftu, 484 & name='buffer1', routine = 'hubbard_term') 485 do ka=1,nua 486 is=isa(ka) 487 spp => species(is) 488 ndftuproj = spp%nprojsdftu 489 if( ndftuproj == 0 ) cycle 490 do ispin=1,nspin 491 call globalize_sum(occu(1:ndftuproj,1:ndftuproj,ka,ispin), 492 & buffer1(1:ndftuproj,1:ndftuproj)) 493 occu(1:ndftuproj,1:ndftuproj,ka,ispin) = 494 & buffer1(1:ndftuproj,1:ndftuproj) 495 end do 496 enddo 497 call de_alloc(buffer1, name='buffer1', routine = 'hubbard_term') 498#endif 499 500 DFTU_dPop = 0._dp 501 do ka=1,nua 502 is=isa(ka) 503 spp => species(is) 504 ndftuproj = spp%nprojsdftu 505 if ( ndftuproj == 0 ) cycle 506 507 oc = 0.0_dp 508 509 if( dftu_shift .and. Node == 0 ) then 510 write(6,*) 'hubbard_term: projector occupations' 511 write(6,*) 'hubbard_term: atom, species: ',ka, is 512 endif 513 do ko = 1, ndftuproj 514 do ikb = 1, ndftuproj 515 if( dftu_shift .and. Node == 0 ) then 516 write(6,'(2i4,2f12.5)') ko, ikb, 517 & (occu(ko,ikb,ka,ispin),ispin=1,nspin) 518 endif 519 520 do ispin = 1, nspin 521 dn = occu(ko,ikb,ka,ispin)-occu_old(ko,ikb,ka,ispin) 522 DFTU_dPop = max(DFTU_dPop,dabs(dn)) 523 enddo 524 enddo 525 if ( dftu_shift ) then 526 do ispin = 1 , nspin 527 oc(ispin) = oc(ispin) + occu(ko,ko,ka,ispin) 528 end do 529 endif 530 enddo 531 if ( dftu_shift .and. Node == 0 ) then 532 write(6,'(a,/,a,3f12.6)') 533 & 'hubbard_term: Total projector shell', 534 & 'Occupations: ', (oc(ispin),ispin=1,nspin) 535 & ,sum(oc) 536 end if 537 enddo 538 539 540#ifdef MPI 541 ! DFTU_dPop need not be globalized. 542 ! The occupations are already globalized 543#endif 544 if ( Node == 0 ) then 545 write(6,'(a,f12.6)') 546 & 'hubbard_term: maximum change in local occup.',DFTU_dPop 547 endif 548 549 pop_conv = ( DFTU_dPop < dtol_dftupop ) 550 if( .not. pop_conv ) occu_old = occu 551 552 recalc_hamilt: if( .not. pop_conv .or. first .or. 553 & .not. matrix_elements_only ) then 554 if ( Node == 0 ) then 555 if ( matrix_elements_only ) then 556 write(6,'(a)')'hubbard_term: recalculating Hamiltonian' 557 else 558 write(6,'(a)') 559 & 'hubbard_term: recalculating Hamiltonian and forces' 560 endif 561 endif 562 563 Hdftu = 0.0_dp 564 565! Initialize neighb subroutine 566 call mneighb( scell, rmax, na, xa, 0, 0, nna ) 567 568 do ka = 1, na 569 kua = indxua(ka) 570 ks = isa(ka) 571 spp => species(ks) 572 ndftuproj = spp%nprojsdftu 573 if( ndftuproj == 0 ) cycle 574 575! Find neighbour atoms 576 call mneighb( scell, rmax, na, xa, ka, 0, nna ) 577 578! Find neighbour orbitals 579 nno = 0 580 do ina = 1, nna 581 ia = iana(ina) 582 is = isa(ia) 583 rki = sqrt(r2ki(ina)) 584 do io = lasto(ia-1)+1, lasto(ia) 585 586! Only calculate if needed locally 587 if (listedall(io)) then 588 ioa = iphorb(io) 589 ig = orb_gindex(is,ioa) 590 rci = rcut(is,ioa) 591 592! Find if orbital is within range 593 within = .false. 594 do ko = 1, ndftuproj 595 dftuidx = spp%pjdftu_index(ko) 596 pp => spp%pjdftu(dftuidx) 597 if ( rci + pp%cutoff > rki ) 598 & within = .true. 599 end do 600 601! Find overlap between neighbour orbitals and DFT+U projectors 602 if (within) then 603! Check maxno - if too small then increase array sizes 604 if (nno.eq.maxno) then 605 maxno = maxno + 10 606 call re_alloc( iano, 1, maxno, name='iano', 607 & copy=.true., routine='hubbard_term' ) 608 call re_alloc( iono, 1, maxno, name='iono', 609 & copy=.true., routine='hubbard_term' ) 610 call re_alloc( xno, 1, 3, 1, maxno, name='xno', 611 & copy=.true., routine='hubbard_term' ) 612 call re_alloc( Ski,1, maxdftu, 1,maxno,name='Ski', 613 & copy=.true., routine='hubbard_term' ) 614 call re_alloc( grSki, 1, 3, 1, maxdftu, 1, maxno, 615 & name='grSki', routine='hubbard_term', 616 & copy=.true. ) 617 endif 618 nno = nno + 1 619 iono(nno) = io 620 iano(nno) = ia 621 do ix = 1,3 622 xno(ix,nno) = xki(ix,ina) 623 enddo 624 do ko = 1, ndftuproj 625 kg = dftu_gindex(ks,ko) 626 call new_matel( 'S', kg, ig, xki(1:3,ina), 627 & Ski(ko,nno), grSki(1:3,ko,nno) ) 628 enddo 629 630 endif ! If on orbitals within range 631 632 endif ! If on orbitals within the local node 633 634 enddo ! end loop on neighbour orbitals 635 636 enddo ! end loop on neighbour atoms 637 638 639! Loop on neighbour orbitals 640 do ino = 1,nno 641 io = iono(ino) 642 ia = iano(ino) 643 644 ! Only atoms in the unit cell 645 if ( ia > nua) cycle 646 647 ! Only local orbitals 648 call GlobalToLocalOrb(io,Node,Nodes,iio) 649 if ( iio < 1 ) cycle 650 651! Scatter filter/density matrix row of orbital io 652 do j = 1,numd(iio) 653 ind = listdptr(iio)+j 654 jo = listd(ind) 655 listed(jo) = .true. 656 do ispin = 1,nspin 657 Di(jo,ispin) = Di(jo,ispin) + Dscf(ind,ispin) 658 enddo 659 enddo 660 661! Find matrix elements with other neighbour orbitals 662 do jno = 1,nno 663 jo = iono(jno) 664 ja = iano(jno) 665 666 if ( listed(jo) ) then 667! Loop on DFT+U projectors 668 do ko = 1, ndftuproj 669 lko = spp%pjdftu_l(ko) 670 nprin_ko = spp%pjdftu_n(ko) 671 dftuidx = spp%pjdftu_index(ko) 672 Ueff = spp%pjdftunl_U(dftuidx) - 673 & spp%pjdftunl_J(dftuidx) 674 if( dftu_shift ) Ueff = 2.0_dp * Ueff 675 Sik = Ski(ko,ino) 676 Sjk = Ski(ko,jno) 677 do ispin=1, nspin 678 Vi(jo,ispin) = Vi(jo,ispin) + 679 & 0.5_dp * Sik * Sjk * Ueff 680 Dij = Di(jo,ispin) 681 Cijk = Ueff * Dij * Sjk 682 if(.not. matrix_elements_only) then 683 do ix=1,3 684 fik = Cijk * grSki(ix,ko,ino) 685 fa(ix,ia) = fa(ix,ia) - fik 686 fa(ix,kua) = fa(ix,kua) + fik 687 do jx = 1, 3 688 stress(jx,ix) = stress(jx,ix) + 689 & xno(jx,ino)*fik/volume 690 enddo 691 enddo 692 endif 693 694 if( .not. dftu_shift ) then 695 do ikb = 1, ndftuproj 696 lkb = spp%pjdftu_l(ikb) 697 nprin_kb = spp%pjdftu_n(ikb) 698 if( lko .eq. lkb .and. 699 & nprin_ko .eq. nprin_kb ) then 700! For the time being we will use the formulation 701! of Dudarev and collaborators 702 Cijk = occu(ko,ikb,kua,ispin) * 703 & Ueff * Ski(ikb,jno) 704 Vi(jo,ispin) = Vi(jo,ispin) - 705 & Sik * Cijk 706 if(.not. matrix_elements_only) then 707 do ix=1,3 708 fik = -2.0_dp * Dij * Cijk * 709 & grSki(ix,ko,ino) 710 fa(ix,ia) = fa(ix,ia) - fik 711 fa(ix,kua) = fa(ix,kua) + fik 712 do jx = 1, 3 713 stress(jx,ix)= stress(jx,ix) + 714 & xno(jx,ino) * fik / volume 715 enddo 716 enddo 717 endif 718 endif 719 enddo 720 endif 721 722 enddo ! Enddo ispin 723 enddo ! Enddo ko, on DFT+U projectors 724 endif ! Endif if the orbital jo is listed 725 enddo ! Enddo on neighbour orbitals jno 726 727! Add to Hdftu and restore Di, Vi and listed 728 do j = 1,numh(iio) 729 ind = listhptr(iio)+j 730 jo = listh(ind) 731 do ispin = 1,nspin 732 Hdftu(ind,ispin) = Hdftu(ind,ispin) + 733 & Vi(jo,ispin) 734 Vi(jo,ispin) = 0.0_dp 735 Di(jo,ispin) = 0.0_dp 736 enddo 737 listed(jo) = .false. 738 enddo 739 740 enddo ! enddo ino = 1,nno 741 enddo ! enddo loop on orbitals with DFT+U projectors(ka=1,na) 742 743 endif recalc_hamilt 744 745 endif occupations 746 747 add_hamilt: if ( DFTU_pop_iter > 0 ) then 748 749 ! Add the Hamiltonian elements and calculate energy 750 do iio = 1 , no_l 751 do j = 1, numh(iio) 752 ind = listhptr(iio) + j 753 do ispin = 1, nspin 754 Edftu = Edftu + 755 & 0.5_dp * Hdftu(ind,ispin) * Dscf(ind,ispin) 756 H(ind,ispin) = H(ind,ispin) + Hdftu(ind,ispin) 757 end do 758 end do 759 end do 760 761 do ia = 1,nua 762 is = isa(ia) 763 spp => species(is) 764 ndftuproj = spp%nprojsdftu 765 if( ndftuproj .gt. 0 ) then 766 do ispin = 1, nspin 767 do ko = 1, ndftuproj 768 dftuidx = spp%pjdftu_index(ko) 769 Ueff = spp%pjdftunl_U(dftuidx) - 770 & spp%pjdftunl_J(dftuidx) 771 if( dftu_shift ) Ueff = 2.0_dp * Ueff 772 DEdftu = DEdftu + 773 & 0.25_dp * (3.0_dp-dble(nspin)) * Ueff * 774 & occu(ko,ko,ia,ispin) 775 enddo 776 enddo 777 endif 778 enddo 779 780 endif add_hamilt 781 782! Deallocate local memory 783 784 call de_alloc( grSki, name='grSki' ) 785 call de_alloc( Ski, name='Ski' ) 786 call de_alloc( xno, name='xno' ) 787 call de_alloc( iono, name='iono' ) 788 call de_alloc( iano, name='iano' ) 789 790 call de_alloc( listedall, name='listedall' ) 791 call de_alloc( listed, name='listed' ) 792 call de_alloc( Vi, name='Vi' ) 793 call de_alloc( Di, name='Di' ) 794 795! Stop time counter 796 call timer( 'hubbard_term', 2 ) 797 798 end subroutine hubbard_term 799! 800 end module m_dftu 801