coulombmatrix.F90 84.2 KB
 Daniel Wortmann committed Jun 21, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 ``````! ! Calculates the Coulomb matrix ! ! v = < M | v | M > ! k,IJ k,I k,J ! ! with the mixed-basis functions M (indices I and J). ! ! Note that ! * ! v = v . ! k,JI k,IJ ! ! In the code: coulomb(IJ,k) = v where only the upper triangle (I<=J) is stored. ! k,IJ ! ! The Coulomb matrix v(IJ,k) diverges at the Gamma-point. Here, we apply the decomposition ! ! (0) (1) * 2-l (0)* (0) (1)* m (1) ! v = v + SUM v * Y (k) / k with v = v , v = (-1) v ! k,IJ IJ lm IJ lm JI IJ JI,lm IJ,l,-m ! ! where a = atom index, R = position vector, T = Wigner-Seitz radius (scalar). ! a 0 ! (0) ! In the code: coulomb(IJ,1) = v where only the upper triangle (I<=J) is stored, ! IJ ! (1) `````` Daniel Wortmann committed Jul 14, 2017 29 ``````! coulfac(IJ,lm) = v IJ,lm `````` Daniel Wortmann committed Jun 21, 2017 30 31 32 33 ``````! ! For the PW contribution we have to construct plane waves within the MT spheres with the help ! of spherical Bessel functions. The value lexp (LEXP in gwinp) is the corresponding cutoff. ! `````` Daniel Wortmann committed Jul 14, 2017 34 ``````MODULE m_coulomb `````` Daniel Wortmann committed Jun 21, 2017 35 `````` `````` Daniel Wortmann committed Jul 14, 2017 36 ``````CONTAINS `````` Daniel Wortmann committed Jun 21, 2017 37 38 `````` `````` Daniel Wortmann committed Jul 14, 2017 39 `````` SUBROUTINE coulombmatrix(mpi,atoms,kpts,cell, sym, hybrid, xcpot,l_restart) `````` Daniel Wortmann committed Jun 21, 2017 40 `````` `````` Daniel Wortmann committed Jul 14, 2017 41 42 43 44 45 46 47 48 `````` USE m_constants , ONLY : pi_const USE m_olap , ONLY : olap_pw,gptnorm USE m_trafo , ONLY : symmetrize,bramat_trafo USE m_util , ONLY : sphbessel,intgrf,intgrf_init, harmonicsr,primitivef USE m_hsefunctional, ONLY : change_coulombmatrix USE m_wrapper USE m_icorrkeys USE m_io_hybrid `````` Daniel Wortmann committed Jun 21, 2017 49 50 `````` USE m_types `````` Daniel Wortmann committed Jul 14, 2017 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 `````` IMPLICIT NONE TYPE(t_xcpot),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_sym),INTENT(IN) :: sym TYPE(t_cell),INTENT(IN) :: cell TYPE(t_kpts),INTENT(IN) :: kpts TYPE(t_atoms),INTENT(IN) :: atoms ! - scalars - LOGICAL , INTENT(IN) :: l_restart ! - local scalars - INTEGER :: maxfac INTEGER :: inviop INTEGER :: nqnrm,iqnrm,iqnrm1,iqnrm2, iqnrmstart,iqnrmstep INTEGER :: itype,l ,ix,iy,iy0,i,j,lm,l1,l2,m1, m2,ineq,idum,ikpt,ikpt0,ikpt1 INTEGER :: lm1,lm2,itype1,itype2,ineq1,ineq2,n, n1,n2,ng INTEGER :: ic,ic1,ic2,ic3,ic4,ic5,ic6,ic7,ic8 INTEGER :: igpt,igpt1,igpt2,igptp,igptp1,igptp2 INTEGER :: isym,isym1,isym2,igpt0 INTEGER :: maxlcut,ok INTEGER :: m INTEGER :: ikptmin,ikptmax,nkminmax LOGICAL :: lsym REAL :: rdum,rdum1,rdum2 REAL :: svol,qnorm,qnorm1,qnorm2,gnorm REAL :: fcoulfac REAL :: time1,time2 COMPLEX :: cdum,cdum1,cexp,csum COMPLEX , PARAMETER :: img = (0d0,1d0) ! - local arrays - INTEGER :: g(3) INTEGER :: nbasm1(kpts%nkptf) INTEGER , ALLOCATABLE :: pqnrm(:,:) INTEGER :: rrot(3,3,sym%nsym),invrrot(3,3,sym%nsym) INTEGER , ALLOCATABLE :: iarr(:),POINTER(:,:,:,:)!,pointer(:,:,:) INTEGER :: igptmin(kpts%nkpt),igptmax(kpts%nkpt) INTEGER , ALLOCATABLE :: nsym_gpt(:,:), sym_gpt(:,:,:) INTEGER :: nsym1(kpts%nkpt+1), sym1(sym%nsym,kpts%nkpt+1) LOGICAL :: calc_mt(kpts%nkpt) REAL :: q(3),q1(3),q2(3) REAL :: integrand(atoms%jmtd),primf1(atoms%jmtd), primf2(atoms%jmtd) REAL :: mat(hybrid%maxindxm1*(hybrid%maxindxm1+1)/2) REAL :: moment(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype), moment2(hybrid%maxindxm1,atoms%ntype) REAL :: sphbes(atoms%jmtd,0:hybrid%maxlcutm1) REAL :: sphbesmoment1(atoms%jmtd,0:hybrid%maxlcutm1) REAL :: rarr(0:hybrid%lexp+1),rarr1(0:hybrid%maxlcutm1) REAL , ALLOCATABLE :: fac(:),sfac(:),facfac(:) REAL , ALLOCATABLE :: gmat(:,:),qnrm(:) REAL , ALLOCATABLE :: sphbesmoment(:,:,:) REAL , ALLOCATABLE :: sphbes0(:,:,:) REAL , ALLOCATABLE :: olap(:,:,:,:),integral(:,:,:,:) REAL , ALLOCATABLE :: gridf(:,:) COMPLEX :: cexp1(atoms%ntype),csumf(9) COMPLEX :: structconst((2*hybrid%lexp+1)**2 ,atoms%nat,atoms%nat, kpts%nkpt) ! nw = 1 COMPLEX :: y((hybrid%lexp+1)**2),y1((hybrid%lexp+1)**2), y2((hybrid%lexp+1)**2) COMPLEX :: dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1,sym%nsym) COMPLEX , ALLOCATABLE :: smat(:,:) COMPLEX , ALLOCATABLE :: coulmat(:,:) COMPLEX , ALLOCATABLE :: carr2(:,:),carr2a(:,:), carr2b(:,:) COMPLEX , ALLOCATABLE :: structconst1(:,:) REAL , ALLOCATABLE :: coulomb_mt1(:,:,:,:,:) !REAL , ALLOCATABLE :: coulomb(:,:) !At the moment we always calculate a complex coulomb matrix REAL , ALLOCATABLE :: coulomb_mt2_r(:,:,:,:,:), coulomb_mt3_r(:,:,:,:) REAL , ALLOCATABLE :: coulomb_mtir_r(:,:,:), coulombp_mtir_r(:,:) COMPLEX , ALLOCATABLE :: coulomb(:,:) COMPLEX , ALLOCATABLE :: coulomb_mt2_c(:,:,:,:,:), coulomb_mt3_c(:,:,:,:) COMPLEX , ALLOCATABLE :: coulomb_mtir_c(:,:,:), coulombp_mtir_c(:,:) INTEGER :: ishift,ishift1 INTEGER :: iatom,iatom1 INTEGER :: indx1,indx2,indx3,indx4 LOGICAL :: l_found,l_warn,l_warned, l_plot = .FALSE.!.true.!.false. TYPE(t_mat) :: olapm,coulhlp `````` Daniel Wortmann committed Jun 21, 2017 136 137 138 `````` `````` Daniel Wortmann committed Jul 14, 2017 139 140 `````` CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,gridf) `````` Daniel Wortmann committed Jun 21, 2017 141 `````` `````` Daniel Wortmann committed Jul 14, 2017 142 143 `````` nbasm1 = hybrid%nbasp + hybrid%ngptm(:) `````` Daniel Wortmann committed Jun 21, 2017 144 `````` `````` Daniel Wortmann committed Jul 14, 2017 145 146 147 148 `````` svol = SQRT(cell%vol) fcoulfac = 4*pi_const/cell%vol maxlcut = MAXVAL( atoms%lmax ) maxfac = MAX(2*maxlcut+hybrid%maxlcutm1+1,4*MAX(hybrid%maxlcutm1,hybrid%lexp)+1) `````` Daniel Wortmann committed Jun 21, 2017 149 `````` `````` Daniel Wortmann committed Jul 14, 2017 150 151 152 153 154 155 156 157 158 `````` ALLOCATE ( fac( 0:maxfac),sfac( 0:maxfac),facfac(-1:maxfac) ) fac(0) = 1 ! sfac(0) = 1 ! Define: facfac(-1:0) = 1 ! fac(i) = i! DO i=1,maxfac ! sfac(i) = sqrt(i!) fac(i) = fac(i-1)*i ! facfac(i) = (2i+1)!! sfac(i) = sfac(i-1)*SQRT(i*1d0) ! facfac(i) = facfac(i-1)*(2*i+1) ! END DO `````` Daniel Wortmann committed Jun 21, 2017 159 160 161 `````` `````` Daniel Wortmann committed Jul 14, 2017 162 163 `````` ! Calculate the structure constant CALL structureconstant(structconst,cell,hybrid, atoms,kpts, mpi) `````` Daniel Wortmann committed Jun 21, 2017 164 165 166 167 `````` `````` Daniel Wortmann committed Jul 14, 2017 168 `````` IF ( mpi%irank == 0 ) WRITE(6,'(//A)') '### subroutine: coulombmatrix ###' `````` Daniel Wortmann committed Jun 21, 2017 169 `````` `````` Daniel Wortmann committed Jul 14, 2017 170 171 172 `````` ! ! Matrix allocation ! `````` Daniel Wortmann committed Jun 21, 2017 173 `````` `````` Daniel Wortmann committed Jul 14, 2017 174 `````` IF(ALLOCATED(coulomb)) DEALLOCATE (coulomb) `````` Daniel Wortmann committed Jun 21, 2017 175 `````` `````` Daniel Wortmann committed Jul 14, 2017 176 177 178 `````` ALLOCATE ( coulomb(hybrid%maxbasm1*(hybrid%maxbasm1+1)/2,kpts%nkpt) , stat = ok ) IF( ok .NE. 0 ) STOP 'coulombmatrix: failure allocation coulomb matrix' coulomb = 0 `````` Daniel Wortmann committed Jun 21, 2017 179 `````` `````` Daniel Wortmann committed Jul 14, 2017 180 `````` IF ( mpi%irank == 0 ) WRITE(6,'(/A,F6.1," MB")') 'Size of coulomb matrix:',16d0/1048576*SIZE(coulomb) `````` Daniel Wortmann committed Jun 21, 2017 181 `````` `````` Daniel Wortmann committed Jul 14, 2017 182 183 184 185 `````` ! ! Generate Symmetry: ! Reduce list of g-Points so that only one of each symm-equivalent is calculated ! `````` Daniel Wortmann committed Jun 21, 2017 186 187 ``````# ifndef CPP_NOCOULSYM `````` Daniel Wortmann committed Jul 14, 2017 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 `````` IF ( mpi%irank == 0 ) WRITE(6,'(/A)',advance='no') 'Setup for symmetry...' CALL cpu_TIME(time1) ! calculate rotations in reciprocal space DO isym = 1,sym%nsym IF( isym .LE. sym%nop ) THEN inviop = sym%invtab(isym) rrot(:,:,isym) = TRANSPOSE(sym%mrot(:,:,inviop)) DO l = 0,hybrid%maxlcutm1 dwgn(:,:,l,isym) = TRANSPOSE(hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1,-hybrid%maxlcutm1:hybrid%maxlcutm1,l,isym) ) END DO ELSE inviop = isym - sym%nop rrot(:,:,isym) = -rrot(:,:,inviop) dwgn(:,:,:,isym) = dwgn(:,:,:,inviop) DO l = 0,hybrid%maxlcutm1 DO m1 = -l,l `````` Daniel Wortmann committed Jun 21, 2017 204 `````` DO m2 = -l,-1 `````` Daniel Wortmann committed Jul 14, 2017 205 206 207 `````` cdum = dwgn(m1, m2,l,isym) dwgn(m1, m2,l,isym) = dwgn(m1,-m2,l,isym) * (-1)**m2 dwgn(m1,-m2,l,isym) = cdum * (-1)**m2 `````` Daniel Wortmann committed Jun 21, 2017 208 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 209 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 210 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 `````` END IF END DO invrrot(:,:,:sym%nop) = rrot(:,:,sym%invtab) IF (sym%nsym > sym%nop) THEN invrrot(:,:,sym%nop+1:) = rrot(:,:,sym%invtab+sym%nop) END IF ! Get symmetry operations that leave bk(:,ikpt) invariant -> sym1 nsym1 = 0 DO ikpt = 1,kpts%nkpt isym1 = 0 DO isym = 1,sym%nsym ! temporary fix until bramat_trafo is correct ! for systems with symmetries including translations IF ( isym > sym%nop ) THEN isym2 = isym-sym%nop ELSE isym2 = isym END IF IF ( ANY(sym%tau(:,isym2) /= 0) ) CYCLE IF(ALL(ABS(MATMUL(rrot(:,:,isym),kpts%bk(:,ikpt)) -kpts%bk(:,ikpt)).LT.1d-12)) THEN isym1 = isym1 + 1 sym1(isym1,ikpt) = isym END IF END DO nsym1(ikpt) = isym1 END DO ! Define reduced lists of G points -> pgptm1(:,ikpt), ikpt=1,..,nkpt !ALLOCATE ( hybrid%pgptm1(hybrid%maxgptm,kpts%nkpt)) !in mixedbasis ALLOCATE (iarr(hybrid%maxgptm), POINTER(kpts%nkpt,& MINVAL(hybrid%gptm(1,:))-1:MAXVAL(hybrid%gptm(1,:))+1,& MINVAL(hybrid%gptm(2,:))-1:MAXVAL(hybrid%gptm(2,:))+1,& MINVAL(hybrid%gptm(3,:))-1:MAXVAL(hybrid%gptm(3,:))+1)) hybrid%pgptm1 = 0 ; iarr = 0 ; POINTER = 0 DO ikpt = 1,kpts%nkpt DO igpt = 1,hybrid%ngptm(ikpt) g = hybrid%gptm(:,hybrid%pgptm(igpt,ikpt)) POINTER(ikpt,g(1),g(2),g(3)) = igpt END DO iarr = 0 j = 0 DO igpt = hybrid%ngptm(ikpt),1,-1 IF (iarr(igpt).EQ.0) THEN j = j + 1 hybrid%pgptm1(j,ikpt) = igpt DO isym1 = 1,nsym1(ikpt) g = MATMUL ( rrot(:,:,sym1(isym1,ikpt)) , hybrid%gptm(:,hybrid%pgptm(igpt,ikpt)) ) i = POINTER(ikpt,g(1),g(2),g(3)) IF(i.EQ.0) STOP 'coulombmatrix: zero pointer (bug?)' `````` Daniel Wortmann committed Jun 21, 2017 261 `````` iarr(i) = 1 `````` Daniel Wortmann committed Jul 14, 2017 262 263 264 265 266 267 `````` END DO END IF END DO hybrid%ngptm1(ikpt) = j END DO DEALLOCATE ( iarr ) `````` Daniel Wortmann committed Jun 21, 2017 268 `````` `````` Daniel Wortmann committed Jul 14, 2017 269 270 271 `````` IF ( mpi%irank == 0 ) WRITE(6,'(12X,A)',advance='no') 'done' CALL cpu_TIME(time2) IF ( mpi%irank == 0 ) WRITE(6,'(2X,A,F8.2,A)') '( Timing:', time2-time1, ' )' `````` Daniel Wortmann committed Jun 21, 2017 272 `````` `````` Daniel Wortmann committed Jul 14, 2017 273 `````` ! no symmetry used `````` Daniel Wortmann committed Jun 21, 2017 274 275 ``````# else `````` Daniel Wortmann committed Jul 14, 2017 276 277 278 279 280 `````` ALLOCATE ( hybrid%pgptm1(hybrid%maxgptm,kpts%nkpt) ) DO ikpt = 1,kpts%nkpt hybrid%pgptm1(:,ikpt) = (/ (igpt0, igpt0 = 1,hybrid%maxgptm) /) hybrid%ngptm1(ikpt) = hybrid%ngptm(ikpt) END DO `````` Daniel Wortmann committed Jun 21, 2017 281 282 283 `````` # endif `````` Daniel Wortmann committed Jul 14, 2017 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 `````` ! Distribute the work as equally as possible over the processes ikptmin = 1 ikptmax = kpts%nkpt igptmin = 1 igptmax = hybrid%ngptm1(:kpts%nkpt) calc_mt = .TRUE. nkminmax = kpts%nkpt IF ( mpi%irank == 0 ) WRITE(6,'(A)',advance='no') 'Preparations...' CALL cpu_TIME(time1) ! Define gmat (symmetric) i = (hybrid%lexp+1)**2 ALLOCATE ( gmat(i,i) ) gmat = 0 lm1 = 0 DO l1=0,hybrid%lexp DO m1=-l1,l1 `````` Daniel Wortmann committed Jun 21, 2017 302 303 `````` lm1 = lm1 + 1 lm2 = 0 `````` Daniel Wortmann committed Jul 14, 2017 304 305 306 307 308 309 310 311 312 `````` lp1:DO l2=0,l1 DO m2=-l2,l2 lm2 = lm2 + 1 IF(lm2.GT.lm1) EXIT lp1 ! Don't cross the diagonal! gmat(lm1,lm2) = sfac(l1+l2+m2-m1)*sfac(l1+l2+m1-m2)/& ( sfac(l1+m1)*sfac(l1-m1)*sfac(l2+m2)*sfac(l2-m2) ) /& SQRT(1d0*(2*l1+1)*(2*l2+1)*(2*(l1+l2)+1))*(4*pi_const)**1.5d0 gmat(lm2,lm1) = gmat(lm1,lm2) END DO `````` Daniel Wortmann committed Jun 21, 2017 313 `````` END DO LP1 `````` Daniel Wortmann committed Jul 14, 2017 314 315 316 317 318 `````` END DO END DO ! Calculate moments of MT functions DO itype=1,atoms%ntype DO l=0,hybrid%lcutm1(itype) `````` Daniel Wortmann committed Jun 21, 2017 319 `````` DO i=1,hybrid%nindxm1(l,itype) `````` Daniel Wortmann committed Jul 14, 2017 320 321 322 `````` ! note that hybrid%basm1 already contains the factor rgrid moment(i,l,itype) = intgrf(atoms%rmsh(:,itype)**(l+1)*hybrid%basm1(:,i,l,itype),& atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) `````` Daniel Wortmann committed Jun 21, 2017 323 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 324 325 `````` END DO DO i =1,hybrid%nindxm1(0,itype) `````` Daniel Wortmann committed Jun 21, 2017 326 `````` moment2(i,itype) = intgrf(atoms%rmsh(:,itype)**3*hybrid%basm1(:,i,0,itype),& `````` Daniel Wortmann committed Jul 14, 2017 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 `````` atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) END DO END DO ! Look for different qnorm = |k+G|, definition of qnrm and pqnrm. CALL getnorm(kpts,hybrid%gptm,hybrid%ngptm,hybrid%pgptm, qnrm,nqnrm,pqnrm,cell) ALLOCATE ( sphbesmoment(0:hybrid%lexp,atoms%ntype,nqnrm),& olap(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype,nqnrm),& integral(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype,nqnrm) ) sphbes = 0 sphbesmoment = 0 sphbesmoment1 = 0 olap = 0 integral = 0 ! Calculate moments of spherical Bessel functions (for (2) and (3)) (->sphbesmoment) ! Calculate overlap of spherical Bessel functions with basis functions (for (2)) (->olap) ! Calculate overlap of sphbesmoment1(r,l) with basis functions (for (2)) (->integral) ! We use sphbes(r,l) = j_l(qr) ! and sphbesmoment1(r,l) = 1/r**(l-1) * INT(0..r) r'**(l+2) * j_l(qr') dr' ! + r**(l+2) * INT(r..S) r'**(1-l) * j_l(qr') dr' . iqnrmstart = mpi%irank + 1 iqnrmstep = mpi%isize DO iqnrm = iqnrmstart,nqnrm,iqnrmstep qnorm = qnrm(iqnrm) DO itype = 1,atoms%ntype `````` Daniel Wortmann committed Jun 21, 2017 354 355 356 357 `````` ng = atoms%jri(itype) rdum = atoms%rmt(itype) sphbes = 0 sphbesmoment1 = 0 `````` Daniel Wortmann committed Jul 14, 2017 358 359 360 361 362 363 `````` IF(qnorm.EQ.0) THEN sphbesmoment(0,itype,iqnrm) = rdum**3 / 3 DO i = 1,ng sphbes(i,0) = 1 sphbesmoment1(i,0) = atoms%rmsh(i,itype)**2 / 3 + ( rdum**2 - atoms%rmsh(i,itype)**2 ) / 2 END DO `````` Daniel Wortmann committed Jun 21, 2017 364 `````` ELSE `````` Daniel Wortmann committed Jul 14, 2017 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 `````` CALL sphbessel(rarr,qnorm*rdum,hybrid%lexp+1) DO l = 0,hybrid%lexp sphbesmoment(l,itype,iqnrm) = rdum**(l+2) * rarr(l+1) / qnorm END DO DO i = ng,1,-1 rdum = atoms%rmsh(i,itype) CALL sphbessel(rarr,qnorm*rdum,hybrid%lcutm1(itype)+1) DO l = 0,hybrid%lcutm1(itype) sphbes(i,l) = rarr(l) IF(l.NE.0) THEN ; rdum1 = -rdum**(1-l) * rarr(l-1) ELSE ; rdum1 = -COS(qnorm*rdum) / qnorm ENDIF IF(i.EQ.ng) rarr1(l) = rdum1 sphbesmoment1(i,l) = (rdum**(l+2)*rarr(l+1)/rdum**(l+1)& + ( rarr1(l) - rdum1 ) * rdum**l ) / qnorm END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 382 383 `````` END IF DO l = 0,hybrid%lcutm1(itype) `````` Daniel Wortmann committed Jul 14, 2017 384 385 386 387 388 `````` DO n = 1,hybrid%nindxm1(l,itype) ! note that hybrid%basm1 already contains one factor rgrid olap(n,l,itype,iqnrm) = & intgrf(atoms%rmsh(:,itype)*hybrid%basm1(:,n,l,itype)*sphbes(:,l),& atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) `````` Daniel Wortmann committed Jun 21, 2017 389 `````` `````` Daniel Wortmann committed Jul 14, 2017 390 391 392 `````` integral(n,l,itype,iqnrm) = & intgrf(atoms%rmsh(:,itype)*hybrid%basm1(:,n,l,itype)*sphbesmoment1(:,l),& atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) `````` Daniel Wortmann committed Jun 21, 2017 393 `````` `````` Daniel Wortmann committed Jul 14, 2017 394 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 395 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 396 397 `````` END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 398 `````` `````` Daniel Wortmann committed Jul 14, 2017 399 400 401 402 403 404 `````` IF ( mpi%irank == 0 ) THEN WRITE(6,'(18X,A)',advance='no') 'done' CALL cpu_TIME(time2) WRITE(6,'(2X,A,F8.2,A)',advance='no') '( Timing:',time2-time1,' )' WRITE(6,*) END IF `````` Daniel Wortmann committed Jun 21, 2017 405 `````` `````` Daniel Wortmann committed Jul 14, 2017 406 407 408 `````` ! ! (1) Case < MT | v | MT > ! `````` Daniel Wortmann committed Jun 21, 2017 409 `````` `````` Daniel Wortmann committed Jul 14, 2017 410 `````` IF( mpi%irank == 0 ) WRITE(6,'(A)',advance='no') '< MT | v | MT > contribution...' `````` Daniel Wortmann committed Jun 21, 2017 411 `````` `````` Daniel Wortmann committed Jul 14, 2017 412 `````` CALL cpu_TIME(time1) `````` Daniel Wortmann committed Jun 21, 2017 413 `````` `````` Daniel Wortmann committed Jul 14, 2017 414 `````` IF ( ANY( calc_mt ) ) THEN `````` Daniel Wortmann committed Jun 21, 2017 415 `````` `````` Daniel Wortmann committed Jul 14, 2017 416 `````` ! (1a) r,r' in same MT `````` Daniel Wortmann committed Jun 21, 2017 417 `````` `````` Daniel Wortmann committed Jul 14, 2017 418 419 420 421 `````` ix = 0 iy = 0 iy0 = 0 DO itype=1,atoms%ntype `````` Daniel Wortmann committed Jun 21, 2017 422 `````` DO ineq=1,atoms%neq(itype) `````` Daniel Wortmann committed Jul 14, 2017 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 `````` ! Here the diagonal block matrices do not depend on ineq. In (1b) they do depend on ineq, though, DO l=0,hybrid%lcutm1(itype) DO n2=1,hybrid%nindxm1(l,itype) ! note that hybrid%basm1 already contains the factor rgrid CALL primitivef(primf1,hybrid%basm1(:,n2,l,itype)& *atoms%rmsh(:,itype)**(l+1),atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,itype,atoms%ntype) ! -itype is to enforce inward integration CALL primitivef(primf2,hybrid%basm1(:,n2,l,itype)& /atoms%rmsh(:,itype)**l,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,-itype,atoms%ntype) primf1 = primf1 / atoms%rmsh(:,itype)**l primf2 = primf2 * atoms%rmsh(:,itype)**(l+1) DO n1=1,n2 integrand = hybrid%basm1(:,n1,l,itype) * (primf1 + primf2) ! call intgr0( (4*pimach())/(2*l+1)*integrand,rmsh(1,itype),dx(itype),jri(itype),mat(n2*(n2-1)/2+n1) ) mat(n2*(n2-1)/2+n1) = (4*pi_const)/(2*l+1)& * intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,& atoms%ntype,itype,gridf) END DO `````` Daniel Wortmann committed Jun 21, 2017 443 444 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 445 446 447 448 449 450 451 452 453 454 455 456 457 `````` ! distribute mat for m=-l,l on coulomb in block-matrix form DO M=-l,l DO n2=1,hybrid%nindxm1(l,itype) ix = ix + 1 iy = iy0 DO n1=1,n2 iy = iy + 1 i = ix*(ix-1)/2+iy j = n2*(n2-1)/2+n1 coulomb(i,kpts%nkpt) = mat(j) END DO END DO iy0 = ix `````` Daniel Wortmann committed Jun 21, 2017 458 459 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 460 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 461 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 462 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 463 `````` `````` Daniel Wortmann committed Jul 14, 2017 464 `````` ! (1b) r,r' in different MT `````` Daniel Wortmann committed Jun 21, 2017 465 `````` `````` Daniel Wortmann committed Jul 14, 2017 466 467 468 `````` ALLOCATE( coulmat(hybrid%nbasp,hybrid%nbasp), stat=ok) IF( ok .NE. 0 ) STOP 'coulombmatrix: failure allocation coulmat' coulmat = 0 `````` Daniel Wortmann committed Jun 21, 2017 469 `````` `````` Daniel Wortmann committed Jul 14, 2017 470 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 471 `````` `````` Daniel Wortmann committed Jul 14, 2017 472 `````` DO ikpt=ikptmin,ikptmax `````` Daniel Wortmann committed Jun 21, 2017 473 `````` `````` Daniel Wortmann committed Jul 14, 2017 474 475 `````` ! only the first rank handles the MT-MT part IF ( calc_mt(ikpt) ) THEN `````` Daniel Wortmann committed Jun 21, 2017 476 477 478 479 `````` ix = 0 ic2 = 0 DO itype2=1,atoms%ntype `````` Daniel Wortmann committed Jul 14, 2017 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 `````` DO ineq2=1,atoms%neq(itype2) ic2 = ic2 + 1 lm2 = 0 DO l2=0,hybrid%lcutm1(itype2) DO m2=-l2,l2 lm2 = lm2 + 1 DO n2=1,hybrid%nindxm1(l2,itype2) ix = ix + 1 iy = 0 ic1 = 0 lp2: DO itype1=1,itype2 DO ineq1=1,atoms%neq(itype1) ic1 = ic1 + 1 lm1 = 0 DO l1=0,hybrid%lcutm1(itype1) DO m1=-l1,l1 lm1 = lm1 + 1 DO n1=1,hybrid%nindxm1(l1,itype1) iy = iy + 1 IF(iy.GT.ix) EXIT lp2 ! Don't cross the diagonal! rdum = (-1)**(l2+m2)* moment(n1,l1,itype1)*moment(n2,l2,itype2)*gmat(lm1,lm2) l = l1 + l2 lm = l**2 + l + m1 - m2 + 1 idum = ix*(ix-1)/2+iy coulmat(iy,ix) = coulomb(idum,kpts%nkpt)& + EXP(img* 2*pi_const * dot_PRODUCT(kpts%bk(:,ikpt),& atoms%taual(:,ic2)-atoms%taual(:,ic1))) *rdum * structconst(lm,ic1,ic2,ikpt) coulmat(ix,iy) = CONJG(coulmat(iy,ix)) END DO END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 513 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 514 `````` END DO lp2 `````` Daniel Wortmann committed Jun 21, 2017 515 `````` `````` Daniel Wortmann committed Jul 14, 2017 516 517 `````` END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 518 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 519 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 520 521 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 522 523 524 525 526 527 `````` IF ( sym%invs) THEN !symmetrize makes the Coulomb matrix real symmetric CALL symmetrize(coulmat,hybrid%nbasp,hybrid%nbasp,3,.FALSE.,& atoms,hybrid%lcutm1,hybrid%maxlcutm1,& hybrid%nindxm1,sym) ENDIF `````` Daniel Wortmann committed Jun 21, 2017 528 `````` `````` Daniel Wortmann committed Jul 14, 2017 529 `````` coulomb(:hybrid%nbasp*(hybrid%nbasp+1)/2,ikpt) = packmat(coulmat) `````` Daniel Wortmann committed Jun 21, 2017 530 `````` `````` Daniel Wortmann committed Jul 14, 2017 531 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 532 `````` `````` Daniel Wortmann committed Jul 14, 2017 533 534 `````` END DO IF ( ANY( calc_mt ) ) DEALLOCATE( coulmat ) `````` Daniel Wortmann committed Jun 21, 2017 535 `````` `````` Daniel Wortmann committed Jul 14, 2017 536 537 538 539 540 541 `````` IF ( mpi%irank == 0 ) THEN WRITE(6,'(2X,A)',advance='no') 'done' CALL cpu_TIME(time2) WRITE(6,'(2X,A,F8.2,A)',advance='no') '( Timing:',time2-time1,' )' WRITE(6,*) END IF `````` Daniel Wortmann committed Jun 21, 2017 542 `````` `````` Daniel Wortmann committed Jul 14, 2017 543 `````` IF(hybrid%maxgptm.EQ.0) GOTO 1 ! skip calculation of plane-wave contribution if mixed basis does not contain plane waves `````` Daniel Wortmann committed Jun 21, 2017 544 `````` `````` Daniel Wortmann committed Jul 14, 2017 545 546 547 `````` ! ! (2) Case < MT | v | PW > ! `````` Daniel Wortmann committed Jun 21, 2017 548 `````` `````` Daniel Wortmann committed Jul 14, 2017 549 `````` IF( mpi%irank == 0 ) WRITE(6,'(A)',advance='no') '< MT | v | PW > contribution...' `````` Daniel Wortmann committed Jun 21, 2017 550 `````` `````` Daniel Wortmann committed Jul 14, 2017 551 `````` CALL cpu_TIME(time1) `````` Daniel Wortmann committed Jun 21, 2017 552 `````` `````` Daniel Wortmann committed Jul 14, 2017 553 554 555 `````` ! (2a) r in MT, r' everywhere ! (2b) r,r' in same MT ! (2c) r,r' in different MT `````` Daniel Wortmann committed Jun 21, 2017 556 `````` `````` Daniel Wortmann committed Jul 14, 2017 557 558 559 `````` ALLOCATE( coulmat(hybrid%nbasp,hybrid%maxgptm), stat=ok ) IF( ok .NE. 0 ) STOP 'coulombmatrix: failure allocation coulmat' coulmat = 0 `````` Daniel Wortmann committed Jun 21, 2017 560 `````` `````` Daniel Wortmann committed Jul 14, 2017 561 `````` DO ikpt = ikptmin,ikptmax !1,kpts%nkpt `````` Daniel Wortmann committed Jun 21, 2017 562 `````` `````` Daniel Wortmann committed Jul 14, 2017 563 `````` coulmat = 0 `````` Daniel Wortmann committed Jun 21, 2017 564 `````` `````` Daniel Wortmann committed Jul 14, 2017 565 566 `````` ! start to loop over interstitial plane waves DO igpt0 = igptmin(ikpt),igptmax(ikpt) !1,hybrid%ngptm1(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 567 568 `````` igpt = hybrid%pgptm1(igpt0,ikpt) igptp = hybrid%pgptm(igpt,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 569 570 571 `````` ix = hybrid%nbasp + igpt q = MATMUL ( kpts%bk(:,ikpt) + hybrid%gptm(:,igptp), cell%bmat ) qnorm = SQRT(SUM(q**2)) `````` Daniel Wortmann committed Jun 21, 2017 572 `````` iqnrm = pqnrm(igpt,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 573 `````` IF(ABS(qnrm(iqnrm)-qnorm).GT.1d-12) STOP 'coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)' ! We shouldn't stop here! `````` Daniel Wortmann committed Jun 21, 2017 574 `````` `````` Daniel Wortmann committed Jul 14, 2017 575 576 `````` CALL harmonicsr(y1,MATMUL(kpts%bk(:,kpts%nkpt),cell%bmat),2) CALL harmonicsr(y2,MATMUL(hybrid%gptm(:,igptp) ,cell%bmat),2) `````` Daniel Wortmann committed Jun 21, 2017 577 `````` CALL harmonicsr(y,q,hybrid%lexp) `````` Daniel Wortmann committed Jul 14, 2017 578 `````` y1 = CONJG(y1) ; y2 = CONJG(y2) ; y = CONJG(y) `````` Daniel Wortmann committed Jun 21, 2017 579 580 581 582 `````` iy = 0 ic = 0 DO itype = 1,atoms%ntype `````` Daniel Wortmann committed Jul 14, 2017 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 `````` DO ineq = 1,atoms%neq(itype) ic = ic + 1 lm = 0 DO l = 0,hybrid%lcutm1(itype) DO M = -l,l lm = lm + 1 ! calculate sum over lm and centers for (2c) -> csum, csumf csum = 0 csumf = 0 ic1 = 0 DO itype1=1,atoms%ntype DO ineq1=1,atoms%neq(itype1) ic1 = ic1 + 1 cexp = 4*pi_const * EXP( img * 2*pi_const& * ( dot_PRODUCT( kpts%bk(:,ikpt)+hybrid%gptm(:,igptp),atoms%taual(:,ic1) )& - dot_PRODUCT( kpts%bk(:,ikpt),atoms%taual(:,ic) ) ) ) lm1 = 0 DO l1=0,hybrid%lexp l2 = l + l1 ! for structconst idum = 1 cdum = sphbesmoment(l1,itype1,iqnrm) * img**(l1) * cexp DO m1=-l1,l1 lm1 = lm1 + 1 m2 = M - m1 ! for structconst lm2 = l2**2 + l2 + m2 + 1 ! csum = csum - idum * gmat(lm1,lm) * y(lm1) * cdum * structconst(lm2,ic,ic1,ikpt) idum = -idum ! factor (-1)*(l1+m1) END DO END DO ! add contribution of (2c) to csum and csumf coming from linear and quadratic orders of Y_lm*(G) / G * j_(l+1)(GS) IF(ikpt.EQ.1.AND.l.LE.2) THEN cexp = EXP(img*2*pi_const * dot_PRODUCT(hybrid%gptm(:,igptp),atoms%taual(:,ic1)))& * gmat(lm,1) * 4*pi_const/cell%vol csumf(lm) = csumf(lm) - cexp * SQRT(4*pi_const) *& img**l * sphbesmoment(0,itype1,iqnrm) / facfac(l-1) IF(l.EQ.0) THEN IF(igpt.NE.1) THEN csum = csum - cexp * ( sphbesmoment(0,itype1,iqnrm)*atoms%rmt(itype1)**2 -& sphbesmoment(2,itype1,iqnrm)*2d0/3 ) / 10 ELSE csum = csum - cexp * atoms%rmt(itype1)**5/30 END IF ELSE IF(l.EQ.1) THEN csum = csum + cexp * img * SQRT(4*pi_const)& * sphbesmoment(1,itype1,iqnrm) * y(lm) / 3 END IF END IF END DO `````` Daniel Wortmann committed Jun 21, 2017 635 636 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 637 638 639 `````` ! add contribution of (2a) to csumf IF(ikpt.EQ.1.AND.igpt.EQ.1.AND.l.LE.2) THEN csumf(lm)=csumf(lm) + (4*pi_const)**2 * img**l / facfac(l) `````` Daniel Wortmann committed Jun 21, 2017 640 641 `````` END IF `````` Daniel Wortmann committed Jul 14, 2017 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 `````` ! finally define coulomb idum = ix*(ix-1)/2 cdum = (4*pi_const)**2 * img**(l) * y(lm) * EXP(img * 2*pi_const& * dot_PRODUCT(hybrid%gptm(:,igptp),atoms%taual(:,ic))) DO n=1,hybrid%nindxm1(l,itype) iy = iy + 1 IF(ikpt.EQ.1.AND.igpt.EQ.1) THEN IF(l.EQ.0) coulmat(iy,ix-hybrid%nbasp) =& - cdum * moment2(n,itype) / 6 / svol ! (2a) coulmat(iy,ix-hybrid%nbasp) = coulmat(iy,ix-hybrid%nbasp)& + ( - cdum / (2*l+1) * integral(n,l,itype,iqnrm)& ! (2b)& + csum * moment(n,l,itype) ) / svol ! (2c) ELSE coulmat(iy,ix-hybrid%nbasp) = & ( cdum * olap(n,l,itype,iqnrm) / qnorm**2 & ! (2a)& - cdum / (2*l+1) * integral(n,l,itype,iqnrm)& ! (2b)& + csum * moment(n,l,itype) ) / svol ! (2c) END IF `````` Daniel Wortmann committed Jun 21, 2017 662 `````` `````` Daniel Wortmann committed Jul 14, 2017 663 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 664 `````` `````` Daniel Wortmann committed Jul 14, 2017 665 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 666 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 667 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 668 669 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 670 671 672 673 674 675 676 677 678 679 680 681 `````` END DO IF ( sym%invs) THEN CALL symmetrize(coulmat,hybrid%nbasp,hybrid%ngptm(ikpt),1,.FALSE.,& atoms,hybrid%lcutm1,hybrid%maxlcutm1, hybrid%nindxm1,sym) ENDIF M = hybrid%nbasp*(hybrid%nbasp+1)/2 DO i=1,hybrid%ngptm(ikpt) DO j=1,hybrid%nbasp+i M = M + 1 IF(j.LE. hybrid%nbasp) coulomb(M,ikpt) = coulmat(j,i) `````` Daniel Wortmann committed Jun 21, 2017 682 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 683 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 684 `````` `````` Daniel Wortmann committed Jul 14, 2017 685 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 686 `````` `````` Daniel Wortmann committed Jul 14, 2017 687 `````` DEALLOCATE( coulmat,olap,integral ) `````` Daniel Wortmann committed Jun 21, 2017 688 `````` `````` Daniel Wortmann committed Jul 14, 2017 689 690 691 692 693 `````` IF ( mpi%irank == 0 ) THEN WRITE(6,'(2X,A)',advance='no') 'done' CALL cpu_TIME(time2) WRITE(6,'(2X,A,F8.2,A)') '( Timing:',time2-time1,' )' END IF `````` Daniel Wortmann committed Jun 21, 2017 694 `````` `````` Daniel Wortmann committed Jul 14, 2017 695 696 697 `````` ! ! (3) Case < PW | v | PW > ! `````` Daniel Wortmann committed Jun 21, 2017 698 `````` `````` Daniel Wortmann committed Jul 14, 2017 699 `````` IF( mpi%irank == 0 ) WRITE(6,'(A)',advance='no') '< PW | v | PW > contribution...' `````` Daniel Wortmann committed Jun 21, 2017 700 `````` `````` Daniel Wortmann committed Jul 14, 2017 701 `````` CALL cpu_TIME(time1) `````` Daniel Wortmann committed Jun 21, 2017 702 `````` `````` Daniel Wortmann committed Jul 14, 2017 703 `````` ! (3a) r,r' everywhere; r everywhere, r' in MT; r in MT, r' everywhere `````` Daniel Wortmann committed Jun 21, 2017 704 `````` `````` Daniel Wortmann committed Jul 14, 2017 705 706 707 708 709 710 `````` CALL cpu_TIME(time1) ! Calculate the hermitian matrix smat(i,j) = sum(a) integral(MT(a)) exp[i(Gj-Gi)r] dr ALLOCATE ( smat(hybrid%gptmd,hybrid%gptmd) ) smat = 0 DO igpt2=1,hybrid%gptmd DO igpt1=1,igpt2 `````` Daniel Wortmann committed Jun 21, 2017 711 712 `````` g = hybrid%gptm(:,igpt2)-hybrid%gptm(:,igpt1) gnorm = gptnorm(g,cell%bmat) `````` Daniel Wortmann committed Jul 14, 2017 713 714 715 716 `````` IF(gnorm.EQ.0) THEN DO itype=1,atoms%ntype smat(igpt1,igpt2) = smat(igpt1,igpt2) + atoms%neq(itype) * 4*pi_const*atoms%rmt(itype)**3/3 END DO `````` Daniel Wortmann committed Jun 21, 2017 717 `````` ELSE `````` Daniel Wortmann committed Jul 14, 2017 718 719 720 721 722 723 724 725 726 727 `````` ic = 0 DO itype=1,atoms%ntype rdum = atoms%rmt(itype) * gnorm rdum = 4*pi_const * ( SIN(rdum) - rdum * COS(rdum) ) / gnorm**3 DO ineq=1,atoms%neq(itype) ic = ic + 1 smat(igpt1,igpt2) = smat(igpt1,igpt2)& + rdum * EXP( img * 2*pi_const * dot_PRODUCT(atoms%taual(:,ic),g) ) END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 728 `````` END IF `````` Daniel Wortmann committed Jul 14, 2017 729 730 731 `````` smat(igpt2,igpt1) = CONJG(smat(igpt1,igpt2)) END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 732 `````` `````` Daniel Wortmann committed Jul 14, 2017 733 734 `````` ! Coulomb matrix, contribution (3a) DO ikpt=ikptmin,ikptmax `````` Daniel Wortmann committed Jun 21, 2017 735 `````` `````` Daniel Wortmann committed Jul 14, 2017 736 `````` DO igpt0=igptmin(ikpt),igptmax(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 737 738 `````` igpt2 = hybrid%pgptm1(igpt0,ikpt) igptp2 = hybrid%pgptm(igpt2,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 739 740 741 742 743 `````` ix = hybrid%nbasp + igpt2 iy = hybrid%nbasp q2 = MATMUL ( kpts%bk(:,ikpt) + hybrid%gptm(:,igptp2) , cell%bmat ) rdum2 = SUM(q2**2) IF( rdum2 .NE. 0 ) rdum2 = 4*pi_const/rdum2 `````` Daniel Wortmann committed Jun 21, 2017 744 745 `````` DO igpt1=1,igpt2 `````` Daniel Wortmann committed Jul 14, 2017 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 `````` igptp1 = hybrid%pgptm(igpt1,ikpt) iy = iy + 1 q1 = MATMUL ( kpts%bk(:,ikpt) + hybrid%gptm(:,igptp1) , cell%bmat ) idum = ix*(ix-1)/2+iy rdum1 = SUM(q1**2) IF( rdum1 .NE. 0 ) rdum1 = 4*pi_const/rdum1 IF(ikpt.EQ.1) THEN IF(igpt1.NE.1) THEN coulomb(idum,1) = - smat(igptp1,igptp2) * rdum1 / cell%vol END IF IF(igpt2.NE.1) THEN coulomb(idum,1) = coulomb(idum,1) - smat(igptp1,igptp2) * rdum2 / cell%vol END IF ELSE coulomb(idum,ikpt) = - smat(igptp1,igptp2) * ( rdum1 + rdum2 ) / cell%vol END IF `````` Daniel Wortmann committed Jun 21, 2017 763 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 764 765 `````` IF(ikpt.NE.1.OR.igpt2.NE.1) THEN ! coulomb(idum,ikpt) = coulomb(idum,ikpt) + rdum2 ! diagonal term `````` Daniel Wortmann committed Jun 21, 2017 766 `````` END IF ! `````` Daniel Wortmann committed Jul 14, 2017 767 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 768 `````` `````` Daniel Wortmann committed Jul 14, 2017 769 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 770 `````` `````` Daniel Wortmann committed Jul 14, 2017 771 `````` ! (3b) r,r' in different MT `````` Daniel Wortmann committed Jun 21, 2017 772 `````` `````` Daniel Wortmann committed Jul 14, 2017 773 `````` DO ikpt=ikptmin,ikptmax!1,kpts%nkpt `````` Daniel Wortmann committed Jun 21, 2017 774 `````` `````` Daniel Wortmann committed Jul 14, 2017 775 776 777 778 `````` ! group together quantities which depend only on l,m and igpt -> carr2a ALLOCATE( carr2a((hybrid%lexp+1)**2,hybrid%maxgptm),carr2b(atoms%nat,hybrid%maxgptm) ) carr2a = 0 ; carr2b = 0 DO igpt=1,hybrid%ngptm(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 779 780 `````` igptp = hybrid%pgptm(igpt,ikpt) iqnrm = pqnrm(igpt,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 781 `````` q = MATMUL ( kpts%bk(:,ikpt) + hybrid%gptm(:,igptp),cell%bmat) `````` Daniel Wortmann committed Jun 21, 2017 782 `````` CALL harmonicsr(y,q,hybrid%lexp) `````` Daniel Wortmann committed Jul 14, 2017 783 `````` y = CONJG(y) `````` Daniel Wortmann committed Jun 21, 2017 784 785 `````` lm = 0 DO l=0,hybrid%lexp `````` Daniel Wortmann committed Jul 14, 2017 786 787 788 789 `````` DO M=-l,l lm = lm + 1 carr2a(lm,igpt) = 4*pi_const * img**(l) * y(lm) END DO `````` Daniel Wortmann committed Jun 21, 2017 790 791 `````` END DO DO ic = 1,atoms%nat `````` Daniel Wortmann committed Jul 14, 2017 792 793 `````` carr2b(ic,igpt) = EXP ( -img * 2*pi_const * & dot_PRODUCT(kpts%bk(:,ikpt)+hybrid%gptm(:,igptp),atoms%taual(:,ic)) ) `````` Daniel Wortmann committed Jun 21, 2017 794 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 795 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 796 `````` `````` Daniel Wortmann committed Jul 14, 2017 797 798 799 800 801 `````` !finally we can loop over the plane waves (G: igpt1,igpt2) ALLOCATE ( carr2(atoms%nat,(hybrid%lexp+1)**2),& structconst1(atoms%nat,(2*hybrid%lexp+1)**2) ) carr2 = 0 ; structconst1 = 0 DO igpt0=igptmin(ikpt),igptmax(ikpt)!1,hybrid%ngptm1(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 802 `````` igpt2 = hybrid%pgptm1(igpt0,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 803 `````` ix = hybrid%nbasp + igpt2 `````` Daniel Wortmann committed Jun 21, 2017 804 805 806 807 808 `````` igptp2 = hybrid%pgptm(igpt2,ikpt) iqnrm2 = pqnrm(igpt2,ikpt) ic2 = 0 carr2 = 0 DO itype2 = 1,atoms%ntype `````` Daniel Wortmann committed Jul 14, 2017 809 810 811 812 813 814 `````` DO ineq2 = 1,atoms%neq(itype2) ic2 = ic2 + 1 cexp = CONJG ( carr2b(ic2,igpt2) ) lm2 = 0 DO ic1 = 1,atoms%nat structconst1(ic1,:) = structconst(:,ic1,ic2,ikpt) `````` Daniel Wortmann committed Jun 21, 2017 815 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 `````` DO l2 = 0,hybrid%lexp idum = 1 DO m2 = -l2,l2 lm2 = lm2 + 1 cdum = idum * sphbesmoment(l2,itype2,iqnrm2) * cexp * carr2a(lm2,igpt2) IF( cdum .NE. 0 ) THEN lm1 = 0 DO l1 = 0,hybrid%lexp l = l1 + l2 M = -l1 - m2 !first loop of m1 lm = l**2 + l + M DO m1 = -l1,l1 lm1 = lm1 + 1 lm = lm + 1 cdum1= cdum * gmat(lm1,lm2) DO ic1 = 1,atoms%nat carr2(ic1,lm1) = carr2(ic1,lm1) + cdum1 * structconst1(ic1,lm) END DO END DO END DO END IF idum = -idum !factor (-1)**(l+M) END DO END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 841 842 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 843 `````` iy = hybrid%nbasp `````` Daniel Wortmann committed Jun 21, 2017 844 `````` DO igpt1=1,igpt2 `````` Daniel Wortmann committed Jul 14, 2017 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 `````` iy = iy + 1 igptp1 = hybrid%pgptm(igpt1,ikpt) iqnrm1 = pqnrm(igpt1,ikpt) csum = 0 ic = 0 DO itype=1,atoms%ntype DO ineq=1,atoms%neq(itype) ic = ic + 1 cexp = carr2b(ic,igpt1) lm = 0 DO l=0,hybrid%lexp cdum = cexp * sphbesmoment(l,itype,iqnrm1) DO M=-l,l lm = lm + 1 csum = csum + cdum * carr2(ic,lm) * CONJG ( carr2a(lm,igpt1) ) ! for coulomb END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 862 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 863 864 865 `````` END DO idum = ix*(ix-1)/2+iy coulomb(idum,ikpt) = coulomb(idum,ikpt) + csum / cell%vol `````` Daniel Wortmann committed Jun 21, 2017 866 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 `````` END DO DEALLOCATE( carr2,carr2a,carr2b,structconst1 ) END DO !ikpt ! Add corrections from higher orders in (3b) to coulomb(:,1) ! (1) igpt1 > 1 , igpt2 > 1 (finite G vectors) rdum = (4*pi_const)**(1.5d0)/cell%vol**2 * gmat(1,1) DO igpt0 = 1,hybrid%ngptm1(1) igpt2 = hybrid%pgptm1(igpt0,1) ; IF ( igpt2 == 1 ) CYCLE ix = hybrid%nbasp + igpt2 iqnrm2 = pqnrm(igpt2,1) igptp2 = hybrid%pgptm(igpt2,1) q2 = MATMUL(hybrid%gptm(:,igptp2),cell%bmat) qnorm2 = SQRT(SUM(q2**2)) iy = hybrid%nbasp + 1 DO igpt1 = 2,igpt2 `````` Daniel Wortmann committed Jun 21, 2017 883 884 885 886 `````` iy = iy + 1 idum = ix*(ix-1)/2+iy iqnrm1 = pqnrm(igpt1,1) igptp1 = hybrid%pgptm(igpt1,1) `````` Daniel Wortmann committed Jul 14, 2017 887 888 889 `````` q1 = MATMUL(hybrid%gptm(:,igptp1),cell%bmat) qnorm1 = SQRT(SUM(q1**2)) rdum1 = dot_PRODUCT(q1,q2) / (qnorm1*qnorm2) `````` Daniel Wortmann committed Jun 21, 2017 890 891 `````` ic1 = 0 DO itype1 = 1,atoms%ntype `````` Daniel Wortmann committed Jul 14, 2017 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 `````` DO ineq1 = 1,atoms%neq(itype1) ic1 = ic1 + 1 ic2 = 0 DO itype2 = 1,atoms%ntype DO ineq2 = 1,atoms%neq(itype2) ic2 = ic2 + 1 cdum = EXP ( img * 2*pi_const *& ( - dot_PRODUCT(hybrid%gptm(:,igptp1),atoms%taual(:,ic1))& + dot_PRODUCT(hybrid%gptm(:,igptp2),atoms%taual(:,ic2)) ) ) coulomb(idum,1) = coulomb(idum,1) + rdum * cdum * (& - sphbesmoment(1,itype1,iqnrm1) & * sphbesmoment(1,itype2,iqnrm2) * rdum1 / 3& - sphbesmoment(0,itype1,iqnrm1)& * sphbesmoment(2,itype2,iqnrm2) / 6& - sphbesmoment(2,itype1,iqnrm1)& * sphbesmoment(0,itype2,iqnrm2) / 6 & + sphbesmoment(0,itype1,iqnrm1)& * sphbesmoment(1,itype2,iqnrm2) / qnorm2 / 2& + sphbesmoment(1,itype1,iqnrm1)& * sphbesmoment(0,itype2,iqnrm2) / qnorm1 / 2 ) END DO `````` Daniel Wortmann committed Jun 21, 2017 913 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 914 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 915 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 916 917 918 919 920 921 922 923 924 925 926 927 928 `````` END DO END DO ! (2) igpt1 = 1 , igpt2 > 1 (first G vector vanishes, second finite) iy = hybrid%nbasp + 1 DO igpt0 = 1,hybrid%ngptm1(1) igpt2 = hybrid%pgptm1(igpt0,1) ; IF ( igpt2 == 1 ) CYCLE ix = hybrid%nbasp + igpt2 iqnrm2 = pqnrm(igpt2,1) igptp2 = hybrid%pgptm(igpt2,1) qnorm2 = qnrm(iqnrm2) idum = ix*(ix-1)/2+iy DO itype1 = 1,atoms%ntype `````` Daniel Wortmann committed Jun 21, 2017 929 `````` DO ineq1 = 1,atoms%neq(itype1) `````` Daniel Wortmann committed Jul 14, 2017 930 931 932 933 934 935 936 937 938 939 940 941 `````` ic2 = 0 DO itype2 = 1,atoms%ntype DO ineq2 = 1,atoms%neq(itype2) ic2 = ic2 + 1 cdum = EXP ( img * 2*pi_const * dot_PRODUCT(hybrid%gptm(:,igptp2),atoms%taual(:,ic2)) ) coulomb(idum,1) = coulomb(idum,1)& + rdum * cdum * atoms%rmt(itype1)**3 * (& + sphbesmoment(0,itype2,iqnrm2) / 30 * atoms%rmt(itype1)**2& - sphbesmoment(2,itype2,iqnrm2) / 18 & + sphbesmoment(1,itype2,iqnrm2) / 6 / qnorm2 ) END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 942 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 943 944 945 946 947 948 949 950 951 `````` END DO END DO ! (2) igpt1 = 1 , igpt2 = 1 (vanishing G vectors) iy = hybrid%nbasp + 1 ix = hybrid%nbasp + 1 idum = ix*(ix-1)/2+iy DO itype1 = 1,atoms%ntype DO ineq1 = 1,atoms%neq(itype1) `````` Daniel Wortmann committed Jun 21, 2017 952 `````` DO itype2 = 1,atoms%ntype `````` Daniel Wortmann committed Jul 14, 2017 953 954 955 956 957 `````` DO ineq2 = 1,atoms%neq(itype2) coulomb(idum,1) = coulomb(idum,1)& + rdum * atoms%rmt(itype1)**3 * atoms%rmt(itype2)**3 *& ( atoms%rmt(itype1)**2 + atoms%rmt(itype2)**2 ) / 90 END DO `````` Daniel Wortmann committed Jun 21, 2017 958 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 959 960 `````` END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 961 962 `````` `````` Daniel Wortmann committed Jul 14, 2017 963 `````` ! (3c) r,r' in same MT `````` Daniel Wortmann committed Jun 21, 2017 964 `````` `````` Daniel Wortmann committed Jul 14, 2017 965 966 967 968 969 970 `````` ! Calculate sphbesintegral ALLOCATE ( sphbes0(-1:hybrid%lexp+2,atoms%ntype,nqnrm),& & carr2((hybrid%lexp+1)**2,hybrid%maxgptm) ) sphbes0 = 0 ; carr2 = 0 DO iqnrm = 1,nqnrm DO itype = 1,atoms%ntype `````` Daniel Wortmann committed Jun 21, 2017 971 972 `````` rdum = qnrm(iqnrm) * atoms%rmt(itype) CALL sphbessel(sphbes0(0,itype,iqnrm),rdum,hybrid%lexp+2) `````` Daniel Wortmann committed Jul 14, 2017 973 974 975 `````` IF( rdum.NE.0 ) sphbes0(-1,itype,iqnrm) = COS(rdum)/rdum END DO END DO `````` Daniel Wortmann committed Jun 21, 2017 976 `````` `````` Daniel Wortmann committed Jul 14, 2017 977 978 `````` l_warn = ( mpi%irank == 0 ) DO ikpt=ikptmin,ikptmax!1,nkpt `````` Daniel Wortmann committed Jun 21, 2017 979 `````` `````` Daniel Wortmann committed Jul 14, 2017 980 `````` DO igpt = 1,hybrid%ngptm(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 981 `````` igptp = hybrid%pgptm(igpt,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 982 `````` q = MATMUL ( kpts%bk(:,ikpt) + hybrid%gptm(:,igptp), cell%bmat ) `````` Daniel Wortmann committed Jun 21, 2017 983 `````` CALL harmonicsr(carr2(:,igpt),q,hybrid%lexp) `````` Daniel Wortmann committed Jul 14, 2017 984 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 985 `````` `````` Daniel Wortmann committed Jul 14, 2017 986 `````` DO igpt0=igptmin(ikpt),igptmax(ikpt)!1,hybrid%ngptm1(ikpt) `````` Daniel Wortmann committed Jun 21, 2017 987 `````` igpt2 = hybrid%pgptm1(igpt0,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 988 `````` ix = hybrid%nbasp + igpt2 `````` Daniel Wortmann committed Jun 21, 2017 989 990 `````` igptp2 = hybrid%pgptm(igpt2,ikpt) iqnrm2 = pqnrm(igpt2,ikpt) `````` Daniel Wortmann committed Jul 14, 2017 991 992 993 `````` q2 = MATMUL (kpts%bk(:,ikpt) + hybrid%gptm(:,igptp2),cell%bmat) y2 = CONJG ( carr2(:,igpt2) ) iy = hybrid%nbasp `````` Daniel Wortmann committed Jun 21, 2017 994 `````` DO igpt1=1,igpt2 `````` Daniel Wortmann committed Jul 14, 2017 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 `````` iy = iy + 1 igptp1 = hybrid%pgptm(igpt1,ikpt) iqnrm1 = pqnrm(igpt1,ikpt) q1 = MATMUL (kpts%bk(:,ikpt) + hybrid%gptm(:,igptp1),cell%bmat) y1 = carr2(:,igpt1) cexp1 = 0 ic = 0 DO itype=1,atoms%ntype DO ineq=1,atoms%neq(itype) ic = ic + 1 cexp1(itype) = cexp1(itype) +& EXP(img * 2*pi_const * dot_PRODUCT(& (hybrid%gptm(:,igptp2)-hybrid%gptm(:,igptp1)),atoms%taual(:,ic)) ) ENDDO ENDDO lm = 0 cdum = 0 DO l=0,hybrid%lexp cdum1 = 0 DO itype=1,atoms%ntype cdum1 = cdum1 + cexp1(itype)*sphbessel_integral(& atoms,itype,qnrm,nqnrm,& iqnrm1,iqnrm2,l,hybrid,& sphbes0,l_warn,l_warned)& / (2*l+1) l_warn = l_warn .AND. .NOT. l_warned ! only warn once END DO DO M=-l,l lm = lm + 1 cdum = cdum + cdum1 * y1(lm) * y2(lm) ENDDO ENDDO idum = ix*(ix-1)/2+iy coulomb(idum,ikpt) = coulomb(idum,ikpt)+(4*pi_const)**3*cdum / cell%vol `````` Daniel Wortmann committed Jun 21, 2017 1029 `````` END DO `````` Daniel Wortmann committed Jul 14, 2017 1030 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 1031 `````` `````` Daniel Wortmann committed Jul 14, 2017 1032 `````` END DO `````` Daniel Wortmann committed Jun 21, 2017 1033 `````` `````` Daniel Wortmann committed Jul 14, 2017 1034 `````` DEALLOCATE( carr2 ) `````` Daniel Wortmann committed Jun 21, 2017 1035 `````` `````` Daniel Wortmann committed Jul 14, 2017 1036 1037 1038 1039 1040 `````` IF ( mpi%irank == 0 ) THEN WRITE(6,'(2X,A)',advance='no') 'done' CALL cpu_TIME(time2) WRITE(6,'(2X,A,F8.2,A)') '( Timing:',time2-time1,' )' END IF `````` Daniel Wortmann committed Jun 21, 2017 1041 `````` `````` Daniel Wortmann committed Jul 14, 2017 1042 1043 1044 `````` ! ! Symmetry-equivalent G vectors ! `````` Daniel Wortmann committed Jun 21, 2017 1045 1046 ``````# ifndef CPP_NOCOULSYM `````` Daniel Wortmann committed Jul 14, 2017 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 `````` IF ( mpi%irank == 0 ) WRITE(6,'(A)',advance='no') 'Symm.-equiv. matrix elements...' CALL cpu_TIME(time1) ! All elements are needed so send all data to all processes treating the ! respective k-points ALLOCATE ( carr2(hybrid%maxbasm1,2),iarr(hybrid%maxgptm) ) ALLOCATE ( nsym_gpt(hybrid%gptmd,kpts%nkpt),& sym_gpt(MAXVAL(nsym1),hybrid%gptmd,kpts%nkpt) ) nsym_gpt = 0 ; sym_gpt = 0 DO ikpt = ikptmin,ikptmax carr2 = 0 ; iarr = 0 iarr(hybrid%pgptm1(:hybrid%ngptm1(ikpt),ikpt)) = 1 DO igpt0 = 1,hybrid%ngptm1(ikpt) !igptmin(ikpt),igptmax(ikpt) lsym = ( ( igptmin(ikpt) <= igpt0 ) .AND.& ( igptmax(ikpt) >= igpt0 ) ) igpt2 = hybrid%pgptm1(igpt0,ikpt) j = (hybrid``````