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 alloc_neighborhood() 11 !--------------------------------------------------------------------- 12 ! 13 ! This routine allocates and assigns values to the neighborhood, at_sc and sc_at 14 ! 15 USE symm_base, ONLY : nsym 16 USE io_global, ONLY : stdout 17 USE ions_base, ONLY : nat, tau, ityp 18 USE cell_base, ONLY : at, alat 19 USE kinds, ONLY : DP 20 USE constants, ONLY : rytoev 21 USE parameters, ONLY : sc_size 22 USE control_flags, ONLY : dfpt_hub 23 USE ldaU, ONLY : num_uc, max_num_neighbors, neighood, & 24 at_sc, sc_at, Hubbard_V, is_hubbard, is_hubbard_back, & 25 dist_s, ityp_s, deallocate_at_center_type, eps_dist 26 ! 27 IMPLICIT NONE 28 ! 29 ! Local variables 30 ! 31 REAL(DP), ALLOCATABLE :: V(:,:,:), tau_sc(:,:) 32 REAL(DP) :: daux, distm(7) 33 REAL(DP), PARAMETER :: eps1 = 1.d-20 ! threshold for Hubbard_V 34 INTEGER :: i, j, k, l, ii, jj, nx, ny, nz, dimn, & 35 viz, atom, nb1, nb2, isym, l1, l2, l3, na, nb 36 ! 37 CALL start_clock( 'alloc_neigh' ) 38 ! 39 ! Number of atoms in the supercell 40 dimn = num_uc * nat 41 ! 42 ALLOCATE (V(nat,dimn,3)) 43 ALLOCATE (tau_sc(3,dimn)) 44 IF (.NOT.ALLOCATED(ityp_s)) ALLOCATE (ityp_s(dimn)) 45 IF (.NOT.ALLOCATED(dist_s)) ALLOCATE (dist_s(nat,dimn)) 46 ! 4713 CONTINUE 48 ! 49 V(:,:,:) = 0.d0 50 ! 51 ! Build V from Hubbard_V (input plus equivalent neighborhood) 52 ! 53 DO na = 1, nat 54 DO nb = 1, dimn 55 IF (ABS(Hubbard_V(na,nb,1)) /= 0.0_dp .OR. & 56 ABS(Hubbard_V(na,nb,2)) /= 0.0_dp .OR. & 57 ABS(Hubbard_V(na,nb,3)) /= 0.0_dp ) THEN 58 ! 59 V(na,nb,1:3) = Hubbard_V(na,nb,1:3) 60 ! 61 ENDIF 62 ENDDO 63 ENDDO 64 ! 65 IF (.NOT.ALLOCATED(neighood)) ALLOCATE (neighood(nat)) 66 IF (.NOT.ALLOCATED(at_sc)) ALLOCATE (at_sc(dimn)) 67 IF (.NOT.ALLOCATED(sc_at)) ALLOCATE (sc_at(1:nat,-sc_size:sc_size,-sc_size:sc_size,-sc_size:sc_size)) 68 ! 69 ! Initialization of various quantities in the supercell 70 ! First initialization in the real unit cell 71 ! 72 atom = 0 73 ! 74 DO na = 1, nat 75 atom = atom + 1 76 at_sc(atom)%n(1) = 0 77 at_sc(atom)%n(2) = 0 78 at_sc(atom)%n(3) = 0 79 at_sc(atom)%at = na ! it is equivalent to itself 80 sc_at(na,0,0,0) = atom 81 tau_sc(:,atom) = tau(:,na) 82 ityp_s(atom) = ityp(na) 83 ENDDO 84 ! 85 ! Now initialization in the virtual cells of the supercell 86 ! (in the Cartesian framework) 87 ! 88 DO nx = -sc_size, sc_size 89 DO ny = -sc_size, sc_size 90 DO nz = -sc_size, sc_size 91 ! 92 IF ( nx.NE.0 .OR. ny.NE.0 .OR. nz.NE.0 ) THEN 93 ! 94 DO na = 1, nat 95 ! 96 atom = atom + 1 97 at_sc(atom)%n(1) = nx 98 at_sc(atom)%n(2) = ny 99 at_sc(atom)%n(3) = nz 100 at_sc(atom)%at = na 101 sc_at(na,nx,ny,nz) = atom 102 tau_sc(:,atom) = tau(:,na) + DBLE(nx)*at(:,1) & 103 + DBLE(ny)*at(:,2) & 104 + DBLE(nz)*at(:,3) 105 ityp_s(atom) = ityp(na) 106 ! 107 ENDDO 108 ! 109 ENDIF 110 ! 111 ENDDO 112 ENDDO 113 ENDDO 114 ! 115 ! Initialize is_hubbard and is_hubbard_back 116 ! 117 DO na = 1, nat 118 ! 119 DO nb = 1, dimn 120 ! 121 IF (ABS(V(na,nb,1)) /= 0.0_dp .OR. & 122 ABS(V(na,nb,2)) /= 0.0_dp .OR. & 123 ABS(V(na,nb,3)) /= 0.0_dp ) THEN 124 ! 125 IF (ABS(V(na,nb,1)) /= 0.0_dp) THEN 126 is_hubbard(ityp(na)) = .TRUE. 127 is_hubbard(ityp_s(nb)) = .TRUE. 128 ENDIF 129 ! 130 IF (ABS(V(na,nb,2)) /= 0.0_dp) THEN 131 is_hubbard(ityp(na)) = .TRUE. 132 is_hubbard_back(ityp_s(nb)) = .TRUE. 133 ENDIF 134 ! 135 IF (ABS(V(na,nb,3)) /= 0.0_dp) THEN 136 is_hubbard_back(ityp(na)) = .TRUE. 137 is_hubbard_back(ityp_s(nb)) = .TRUE. 138 ENDIF 139 ! 140 ENDIF 141 ! 142 ENDDO 143 ! 144 ENDDO 145 ! 146 ! Compute the distances between atom na and all other atoms 147 ! in the supercell 148 ! 149 dist_s(:,:) = 0.d0 150 ! 151 DO na = 1, nat 152 ! 153 DO nb = 1, dimn 154 ! 155 DO k = 1, 3 156 dist_s(na,nb) = dist_s(na,nb) + (tau_sc(k,na) - tau_sc(k,nb))**2 157 ENDDO 158 dist_s(na,nb) = DSQRT(dist_s(na,nb)) * alat 159 ! 160 ! Uncomment the following lines if you want that the 161 ! distances are printed in the output file 162 ! 163 !IF (.NOT.dfpt_hub) WRITE(stdout,*) na, nb, ityp_s(na), ityp_s(nb), dist_s(na,nb) 164 ! 165 ENDDO 166 ! 167 ENDDO 168 ! 169 ! Order distances 170 ! 171 IF (.NOT.dfpt_hub) WRITE(stdout,'(5x,"First shells distances (in Bohr):")') 172 ! 173 ! Here we determine 7 smallest distances from the first atom 174 ! to its neighbors (in the increasing order). If larger distances 175 ! are needed, then 7 has to be changed to a larger number. 176 ! 177 l = 0 178 distm(:) = 10000.d0 179 ! 180 DO na = 1, 7 181 DO nb = 1, dimn 182 DO k = 1, l 183 IF ( ABS(dist_s(1,nb)-distm(k)).LE.1.d-6 ) GO TO 49 184 ENDDO 185 distm(l+1) = MIN(distm(l+1), dist_s(1,nb)) 18649 CONTINUE 187 ENDDO 188 l = l + 1 189 IF (.NOT.dfpt_hub) WRITE(stdout,'(5x,"shell:",2x,i2,2x,f10.6)') l, distm(l) 190 ENDDO 191 ! 192 ! Find the missing Hubbard_V(i,j,2) i.e. standard-background interactions. 193 ! This loop is needed only when the background is included. 194 ! 195 IF ( ANY(V(:,:,2).GE.eps1) ) THEN 196 ! 197 DO i = 1, nat 198 DO j = 1, dimn 199 DO l = 1, nat 200 DO k = 1, dimn 201 ! 202 IF ( ityp_s(k).EQ.ityp_s(j) .AND. & 203 ityp(l).EQ.ityp(i) .AND. & 204 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 205 ABS(V(l,k,2)).GE.eps1 ) THEN 206 ! 207 IF (Hubbard_V(i,j,2).EQ.0.d0) Hubbard_V(i,j,2) = V(l,k,2) 208 GO TO 22 209 ! 210 ENDIF 211 ! 212 ENDDO ! k 213 ENDDO ! l 21422 CONTINUE 215 ENDDO ! j 216 ENDDO ! i 217 ! 218 DO i = 1, nat 219 DO j = 1, dimn 220 DO l = 1, nat 221 DO k = 1, dimn 222 ! 223 IF ( ityp_s(k).EQ.ityp(i) .AND. & 224 ityp(l).EQ.ityp_s(j) .AND. & 225 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 226 ABS(V(l,k,2)).GE.eps1 ) THEN 227 ! 228 IF (Hubbard_V(i,j,4).EQ.0.d0) Hubbard_V(i,j,4) = V(l,k,2) 229 GO TO 23 230 ! 231 ENDIF 232 ! 233 ENDDO ! k 234 ENDDO ! l 23523 CONTINUE 236 ENDDO ! j 237 ENDDO ! i 238 ENDIF 239 ! 240 ! Find the missing Hubbard_V(i,j,1) i.e. standard-standard interactions. 241 ! 242 DO i = 1, nat 243 DO j = 1, dimn 244 DO l = 1, nat 245 DO k = 1, dimn 246 ! 247 IF ( ityp_s(k).EQ.ityp_s(j) .AND. & 248 ityp(l).EQ.ityp(i) .AND. & 249 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 250 ( ABS(V(l,k,1)).GE.eps1 ) ) THEN 251 ! 252 IF (Hubbard_V(i,j,1).EQ.0.d0) Hubbard_V(i,j,1) = V(l,k,1) 253 GO TO 32 254 ! 255 ENDIF 256 ! 257 IF ( ityp_s(k).EQ.ityp(i) .AND. & 258 ityp(l).EQ.ityp_s(j) .AND. & 259 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 260 ( ABS(V(l,k,1)).GE.eps1 ) ) THEN 261 ! 262 IF (Hubbard_V(i,j,1).EQ.0.d0) Hubbard_V(i,j,1) = V(l,k,1) 263 GO TO 32 264 ! 265 ENDIF 266 ! 267 ENDDO 268 ENDDO 26932 CONTINUE 270 ENDDO 271 ENDDO 272 ! 273 ! Find the missing Hubbard_V(i,j,3) i.e. background-background interactions. 274 ! This loop is needed only when the background is included. 275 ! 276 IF ( ANY(V(:,:,3).GE.eps1) ) THEN 277 DO i = 1, nat 278 DO j = 1, dimn 279 DO l = 1, nat 280 DO k = 1, dimn 281 ! 282 IF ( ityp_s(k).EQ.ityp_s(j) .AND. & 283 ityp(l).EQ.ityp(i) .AND. & 284 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 285 ( ABS(V(l,k,3)).GE.eps1 ) ) THEN 286 ! 287 IF (Hubbard_V(i,j,3).EQ.0.d0) Hubbard_V(i,j,3) = V(l,k,3) 288 GO TO 33 289 ! 290 ENDIF 291 ! 292 IF ( ityp_s(k).EQ.ityp(i) .AND. & 293 ityp(l).EQ.ityp_s(j) .AND. & 294 ABS(dist_s(l,k)-dist_s(i,j)).LE.eps_dist .AND. & 295 ( ABS(V(l,k,3)).GE.eps1) ) THEN 296 ! 297 IF (Hubbard_V(i,j,3).EQ.0.d0) Hubbard_V(i,j,3) = V(l,k,3) 298 GO TO 33 299 ! 300 ENDIF 301 ! 302 ENDDO ! k 303 ENDDO ! l 30433 CONTINUE 305 ENDDO ! j 306 ENDDO ! i 307 ENDIF 308 ! 309 IF (.NOT.dfpt_hub) WRITE(stdout,'(/5x,"i",4x,"j",2x,"dist (Bohr)", & 310 & 7x,"stan-stan",x,"stan-bac",x,"bac-bac",x,"bac-stan")') 311 ! stan-stan = standard-standard 312 ! stan-bac = standard-background 313 ! bac-bac = background-background 314 ! bac-stan = background-standard 315 ! 316 ! Determine how many neighbors there are and what are their indices. 317 ! 318 max_num_neighbors = 0 319 ! 320 DO i = 1, nat 321 ! 322 neighood(i)%num_neigh = 0 323 ! 324 DO j = 1, dimn 325 IF ( ABS(Hubbard_V(i,j,1)).GE.eps1 .OR. & 326 ABS(Hubbard_V(i,j,2)).GE.eps1 .OR. & 327 ABS(Hubbard_V(i,j,3)).GE.eps1 .OR. & 328 i==j ) THEN 329 ! 330 ! Diagonal term i=j needed in init_nsg, nsg_adj and mix_rho 331 ! 332 neighood(i)%num_neigh = neighood(i)%num_neigh + 1 333 ! 334 ENDIF 335 ENDDO 336 ! 337 IF (max_num_neighbors .LT. neighood(i)%num_neigh) THEN 338 max_num_neighbors = neighood(i)%num_neigh 339 ENDIF 340 ! 341 ALLOCATE (neighood(i)%neigh(1:neighood(i)%num_neigh)) 342 ! 343 viz = 0 344 ! 345 DO j = 1, dimn 346 IF ( ABS(Hubbard_V(i,j,1)).GE.eps1 .OR. & 347 ABS(Hubbard_V(i,j,2)).GE.eps1 .OR. & 348 ABS(Hubbard_V(i,j,3)).GE.eps1 .OR. & 349 i==j ) THEN 350 ! 351 ! Diagonal term i=j needed in init_nsg, nsg_adj and mix_rho 352 ! 353 IF (.NOT.dfpt_hub) WRITE(stdout,'(2x,i4,x,i4,x,f12.8,x,a5,4(x,f8.4))') & 354 i, j, dist_s(i,j), ' V = ', (Hubbard_V(i,j,k)*rytoev, k=1,4) 355 ! 356 viz = viz + 1 357 ! 358 neighood(i)%neigh(viz) = j 359 ! 360 ENDIF 361 ENDDO 362 ! 363 ENDDO 364 ! 365 nb2 = 0 366 ! 367 DO i = 1, nat 368 ! 369 DO nb1 = 1, neighood(i)%num_neigh 370 ! 371 j = neighood(i)%neigh(nb1) 372 ! 373 DO isym = 1, nsym 374 ! 375 IF ( ABS(Hubbard_V(i,j,1)) /= 0.0_dp .OR. & 376 ABS(Hubbard_V(i,j,2)) /= 0.0_dp .OR. & 377 ABS(Hubbard_V(i,j,3)) /= 0.0_dp ) THEN 378 ! 379 CALL symonpair(i,j,isym,ii,jj) 380 ! 381 IF (ABS(dist_s(i,j)-dist_s(ii,jj)).GT.eps_dist) THEN 382 WRITE(stdout,'(/2x,"Different distances between couples!")') 383 WRITE(stdout,'(2x,"Original couple: ",2x,i4,2x,i4,2x,"dist =",1x,f8.4,1x,"(Bohr)")') & 384 i, j, dist_s(i,j) 385 WRITE(stdout,'(2x,"New additional couple:",2x,i4,2x,i4,2x,"dist =",1x,f8.4,1x,"(Bohr)")') & 386 ii, jj, dist_s(ii,jj) 387 CALL errore('alloc_neighborhood', 'Probably a larger sc_size is needed',1) 388 ENDIF 389 ! 390 IF ( ABS(Hubbard_V(ii,jj,1)).LT.eps1 .AND. & 391 ABS(Hubbard_V(ii,jj,2)).LT.eps1 .AND. & 392 ABS(Hubbard_V(ii,jj,3)).LT.eps1 ) THEN 393 ! 394 Hubbard_V(ii,jj,:) = Hubbard_V(i,j,:) 395 WRITE(stdout,'(2x,"Found a new Hubbard_V element from the symmetry analysis!")') 396 WRITE(stdout,'(2x,"Original couple: ",2x,i4,2x,i4,2x,"dist =",1x,f8.4,1x,"(Bohr)")') & 397 i, j, dist_s(i,j) 398 WRITE(stdout,'(2x,"New additional couple:",2x,i4,2x,i4,2x,"dist =",1x,f8.4,1x,"(Bohr)")') & 399 ii, jj, dist_s(ii,jj) 400 ! 401 nb2 = nb2 + 1 402 ! 403 ELSEIF ( ABS(Hubbard_V(ii,jj,1)-Hubbard_V(i,j,1)) + & 404 ABS(Hubbard_V(ii,jj,2)-Hubbard_V(i,j,2)) + & 405 ABS(Hubbard_V(ii,jj,3)-Hubbard_V(i,j,3)).GT.eps1 ) THEN 406 ! 407 ! This V element has been set to something different before 408 ! possibly, wrong assignement/symm ops... 409 ! 410 ! Some more output for the standard-standard case 411 ! 412 IF (ABS(Hubbard_V(ii,jj,1)-Hubbard_V(i,j,1)).GT.eps1) THEN 413 ! 414 WRITE(stdout,'(/4x,"WARNING! Equivalent couples with some small variations in V :")') 415 WRITE(stdout,'(4x,"V (",i3,",",i5,") =",1x,f10.6)') & 416 i, j, Hubbard_V(i,j,1)*rytoev 417 WRITE(stdout,'(4x,"V (",i3,",",i5,") =",1x,f10.6)') & 418 ii, jj, Hubbard_V(ii,jj,1)*rytoev 419 WRITE(stdout,'(4x,"delta(V) =",1x,f10.6)') & 420 ABS(Hubbard_V(ii,jj,1)-Hubbard_V(i,j,1))*rytoev 421 ! 422 IF ( ABS(Hubbard_V(ii,jj,1)-Hubbard_V(i,j,1))*rytoev .GE. 5.0d-3 ) THEN 423 WRITE(stdout,*) ' i, j, V:', i, j, Hubbard_V(i,j,1)*rytoev 424 WRITE(stdout,*) 'ii, jj, V:', ii, jj, Hubbard_V(ii,jj,1)*rytoev 425 CALL errore('alloc_neighborhood', & 426 & 'Problem with Hubbard_V after a symmetry analysis',1) 427 ELSE 428 WRITE(stdout,'(4x,"Setting them to be equal...")') 429 Hubbard_V(ii,jj,1) = Hubbard_V(i,j,1) 430 ENDIF 431 ! 432 ENDIF 433 ! 434 IF ( ABS(Hubbard_V(ii,jj,2)-Hubbard_V(i,j,2))*rytoev .GE. 5.0d-3 .OR. & 435 ABS(Hubbard_V(ii,jj,3)-Hubbard_V(i,j,3))*rytoev .GE. 5.0d-3 ) THEN 436 CALL errore('alloc_neighborhood', & 437 & 'Problem with Hubbard_V after a symmetry analysis',1) 438 ENDIF 439 ! 440 ENDIF 441 ! 442 ENDIF 443 ! 444 ENDDO ! isym 445 ! 446 ENDDO ! nb1 447 ! 448 ENDDO ! i 449 ! 450 IF ( nb2.GT.0 ) THEN 451 DO na = 1, nat 452 CALL deallocate_at_center_type ( neighood(na) ) 453 ENDDO 454 WRITE(stdout,'(2x,"Check all Hubbard_V again based on the distances...")') 455 GO TO 13 456 ENDIF 457 ! 458 DEALLOCATE (V) 459 DEALLOCATE (tau_sc) 460 ! 461 CALL stop_clock( 'alloc_neigh' ) 462 ! 463 RETURN 464 ! 465END SUBROUTINE alloc_neighborhood 466!----------------------------------------------------------------------- 467 468!----------------------------------------------------------------------- 469FUNCTION type_interaction (na1, m1, na2, m2) 470 !----------------------------------------------------------------------- 471 ! 472 ! This function determines type of the interaction 473 ! 474 USE ions_base, ONLY : ityp 475 USE ldaU, ONLY : Hubbard_l 476 ! 477 IMPLICIT NONE 478 INTEGER :: na1, na2, & ! atom number (in the unit cell) 479 m1, m2 ! magnetic quantum number 480 INTEGER :: type_interaction, nt1, nt2 481 ! 482 nt1 = ityp(na1) 483 nt2 = ityp(na2) 484 ! 485 IF ( m1 .LE. 2*Hubbard_l(nt1)+1 .AND. & 486 m2 .LE. 2*Hubbard_l(nt2)+1) THEN 487 ! 488 ! standard-standard term 489 type_interaction = 1 490 ! 491 ELSEIF ( m1 .LE. 2*Hubbard_l(nt1)+1 .AND. & 492 m2 .GT. 2*Hubbard_l(nt2)+1 ) THEN 493 ! 494 ! standard-background term 495 type_interaction = 2 496 ! 497 ELSEIF ( m1 .GT. 2*Hubbard_l(nt1)+1 .AND. & 498 m2 .GT. 2*Hubbard_l(nt2)+1 ) THEN 499 ! 500 ! background-background term 501 type_interaction = 3 502 ! 503 ELSE 504 ! 505 ! background-standard term 506 type_interaction = 4 507 ! 508 ENDIF 509 ! 510 RETURN 511 ! 512END FUNCTION type_interaction 513!----------------------------------------------------------------------- 514 515!----------------------------------------------------------------------- 516SUBROUTINE symonpair (at1, at2, p_sym, rat1, rat2) 517 !----------------------------------------------------------------------- 518 ! 519 ! This subroutine acts the symmetry operatorion number p_sym over the 520 ! atomic position of atom at1 in the reference unit cell and over the 521 ! atomic position of atom at2 in the supercell. 522 ! Then it identifies, in the reference unit cell, 523 ! the equivalent atom to the first 'rotation' above, storing the result 524 ! in rat1 and the equivalent atom to the second 'rotation', storing the 525 ! result in rat2 (temporarily). 526 ! The difference dx = O(p_sym)at1 - rat1 is also subtracted from 527 ! O(p_sym)at2 to find the final position of the second atom. 528 ! To determine in which unit cell the second atom ended up, we 529 ! compute the difference between its final position and the position 530 ! of its equivalent in the reference unit cell. 531 ! Finally, we store in rat2 the position, in the full list of atoms in 532 ! the supercell, the rotated and translated second atom corresponds to. 533 ! 534 USE symm_base, ONLY : s,ft,nsym,irt 535 USE ions_base, ONLY : nat,ityp 536 USE cell_base, ONLY : bg 537 USE fft_base, ONLY : dfftp 538 USE ldaU, ONLY : atom_pos, at_sc, sc_at, num_uc 539 USE kinds 540 USE io_global, ONLY : stdout 541 ! 542 IMPLICIT NONE 543 INTEGER, INTENT(IN) :: at1, at2, p_sym 544 INTEGER, INTENT(OUT) :: rat1, rat2 545 ! rat1 = atom equiv to at1 after 'rotation' and eventual integer translation 546 ! rat2 = atom equiv to at2 after 'rotation' and the same translation above. 547 ! 548 ! Local variables 549 ! 550 INTEGER :: i, j, at, dr(3), equiv_2, dimn 551 ! dr(1), dr(2), dr(3) = location of the unit cell where at2 goes after sym. operation 552 ! 553 REAL(DP) :: diff, x2(3), r1(3), r2(3), dx(3), ss(3,3) 554 REAL(DP), PARAMETER :: eps = 5.d-6 555 ! 556 ! Number of atoms in the supercell 557 dimn = num_uc * nat 558 ! 559 ! Convert the symmetry matrix from integer to real type 560 ! for a given symmetry operation p_sym 561 ! 562 DO i = 1, 3 563 DO j = 1, 3 564 ss(i,j) = DBLE(s(i,j,p_sym)) 565 ENDDO 566 ENDDO 567 ! 568 ! Determine the coordinates of at2 in the supercell 569 ! 570 DO i = 1, 3 571 x2(i) = atom_pos(at_sc(at2)%at,i) + DBLE(at_sc(at2)%n(i)) 572 ENDDO 573 ! 574 ! Apply the symmetry operation: S*at - f 575 ! 576 DO i = 1, 3 577 ! 578 r1(i) = 0.d0 579 r2(i) = 0.d0 580 ! 581 ! Rotate at1 and at2 582 ! 583 DO j = 1, 3 584 r1(i) = r1(i) + atom_pos(at1,j) * ss(j,i) 585 r2(i) = r2(i) + x2(j) * ss(j,i) 586 ENDDO 587 ! 588 ! Subtract vectors of fractional translations 589 ! 590 r1(i) = r1(i) - ft(i, p_sym) 591 r2(i) = r2(i) - ft(i, p_sym) 592 ! 593 ENDDO 594 ! 595 ! Atom equivalent to r2 : irt(p_sym,at2) 596 ! 597 equiv_2 = at_sc(at2)%at 598 ! 599 diff = 1.d0 600 at = 1 601 ! 602 DO WHILE ( at.LE.nat .AND. diff > eps ) 603 ! 604 IF ( ityp(at).EQ.ityp(equiv_2) ) THEN 605 diff = 0.d0 606 DO i = 1, 3 607 dx(i) = r2(i) - atom_pos(at,i) 608 diff = diff + ABS(dx(i)-NINT(dx(i))) 609 ENDDO 610 ELSE 611 diff = 1.d0 612 ENDIF 613 ! 614 at = at + 1 615 ! 616 ENDDO 617 ! 618 IF ( diff > eps ) THEN 619 WRITE(stdout,*) "diff > 0, diff= ", diff, "at1= ", at1, "at2= ", at2 620 CALL errore('symonpair', 'No atom equivalent to r2', 1) 621 ENDIF 622 ! 623 rat2 = at - 1 624 ! 625 ! Atom equivalent to r1 : irt(p_sym,at1) 626 ! 627 diff = 1.d0 628 at = 1 629 ! 630 DO WHILE ( at.LE.nat .AND. diff > eps ) 631 ! 632 IF ( ityp(at).EQ.ityp(at1) ) THEN 633 diff = 0.d0 634 DO i = 1, 3 635 dx(i) = r1(i) - atom_pos(at,i) 636 diff = diff + ABS(dx(i)-NINT(dx(i))) 637 ENDDO 638 ELSE 639 diff = 1.d0 640 ENDIF 641 ! 642 at = at + 1 643 ! 644 ENDDO 645 ! 646 IF ( diff > eps ) THEN 647 WRITE(stdout,*) "diff > 0, diff= ", diff, "at1= ", at1, "at2= ", at2 648 CALL errore('symonpair', 'No atom equivalent to r1', 1) 649 ENDIF 650 ! 651 rat1 = at - 1 652 ! 653 IF (rat1 > nat .OR. rat1 < 1) THEN 654 WRITE(stdout,*) "Index of the first rotated atom=", rat1 655 WRITE(stdout,*) "Number of atoms in the original unit cell=", nat 656 CALL errore('symonpair', 'Out of bounds', 1) 657 ENDIF 658 ! 659 DO i = 1, 3 660 r2(i) = r2(i) - dx(i) 661 dr(i) = NINT(r2(i) - atom_pos(rat2,i)) 662 ENDDO 663 ! 664 rat2 = sc_at(rat2, dr(1), dr(2), dr(3)) 665 ! 666 IF (rat2 > dimn) THEN 667 WRITE(stdout,*) "Index of the second rotated atom=", rat2 668 WRITE(stdout,*) "Number of atoms in the supercell=", dimn 669 WRITE(stdout,*) "Probably a larger sc_size is needed" 670 CALL errore('symonpair', 'Out of bounds', 1) 671 ELSEIF (rat2 < 1) THEN 672 WRITE(stdout,*) "Index of the second rotated atom=", rat2 673 CALL errore('symonpair', 'Out of bounds', 1) 674 ENDIF 675 ! 676 RETURN 677 ! 678END SUBROUTINE symonpair 679!----------------------------------------------------------------------- 680 681!----------------------------------------------------------------------- 682FUNCTION find_viz (center, atom) 683 !------------------------------------------------------------------------ 684 ! 685 ! This function finds by which number the atom 'atom' (in the full supercell 686 ! list) is indexed amongst the neighbour atoms of the atom 'center' (in the 687 ! reference unit cell) 688 ! 689 USE ldaU, ONLY : neighood 690 ! 691 IMPLICIT NONE 692 INTEGER :: center, atom 693 INTEGER :: i, find_viz 694 ! 695 i = 1 696 ! 697 DO WHILE ( i.LE.neighood(center)%num_neigh .AND. & 698 neighood(center)%neigh(i).NE.atom ) 699 i = i + 1 700 ENDDO 701 ! 702 IF ( i.LE.neighood(center)%num_neigh ) THEN 703 find_viz = i 704 ELSE 705 find_viz = -1 706 WRITE(*,*) "find_viz(",center,atom,")", neighood(center)%num_neigh, i 707 CALL errore('find_viz', 'atom is not neighbour of center', 1) 708 ENDIF 709 ! 710 RETURN 711 ! 712END FUNCTION find_viz 713!----------------------------------------------------------------------- 714 715!----------------------------------------------------------------------- 716SUBROUTINE phase_factor (ik) 717 !------------------------------------------------------------------------ 718 ! 719 ! This routine computes the phases for each atom at a given k point 'ik' 720 ! 721 USE kinds, ONLY : DP 722 USE klist, ONLY : xk 723 USE ions_base, ONLY : nat, ityp 724 USE cell_base, ONLY : at, tpiba 725 USE ldaU, ONLY : at_sc, ldim_u, neighood, phase_fac, num_uc 726 USE constants, ONLY : tpi 727 ! 728 IMPLICIT NONE 729 INTEGER, INTENT(IN) :: ik 730 ! 731 ! Local variables 732 ! 733 INTEGER :: na1, viz, na2, nt1, nt2, i, j, dimn 734 REAL(DP) :: angle, sum_j 735 ! 736 ! Number of atoms in the supercell 737 dimn = num_uc * nat 738 ! 739 IF (.NOT.ALLOCATED(phase_fac)) ALLOCATE(phase_fac(dimn)) 740 ! 741 DO na1 = 1, nat 742 ! 743 nt1 = ityp(na1) 744 ! 745 IF ( ldim_u(nt1).GT.0 ) THEN 746 ! 747 DO viz = 1, neighood(na1)%num_neigh 748 ! 749 na2 = neighood(na1)%neigh(viz) 750 angle = 0.d0 751 ! 752 DO i = 1, 3 753 sum_j = 0.d0 754 DO j = 1, 3 755 sum_j = sum_j + DBLE(at_sc(na2)%n(j)) * at(i,j) 756 ENDDO 757 angle = angle + sum_j * xk(i,ik) 758 ENDDO 759 ! 760 angle = tpi * angle 761 phase_fac(na2) = CMPLX(COS(angle),SIN(angle),kind=DP) 762 ! 763 ENDDO 764 ! 765 ENDIF 766 ! 767 ENDDO 768 ! 769 RETURN 770 ! 771END SUBROUTINE phase_factor 772!----------------------------------------------------------------------- 773 774!------------------------------------------------------------------------- 775SUBROUTINE alloc_atom_pos() 776 !----------------------------------------------------------------------- 777 ! 778 ! This routine computes the atomic position in the primitive basis coordinates 779 ! 780 USE ions_base, ONLY : nat, tau 781 USE cell_base, ONLY : bg 782 USE ldaU, ONLY : atom_pos 783 ! 784 IMPLICIT NONE 785 INTEGER :: na, ipol 786 ! counters on atoms and coordinates 787 ! 788 ALLOCATE(atom_pos(nat,3)) 789 ! 790 DO na = 1, nat 791 DO ipol = 1, 3 792 atom_pos(na,ipol) = bg(1,ipol)*tau(1,na) + & 793 bg(2,ipol)*tau(2,na) + & 794 bg(3,ipol)*tau(3,na) 795 ENDDO 796 ENDDO 797 ! 798 RETURN 799 ! 800END SUBROUTINE alloc_atom_pos 801!------------------------------------------------------------------------- 802 803!------------------------------------------------------------------------- 804SUBROUTINE write_V 805 !----------------------------------------------------------------------- 806 ! 807 ! This routine writes Hubbard_V to file 808 ! 809 USE io_files, ONLY : seqopn 810 USE io_global, ONLY : ionode 811 USE ldaU, ONLY : Hubbard_V 812 ! 813 IMPLICIT NONE 814 INTEGER :: i, j, k, iunit 815 LOGICAL :: exst 816 INTEGER, EXTERNAL :: find_free_unit 817 ! 818 iunit = find_free_unit() 819 ! 820 IF (ionode) THEN 821 CALL seqopn( iunit, 'HubbardV.txt', 'FORMATTED', exst ) 822 DO i = 1, SIZE(Hubbard_V,1) 823 DO j = 1, SIZE(Hubbard_V,2) 824 DO k = 1, SIZE(Hubbard_V,3) 825 IF (Hubbard_V(i,j,k) > 1.d-20) WRITE(iunit,*) i, j, k, Hubbard_V(i,j,k) 826 ENDDO 827 ENDDO 828 ENDDO 829 CLOSE(UNIT=iunit, STATUS='KEEP') 830 ENDIF 831 ! 832 RETURN 833 ! 834END SUBROUTINE write_V 835!------------------------------------------------------------------------- 836 837!------------------------------------------------------------------------- 838SUBROUTINE read_V 839 !----------------------------------------------------------------------- 840 ! 841 ! This routine reads Hubbard_V from file 842 ! 843 USE kinds, ONLY : DP 844 USE io_files, ONLY : seqopn 845 USE io_global, ONLY : ionode, ionode_id 846 USE ldaU, ONLY : Hubbard_V 847 USE mp_images, ONLY : intra_image_comm 848 USE mp, ONLY : mp_bcast 849 ! 850 IMPLICIT NONE 851 REAL(DP) :: V 852 INTEGER :: i, j, k, iunit, ierr 853 LOGICAL :: exst 854 INTEGER, EXTERNAL :: find_free_unit 855 ! 856 iunit = find_free_unit() 857 ! 858 Hubbard_V(:,:,:) = 0.0d0 859 ! 860 IF (ionode) THEN 861 CALL seqopn( iunit, 'HubbardV.txt', 'FORMATTED', exst ) 862 IF (exst) THEN 86310 READ(iunit,*,END=11,IOSTAT=ierr) i, j, k, V 864 IF ( ierr/=0 ) THEN 865 CALL mp_bcast( ierr, ionode_id, intra_image_comm ) 866 CALL errore('read_V', 'Reading Hubbard_V', 1) 867 ENDIF 868 Hubbard_V(i,j,k) = V 869 GO TO 10 870 ELSE 871 CALL errore('read_V','File HubbardV.txt was not found...',1) 872 ENDIF 87311 CLOSE( UNIT=iunit, STATUS='KEEP' ) 874 ENDIF 875 ! Broadcast Hubbard_V across all processors 876 CALL mp_bcast( Hubbard_V, ionode_id, intra_image_comm ) 877 ! 878 RETURN 879 ! 880END SUBROUTINE read_V 881!------------------------------------------------------------------------- 882