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!-------------------------------------------------------------------------- 9MODULE ldaU 10 !-------------------------------------------------------------------------- 11 ! 12 ! The quantities needed in DFT+U and extended DFT+U calculations. 13 ! 14 USE kinds, ONLY : DP 15 USE upf_params, ONLY : lqmax 16 ! FIXME: lqmax should not be used (see starting_ns* below) 17 USE parameters, ONLY : ntypx, natx, sc_size 18 USE basis, ONLY : natomwfc 19 USE ions_base, ONLY : nat, ntyp => nsp, ityp 20 USE control_flags, ONLY : dfpt_hub 21 ! 22 SAVE 23 ! 24 COMPLEX(DP), ALLOCATABLE :: wfcU(:,:) 25 !! atomic wfcs with U term 26 COMPLEX(DP), ALLOCATABLE :: d_spin_ldau(:,:,:) 27 !! the rotations in spin space for all symmetries 28 REAL(DP) :: eth 29 !! the Hubbard contribution to the energy 30 REAL(DP) :: Hubbard_U(ntypx) 31 !! the Hubbard U 32 REAL(DP) :: Hubbard_U_back(ntypx) 33 !! the Hubbard U on background states 34 REAL(DP) :: Hubbard_J0(ntypx) 35 !! the Hubbard J, in simplified DFT+U 36 REAL(DP) :: Hubbard_J(3,ntypx) 37 !! extra Hubbard parameters for full DFT+U: 38 !! * p: J(1)=J 39 !! * d: J(1)=J, J(2)=B 40 !! * f: J(1)=J, J(2)=E2, J(3)=E3 41 REAL(DP) :: Hubbard_alpha(ntypx) 42 !! the Hubbard alpha (used to calculate U) 43 REAL(DP) :: Hubbard_alpha_back(ntypx) 44 !! the Hubbard alpha (used to calculate U on background states) 45 REAL(DP) :: Hubbard_beta(ntypx) 46 !! the Hubbard beta (used to calculate J0) 47 REAL(DP) :: starting_ns(lqmax,2,ntypx) 48 !! starting ns 49 !! FIXME: allocate dynamically 50 REAL(DP) :: starting_ns_back(lqmax,2,ntypx) 51 !! starting ns on background states 52 !! FIXME: allocate dynamically, or better, remove 53 INTEGER :: nwfcU 54 !! total no. of atomic wavefunctions having U term 55 INTEGER :: niter_with_fixed_ns 56 !! no. of iterations with fixed ns 57 INTEGER :: lda_plus_u_kind 58 !! 0 --> Simplified rotationally-invariant formulation of DFT+U 59 !! 1 --> Full formulation of DFT+U 60 !! 2 --> Simplified rotationally-invariant formulation of DFT+U+V 61 INTEGER :: Hubbard_l(ntypx) 62 !! the angular momentum of Hubbard states 63 INTEGER :: Hubbard_l_back(ntypx) 64 !! the angular momentum of Hubbard background states 65 INTEGER :: Hubbard_l1_back(ntypx) 66 !! the angular momentum of the second channel of iHubbard background states 67 INTEGER :: Hubbard_lmax = 0 68 !! maximum angular momentum of Hubbard states 69 INTEGER :: Hubbard_lmax_back = 0 70 !! maximum angular momentum of Hubbard background states 71 INTEGER :: lback(ntypx) 72 !! the angular momentum of background states 73 INTEGER :: l1back(ntypx) 74 !! like above for the second background state 75 !! (there is a possibility to have two background states at the same U) 76 INTEGER :: ldmx = -1 77 !! max dimension of the manifold where the Hubbard correction 78 !! is applied (max of 2*Hubbard_l+1 over all atoms) 79 INTEGER :: ldmx_b = -1 80 !! like above for background states 81 INTEGER :: ldmx_tot 82 !! max value of ldim_u(nt) = ldim_u(nt) + ldim_back(nt) over all ntyp 83 INTEGER :: l0 84 !! index in the array of Hubbard states. 85 !! for every atom one has 2*Hubbard_l+1 + 2*Hubbard_l_back +1 = ldim_u states for example. 86 INTEGER :: l0b 87 !! like above for background states 88 LOGICAL :: is_hubbard(ntypx) 89 !! .TRUE. if this atom species has U correction 90 LOGICAL :: is_hubbard_back(ntypx) 91 !! .TRUE. if this atom species has U correction for background states 92 LOGICAL :: lda_plus_u 93 !! .TRUE. if DFT+U (or extended) calculation is performed 94 LOGICAL :: conv_ns 95 !! .TRUE. if ns are converged 96 LOGICAL :: reserv(ntypx), reserv_back(ntypx) 97 !! reservoir states 98 LOGICAL :: backall(ntypx) 99 !! if .TRUE. two l channels can be used in the background (lback and l1back should be 100 !! specified in input for all the types for which backall is .true.) 101 LOGICAL :: hub_back 102 !! .TRUE. if at least one species has Hubbard_U in the background states 103 LOGICAL :: hub_pot_fix 104 !! if .TRUE. do not include into account the change of the Hubbard potential 105 !! during the SCF cycle (needed to compute U self-consistently with supercells) 106 LOGICAL :: iso_sys 107 !! .TRUE. if the system is isolated (the code diagonalizes 108 !! and prints the full occupation matrix) 109 CHARACTER(len=30) :: U_projection 110 !! 'atomic', 'ortho-atomic', 'file' 111 CHARACTER(len=80) :: Hubbard_parameters 112 !! if 'input' then read Hubbard_V from input, 113 !! if 'file' read them from the file called 'parameters.in' (used only with lda_plus_u_kind = 2) 114 INTEGER, ALLOCATABLE :: oatwfc(:) 115 !! specifies how input coordinates are given 116 INTEGER, ALLOCATABLE :: oatwfc_back(:), oatwfc_back1(:) 117 !! specifies how input coordinates are given for background states 118 INTEGER, ALLOCATABLE :: offsetU(:) 119 !! offset of atomic wfcs used for projections 120 INTEGER, ALLOCATABLE :: offsetU_back(:), offsetU_back1(:) 121 !! offset of atomic wfcs used for projections (background states) 122 INTEGER, ALLOCATABLE :: ldim_u(:) 123 !! number of U states for each atom 124 INTEGER, ALLOCATABLE :: ldim_back(:) 125 !! number of U background states for each atom 126 INTEGER, ALLOCATABLE :: ll(:,:) 127 !! ll(i=1:ldim_u) = Hubbard_l if 0 <= i <= l0 128 !! = Hubbard_l_back if l0+1 < i <= ldim_u 129 REAL(DP), ALLOCATABLE :: q_ae(:,:,:) 130 !! coefficients for projecting onto beta functions 131 REAL(DP), ALLOCATABLE :: q_ps(:,:,:) 132 !! (matrix elements on AE and PS atomic wfcs) 133 !! 134 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 135 !!!!!!!!!!!!!!!!!!!!! Hubbard V part !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 136 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 137 ! 138 ! Inter atomic interaction should be cut off at some distance 139 ! that is the reason of having so many unitcell information. 140 ! 141 REAL(DP) :: Hubbard_V(natx,natx*(2*sc_size+1)**3,4) 142 !! The Hubbard_V(I,J,int_type) gives the interaction between atom I (in the unit cell) 143 !! with atom J (in the supercell). 144 !! If int_type=1, the interaction is between standard orbitals, 145 !! If int_type=2, the interaction is between standard (on I) and background (on J) orbitals, 146 !! If int_type=3, the interaction is between background orbitals, 147 !! If int_type=4, the interaction is between background (on I) and standard (on J) orbitals. 148 !! Hubbard_V(I,J,4) is equal to Hubbard_V(J,I,2). It is useful 149 !! in cases where Hubbard_V(I,J,2) /= 0 but I is outside the unit cell, J inside. 150 INTEGER :: num_uc 151 !! Number of unit cells in the supercell = (2*sc_size+1)**3 152 INTEGER :: max_num_neighbors 153 !! the maximum number of neighbors 154 REAL(DP), ALLOCATABLE :: atom_pos(:,:) 155 !! matrix with dimensions nat x 3 with the atomic coordinates in 156 !! the primitive basis. 157 REAL(DP), ALLOCATABLE :: dist_s(:,:) 158 !! Distance between atoms in the 3x3x3 supercell (if sc_size = 1) 159 INTEGER, ALLOCATABLE :: ityp_s(:) 160 !! Type of atoms in the 3x3x3 supercell (if sc_size = 1) 161 COMPLEX(DP), ALLOCATABLE :: nsg(:,:,:,:,:), nsgnew(:,:,:,:,:) 162 !! Generalized occupation matrices, which depend on two atomic sites. 163 !! These matrices nsg(at1,m1,viz,m2,sp) store the expectation value: 164 !! <C^\dagger_{at1,m1,sp}C_{viz,m2,sp}>, where sp = spin and 165 !! viz identifies the atom in the neighborhood of at1. 166 COMPLEX(DP), ALLOCATABLE :: v_nsg(:,:,:,:,:) 167 !! The kernel of the Hubbard potential (see above for the meaning of the 168 !! size of the array) 169 COMPLEX(DP), ALLOCATABLE :: phase_fac(:) 170 !! Phase factor (it is 1 if we have only Hubbard U) 171 INTEGER, ALLOCATABLE :: sc_at(:,:,:,:) 172 !! Matrix with ranges [1:nat], gives the corresponding atom in the supercell ordering 173 REAL(DP), PARAMETER :: eps_dist = 6.d-4 174 !! Threshold for comparing inter-atomic distances 175 ! 176 TYPE position 177 INTEGER :: at, n(3) 178 !! Identifies the atom: it is equivalent to atom 'at' in the unit cell, 179 !! and it is located in the unit cell (n(1),n(2),n(3)) 180 ENDTYPE position 181 ! 182 TYPE at_center 183 integer :: num_neigh 184 !! Number of neighbours (i.e., other V-interacting atoms) 185 !! of the considered atom 186 integer, ALLOCATABLE :: neigh(:) 187 !! Vector indicating which are the neigbours 188 ENDTYPE at_center 189 ! 190 TYPE(position), ALLOCATABLE :: at_sc(:) 191 !! Vector with all the atoms in the supercell 192 ! 193 TYPE(at_center), ALLOCATABLE :: neighood(:) 194 !! Vector with the information about the neighbours 195 !! for all the atoms in the unit cell 196 ! 197CONTAINS 198 ! 199 SUBROUTINE init_lda_plus_u ( psd, nspin, noncolin ) 200 ! 201 !! NOTE: average_pp must be called before init_lda_plus_u 202 ! 203 IMPLICIT NONE 204 ! 205 CHARACTER (LEN=2), INTENT(IN) :: psd(:) 206 INTEGER, INTENT(IN) :: nspin 207 LOGICAL, INTENT(IN) :: noncolin 208 ! 209 INTEGER, EXTERNAL :: set_Hubbard_l, set_Hubbard_l_back 210 INTEGER :: na, nt 211 LOGICAL :: lba, lb 212 ! 213 lba = .FALSE. 214 lb = .FALSE. 215 hub_back = .FALSE. 216 ! 217 is_hubbard(:) = .FALSE. 218 is_hubbard_back(:) = .FALSE. 219 ! 220 IF ( .NOT. lda_plus_u ) THEN 221 Hubbard_lmax = 0 222 RETURN 223 ENDIF 224 ! 225 Hubbard_lmax = -1 226 ! Set the default of Hubbard_l for the species which have 227 ! Hubbard_U=0 (in that case set_Hubbard_l will not be called) 228 Hubbard_l(:) = -1 229 ! 230 ! Background part 231 ! 232 Hubbard_lmax_back = -1 233 Hubbard_l_back(:) = -1 234 Hubbard_l1_back(:) = -1 235 ldmx_tot = -1 236 ! 237 IF (.NOT.ALLOCATED (ldim_u) ) ALLOCATE(ldim_u(ntyp)) 238 ldim_u(:)=-1 239 IF (.NOT.ALLOCATED (ldim_back) ) ALLOCATE(ldim_back(ntyp)) 240 ldim_back(:)=-1 241 ! 242 IF ( lda_plus_u_kind == 0 ) THEN 243 ! 244 ! DFT+U (simplified) 245 ! 246 DO nt = 1, ntyp 247 ! 248 is_hubbard(nt) = Hubbard_U(nt) /= 0.0_DP .OR. & 249 Hubbard_U_back(nt)/= 0.0_dp .OR. & 250 Hubbard_alpha(nt) /= 0.0_DP .OR. & 251 Hubbard_alpha_back(nt) /= 0.0_dp .OR. & 252 Hubbard_J0(nt) /= 0.0_DP .OR. & 253 Hubbard_beta(nt) /= 0.0_DP 254 ! 255 is_hubbard_back(nt) = Hubbard_U_back(nt)/= 0.0_dp .OR. & 256 Hubbard_alpha_back(nt) /= 0.0_dp 257 ! 258 IF ( is_hubbard(nt) ) THEN 259 Hubbard_l(nt) = set_Hubbard_l( psd(nt) ) 260 Hubbard_lmax = MAX( Hubbard_lmax, Hubbard_l(nt) ) 261 ldmx = MAX( ldmx, 2*Hubbard_l(nt)+1 ) 262 ldim_u(nt) = 2*Hubbard_l(nt)+1 263 ENDIF 264 ! 265 IF ( is_hubbard_back(nt) ) THEN 266 lb = .TRUE. 267 hub_back = .TRUE. 268 ! 269 ! if .not.backall set_hubbard_l_back determines the Hubbard_l_back; otherwise 270 ! Hubbard_l_back and Hubbard_l1_back are set from input (by lback and l1back). 271 ! 272 IF (.NOT.backall(nt)) THEN 273 ! In this case there is only one Hubbard channel for background states 274 Hubbard_l_back(nt) = set_Hubbard_l_back( psd(nt) ) 275 Hubbard_lmax_back = MAX( Hubbard_lmax_back, Hubbard_l_back(nt) ) 276 ldmx_b = MAX( ldmx_b, 2*Hubbard_l_back(nt)+1) 277 ldim_back(nt) = 2 * Hubbard_l_back(nt) + 1 278 ELSE 279 ! In this case there are two Hubbard channels for background states. 280 ! Note: the same U_back is used for these two background channels. 281 lba = .TRUE. 282 Hubbard_l_back(nt) = lback(nt) 283 Hubbard_l1_back(nt) = l1back(nt) 284 Hubbard_lmax_back = MAX( Hubbard_lmax_back, Hubbard_l1_back(nt) ) 285 ldmx_b = MAX( ldmx_b, 2*Hubbard_l_back(nt)+2*Hubbard_l1_back(nt)+2 ) 286 ldim_back(nt) = 2 * (Hubbard_l_back(nt) + Hubbard_l1_back(nt) + 1) 287 ENDIF 288 ldim_u(nt) = ldim_u(nt) + ldim_back(nt) !2 * Hubbard_l1_back(nt) + 1 289 Hubbard_lmax_back = MAX( Hubbard_lmax_back, Hubbard_l_back(nt) ) 290 ELSE 291 backall(nt) = .FALSE. 292 ENDIF 293 ! 294 ldmx_tot = MAX( ldmx_tot, ldim_u(nt) ) 295 ! 296 ENDDO !nt 297 ! 298 IF (ALLOCATED(ll)) DEALLOCATE (ll) 299 ALLOCATE(ll(ldmx_tot,ntyp)) 300 ! 301 ! ll is a label of all the Hubbard states telling the l of that states. 302 ! It is equal to Hubbard_l for the first 2*Hubbard_l+1 states, 303 ! lback for the next 2*Hubbard_l_back+1, 304 ! l1back for the next 2*Hubbard_l1_back+1 305 ! (if there are two channels in the background). 306 ! 307 ll(:,:) = -1 308 DO nt = 1,ntyp 309 IF (Hubbard_l(nt).GE.0) THEN 310 ll(1:2*Hubbard_l(nt)+1,nt) = Hubbard_l(nt) 311 l0 = 2*Hubbard_l(nt)+1 312 ELSE 313 l0 = 0 314 ENDIF 315 IF (Hubbard_l_back(nt).GE.0) THEN 316 ll(l0+1:l0+2*Hubbard_l_back(nt)+1,nt) = & 317 Hubbard_l_back(nt) 318 l0b = l0 + 2*Hubbard_l_back(nt)+1 319 ELSE 320 l0b = 0 321 ENDIF 322 IF (backall(nt) .AND. Hubbard_l1_back(nt).GE.0) THEN 323 ll(l0b+1:l0b+2*Hubbard_l1_back(nt)+1,nt) = & 324 Hubbard_l1_back(nt) 325 ENDIF 326 ENDDO 327 ! 328 ELSEIF ( lda_plus_u_kind == 1 ) THEN 329 ! 330 ! DFT+U (full) 331 ! 332 IF ( U_projection == 'pseudo' ) CALL errore( 'init_lda_plus_u', & 333 & 'full DFT+U not implemented with pseudo projection type', 1 ) 334 ! 335 IF (noncolin) THEN 336 IF ( .NOT. ALLOCATED (d_spin_ldau) ) ALLOCATE( d_spin_ldau(2,2,48) ) 337 CALL comp_dspinldau() 338 ENDIF 339 ! 340 DO nt = 1, ntyp 341 IF (Hubbard_alpha(nt)/=0.d0 ) CALL errore( 'init_lda_plus_u', & 342 'full DFT+U does not support Hubbard_alpha calculation', 1 ) 343 344 is_hubbard(nt) = Hubbard_U(nt)/= 0.0_dp .OR. & 345 ANY( Hubbard_J(:,nt) /= 0.0_dp ) 346 347 IF ( is_hubbard(nt) ) THEN 348 ! 349 Hubbard_l(nt) = set_Hubbard_l( psd(nt) ) 350 Hubbard_lmax = MAX( Hubbard_lmax, Hubbard_l(nt) ) 351 ldmx = MAX( ldmx, 2*Hubbard_l(nt)+1 ) 352 ldim_u(nt) = 2*Hubbard_l(nt)+1 353 ! 354 IF (Hubbard_U(nt) == 0.0_dp) Hubbard_U(nt) = 1.d-14 355 ! 356 IF ( Hubbard_l(nt) == 2 ) THEN 357 IF ( Hubbard_J(2,nt) == 0.d0 ) & 358 Hubbard_J(2,nt) = 0.114774114774d0 * Hubbard_J(1,nt) 359 ELSEIF ( Hubbard_l(nt) == 3 ) THEN 360 IF ( Hubbard_J(2,nt) == 0.d0 ) & 361 Hubbard_J(2,nt) = 0.002268d0 * Hubbard_J(1,nt) 362 IF ( Hubbard_J(3,nt)==0.d0 ) & 363 Hubbard_J(3,nt) = 0.0438d0 * Hubbard_J(1,nt) 364 ENDIF 365 ENDIF 366 ! 367 ENDDO 368 ! 369 ELSEIF ( lda_plus_u_kind == 2 ) THEN 370 ! 371 ! DFT+U+V (simplified) 372 ! 373 ! Number of cells in the supercell 374 num_uc = (2*sc_size+1)**3 375 ! 376 ! Setup atomic positions in the primitive basis coordinates 377 ! 378 CALL alloc_atom_pos() 379 ! 380 ! In this case is_hubbard is set inside alloc_neighborhood 381 ! 382 CALL alloc_neighborhood() 383 ! 384 ldmx = 0 385 ldmx_b = 0 386 lba = .FALSE. 387 ! 388 DO nt = 1, ntyp 389 ! 390 ! Here we account for the remaining cases when we need to 391 ! setup is_hubbard 392 ! 393 is_hubbard(nt) = is_hubbard(nt) .OR. & 394 Hubbard_alpha(nt) /= 0.0_dp .OR. & 395 Hubbard_beta(nt) /= 0.0_dp .OR. & 396 Hubbard_J0(nt) /= 0.0_dp 397 ! 398 is_hubbard_back(nt) = is_hubbard_back(nt) .OR. & 399 Hubbard_alpha_back(nt) /= 0.0_dp 400 ! 401 IF ( is_hubbard(nt) ) THEN 402 ! 403 IF (Hubbard_l(nt).EQ.-1) Hubbard_l(nt) = set_Hubbard_l( psd(nt) ) 404 ! 405 IF (Hubbard_l(nt).GE.0) THEN 406 ldim_u(nt) = 2 * Hubbard_l(nt) + 1 407 Hubbard_lmax = MAX( Hubbard_lmax, Hubbard_l(nt) ) 408 ldmx = MAX( ldmx, ldim_u(nt) ) 409 ENDIF 410 ! 411 IF ( is_hubbard_back(nt) ) THEN 412 ! 413 lb = .TRUE. 414 hub_back = .TRUE. 415 ! 416 ! if .not.backall set_hubbard_l_back determines the Hubbard_l_back; otherwise 417 ! Hubbard_l_back and Hubbard_l1_back are set from input (by lback and l1back). 418 ! 419 IF (.NOT.backall(nt)) THEN 420 ! In this case there is only one Hubbard channel for background states 421 IF (Hubbard_l_back(nt).EQ.-1) & 422 Hubbard_l_back(nt) = set_Hubbard_l_back( psd(nt) ) 423 IF (Hubbard_l_back(nt).GE.0) & 424 ldim_back(nt) = 2 * Hubbard_l_back(nt) + 1 425 ELSE 426 ! In this case there are two Hubbard channels for background states. 427 ! Note: the same Hubbard parameter is used for these two background channels. 428 lba = .TRUE. 429 IF (Hubbard_l1_back(nt).EQ.-1) THEN 430 Hubbard_l_back(nt) = lback(nt) 431 Hubbard_l1_back(nt) = l1back(nt) 432 ENDIF 433 ldim_back(nt) = 2 * Hubbard_l_back(nt) + 2 * Hubbard_l1_back(nt) + 2 434 Hubbard_lmax_back = MAX( Hubbard_lmax_back, Hubbard_l1_back(nt) ) 435 ENDIF 436 ! 437 ldim_u(nt) = ldim_u(nt) + ldim_back(nt) ! 2 * Hubbard_l1_back(nt) + 1 438 ldmx_b = MAX( ldmx_b, ldim_back(nt) ) 439 Hubbard_lmax_back = MAX( Hubbard_lmax_back, Hubbard_l_back(nt) ) 440 ! 441 ENDIF 442 ! 443 ELSE 444 ! 445 ldim_u(nt) = 0 446 ldim_back(nt) = 0 447 ! 448 ENDIF 449 ! 450 ldmx_tot = MAX( ldmx_tot, ldim_u(nt) ) 451 ! 452 ENDDO 453 ! 454 ! The allocation should be moved into scf_mod ? 455 ! 456 ALLOCATE ( v_nsg ( ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin ) ) 457 ALLOCATE ( nsg ( ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin ) ) 458 ALLOCATE ( nsgnew( ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin ) ) 459 ALLOCATE ( phase_fac(nat*num_uc)) 460 ALLOCATE ( ll(ldmx_tot, ntyp)) 461 ! 462 ! ll is a label of all the Hubbard states telling the l of that states. 463 ! It is equal to Hubbard_l for the first 2*Hubbard_l+1 states, 464 ! lback for the next 2*Hubbard_l_back+1, 465 ! l1back for the next 2*Hubbard_l1_back+1 466 ! (if there are two channels in the background). 467 ! 468 ll(:,:) = -1 469 ! 470 DO nt = 1, ntyp 471 ! 472 IF (Hubbard_l(nt).GE.0) THEN 473 ll(1:2*Hubbard_l(nt)+1,nt) = Hubbard_l(nt) 474 l0 = 2*Hubbard_l(nt)+1 475 ELSE 476 l0 = 0 477 ENDIF 478 ! 479 IF (Hubbard_l_back(nt).GE.0) THEN 480 ll(l0+1:l0+2*Hubbard_l_back(nt)+1,nt) = Hubbard_l_back(nt) 481 l0b = l0 + 2*Hubbard_l_back(nt)+1 482 ELSE 483 l0b = 0 484 ENDIF 485 ! 486 IF (backall(nt) .AND. Hubbard_l1_back(nt).GE.0) & 487 & ll(l0b+1:l0b+2*Hubbard_l1_back(nt)+1,nt) = Hubbard_l1_back(nt) 488 ! 489 ENDDO 490 ! 491 ELSE 492 ! 493 CALL errore( 'init_lda_plus_u', 'Not allowed value of lda_plus_u_kind', 1 ) 494 ! 495 ENDIF 496 ! 497 IF ( Hubbard_lmax == -1 ) CALL errore( 'init_lda_plus_u', & 498 'lda_plus_u calculation but Hubbard_l not set', 1 ) 499 ! 500 IF ( Hubbard_lmax > 3 ) & 501 CALL errore( 'init_lda_plus_u', 'Hubbard_l should not be > 3 ', 1 ) 502 ! 503 ! Compute the index of atomic wfcs used as projectors 504 ! 505 IF ( .NOT.ALLOCATED(oatwfc)) ALLOCATE( oatwfc(nat) ) 506 CALL offset_atom_wfc ( .FALSE., 1, oatwfc, nwfcU ) 507 ! 508 IF ( lb .AND. .NOT.ALLOCATED(oatwfc_back)) THEN 509 ALLOCATE ( oatwfc_back(nat) ) 510 CALL offset_atom_wfc ( .FALSE., 2, oatwfc_back, nwfcU ) 511 ENDIF 512 ! 513 IF ( lba .AND. .NOT.ALLOCATED(oatwfc_back1)) THEN 514 ALLOCATE ( oatwfc_back1(nat) ) 515 CALL offset_atom_wfc ( .FALSE., 3, oatwfc_back1, nwfcU ) 516 ENDIF 517 ! 518 ! nwfcU is set to natomwfc by the routine above 519 IF ( nwfcU /= natomwfc ) & 520 CALL errore( 'offset_atom_wfc', 'wrong number of wavefunctions', 1 ) 521 ! 522 ! For each atom, compute the index of its projectors (among projectors only) 523 ! 524 IF ( .NOT.ALLOCATED(offsetU)) ALLOCATE( offsetU(nat) ) 525 CALL offset_atom_wfc( .TRUE., 1, offsetU, nwfcU ) 526 ! 527 IF ( lb .AND. .NOT.ALLOCATED(offsetU_back)) THEN 528 ALLOCATE ( offsetU_back(nat) ) 529 CALL offset_atom_wfc ( .TRUE., 2, offsetU_back, nwfcU ) 530 ENDIF 531 ! 532 IF ( lba .AND. .NOT.ALLOCATED(offsetU_back1)) THEN 533 ALLOCATE ( offsetU_back1(nat) ) 534 CALL offset_atom_wfc ( .TRUE., 3, offsetU_back1, nwfcU ) 535 ENDIF 536 ! nwfcU is set to natomwfc by the routine above 537 ! 538 RETURN 539 ! 540 END SUBROUTINE init_lda_plus_u 541 ! 542 SUBROUTINE deallocate_ldaU( flag ) 543 ! 544 LOGICAL, INTENT(IN) :: flag 545 INTEGER :: na 546 ! 547 IF ( flag ) THEN 548 IF ( ALLOCATED( oatwfc ) ) DEALLOCATE( oatwfc ) 549 IF ( ALLOCATED( oatwfc_back ) ) DEALLOCATE( oatwfc_back ) 550 IF ( ALLOCATED( oatwfc_back1 ) ) DEALLOCATE( oatwfc_back1 ) 551 IF ( ALLOCATED( offsetU ) ) DEALLOCATE( offsetU ) 552 IF ( ALLOCATED( offsetU_back ) ) DEALLOCATE( offsetU_back ) 553 IF ( ALLOCATED( offsetU_back1 ) ) DEALLOCATE( offsetU_back1 ) 554 IF ( ALLOCATED( q_ae ) ) DEALLOCATE( q_ae ) 555 IF ( ALLOCATED( q_ps ) ) DEALLOCATE( q_ps ) 556 IF ( ALLOCATED( d_spin_ldau )) DEALLOCATE( d_spin_ldau ) 557 IF ( ALLOCATED( ll ) ) DEALLOCATE( ll ) 558 IF ( ALLOCATED( v_nsg ) ) DEALLOCATE( v_nsg ) 559 IF ( ALLOCATED( nsg ) ) DEALLOCATE( nsg ) 560 IF ( ALLOCATED( nsgnew ) ) DEALLOCATE( nsgnew ) 561 IF ( ALLOCATED( phase_fac ) ) DEALLOCATE( phase_fac ) 562 IF ( ALLOCATED( atom_pos ) ) DEALLOCATE( atom_pos ) 563 IF ( ALLOCATED( at_sc ) ) DEALLOCATE( at_sc ) 564 IF ( ALLOCATED( sc_at ) ) DEALLOCATE( sc_at ) 565 IF ( ALLOCATED( neighood ) ) THEN 566 DO na = 1, nat 567 CALL deallocate_at_center_type ( neighood(na) ) 568 ENDDO 569 DEALLOCATE( neighood ) 570 ENDIF 571 IF ( ALLOCATED( ldim_u ) ) DEALLOCATE( ldim_u ) 572 IF ( ALLOCATED( ldim_back ) ) DEALLOCATE( ldim_back ) 573 END IF 574 ! 575 IF ( ALLOCATED( wfcU ) ) DEALLOCATE( wfcU ) 576 ! 577 IF (.NOT.dfpt_hub) THEN 578 IF ( ALLOCATED( dist_s ) ) DEALLOCATE( dist_s ) 579 IF ( ALLOCATED( ityp_s ) ) DEALLOCATE( ityp_s ) 580 ENDIF 581 ! 582 RETURN 583 ! 584 END SUBROUTINE deallocate_ldaU 585 ! 586 SUBROUTINE deallocate_at_center_type (neighood_) 587 ! 588 IMPLICIT NONE 589 TYPE (at_center) :: neighood_ 590 ! 591 neighood_%num_neigh = 0 592 ! 593 IF ( ALLOCATED( neighood_%neigh ) ) DEALLOCATE( neighood_%neigh ) 594 ! 595 RETURN 596 ! 597 END SUBROUTINE deallocate_at_center_type 598 ! 599 SUBROUTINE copy_U_wfc( swfcatom, noncolin ) 600 ! 601 ! Copy (orthogonalized) atomic wavefunctions "swfcatom" 602 ! having a Hubbard U correction to array "wfcU" 603 ! 604 IMPLICIT NONE 605 COMPLEX(KIND=DP), INTENT(IN) :: swfcatom(:,:) 606 LOGICAL, INTENT(IN), OPTIONAL :: noncolin 607 LOGICAL :: twice 608 INTEGER :: na, nt, m1, m2 609 610 IF ( PRESENT(noncolin) ) THEN 611 twice = noncolin 612 ELSE 613 twice = .FALSE. 614 ENDIF 615 ! 616 DO na = 1, nat 617 nt = ityp(na) 618 IF ( is_hubbard(nt) ) THEN 619 m1 = 1 620 m2 = 2*hubbard_l(nt)+1 621 IF ( twice ) m2 = 2*m2 622 wfcU(:,offsetU(na)+m1:offsetU(na)+m2) = swfcatom(:,oatwfc(na)+m1:oatwfc(na)+m2) 623 ENDIF 624 IF (is_hubbard_back(nt)) THEN 625 m1 = 1 626 m2 = 2*Hubbard_l_back(nt)+1 627 wfcU(:,offsetU_back(na)+m1:offsetU_back(na)+m2) = & 628 swfcatom(:,oatwfc_back(na)+m1:oatwfc_back(na)+m2) 629 IF (backall(nt)) THEN 630 m1 = 1 631 m2 = 2*Hubbard_l1_back(nt)+1 632 wfcU(:,offsetU_back1(na)+m1:offsetU_back1(na)+m2) = & 633 swfcatom(:,oatwfc_back1(na)+m1:oatwfc_back1(na)+m2) 634 ENDIF 635 ENDIF 636 ENDDO 637 ! 638 RETURN 639 ! 640 END SUBROUTINE copy_U_wfc 641 ! 642END MODULE ldaU 643! 644