1! 2! Copyright (C) 2001-2020 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!----------------------------------------------------------------------- 10SUBROUTINE new_nsg() 11 !----------------------------------------------------------------------- 12 !! This routine computes the new value for nsgnew (the occupation matrices) 13 !! of the DFT+U+V approach. 14 !! These quantities are defined as follows: 15 !! 16 !! nsg_{I,m1,J,m2,s} = \sum_{k,v} 17 !! f_{kv} <\fi^{at}_{I,m1}|\psi_{k,v,s}><\psi_{k,v,s}|\fi^{at}_{J,m2}> 18 ! 19 USE io_global, ONLY : stdout 20 USE kinds, ONLY : DP 21 USE ions_base, ONLY : nat, ityp, ntyp => nsp 22 USE klist, ONLY : nks, ngk 23 USE ldaU, ONLY : Hubbard_l, q_ae, U_projection, wfcU, nwfcU, & 24 ldim_u, ll, neighood, at_sc, nsgnew, phase_fac, & 25 max_num_neighbors, Hubbard_l_back, backall, & 26 offsetU, offsetU_back, offsetU_back1, is_hubbard_back 27 USE symm_base, ONLY : d1, d2, d3 28 USE lsda_mod, ONLY : lsda, current_spin, nspin, isk 29 USE symm_base, ONLY : nsym, irt 30 USE wvfct, ONLY : nbnd, npw, npwx, wg 31 USE control_flags, ONLY : gamma_only 32 USE wavefunctions, ONLY : evc 33 USE io_files, ONLY : nwordwfc, iunwfc, nwordwfcU, iunsat, iunhub 34 USE buffers, ONLY : get_buffer 35 USE mp_global, ONLY : inter_pool_comm 36 USE mp, ONLY : mp_sum 37 USE becmod, ONLY : bec_type, calbec, & 38 allocate_bec_type, deallocate_bec_type 39 ! 40 IMPLICIT NONE 41 ! 42 TYPE (bec_type) :: proj 43 ! proj(nwfcU,nbnd) 44 INTEGER :: ik, ibnd, is, i, na, nb, nt, isym, m1, m2, m11, m22, m0, m00, ldim,i_type 45 INTEGER :: na1, na2, viz, nt1, nt2, ldim1, ldim2, ldimb, nb1, nb2, viz_b 46 INTEGER :: off, off1, off2, off3, eq_na2 47 INTEGER :: ldim_std1, ldim_std2 48 ! in the nwfcU ordering 49 COMPLEX(DP), ALLOCATABLE :: nrg(:,:,:,:,:) 50 COMPLEX(DP) :: phase 51 INTEGER, EXTERNAL :: find_viz, type_interaction 52 ! 53 CALL start_clock('new_nsg') 54 ! 55 ldim = 0 56 DO nt = 1, ntyp 57 ldim = MAX(ldim,ldim_u(nt)) 58 ENDDO 59 ! 60 ALLOCATE ( nrg (ldim, ldim, max_num_neighbors, nat, nspin) ) 61 nrg (:,:,:,:,:) = (0.d0, 0.d0) 62 nsgnew (:,:,:,:,:) = (0.d0, 0.d0) 63 ! 64 CALL allocate_bec_type ( nwfcU, nbnd, proj ) 65 ! 66 ! D_Sl for l=1, l=2 and l=3 are already initialized, for l=0 D_S0 is 1 67 ! 68 ! Offset of atomic wavefunctions initialized in setup and stored in offsetU 69 ! 70 ! we start a loop on k points 71 ! 72 DO ik = 1, nks 73 ! 74 IF (lsda) current_spin = isk(ik) 75 ! 76 npw = ngk(ik) 77 ! 78 IF (nks > 1) CALL get_buffer (evc, nwordwfc, iunwfc, ik) 79 ! 80 ! make the projection 81 ! 82 IF ( U_projection == 'pseudo' ) THEN 83 CALL errore('new_nsg', 'U_projection = pseudo is not supported',1) 84 !CALL compute_pproj( ik, q_ae, proj ) 85 ELSE 86 IF (nks > 1) CALL get_buffer (wfcU, nwordwfcU, iunhub, ik) 87 CALL calbec ( npw, wfcU, evc, proj ) 88 END IF 89 ! 90 ! compute the phase factors for this k and put the result 91 ! in the phase_fac array 92 ! 93 CALL phase_factor(ik) 94 ! 95 ! compute the occupation matrix (nsg_{I,m1,J,m2,s}) of the 96 ! atomic orbitals 97 ! 98 DO ibnd = 1, nbnd 99 ! 100 DO na1 = 1, nat 101 ! 102 nt1 = ityp(na1) 103 ! 104 IF (ldim_u(nt1).GT.0) THEN 105 ! 106 ldim1 = ldim_u(nt1) 107 ! 108 DO viz =1, neighood(na1)%num_neigh 109 ! 110 na2 = neighood(na1)%neigh(viz) 111 eq_na2 = at_sc(na2)%at 112 nt2 = ityp(eq_na2) 113 ldim2 = ldim_u(nt2) 114 ! 115 IF (na1.GT.na2) THEN 116 ! 117 DO m1 = 1, ldim1 118 DO m2 = 1, ldim2 119 nrg(m2,m1,viz,na1,current_spin) = & 120 CONJG(nrg(m1,m2,find_viz(na2,na1),na2,current_spin)) 121 ENDDO 122 ENDDO 123 ! 124 ELSE 125 ! 126 phase = phase_fac(na2) 127 ! 128 DO m1 = 1, ldim1 129 ! 130 off1 = offsetU(na1) + m1 131 ! 132 IF ( is_hubbard_back(nt1) .AND. m1.GT.2*Hubbard_l(nt1)+1 ) THEN 133 ! 134 off1 = offsetU_back(na1) + m1 - 2*Hubbard_l(nt1) - 1 135 ! 136 IF ( backall(nt1) .AND. & 137 m1.GT.2*(Hubbard_l(nt1)+Hubbard_l_back(nt1)+1)) & 138 off1 = offsetU_back1(na1) + m1 - & 139 2*(Hubbard_l(nt1)+Hubbard_l_back(nt1)+1) 140 ! 141 ENDIF 142 ! 143 DO m2 = 1, ldim2 144 ! 145 off2 = offsetU(eq_na2) + m2 146 ! 147 IF ( is_hubbard_back(nt1) .AND. m2.GT.2*Hubbard_l(nt2)+1 ) THEN 148 ! 149 off2 = offsetU_back(eq_na2) + m2 - 2*Hubbard_l(nt2) - 1 150 ! 151 IF ( backall(nt2) .AND. & 152 m2.GT.2*(Hubbard_l(nt2)+Hubbard_l_back(nt2)+1)) & 153 off2 = offsetU_back1(eq_na2) + m2 - & 154 2*(Hubbard_l(nt2)+Hubbard_l_back(nt2)+1) 155 ! 156 ENDIF 157 ! 158 IF (gamma_only) THEN 159 nrg(m2,m1,viz,na1,current_spin) = & 160 nrg(m2,m1,viz,na1,current_spin) + & 161 DBLE( wg(ibnd,ik) * ( proj%r(off1,ibnd) * & 162 CONJG(proj%r(off2,ibnd) * phase) ) ) 163 ELSE 164 nrg(m2,m1,viz,na1,current_spin) = & 165 nrg(m2,m1,viz,na1,current_spin) + & 166 DBLE( wg(ibnd,ik) * ( proj%k(off1,ibnd) * & 167 CONJG(proj%k(off2,ibnd) * phase) ) ) 168 ENDIF 169 ! 170 ENDDO ! m2 171 ENDDO ! m1 172 ! 173 ENDIF 174 ! 175 ENDDO 176 ENDIF 177 ENDDO 178 ENDDO 179 ! 180 ENDDO ! ik 181 ! 182 CALL deallocate_bec_type (proj) 183 ! 184 CALL mp_sum( nrg, inter_pool_comm ) 185 ! 186 IF (nspin.EQ.1) nrg = 0.5d0 * nrg 187 ! 188 ! impose hermiticity of n_{m1,m2} 189 ! 190 DO ibnd = 1,nbnd 191 DO na1 = 1, nat 192 nt1 = ityp(na1) 193 IF (ldim_u(nt1).GT.0) THEN 194 ldim1 = ldim_u(nt1) 195 DO viz = 1, neighood(na1)%num_neigh 196 na2 = neighood(na1)%neigh(viz) 197 IF (na1.GT.na2) THEN 198 eq_na2 = at_sc(na2)%at 199 nt2 = ityp (eq_na2) 200 ldim2 = ldim_u(nt2) 201 DO m1 = 1, ldim1 202 DO m2 = 1, ldim2 203 nrg(m2,m1,viz,na1,current_spin) = & 204 ( nrg(m2,m1,viz,na1,current_spin) + & 205 CONJG(nrg(m1,m2,find_viz(na2,na1),na2,current_spin)) )*0.5d0 206 nrg(m1,m2,find_viz(na2,na1),na2,current_spin) = & 207 CONJG(nrg(m2,m1,viz,na1,current_spin)) 208 ENDDO 209 ENDDO 210 ENDIF 211 ENDDO 212 ENDIF 213 ENDDO 214 ENDDO 215 ! 216 ! symmetrize the quantities nr -> ns 217 ! 218 DO is = 1, nspin 219 ! 220 DO na1 = 1, nat 221 ! 222 nt1 = ityp(na1) 223 ! 224 IF (ldim_u(nt1).GT.0) THEN 225 ! 226 ldim_std1 = 2*Hubbard_l(nt1)+1 227 ! 228 DO viz = 1, neighood(na1)%num_neigh 229 ! 230 na2 = neighood(na1)%neigh(viz) 231 eq_na2 = at_sc(na2)%at 232 nt2 = ityp(eq_na2) 233 ldim_std2 = 2*Hubbard_l(nt2)+1 234 ! 235 IF (na1.GT.na2) THEN 236 ! 237 ! we don't need to compute again 238 DO m1 = 1, ldim_u(nt1) 239 DO m2 = 1, ldim_u(nt2) 240 nsgnew(m2,m1,viz,na1,is) = & 241 CONJG(nsgnew(m1,m2,find_viz(na2,na1),na2,is)) 242 ENDDO 243 ENDDO 244 ! 245 ELSE 246 ! 247 DO m1 = 1, ldim_u(nt1) 248 ! 249 off = 1 250 off2 = 2*Hubbard_l(nt1)+1 251 ! 252 IF ( is_hubbard_back(nt1) .AND. m1.GT.ldim_std1 ) THEN 253 ! 254 off = ldim_std1 + 1 255 off2 = ldim_std1 + 2*Hubbard_l_back(nt1) + 1 256 ! 257 IF ( backall(nt1) .AND. m1.GT.(ldim_std1+2*Hubbard_l_back(nt1)+1) ) THEN 258 off = ldim_std1 + 2*Hubbard_l_back(nt1) + 2 259 off2 = ldim_u(nt1) 260 ENDIF 261 ! 262 ENDIF 263 ! 264 DO m2 = 1, ldim_u(nt2) 265 ! 266 off1 = 1 267 off3 = 2*Hubbard_l(nt2)+1 268 ! 269 IF ( is_hubbard_back(nt2) .AND. m2.GT.ldim_std2 ) THEN 270 ! 271 off1 = ldim_std2 + 1 272 off3 = ldim_std2 + 2*Hubbard_l_back(nt2) + 1 273 ! 274 IF ( backall(nt2) .AND. m2.GT.(ldim_std2+2*Hubbard_l_back(nt2)+1) ) THEN 275 off1 = ldim_std2 + 2*Hubbard_l_back(nt2) + 2 276 off3 = ldim_u(nt2) 277 ENDIF 278 ! 279 ENDIF 280 ! 281 ! Perform symmetrization using all available symmetries 282 ! 283 DO isym = 1, nsym 284 ! 285 CALL symonpair(na1,na2,isym,nb1,nb2) 286 ! 287 viz_b = find_viz(nb1,nb2) 288 ! 289 DO m0 = off, off2 290 DO m00 = off1, off3 291 IF (ll(m1,nt1).EQ.0 .AND. & 292 m1.ge.off.and.m1.le.off2.and. & 293 m2.ge.off1.and.m2.le.off3.and. & 294 ll(m2,nt2).EQ.0) THEN 295 nsgnew(m2,m1,viz,na1,is) = & 296 nsgnew(m2,m1,viz,na1,is) + & 297 nrg(m00,m0,viz_b,nb1,is) / nsym 298 ELSE IF (ll(m1,nt1).EQ.1 .AND. & 299 m1.ge.off.and.m1.le.off2.and. & 300 m2.ge.off1.and.m2.le.off3.and. & 301 ll(m2,nt2).EQ.0) THEN 302 nsgnew(m2,m1,viz,na1,is) = & 303 nsgnew(m2,m1,viz,na1,is) + & 304 d1(m0-off+1,m1-off+1,isym) * & 305 nrg(m00,m0,viz_b,nb1,is) * & 306 1.d0 / nsym 307 ELSE IF (ll(m1,nt1).EQ.2 .AND. & 308 m1.ge.off.and.m1.le.off2.and. & 309 m2.ge.off1.and.m2.le.off3.and. & 310 ll(m2,nt2).EQ.0) THEN 311 nsgnew(m2,m1,viz,na1,is) = & 312 nsgnew(m2,m1,viz,na1,is) + & 313 d2(m0-off+1,m1-off+1,isym) * & 314 nrg(m00,m0,viz_b,nb1,is) * & 315 1.d0 / nsym 316 ELSE IF (ll(m1,nt1).EQ.3 .AND. & 317 m1.ge.off.and.m1.le.off2.and. & 318 m2.ge.off1.and.m2.le.off3.and. & 319 ll(m2,nt2).EQ.0) THEN 320 nsgnew(m2,m1,viz,na1,is) = & 321 nsgnew(m2,m1,viz,na1,is) + & 322 d3(m0-off+1,m1-off+1,isym) * & 323 nrg(m00,m0,viz_b,nb1,is) * & 324 1.d0 / nsym 325 326 ELSE IF (ll(m1,nt1).EQ.0 .AND. & 327 m1.ge.off.and.m1.le.off2.and. & 328 m2.ge.off1.and.m2.le.off3.and. & 329 ll(m2,nt2).EQ.1) THEN 330 nsgnew(m2,m1,viz,na1,is) = & 331 nsgnew(m2,m1,viz,na1,is) + & 332 nrg(m00,m0,viz_b,nb1,is) * & 333 d1(m00-off1+1,m2-off1+1,isym) / nsym 334 ELSE IF (ll(m1,nt1).EQ.1 .AND. & 335 m1.ge.off.and.m1.le.off2.and. & 336 m2.ge.off1.and.m2.le.off3.and. & 337 ll(m2,nt2).EQ.1) THEN 338 nsgnew(m2,m1,viz,na1,is) = & 339 nsgnew(m2,m1,viz,na1,is) + & 340 d1(m0-off+1,m1-off+1,isym) * & 341 nrg(m00,m0,viz_b,nb1,is) * & 342 d1(m00-off1+1,m2-off1+1,isym) / nsym 343 ELSE IF (ll(m1,nt1).EQ.2 .AND. & 344 m1.ge.off.and.m1.le.off2.and. & 345 m2.ge.off1.and.m2.le.off3.and. & 346 ll(m2,nt2).EQ.1) THEN 347 nsgnew(m2,m1,viz,na1,is) = & 348 nsgnew(m2,m1,viz,na1,is) + & 349 d2(m0-off+1,m1-off+1,isym) * & 350 nrg(m00,m0,viz_b,nb1,is) * & 351 d1(m00-off1+1,m2-off1+1,isym) / nsym 352 ELSE IF (ll(m1,nt1).EQ.3 .AND. & 353 m1.ge.off.and.m1.le.off2.and. & 354 m2.ge.off1.and.m2.le.off3.and. & 355 ll(m2,nt2).EQ.1) THEN 356 nsgnew(m2,m1,viz,na1,is) = & 357 nsgnew(m2,m1,viz,na1,is) + & 358 d3(m0-off+1,m1-off+1,isym) * & 359 nrg(m00,m0,viz_b,nb1,is) * & 360 d1(m00-off1+1,m2-off1+1,isym) / nsym 361 362 ELSE IF (ll(m1,nt1).EQ.0 .AND. & 363 m1.ge.off.and.m1.le.off2.and. & 364 m2.ge.off1.and.m2.le.off3.and. & 365 ll(m2,nt2).EQ.2) THEN 366 nsgnew(m2,m1,viz,na1,is) = & 367 nsgnew(m2,m1,viz,na1,is) + & 368 nrg(m00,m0,viz_b,nb1,is) * & 369 d2(m00-off1+1,m2-off1+1,isym) / nsym 370 ELSE IF (ll(m1,nt1).EQ.1 .AND. & 371 m1.ge.off.and.m1.le.off2.and. & 372 m2.ge.off1.and.m2.le.off3.and. & 373 ll(m2,nt2).EQ.2) THEN 374 nsgnew(m2,m1,viz,na1,is) = & 375 nsgnew(m2,m1,viz,na1,is) + & 376 d1(m0-off+1,m1-off+1,isym) * & 377 nrg(m00,m0,viz_b,nb1,is) * & 378 d2(m00-off1+1,m2-off1+1,isym) / nsym 379 ELSE IF (ll(m1,nt1).EQ.2 .AND. & 380 m1.ge.off.and.m1.le.off2.and. & 381 m2.ge.off1.and.m2.le.off3.and. & 382 ll(m2,nt2).EQ.2) THEN 383 nsgnew(m2,m1,viz,na1,is) = & 384 nsgnew(m2,m1,viz,na1,is) + & 385 d2(m0-off+1,m1-off+1,isym) * & 386 nrg(m00,m0,viz_b,nb1,is) * & 387 d2(m00-off1+1,m2-off1+1,isym) / nsym 388 ELSE IF (ll(m1,nt1).EQ.3 .AND. & 389 m1.ge.off.and.m1.le.off2.and. & 390 m2.ge.off1.and.m2.le.off3.and. & 391 ll(m2,nt2).EQ.2) THEN 392 nsgnew(m2,m1,viz,na1,is) = & 393 nsgnew(m2,m1,viz,na1,is) + & 394 d3(m0-off+1,m1-off+1,isym) * & 395 nrg(m00,m0,viz_b,nb1,is) * & 396 d2(m00-off1+1,m2-off1+1,isym) / nsym 397 398 ELSE IF (ll(m1,nt1).EQ.0 .AND. & 399 m1.ge.off.and.m1.le.off2.and. & 400 m2.ge.off1.and.m2.le.off3.and. & 401 ll(m2,nt2).EQ.3) THEN 402 nsgnew(m2,m1,viz,na1,is) = & 403 nsgnew(m2,m1,viz,na1,is) + & 404 nrg(m00,m0,viz_b,nb1,is) * & 405 d3(m00-off1+1,m2-off1+1,isym) / nsym 406 ELSE IF (ll(m1,nt1).EQ.1 .AND. & 407 m1.ge.off.and.m1.le.off2.and. & 408 m2.ge.off1.and.m2.le.off3.and. & 409 ll(m2,nt2).EQ.3) THEN 410 nsgnew(m2,m1,viz,na1,is) = & 411 nsgnew(m2,m1,viz,na1,is) + & 412 d1(m0-off+1,m1-off+1,isym) * & 413 nrg(m00,m0,viz_b,nb1,is) * & 414 d3(m00-off1+1,m2-off1+1,isym) / nsym 415 ELSE IF (ll(m1,nt1).EQ.2 .AND. & 416 m1.ge.off.and.m1.le.off2.and. & 417 m2.ge.off1.and.m2.le.off3.and. & 418 ll(m2,nt2).EQ.3) THEN 419 nsgnew(m2,m1,viz,na1,is) = & 420 nsgnew(m2,m1,viz,na1,is) + & 421 d2(m0-off+1,m1-off+1,isym) * & 422 nrg(m00,m0,viz_b,nb1,is) * & 423 d3(m00-off1+1,m2-off1+1,isym) / nsym 424 ELSE IF (ll(m1,nt1).EQ.3 .AND. & 425 m1.ge.off.and.m1.le.off2.and. & 426 m2.ge.off1.and.m2.le.off3.and. & 427 ll(m2,nt2).EQ.3) THEN 428 nsgnew(m2,m1,viz,na1,is) = & 429 nsgnew(m2,m1,viz,na1,is) + & 430 d3(m0-off+1,m1-off+1,isym) * & 431 nrg(m00,m0,viz_b,nb1,is) * & 432 d3(m00-off1+1,m2-off1+1,isym) / nsym 433 ELSE 434 CALL errore ('new_nsg', & 435 'angular momentum not implemented for at least one type', & 436 ABS(Hubbard_l(nt1)) ) 437 END IF 438 ENDDO !m00 439 ENDDO !m0 440 ! 441 ENDDO !isym 442 ! 443 ENDDO !m2 444 ! 445 ENDDO !m1 446 ! 447 ENDIF !na1 > na2 448 ! 449 ENDDO !viz 450 ! 451 ENDIF !ldim_u > 0 452 ! 453 ENDDO !na1 454 ! 455 ENDDO !is 456 ! 457 DEALLOCATE(nrg) 458 ! 459 CALL stop_clock('new_nsg') 460 ! 461 RETURN 462 ! 463END SUBROUTINE new_nsg 464