1c $Id$ 2 3* ********************************* 4* * * 5* * nwpw_xc_init * 6* * * 7* ********************************* 8 subroutine nwpw_xc_init(nion0,nkatm0, 9 > nprj,nbasis,n1dgrid,psp_type,lmax0, 10 > nprj_max0,l_prj,m_prj,b_prj) 11 implicit none 12 integer nion0,nkatm0 13 integer nprj(*),nbasis(*),n1dgrid(*),psp_type(*),lmax0(*) 14 integer nprj_max0 15 integer l_prj(nprj_max0,*),m_prj(nprj_max0,*),b_prj(nprj_max0,*) 16 17 18#include "bafdecls.fh" 19#include "errquit.fh" 20#include "nwpw_xc.fh" 21 22c **** local variables **** 23 logical ok 24 integer ii,ia,nsize 25 integer ic 26 integer l,m,li,mi,lj,mj,bi,bj,iprj,jprj,lm 27 integer i_p,i_t 28 integer i_tlm,i_plm,i_rlm 29 integer mtr_size 30 integer nb 31 real*8 tmp_theta,cs_theta,tmp_phi,tmp_gaunt,tmp_gaunt2,tmp_gaunt3 32 real*8 angle_phi 33 integer paw_vxc_size 34 35c **** external functions **** 36 integer control_gga,control_lmax_multipole 37 external control_gga,control_lmax_multipole 38 real*8 rtheta_lm,drtheta_lm,nwpw_gaunt,nwpw_gaunt2,nwpw_gaunt3 39 real*8 rtheta_lm_div 40 external rtheta_lm,drtheta_lm,nwpw_gaunt,nwpw_gaunt2,nwpw_gaunt3 41 external rtheta_lm_div 42 43 44 call nwpw_timing_start(4) 45 ok = .true. 46 47 nion = nion0 48 nkatm = nkatm0 49 gga = control_gga() 50 51 nprj_max = 0 52 nbasis_max = 0 53 n1dgrid_max = 0 54 do ia=1,nkatm 55 if (nprj(ia) .gt.nprj_max) nprj_max = nprj(ia) 56 if (nbasis(ia) .gt.nbasis_max) nbasis_max = nbasis(ia) 57 if (n1dgrid(ia).gt.n1dgrid_max) n1dgrid_max = n1dgrid(ia) 58 end do 59 60 mult_l_max = control_lmax_multipole() 61 if (mult_l_max.lt.0) then 62 mult_l_max = 0 63 do ia=1,nkatm 64 if (psp_type(ia).eq.4) then 65 if (mult_l_max.lt.(2*lmax0(ia))) mult_l_max = 2*lmax0(ia) 66 end if 67 end do 68 end if 69 !call nwpw_gaunt2_init(.false.,mult_l_max) 70 !call nwpw_gaunt3_init(.false.,mult_l_max) 71 72* ***paw_xc energies *** 73 ok = BA_alloc_get(mt_dbl,nion,"paw_xc_e",paw_xc_e(2),paw_xc_e(1)) 74 call dcopy(nion,0.0d0,0,dbl_mb(paw_xc_e(1)),1) 75 76 77* *** xc matrix arrays - nbasis_max**2 * (mult_l_max+1)**2 *** 78 mtr_size = 2*(nbasis_max**2)*((mult_l_max+1)**2) 79 ok = ok.and. 80 > BA_alloc_get(mt_dbl,mtr_size,"paw_xc_matr", 81 > paw_xc_matr(2),paw_xc_matr(1)) 82 83 if (gga.ge.10) then 84 ok = ok.and. 85 > BA_alloc_get(mt_dbl,3*mtr_size,"paw_xc_dmatr_u", 86 > paw_xc_dmatr(2),paw_xc_dmatr(1)) 87 end if 88 if (.not.ok) 89 > call errquit("init_paw_vxc: error allocating paw_xc_matr", 90 > mtr_size,0) 91 92 93 mtr_size = 2*nprj_max*nprj_max 94 ok = BA_alloc_get(mt_dbl,mtr_size,"paw_xc_pot", 95 > paw_xc_pot(2),paw_xc_pot(1)) 96 if (.not.ok) 97 > call errquit("init_paw_vxc: error allocating paw_xc_pot", 98 > mtr_size,0) 99 100 101 102* *** spherical grid arrays *** 103 if(mult_l_max .eq. 0 ) then 104 paw_xc_nphi = 1 105 paw_xc_ntheta = 1 106 else 107 paw_xc_nphi = 3*mult_l_max 108 paw_xc_ntheta = 3*mult_l_max 109 end if 110 111 112 ok = BA_alloc_get(mt_dbl,paw_xc_nphi,"paw_xc_angle_phi", 113 > paw_xc_angle_phi(2),paw_xc_angle_phi(1)) 114 115 ok = ok.and. 116 > BA_alloc_get(mt_dbl,paw_xc_ntheta,"paw_xc_cos_theta", 117 > paw_xc_cos_theta(2),paw_xc_cos_theta(1)) 118 119 ok = ok.and. 120 > BA_alloc_get(mt_dbl,paw_xc_nphi,"paw_xc_w_phi", 121 > paw_xc_w_phi(2),paw_xc_w_phi(1)) 122 123 ok = ok.and. 124 > BA_alloc_get(mt_dbl,paw_xc_ntheta,"paw_xc_w_theta", 125 > paw_xc_w_theta(2),paw_xc_w_theta(1)) 126 127 ok = ok.and. 128 > BA_alloc_get(mt_dbl, 129 > paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2, 130 > "paw_xc_tlm", 131 > paw_xc_tlm(2),paw_xc_tlm(1)) 132 133 134* **** used for generating derivatives of tlm's **** 135 if (gga.ge.10) then 136 137c **** derivatives wrt to theta **** 138 ok = ok.and. 139 > BA_alloc_get(mt_dbl, 140 > paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2, 141 > "paw_xc_dtlm_theta", 142 > paw_xc_dtlm_theta(2),paw_xc_dtlm_theta(1)) 143 144c **** derivatives wrt to phi **** 145 ok = ok.and. 146 > BA_alloc_get(mt_dbl, 147 > paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2, 148 > "paw_xc_dtlm_phi", 149 > paw_xc_dtlm_phi(2),paw_xc_dtlm_phi(1)) 150 151 end if 152 153 call nwpw_get_spher_grid(paw_xc_ntheta,paw_xc_nphi, 154 > dbl_mb(paw_xc_angle_phi(1)), 155 > dbl_mb(paw_xc_cos_theta(1)), 156 > dbl_mb(paw_xc_w_theta(1)), 157 > dbl_mb(paw_xc_w_phi(1))) 158 159c **** define tlm's **** 160 i_tlm = 0 161 do i_t=1,paw_xc_ntheta 162 do i_p=1,paw_xc_nphi 163 do l=0,mult_l_max 164 do m=-l,l 165 tmp_theta = rtheta_lm(l,m, 166 > dbl_mb(paw_xc_cos_theta(1)+i_t-1)) 167 angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1) 168 if (m.lt.0) then 169 tmp_phi = dsin(abs(m)*angle_phi) 170 else if (m.gt.0) then 171 tmp_phi = dcos(abs(m)*angle_phi) 172 else 173 tmp_phi = 1.0d0 174 end if 175 dbl_mb(paw_xc_tlm(1)+i_tlm) = tmp_theta*tmp_phi 176 i_tlm = i_tlm + 1 177 end do 178 end do 179 end do 180 end do 181 182 if (gga.ge.10) then 183 184c **** define derivative wrt to theta **** 185 i_tlm = 0 186 do i_t=1,paw_xc_ntheta 187 do i_p=1,paw_xc_nphi 188 do l=0,mult_l_max 189 do m=-l,l 190 tmp_theta = drtheta_lm(l,m, 191 > dbl_mb(paw_xc_cos_theta(1)+i_t-1)) 192 angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1) 193 if (m.lt.0) then 194 tmp_phi = dsin(abs(m)*angle_phi) 195 else if (m.gt.0) then 196 tmp_phi = dcos(abs(m)*angle_phi) 197 else 198 tmp_phi = 1.0d0 199 end if 200 dbl_mb(paw_xc_dtlm_theta(1)+i_tlm) = tmp_theta*tmp_phi 201 i_tlm = i_tlm + 1 202 end do 203 end do 204 end do 205 end do 206 207c **** define derivative wrt to phi **** 208 i_tlm = 0 209 do i_t=1,paw_xc_ntheta 210 do i_p=1,paw_xc_nphi 211 do l=0,mult_l_max 212 do m=-l,l 213 if (m.eq.0) then 214 dbl_mb(paw_xc_dtlm_phi(1)+ i_tlm) = 0.0d0 215 else 216 tmp_theta = rtheta_lm_div(l,m, 217 > dbl_mb(paw_xc_cos_theta(1)+i_t-1)) 218 219 angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1) 220 if (m.lt.0) then 221 tmp_phi = abs(m)*dcos(abs(m)*angle_phi) 222 else if (m.gt.0) then 223 tmp_phi = -abs(m)*dsin(abs(m)*angle_phi) 224 else 225 tmp_phi = 0.0d0 226 end if 227 dbl_mb(paw_xc_dtlm_phi(1)+ i_tlm) = tmp_theta*tmp_phi 228 end if 229 i_tlm = i_tlm + 1 230 end do 231 end do 232 end do 233 end do 234 235 end if 236 237* *** temp arrays *** 238 nsize = 2*n1dgrid_max 239 ok = ok.and. 240 > BA_alloc_get(mt_dbl,nsize,"rho_ae",rho_ae(2),rho_ae(1)) 241 ok = ok.and. 242 > BA_alloc_get(mt_dbl,nsize,"rho_ps",rho_ps(2),rho_ps(1)) 243 244* **** allocate gradient's and agr's **** 245 if (gga.ge.10) then 246 nsize = 2*3*n1dgrid_max 247 ok = ok.and. 248 > BA_alloc_get(mt_dbl,nsize, 249 > "paw_tmp_drho_ae", 250 > rho_ae_prime(2),rho_ae_prime(1)) 251 ok = ok.and. 252 > BA_alloc_get(mt_dbl,nsize, 253 > "paw_tmp_drho_ps", 254 > rho_ps_prime(2),rho_ps_prime(1)) 255 nsize = (2*2-1)*n1dgrid_max 256 ok = ok.and. 257 > BA_alloc_get(mt_dbl,nsize,"agr_ae",agr_ae(2),agr_ae(1)) 258 ok = ok.and. 259 > BA_alloc_get(mt_dbl,nsize,"agr_ps",agr_ps(2),agr_ps(1)) 260 ok = ok.and. 261 > BA_alloc_get(mt_dbl,nsize,"fdn_ae",fdn_ae(2),fdn_ae(1)) 262 ok = ok.and. 263 > BA_alloc_get(mt_dbl,nsize,"fdn_ps",fdn_ps(2),fdn_ps(1)) 264 end if 265 266 nsize = 2*n1dgrid_max 267 ok = ok.and. 268 > BA_alloc_get(mt_dbl,nsize,"vxc_ae",vxc_ae(2),vxc_ae(1)) 269 ok = ok.and. 270 > BA_alloc_get(mt_dbl,nsize,"vxc_ps",vxc_ps(2),vxc_ps(1)) 271 ok = ok.and. 272 > BA_alloc_get(mt_dbl,nsize,"exc_ae",exc_ae(2),exc_ae(1)) 273 ok = ok.and. 274 > BA_alloc_get(mt_dbl,nsize,"exc_ps",exc_ps(2),exc_ps(1)) 275 ok = ok.and. 276 > BA_alloc_get(mt_dbl,nsize,"xc_temp", 277 > xc_temp(2),xc_temp(1)) 278 ok = ok.and. 279 > BA_alloc_get(mt_dbl,nsize,"xc_temp2", 280 > xc_temp2(2),xc_temp2(1)) 281 282 283* *** allocate vxclm multipole expansion arrays **** 284 nsize = 2*n1dgrid_max*(mult_l_max+1)**2 285 286 ok = ok.and. 287 > BA_alloc_get(mt_dbl,nsize,"paw_vxc_ae", 288 > paw_vxc_ae(2),paw_vxc_ae(1)) 289 ok = ok.and. 290 > BA_alloc_get(mt_dbl,nsize,"paw_vxc_ps", 291 > paw_vxc_ps(2),paw_vxc_ps(1)) 292 293* *** allocate dvxclm multipole expansion arrays **** 294 if (gga.ge.10) then 295 ok = ok.and. 296 > BA_alloc_get(mt_dbl,3*nsize, 297 > "paw_dvxc_u_ae", 298 > paw_dvxc_ae(2),paw_dvxc_ae(1)) 299 ok = ok.and. 300 > BA_alloc_get(mt_dbl,3*nsize, 301 > "paw_dvxc_u_ps", 302 > paw_dvxc_ps(2),paw_dvxc_ps(1)) 303 end if 304 305* *** allocate rholm multipole expansion arrays **** 306 ok = ok.and. 307 > BA_alloc_get(mt_dbl,nsize,"paw_rho2_ae", 308 > paw_rho2_ae(2),paw_rho2_ae(1)) 309 ok = ok.and. 310 > BA_alloc_get(mt_dbl,nsize,"paw_rho2_ps", 311 > paw_rho2_ps(2),paw_rho2_ps(1)) 312 313* *** allocate rholm_prime multipole expansion arrays - need for GGA's**** 314 if (gga.ge.10) then 315 ok = ok.and. 316 > BA_alloc_get(mt_dbl,nsize, 317 > "paw_rho2_ae_prime", 318 > paw_rho2_ae_prime(2),paw_rho2_ae_prime(1)) 319 ok = ok.and. 320 > BA_alloc_get(mt_dbl,nsize, 321 > "paw_rho2_ps_prime", 322 > paw_rho2_ps_prime(2),paw_rho2_ps_prime(1)) 323 end if 324 if (.not.ok) 325 > call errquit("nwpw_xc_init: error allocating work arrays", 326 > 1,MA_ERR) 327 328c call nwpw_init_gntxc(mult_l_max) 329c if ((gga.ge.10).and.(mult_l_max.ge.1)) then 330c call nwpw_init_gntxc2(mult_l_max) 331c call nwpw_init_gntxc3(mult_l_max) 332c end if 333 334 335* ***** set up index arrays for nwpw_xc_density_solve2 ***** 336 ok = BA_alloc_get(mt_int,nkatm,"nindx_rholm", 337 > nindx_rholm(2),nindx_rholm(1)) 338 ok = ok.and. 339 > BA_alloc_get(mt_int,nkatm,"shift_rholm", 340 > shift_rholm(2),shift_rholm(1)) 341 if (.not.ok) 342 > call errquit("nwpw_xc_init: error allocating work arrays",2,0) 343 344 nsize = 0 345 do ia=1,nkatm 346 int_mb(shift_rholm(1)+ia-1) = nsize 347 do jprj = 1,nprj(ia) 348 lj = l_prj(jprj,ia) 349 mj = m_prj(jprj,ia) 350 do iprj = 1,nprj(ia) 351 li = l_prj(iprj,ia) 352 mi = m_prj(iprj,ia) 353 do l=0,mult_l_max 354 do m=-l,l 355 if ((l.le.(li+lj)) .and. (l.ge.abs(li-lj))) then 356 tmp_gaunt = nwpw_gaunt(.false.,l,m,li,mi,lj,mj) 357 if (dabs(tmp_gaunt).gt.1.0d-11) then 358 nsize = nsize + 1 359 end if 360 end if 361 end do 362 end do 363 364 end do 365 end do 366 int_mb(nindx_rholm(1)+ia-1) = nsize-int_mb(shift_rholm(1)+ia-1) 367 end do 368 369 ok = BA_alloc_get(mt_int,nsize,"bi_rholm",bi_rholm(2),bi_rholm(1)) 370 ok = ok.and. 371 > BA_alloc_get(mt_int,nsize,"bj_rholm",bj_rholm(2),bj_rholm(1)) 372 ok = ok.and. 373 > BA_alloc_get(mt_int,nsize,"lm_rholm",lm_rholm(2),lm_rholm(1)) 374 ok = ok.and. 375 > BA_alloc_get(mt_int,nsize,"iprj_rholm", 376 > iprj_rholm(2),iprj_rholm(1)) 377 ok = ok.and. 378 > BA_alloc_get(mt_int,nsize,"jprj_rholm", 379 > jprj_rholm(2),jprj_rholm(1)) 380 ok = ok.and. 381 > BA_alloc_get(mt_dbl,nsize,"coeff_rholm", 382 > coeff_rholm(2),coeff_rholm(1)) 383 ok = ok.and. 384 > BA_alloc_get(mt_dbl,nsize,"coeff_rholm2", 385 > coeff_rholm2(2),coeff_rholm2(1)) 386 ok = ok.and. 387 > BA_alloc_get(mt_dbl,nsize,"coeff_rholm3", 388 > coeff_rholm3(2),coeff_rholm3(1)) 389 if (.not.ok) 390 > call errquit("nwpw_xc_init: error allocating work arrays",3,0) 391 392 nsize = 0 393 do ia=1,nkatm 394 do jprj = 1,nprj(ia) 395 lj = l_prj(jprj,ia) 396 mj = m_prj(jprj,ia) 397 bj = b_prj(jprj,ia) 398 do iprj = 1,nprj(ia) 399 li = l_prj(iprj,ia) 400 mi = m_prj(iprj,ia) 401 bi = b_prj(iprj,ia) 402 lm = 0 403 do l=0,mult_l_max 404 do m=-l,l 405 if ((l.le.(li+lj)) .and. (l.ge.abs(li-lj))) then 406 tmp_gaunt = nwpw_gaunt(.false.,l,m,li,mi,lj,mj) 407 if ((gga.ge.10).and.(mult_l_max.gt.0)) then 408 tmp_gaunt2=nwpw_gaunt2(.false.,l,m,li,mi,lj,mj) 409 tmp_gaunt3=nwpw_gaunt3(.false.,l,m,li,mi,lj,mj) 410 else 411 tmp_gaunt2 = 0.0d0 412 tmp_gaunt3 = 0.0d0 413 end if 414 if ((dabs(tmp_gaunt) .gt.1.0d-11).or. 415 > (dabs(tmp_gaunt2).gt.1.0d-11).or. 416 > (dabs(tmp_gaunt3).gt.1.0d-11)) then 417 dbl_mb(coeff_rholm(1) +nsize) = tmp_gaunt 418 dbl_mb(coeff_rholm2(1)+nsize) = tmp_gaunt2 419 dbl_mb(coeff_rholm3(1)+nsize) = tmp_gaunt3 420 int_mb(lm_rholm(1)+nsize) = lm 421 int_mb(iprj_rholm(1)+nsize) = iprj 422 int_mb(jprj_rholm(1)+nsize) = jprj 423 int_mb(bi_rholm(1)+nsize) = bi 424 int_mb(bj_rholm(1)+nsize) = bj 425 nsize = nsize + 1 426 end if 427 end if 428 lm = lm + 1 429 end do 430 end do 431 end do 432 end do 433 end do 434 435 436 call nwpw_timing_end(4) 437 return 438 end 439 440* ******************************************** 441* * * 442* * nwpw_xc_solve * 443* * * 444* ******************************************** 445 subroutine nwpw_xc_solve(ii,ia,n1dgrid,nbasis, 446 > phi_ae,phi_ps,phi_ae_prime,phi_ps_prime, 447 > core_ae,core_ps,core_ae_prime,core_ps_prime, 448 > rgrid,log_amesh, 449 > ispin,ne,nprj,sw1,sw2) 450 implicit none 451 integer ii,ia 452 integer n1dgrid,nbasis 453 real*8 phi_ae(n1dgrid,nbasis) 454 real*8 phi_ps(n1dgrid,nbasis) 455 real*8 phi_ae_prime(n1dgrid,nbasis) 456 real*8 phi_ps_prime(n1dgrid,nbasis) 457 real*8 core_ae(n1dgrid) 458 real*8 core_ps(n1dgrid) 459 real*8 core_ae_prime(n1dgrid) 460 real*8 core_ps_prime(n1dgrid) 461 real*8 rgrid(n1dgrid) 462 real*8 log_amesh 463 464 integer ispin,ne(2),nprj 465 real*8 sw1(ne(1)+ne(2),nprj) 466 real*8 sw2(ne(1)+ne(2),nprj) 467 468#include "bafdecls.fh" 469#include "errquit.fh" 470#include "nwpw_xc.fh" 471 472 logical value 473 integer i,j,np,i_tlm,i_tlm1,lmax2,nsize,shift,ig 474 integer i_t,i_p,lm 475 real*8 exc_tmp 476 477* **** external functions **** 478 real*8 log_integrate_def_destroyf 479 external log_integrate_def_destroyf 480 481 call nwpw_timing_start(36) 482 483* *** index for spin down in temp arrays *** 484 lmax2 = (mult_l_max+1)**2 485!$OMP MASTER 486 dbl_mb(paw_xc_e(1)+ii-1)=0.0d0 487!$OMP END MASTER 488 489* *** zero out vxclm multipole arrays *** 490 nsize = ispin*n1dgrid_max*lmax2 491 !call dcopy(nsize,0.0d0,0,dbl_mb(paw_vxc_ae(1)),1) 492 !call dcopy(nsize,0.0d0,0,dbl_mb(paw_vxc_ps(1)),1) 493 call Parallel_shared_vector_zero(.false.,nsize, 494 > dbl_mb(paw_vxc_ae(1))) 495 call Parallel_shared_vector_zero(.true.,nsize, 496 > dbl_mb(paw_vxc_ps(1))) 497 498* *** zero out dvxclm multipole arrays *** 499 if (gga.ge.10) then 500 !call dcopy(3*nsize,0.0d0,0,dbl_mb(paw_dvxc_ae(1)),1) 501 !call dcopy(3*nsize,0.0d0,0,dbl_mb(paw_dvxc_ps(1)),1) 502 call Parallel_shared_vector_zero(.false.,3*nsize, 503 > dbl_mb(paw_dvxc_ae(1))) 504 call Parallel_shared_vector_zero(.true.,3*nsize, 505 > dbl_mb(paw_dvxc_ps(1))) 506 end if 507 508* *** find rholm multipole expansion on spherical grid *** 509 shift = int_mb(shift_rholm(1)+ia-1) 510 call nwpw_xc_density_solve2(ispin,ne,nprj,sw1, 511 > n1dgrid,nbasis,lmax2, 512 > int_mb(nindx_rholm(1)+ia-1), 513 > int_mb(lm_rholm(1)+shift), 514 > int_mb(iprj_rholm(1)+shift), 515 > int_mb(jprj_rholm(1)+shift), 516 > int_mb(bi_rholm(1)+shift), 517 > int_mb(bj_rholm(1)+shift), 518 > dbl_mb(coeff_rholm(1)+shift), 519 > phi_ae,phi_ps,rgrid, 520 > dbl_mb(paw_rho2_ae(1)), 521 > dbl_mb(paw_rho2_ps(1))) 522 523* *** find rholm_prime multipole expansion on spherical grid *** 524 if (gga.ge.10) then 525 shift = int_mb(shift_rholm(1)+ia-1) 526 call nwpw_xc_density_prime_solve2(ispin,ne,nprj,sw1, 527 > n1dgrid,nbasis,lmax2, 528 > int_mb(nindx_rholm(1)+ia-1), 529 > int_mb(lm_rholm(1)+shift), 530 > int_mb(iprj_rholm(1)+shift), 531 > int_mb(jprj_rholm(1)+shift), 532 > int_mb(bi_rholm(1)+shift), 533 > int_mb(bj_rholm(1)+shift), 534 > dbl_mb(coeff_rholm(1)+shift), 535 > phi_ae,phi_ps, 536 > phi_ae_prime,phi_ps_prime,rgrid, 537 > dbl_mb(paw_rho2_ae_prime(1)), 538 > dbl_mb(paw_rho2_ps_prime(1))) 539 end if 540 541 542 543 i_tlm = 0 544 i_tlm1 = 0 545 do i_t = 1,paw_xc_ntheta 546 do i_p = 1,paw_xc_nphi 547 548* *** zero out temp arrays *** 549 nsize = ispin*n1dgrid_max 550 !call dcopy(nsize,0.0d0,0,dbl_mb(rho_ae(1)),1) 551 !call dcopy(nsize,0.0d0,0,dbl_mb(rho_ps(1)),1) 552 !call dcopy(nsize,0.0d0,0,dbl_mb(vxc_ae(1)),1) 553 !call dcopy(nsize,0.0d0,0,dbl_mb(vxc_ps(1)),1) 554 !call dcopy(nsize,0.0d0,0,dbl_mb(exc_ae(1)),1) 555 !call dcopy(nsize,0.0d0,0,dbl_mb(exc_ps(1)),1) 556 call Parallel_shared_vector_zero(.false.,nsize, 557 > dbl_mb(rho_ae(1))) 558 call Parallel_shared_vector_zero(.false.,nsize, 559 > dbl_mb(rho_ps(1))) 560 call Parallel_shared_vector_zero(.false.,nsize, 561 > dbl_mb(vxc_ae(1))) 562 call Parallel_shared_vector_zero(.false.,nsize, 563 > dbl_mb(vxc_ps(1))) 564 call Parallel_shared_vector_zero(.false.,nsize, 565 > dbl_mb(exc_ae(1))) 566 call Parallel_shared_vector_zero(.true., nsize, 567 > dbl_mb(exc_ps(1))) 568 569* *** generate atomic densities on spherical grid *** 570 call nwpw_xc_gen_atomic_densities(n1dgrid,lmax2,ispin, 571 > dbl_mb(paw_rho2_ae(1)), 572 > dbl_mb(paw_xc_tlm(1)+i_tlm), 573 > core_ae, 574 > dbl_mb(rho_ae(1))) 575 call nwpw_xc_gen_atomic_densities(n1dgrid,lmax2,ispin, 576 > dbl_mb(paw_rho2_ps(1)), 577 > dbl_mb(paw_xc_tlm(1)+i_tlm), 578 > core_ps, 579 > dbl_mb(rho_ps(1))) 580 581* *************************************************** 582* *** find exchange-correlation on spherical grid *** 583* *************************************************** 584 585* **** LDA functionals on spherical grid **** 586 if (gga.lt.10) then 587 588* **** LDA functionals of ae and ps atomic densities **** 589 call nwpw_vosko(n1dgrid,n1dgrid,ispin, 590 > dbl_mb(rho_ae(1)), 591 > dbl_mb(exc_ae(1)), 592 > dbl_mb(vxc_ae(1)), 593 > dbl_mb(xc_temp(1))) 594 call nwpw_vosko(n1dgrid,n1dgrid,ispin, 595 > dbl_mb(rho_ps(1)), 596 > dbl_mb(exc_ps(1)), 597 > dbl_mb(vxc_ps(1)), 598 > dbl_mb(xc_temp(1))) 599 600 601* **** GGA functionals on spherical grid **** 602 else 603 604* *** generate the gradient of atomic densities in spherical **** 605* **** coordinates on spherical grid **** 606 call nwpw_xc_gen_atomic_gradients(n1dgrid,lmax2,ispin, 607 > dbl_mb(paw_rho2_ae(1)), 608 > dbl_mb(paw_rho2_ae_prime(1)), 609 > dbl_mb(paw_xc_tlm(1)+i_tlm), 610 > dbl_mb(paw_xc_dtlm_theta(1)+i_tlm), 611 > dbl_mb(paw_xc_dtlm_phi(1)+i_tlm), 612 > rgrid,core_ae,core_ae_prime, 613 > dbl_mb(rho_ae_prime(1))) 614 615 call nwpw_xc_gen_atomic_gradients(n1dgrid,lmax2,ispin, 616 > dbl_mb(paw_rho2_ps(1)), 617 > dbl_mb(paw_rho2_ps_prime(1)), 618 > dbl_mb(paw_xc_tlm(1)+i_tlm), 619 > dbl_mb(paw_xc_dtlm_theta(1)+i_tlm), 620 > dbl_mb(paw_xc_dtlm_phi(1)+i_tlm), 621 > rgrid,core_ps,core_ps_prime, 622 > dbl_mb(rho_ps_prime(1))) 623 624* **** generate agr **** 625 call nwpw_xc_gen_atomic_agr(n1dgrid,lmax2,ispin, 626 > dbl_mb(rho_ae_prime(1)), 627 > dbl_mb(agr_ae(1))) 628 call nwpw_xc_gen_atomic_agr(n1dgrid,lmax2,ispin, 629 > dbl_mb(rho_ps_prime(1)), 630 > dbl_mb(agr_ps(1))) 631 632* **** GGA functionals of ae and ps atomic densities **** 633 call nwpw_gga(gga,n1dgrid,ispin, 634 > dbl_mb(rho_ae(1)), 635 > dbl_mb(agr_ae(1)), 636 > dbl_mb(exc_ae(1)), 637 > dbl_mb(vxc_ae(1)), !* df/dnup, df/dndn 638 > dbl_mb(fdn_ae(1)), !* df/d|grad nup|,df/d|grad ndn|, df/d|grad n| 639 > dbl_mb(xc_temp(1))) 640 call nwpw_gga(gga,n1dgrid,ispin, 641 > dbl_mb(rho_ps(1)), 642 > dbl_mb(agr_ps(1)), 643 > dbl_mb(exc_ps(1)), 644 > dbl_mb(vxc_ps(1)), !* df/dnup, df/dndn 645 > dbl_mb(fdn_ps(1)), !* df/d|grad nup|,df/d|grad ndn|, df/d|grad n| 646 > dbl_mb(xc_temp(1))) 647 648 649c **** if restricted calculate df/d|grad n| *(grad n)/|grad n| **** 650* ***** and put it in rho_prime **** 651 652* **** if unrestricted calculate (df/d|grad nup|* (grad nup)/|grad nup|) **** 653* **** + (df/d|grad n| * (grad n)/|grad n|) **** 654* **** and calculate (df/d|grad ndn|* (grad ndn)/|grad ndn|) **** 655* **** + (df/d|grad n| * (grad n)/|grad n|) **** 656* ***** and put it in rho_prime **** 657 call nwpw_xc_gen_dvxc(n1dgrid,lmax2,ispin, 658 > dbl_mb(fdn_ae(1)), 659 > dbl_mb(agr_ae(1)), 660 > dbl_mb(rho_ae_prime(1))) 661 call nwpw_xc_gen_dvxc(n1dgrid,lmax2,ispin, 662 > dbl_mb(fdn_ps(1)), 663 > dbl_mb(agr_ps(1)), 664 > dbl_mb(rho_ps_prime(1))) 665 666c **** find dvxclm multipole expansion **** 667 call nwpw_xc_addto_dvxclm(n1dgrid,lmax2,ispin, 668 > (dbl_mb(paw_xc_w_theta(1)+i_t-1) 669 > *dbl_mb(paw_xc_w_phi(1) +i_p-1)), 670 > dbl_mb(paw_xc_tlm(1)+i_tlm), 671 > dbl_mb(rho_ae_prime(1)), 672 > dbl_mb(rho_ps_prime(1)), 673 > dbl_mb(paw_dvxc_ae(1)), 674 > dbl_mb(paw_dvxc_ps(1))) 675 676 end if 677 i_tlm = i_tlm + lmax2 678 679 680* ************************************************************************************ 681* *** find vxclm the multipole expansion of atomic vxc=df/dn or (df/dnup,df/dndn) *** 682* ************************************************************************************ 683 call nwpw_xc_addto_vxclm(n1dgrid,lmax2,ispin, 684 > (dbl_mb(paw_xc_w_theta(1)+i_t-1) 685 > *dbl_mb(paw_xc_w_phi(1) +i_p-1)), 686 > dbl_mb(paw_xc_tlm(1)+i_tlm1), 687 > dbl_mb(vxc_ae(1)), 688 > dbl_mb(vxc_ps(1)), 689 > dbl_mb(paw_vxc_ae(1)), 690 > dbl_mb(paw_vxc_ps(1))) 691 i_tlm1 = i_tlm1 + lmax2 692 693* ****************************************************** 694* *** compute the atomic exchange correlation energy *** 695* ****************************************************** 696!$OMP DO 697 do ig=1,n1dgrid 698 dbl_mb(xc_temp(1)+ig-1)= 699 > ( 700 > dbl_mb(rho_ae(1)+ig-1) + 701 > dbl_mb(rho_ae(1)+(ispin-1)*n1dgrid+ig-1) 702 > )*dbl_mb(exc_ae(1)+ig-1) 703 > - 704 > ( 705 > dbl_mb(rho_ps(1)+ig-1)+ 706 > dbl_mb(rho_ps(1)+(ispin-1)*n1dgrid+ig-1) 707 > )* 708 > dbl_mb(exc_ps(1)+ig-1) 709 end do 710!$OMP END DO 711!$OMP MASTER 712 exc_tmp = log_integrate_def_destroyf(0,dbl_mb(xc_temp(1)), 713 > 2,rgrid,log_amesh,n1dgrid) 714 dbl_mb(paw_xc_e(1)+ii-1)=dbl_mb(paw_xc_e(1)+ii-1)+ 715 > exc_tmp* 716 > dbl_mb(paw_xc_w_theta(1)+i_t-1)* 717 > dbl_mb(paw_xc_w_phi(1)+i_p-1) 718!$OMP END MASTER 719!$OMP BARRIER 720 721 end do !i_p 722 end do !i_t 723 724 725* ********************************************************* 726* **** non-local operator computation - LDA part **** 727* ********************************************************* 728 729* **** compute (vxc^a)_nln'l'^lm radial integrals ***** 730 call nwpw_xc_gen_matr(n1dgrid,nbasis,lmax2,ispin, 731 > log_amesh, 732 > phi_ae,phi_ps,rgrid, 733 > dbl_mb(paw_vxc_ae(1)), 734 > dbl_mb(paw_vxc_ps(1)), 735 > dbl_mb(xc_temp(1)), 736 > dbl_mb(paw_xc_matr(1))) 737 738 739* **** xc potential non-local matrix elements **** 740 shift = int_mb(shift_rholm(1)+ia-1) 741 call nwpw_xc_sw1tosw2(ispin,ne,nprj,sw1,sw2, 742 > n1dgrid,nbasis,lmax2, 743 > int_mb(nindx_rholm(1)+ia-1), 744 > int_mb(lm_rholm(1)+shift), 745 > int_mb(iprj_rholm(1)+shift), 746 > int_mb(jprj_rholm(1)+shift), 747 > int_mb(bi_rholm(1)+shift), 748 > int_mb(bj_rholm(1)+shift), 749 > dbl_mb(coeff_rholm(1)+shift), 750 > dbl_mb(paw_xc_matr(1))) 751 752 753* ********************************************************* 754* **** non-local operator computation - GGA part **** 755* ********************************************************* 756 if (gga.ge.10) then 757 758* **** compute (dvxc^a)_nln'l'^lm radial integrals ***** 759 call nwpw_xc_gen_dmatr(n1dgrid,nbasis,lmax2,ispin, 760 > log_amesh, 761 > phi_ae,phi_ps, 762 > phi_ae_prime,phi_ps_prime, 763 > rgrid, 764 > dbl_mb(paw_dvxc_ae(1)), 765 > dbl_mb(paw_dvxc_ps(1)), 766 > dbl_mb(xc_temp(1)), 767 > dbl_mb(paw_xc_dmatr(1))) 768 769 shift = int_mb(shift_rholm(1)+ia-1) 770 call nwpw_xc_sw1tosw2_dmatr(ispin,ne,nprj,sw1,sw2, 771 > n1dgrid,nbasis,lmax2, 772 > int_mb(nindx_rholm(1)+ia-1), 773 > int_mb(lm_rholm(1)+shift), 774 > int_mb(iprj_rholm(1)+shift), 775 > int_mb(jprj_rholm(1)+shift), 776 > int_mb(bi_rholm(1)+shift), 777 > int_mb(bj_rholm(1)+shift), 778 > dbl_mb(coeff_rholm(1)+shift), 779 > dbl_mb(coeff_rholm2(1)+shift), 780 > dbl_mb(coeff_rholm3(1)+shift), 781 > dbl_mb(paw_xc_dmatr(1))) 782 end if !**gga** 783 784c if (np.gt.1) then 785c 786cc **** xc non-local matrix elements **** 787c call D3dB_Vector_Sumall(2*paw_xc_pot(3), 788c > dbl_mb(paw_xc_pot(1))) 789c 790cc **** atomic exchange-correlation energies **** 791c call D3dB_Vector_Sumall(paw_xc_e(3),dbl_mb(paw_xc_e(1))) 792c end if 793c 794c 795 call nwpw_timing_end(36) 796 797 return 798 end 799 800 801* ***************************************** 802* * * 803* * nwpw_xc_end * 804* * * 805* ***************************************** 806 subroutine nwpw_xc_end() 807 implicit none 808 809#include "bafdecls.fh" 810#include "errquit.fh" 811#include "nwpw_xc.fh" 812 813 logical ok 814 integer ms,i 815 816 call nwpw_timing_start(4) 817 818 !call nwpw_gaunt2_end() 819 !call nwpw_gaunt3_end() 820 821 ok = .true. 822 ok = ok.and.BA_free_heap(coeff_rholm(2)) 823 ok = ok.and.BA_free_heap(coeff_rholm2(2)) 824 ok = ok.and.BA_free_heap(coeff_rholm3(2)) 825 ok = ok.and.BA_free_heap(jprj_rholm(2)) 826 ok = ok.and.BA_free_heap(iprj_rholm(2)) 827 ok = ok.and.BA_free_heap(lm_rholm(2)) 828 ok = ok.and.BA_free_heap(bj_rholm(2)) 829 ok = ok.and.BA_free_heap(bi_rholm(2)) 830 ok = ok.and.BA_free_heap(shift_rholm(2)) 831 ok = ok.and.BA_free_heap(nindx_rholm(2)) 832 833 834 ok = ok.and.BA_free_heap(exc_ps(2)) 835 ok = ok.and.BA_free_heap(exc_ae(2)) 836 ok = ok.and.BA_free_heap(vxc_ps(2)) 837 ok = ok.and.BA_free_heap(vxc_ae(2)) 838 ok = ok.and.BA_free_heap(rho_ps(2)) 839 ok = ok.and.BA_free_heap(rho_ae(2)) 840 ok = ok.and.BA_free_heap(xc_temp(2)) 841 ok = ok.and.BA_free_heap(xc_temp2(2)) 842 843 ok = ok.and.BA_free_heap(paw_xc_e(2)) 844 ok = ok.and.BA_free_heap(paw_xc_tlm(2)) 845 ok = ok.and.BA_free_heap(paw_xc_w_theta(2) ) 846 ok = ok.and.BA_free_heap(paw_xc_w_phi(2)) 847 ok = ok.and.BA_free_heap(paw_xc_cos_theta(2)) 848 ok = ok.and.BA_free_heap(paw_xc_angle_phi(2)) 849 850 ok = ok.and.BA_free_heap(paw_xc_pot(2)) 851 ok = ok.and.BA_free_heap(paw_xc_matr(2)) 852 ok = ok.and.BA_free_heap(paw_vxc_ae(2)) 853 ok = ok.and.BA_free_heap(paw_vxc_ps(2)) 854 ok = ok.and.BA_free_heap(paw_rho2_ae(2)) 855 ok = ok.and.BA_free_heap(paw_rho2_ps(2)) 856 857 if (gga.ge.10) then 858 ok = ok.and.BA_free_heap(paw_rho2_ae_prime(2)) 859 ok = ok.and.BA_free_heap(paw_rho2_ps_prime(2)) 860 ok = ok.and.BA_free_heap(rho_ps_prime(2)) 861 ok = ok.and.BA_free_heap(rho_ae_prime(2)) 862 ok = ok.and.BA_free_heap(agr_ae(2)) 863 ok = ok.and.BA_free_heap(agr_ps(2)) 864 ok = ok.and.BA_free_heap(fdn_ae(2)) 865 ok = ok.and.BA_free_heap(fdn_ps(2)) 866 ok = ok.and.BA_free_heap(paw_xc_dtlm_phi(2)) 867 ok = ok.and.BA_free_heap(paw_xc_dtlm_theta(2)) 868 ok = ok.and.BA_free_heap(paw_dvxc_ae(2)) 869 ok = ok.and.BA_free_heap(paw_dvxc_ps(2)) 870 ok = ok.and.BA_free_heap(paw_xc_dmatr(2)) 871 end if 872 873 874 if (.not.ok) 875 > call errquit("nwpw_xc_end: error freeing heap",0,MA_ERR) 876 877 878 call nwpw_timing_end(4) 879 return 880 end 881 882* ****************************************** 883* * * 884* * nwpw_xc_mult_l_max * 885* * * 886* ****************************************** 887 integer function nwpw_xc_mult_l_max() 888 implicit none 889 890#include "bafdecls.fh" 891#include "errquit.fh" 892#include "nwpw_xc.fh" 893 894 nwpw_xc_mult_l_max = mult_l_max 895 return 896 end 897 898 899 900* ****************************************** 901* * * 902* * nwpw_xc_energy_atom * 903* * * 904* ****************************************** 905 real*8 function nwpw_xc_energy_atom() 906 implicit none 907 908#include "bafdecls.fh" 909#include "errquit.fh" 910#include "nwpw_xc.fh" 911 912 integer ii 913 real*8 e 914 915 e = 0.0d0 916 do ii=1,nion 917 e = e + dbl_mb(paw_xc_e(1)+ii-1) 918 end do 919 nwpw_xc_energy_atom = e 920 return 921 end 922 923 924 925* ****************************************** 926* * * 927* * nwpw_xc_gen_sphere_rho * 928* * * 929* ****************************************** 930 subroutine nwpw_xc_gen_sphere_rho(ng,lmax2,rholm,Tlm,rho) 931 implicit none 932 integer ng,lmax2 933 real*8 rholm(ng,lmax2) 934 real*8 Tlm(lmax2) 935 real*8 rho(ng) 936 937 integer lm,ig 938 939!$OMP DO 940 do ig=1,ng 941 do lm=1,lmax2 942 rho(ig) = rho(ig) + rholm(ig,lm)*Tlm(lm) 943 end do 944 end do 945!$OMP END DO 946 947 return 948 end 949 950 951* ****************************************** 952* * * 953* * nwpw_xc_rho_div_r * 954* * * 955* ****************************************** 956 subroutine nwpw_xc_rho_div_r(ic,r,rho) 957 implicit none 958 integer ic 959 real*8 r(*) 960 real*8 rho(*) 961 962 integer ig 963!$OMP DO 964 do ig=1,ic 965 rho(ig) = rho(ig)/r(ig) 966 end do 967!$OMP END DO 968 return 969 end 970 971 972 973* ******************************************************** 974* * * 975* * nwpw_xc_addto_vxclm * 976* * * 977* ******************************************************** 978 subroutine nwpw_xc_addto_vxclm(ic,lmax2,ispin, 979 > alpha, 980 > tlm, 981 > vxc_ae,vxc_ps, 982 > vxclm_ae,vxclm_ps) 983 implicit none 984 integer ic,lmax2,ispin 985 real*8 alpha 986 real*8 tlm(*) 987 real*8 vxc_ae(ic,ispin) 988 real*8 vxc_ps(ic,ispin) 989 real*8 vxclm_ae(ic,lmax2,ispin) 990 real*8 vxclm_ps(ic,lmax2,ispin) 991 992* **** local variables **** 993 integer i,lm,ms 994 995 do ms = 1,ispin 996 do lm = 1,lmax2 997!$OMP DO 998 do i = 1,ic 999 vxclm_ae(i,lm,ms) = vxclm_ae(i,lm,ms) 1000 > + vxc_ae(i,ms)*(tlm(lm)*alpha) 1001 vxclm_ps(i,lm,ms) = vxclm_ps(i,lm,ms) 1002 > + vxc_ps(i,ms)*(tlm(lm)*alpha) 1003 end do 1004!$OMP END DO 1005 end do 1006 end do 1007 return 1008 end 1009 1010* ******************************************************** 1011* * * 1012* * nwpw_xc_gen_matr * 1013* * * 1014* ******************************************************** 1015 subroutine nwpw_xc_gen_matr(ng,nb,lmax2,ispin,log_amesh, 1016 > phi_ae,phi_ps,r, 1017 > vxclm_ae,vxclm_ps, 1018 > tmpC, 1019 > matr) 1020 implicit none 1021 integer ng,nb,lmax2,ispin 1022 real*8 log_amesh 1023 real*8 phi_ae(ng,nb) 1024 real*8 phi_ps(ng,nb) 1025 real*8 r(ng) 1026 real*8 vxclm_ae(ng,lmax2,ispin) 1027 real*8 vxclm_ps(ng,lmax2,ispin) 1028 real*8 tmpC(ng) 1029 real*8 matr(nb,nb,lmax2,ispin) 1030 1031* **** local varialbles **** 1032 integer ig,i,j,lm,ms 1033 real*8 tmp_ae,tmp_ps 1034 1035* **** external functions **** 1036 real*8 log_integrate_def_destroyf 1037 external log_integrate_def_destroyf 1038 1039 do ms=1,ispin 1040 do lm=1,lmax2 1041 do i=1,nb 1042 do j=i,nb 1043!$OMP DO 1044 do ig=1,ng 1045 tmp_ae = phi_ae(ig,i) 1046 > *phi_ae(ig,j)/(r(ig)**2) 1047 tmp_ps = phi_ps(ig,i) 1048 > *phi_ps(ig,j)/(r(ig)**2) 1049 tmpC(ig) = vxclm_ae(ig,lm,ms)*tmp_ae 1050 > - vxclm_ps(ig,lm,ms)*tmp_ps 1051 end do 1052!$OMP END DO 1053!$OMP MASTER 1054 matr(i,j,lm,ms) 1055 > = log_integrate_def_destroyf(0,tmpC,2,r,log_amesh,ng) 1056 if (i.ne.j) matr(j,i,lm,ms) = matr(i,j,lm,ms) 1057!$OMP END MASTER 1058!$OMP BARRIER 1059 end do 1060 end do 1061 end do 1062 end do 1063 return 1064 end 1065 1066 1067* ******************************************************** 1068* * * 1069* * nwpw_xc_gen_atomic_densities * 1070* * * 1071* ******************************************************** 1072 subroutine nwpw_xc_gen_atomic_densities(ng,lmax2,ispin, 1073 > rholm, 1074 > tlm, 1075 > rhocore, 1076 > rho) 1077 implicit none 1078 integer ng,lmax2,ispin 1079 real*8 rholm(ng,lmax2,ispin) 1080 real*8 tlm(lmax2) 1081 real*8 rhocore(ng) 1082 real*8 rho(ng,ispin) 1083 1084 integer ms,i 1085 1086 do ms=1,ispin 1087 call nwpw_xc_gen_sphere_rho(ng,lmax2, 1088 > rholm(1,1,ms), 1089 > tlm, 1090 > rho(1,ms)) 1091 !call daxpy(ng,0.5d0,rhocore,1,rho(1,ms),1) 1092!$OMP DO 1093 do i=1,ng 1094 rho(i,ms) = dabs(rho(i,ms)+0.5d0*rhocore(i)) 1095 end do 1096!$OMP END DO 1097 end do 1098 1099 return 1100 end 1101 1102 1103 1104* ******************************************************** 1105* * * 1106* * nwpw_xc_gen_atomic_gradients * 1107* * * 1108* ******************************************************** 1109 subroutine nwpw_xc_gen_atomic_gradients(ic,lmax2,ispin, 1110 > rholm,rholm_prime, 1111 > tlm,dtlm_theta,dtlm_phi, 1112 > r,rhocore,rhocore_prime, 1113 > rho_prime) 1114 implicit none 1115 integer ic,lmax2,ispin 1116 real*8 rholm(ic,lmax2,ispin) 1117 real*8 rholm_prime(ic,lmax2,ispin) 1118 real*8 tlm(lmax2) 1119 real*8 dtlm_theta(lmax2) 1120 real*8 dtlm_phi(lmax2) 1121 real*8 r(ic) 1122 real*8 rhocore(ic) 1123 real*8 rhocore_prime(ic) 1124 1125 real*8 rho_prime(ic,3,ispin) 1126 1127 1128 integer ms,i 1129 1130 !call dcopy(3*ic*ispin,0.0d0,0,rho_prime,1) 1131 call Parallel_shared_vector_zero(.true.,3*ic*ispin,rho_prime) 1132 1133 do ms=1,ispin 1134 1135* *** find [d/dr paw_rho_ae] on spherical grid *** 1136 call nwpw_xc_gen_sphere_rho(ic,lmax2, 1137 > rholm_prime(1,1,ms), 1138 > tlm, 1139 > rho_prime(1,1,ms)) 1140 1141* *** add core d/dr densities*** 1142 !call daxpy(ic,0.5d0,rhocore_prime,1,rho_prime(1,1,ms),1) 1143!$OMP DO 1144 do i=1,ic 1145 rho_prime(i,1,ms) = rho_prime(i,1,ms) 1146 > + 0.5d0*rhocore_prime(i) 1147 end do 1148!$OMP END DO 1149 1150 1151 !*** only computate radial derivatives if lmax2==1 **** 1152 if (lmax2.gt.1) then 1153 1154* *** find (1/r) * d/dtheta paw_rho_ae on spherical grid *** 1155 call nwpw_xc_gen_sphere_rho(ic,lmax2, 1156 > rholm(1,1,ms), 1157 > dtlm_theta, 1158 > rho_prime(1,2,ms)) 1159 call nwpw_xc_rho_div_r(ic,r,rho_prime(1,2,ms)) 1160 1161* *** find [(r*sin(theta)) * d/dphi paw_rho_ae] on spherical grid *** 1162 call nwpw_xc_gen_sphere_rho(ic,lmax2, 1163 > rholm(1,1,ms), 1164 > dtlm_phi, 1165 > rho_prime(1,3,ms)) 1166 call nwpw_xc_rho_div_r(ic,r,rho_prime(1,3,ms)) 1167 1168 end if 1169 1170 1171 end do 1172 1173 return 1174 end 1175 1176 1177* ************************************ 1178* * * 1179* * nwpw_xc_gen_atomic_agr * 1180* * * 1181* ************************************ 1182* 1183* This function returns the absolute values of the gradient. 1184* 1185* Entry - ic : number of grid points 1186* ispin : restricted/unrestricted 1187* rho_prime(ic,3,ispin) : gradient in spherical coordinates 1188* of atomic spin densites nup and ndn 1189* 1190* Exit - agr_in(*,1): |grad n| if restricted 1191* Exit - agr_in(*,3): |grad nup|, |grad ndn|, and |grad n| if unrestricted 1192 1193 subroutine nwpw_xc_gen_atomic_agr(ic,lmax2,ispin, 1194 > rho_prime, 1195 > agr) 1196 implicit none 1197 integer ic,lmax2,ispin 1198 real*8 rho_prime(ic,3,ispin) 1199 real*8 agr(ic,*) !*(ic,2*ispin-1) 1200 1201 integer ig,ms 1202 1203 1204 !*** only computate radial derivatives if lmax2==1 **** 1205 if (lmax2.eq.1) then 1206 1207c **** compute |grad n| **** 1208!$OMP DO 1209 do ig=1,ic 1210 agr(ig,2*ispin-1) 1211 > = dsqrt( (rho_prime(ig,1,1)+rho_prime(ig,1,ispin))**2) 1212 end do 1213!$OMP END DO 1214 1215c **** compute |grad nup| and |grad ndn| **** 1216 if (ispin.eq.2) then 1217 do ms=1,ispin 1218!$OMP DO 1219 do ig=1,ic 1220 agr(ig,ms) = dsqrt(rho_prime(ig,1,ms)**2) 1221 end do 1222!$OMP END DO 1223 end do 1224 end if 1225 1226 1227 1228 !*** computate theta and phi derivatives if lmax2>1 **** 1229 else 1230c **** compute |grad n| **** 1231!$OMP DO 1232 do ig=1,ic 1233 agr(ig,2*ispin-1) 1234 > = dsqrt( (rho_prime(ig,1,1)+rho_prime(ig,1,ispin))**2 1235 > + (rho_prime(ig,2,1)+rho_prime(ig,2,ispin))**2 1236 > + (rho_prime(ig,3,1)+rho_prime(ig,3,ispin))**2) 1237 end do 1238!$OMP END DO 1239 1240c **** compute |grad nup| and |grad ndn| **** 1241 if (ispin.eq.2) then 1242 do ms=1,ispin 1243!$OMP DO 1244 do ig=1,ic 1245 agr(ig,ms) = dsqrt( rho_prime(ig,1,ms)**2 1246 > + rho_prime(ig,2,ms)**2 1247 > + rho_prime(ig,3,ms)**2) 1248 end do 1249!$OMP END DO 1250 end do 1251 end if 1252 1253 end if 1254 return 1255 end 1256 1257 1258* ************************************ 1259* * * 1260* * nwpw_xc_gen_dvxc * 1261* * * 1262* ************************************ 1263 subroutine nwpw_xc_gen_dvxc(ic,lmax2,ispin, 1264 > fdn,agr,rho_prime) 1265 implicit none 1266 integer ic,lmax2,ispin 1267 real*8 fdn(ic,*) 1268 real*8 agr(ic,*) 1269 real*8 rho_prime(ic,3,ispin) 1270 1271 integer i,j,lm,ms,jmax 1272 real*8 drho1,drho2,drhoa 1273 real*8 eta 1274 parameter (eta=1.0d-20) 1275 1276 !*** only computate radial derivatives if lmax2=1 **** 1277 if (lmax2.eq.1) then 1278 jmax = 1 1279 else 1280 jmax = 3 1281 end if 1282 1283* *** restricted **** 1284 if (ispin.eq.1) then 1285 do j=1,jmax 1286!$OMP DO 1287 do i=1,ic 1288 rho_prime(i,j,1) = (rho_prime(i,j,1)+rho_prime(i,j,1)) 1289 > *(fdn(i,1)/(agr(i,1)+eta)) 1290 end do 1291!$OMP END DO 1292 end do 1293 1294* *** unrestricted **** 1295 else 1296 do j=1,jmax 1297!$OMP DO 1298 do i=1,ic 1299 drho1 = rho_prime(i,j,1) 1300 drho2 = rho_prime(i,j,2) 1301 drhoa = drho1+drho2 1302 rho_prime(i,j,1) = (drho1/(agr(i,1)+eta))*fdn(i,1) 1303 > + (drhoa/(agr(i,3)+eta))*fdn(i,3) 1304 rho_prime(i,j,2) = (drho2/(agr(i,2)+eta))*fdn(i,2) 1305 > + (drhoa/(agr(i,3)+eta))*fdn(i,3) 1306 end do 1307!$OMP END DO 1308 end do 1309 end if 1310 return 1311 end 1312 1313* ************************************ 1314* * * 1315* * nwpw_xc_addto_dvxclm * 1316* * * 1317* ************************************ 1318 subroutine nwpw_xc_addto_dvxclm(ic,lmax2,ispin, 1319 > alpha, 1320 > tlm, 1321 > dvxc_ae,dvxc_ps, 1322 > dvxclm_ae,dvxclm_ps) 1323 implicit none 1324 integer ic,lmax2,ispin 1325 real*8 alpha 1326 real*8 tlm(*) 1327 real*8 dvxc_ae(ic,3,ispin) 1328 real*8 dvxc_ps(ic,3,ispin) 1329 1330 real*8 dvxclm_ae(ic,lmax2,3,ispin) 1331 real*8 dvxclm_ps(ic,lmax2,3,ispin) 1332 1333 integer i,j,lm,ms,jmax 1334 1335 !*** only computate radial derivatives if lmax2=1 **** 1336 if (lmax2.eq.1) then 1337 jmax = 1 1338 else 1339 jmax = 3 1340 end if 1341 1342 do ms = 1,ispin 1343 do j = 1,jmax 1344 do lm = 1,lmax2 1345!$OMP DO 1346 do i = 1,ic 1347 dvxclm_ae(i,lm,j,ms) = dvxclm_ae(i,lm,j,ms) 1348 > + dvxc_ae(i,j,ms)*(tlm(lm)*alpha) 1349 dvxclm_ps(i,lm,j,ms) = dvxclm_ps(i,lm,j,ms) 1350 > + dvxc_ps(i,j,ms)*(tlm(lm)*alpha) 1351 end do 1352!$OMP END DO 1353 end do 1354 end do 1355 end do 1356 return 1357 end 1358 1359* ************************************ 1360* * * 1361* * nwpw_xc_gen_dmatr * 1362* * * 1363* ************************************ 1364 1365 subroutine nwpw_xc_gen_dmatr(ng,nb,lmax2, 1366 > ispin,log_amesh, 1367 > phi_ae,phi_ps, 1368 > phi_ae_prime,phi_ps_prime, 1369 > r, 1370 > dvxclm_ae,dvxclm_ps, 1371 > tmpC, 1372 > dmatr) 1373 implicit none 1374 integer ng,nb,lmax2,ispin 1375 real*8 log_amesh 1376 1377 real*8 phi_ae(ng,nb) 1378 real*8 phi_ps(ng,nb) 1379 real*8 phi_ae_prime(ng,nb) 1380 real*8 phi_ps_prime(ng,nb) 1381 real*8 r(ng) 1382 real*8 dvxclm_ae(ng,lmax2,3,ispin) 1383 real*8 dvxclm_ps(ng,lmax2,3,ispin) 1384 real*8 tmpC(ng) 1385 1386 real*8 dmatr(nb,nb,lmax2,3,ispin) 1387 1388* **** local varialbles **** 1389 integer ig,i,j,lm,ms 1390 real*8 tmp_ae,tmp_ps 1391 1392* **** external functions **** 1393 real*8 log_integrate_def_destroyf 1394 external log_integrate_def_destroyf 1395 1396 !call dcopy(nb*nb*lmax2*3*ispin,0.0d0,0,dmatr,1) 1397 call Parallel_shared_vector_zero(.true.,nb*nb*lmax2*3*ispin,dmatr) 1398 1399* **** radial derivative integral **** 1400 do ms=1,ispin 1401 do lm=1,lmax2 1402 do i=1,nb 1403 do j=i,nb 1404!$OMP DO 1405 do ig=1,ng 1406 tmp_ae 1407 > = ( phi_ae_prime(ig,i)*phi_ae(ig,j) 1408 > + phi_ae(ig,i)*phi_ae_prime(ig,j)) 1409 > /(r(ig)**2) 1410 > - 2.0d0*phi_ae(ig,i)*phi_ae(ig,j) 1411 > /(r(ig)**3) 1412 tmp_ps 1413 > = ( phi_ps_prime(ig,i)*phi_ps(ig,j) 1414 > + phi_ps(ig,i)*phi_ps_prime(ig,j)) 1415 > /(r(ig)**2) 1416 > - 2.0d0*phi_ps(ig,i)*phi_ps(ig,j) 1417 > /(r(ig)**3) 1418 1419 tmpC(ig) = dvxclm_ae(ig,lm,1,ms)*tmp_ae 1420 > - dvxclm_ps(ig,lm,1,ms)*tmp_ps 1421 end do 1422!$OMP END DO 1423!$OMP MASTER 1424 dmatr(i,j,lm,1,ms) 1425 > = log_integrate_def_destroyf(0,tmpC,2,r,log_amesh,ng) 1426 if (i.ne.j) dmatr(j,i,lm,1,ms) = dmatr(i,j,lm,1,ms) 1427!$OMP END MASTER 1428!$OMP BARRIER 1429 end do 1430 end do 1431 end do 1432 end do 1433 1434 1435* *** only comnpute theta and phi integrals if lmax2>1 **** 1436 if (lmax2.gt.1) then 1437 1438* **** theta derivative integral **** 1439 do ms=1,ispin 1440 do lm=1,lmax2 1441 do i=1,nb 1442 do j=i,nb 1443!$OMP DO 1444 do ig=1,ng 1445 tmp_ae = phi_ae(ig,i) 1446 > *phi_ae(ig,j) 1447 > /(r(ig)**3) 1448 tmp_ps = phi_ps(ig,i) 1449 > *phi_ps(ig-1,j) 1450 > /(r(ig)**3) 1451 tmpC(ig) = dvxclm_ae(ig,lm,2,ms)*tmp_ae 1452 > - dvxclm_ps(ig,lm,2,ms)*tmp_ps 1453 end do 1454!$OMP END DO 1455!$OMP MASTER 1456 dmatr(i,j,lm,2,ms) 1457 > = log_integrate_def_destroyf(0,tmpC,2,r,log_amesh,ng) 1458 if (i.ne.j) dmatr(j,i,lm,2,ms) = dmatr(i,j,lm,2,ms) 1459!$OMP END MASTER 1460!$OMP BARRIER 1461 end do 1462 end do 1463 end do 1464 end do 1465 1466* **** phi derivative integral **** 1467 do ms=1,ispin 1468 do lm=1,lmax2 1469 do i=1,nb 1470 do j=i,nb 1471 1472!$OMP DO 1473 do ig=1,ng 1474 tmp_ae = phi_ae(ig,i) 1475 > *phi_ae(ig,j) 1476 > /(r(ig)**3) 1477 tmp_ps = phi_ps(ig,i) 1478 > *phi_ps(ig,j) 1479 > /(r(ig)**3) 1480 tmpC(ig) = dvxclm_ae(ig,lm,3,ms)*tmp_ae 1481 > - dvxclm_ps(ig,lm,3,ms)*tmp_ps 1482 end do 1483!$OMP END DO 1484!$OMP MASTER 1485 dmatr(i,j,lm,3,ms) 1486 > = log_integrate_def_destroyf(0,tmpC,2,r,log_amesh,ng) 1487 if (i.ne.j) dmatr(j,i,lm,3,ms) = dmatr(i,j,lm,3,ms) 1488!$OMP END MASTER 1489!$OMP BARRIER 1490 end do 1491 end do 1492 end do 1493 end do 1494 1495 end if 1496 1497 return 1498 end 1499 1500 1501* ***************************************************** 1502* * * 1503* * nwpw_get_spher_grid * 1504* * * 1505* ***************************************************** 1506 subroutine nwpw_get_spher_grid(ntheta,nphi,angle_phi, 1507 > cos_theta,w_theta,w_phi) 1508 implicit none 1509 1510 integer ntheta 1511 integer nphi 1512 double precision cos_theta(ntheta) 1513 double precision angle_phi(nphi) 1514 double precision w_theta(ntheta) 1515 double precision w_phi(nphi) 1516 1517* *** local variables *** 1518 integer i 1519 real*8 pi 1520 1521 pi = 4.0d0*datan(1.0d0) 1522 1523* *** gaussian quadrature angular grid for cos_theta *** 1524 call nwpw_gauss_weights(-1.0d0,1.0d0,cos_theta,w_theta,ntheta) 1525 1526 if (nphi.gt.1) then 1527* *** linear angular grid for angle_phi*** 1528 do i=1,nphi 1529 angle_phi(i) = 2.0d0*pi*(i-1)/dble(nphi-1) 1530 w_phi(i) = 2.0d0*pi/dble(nphi-1) 1531 end do 1532 w_phi(1) = 0.5d0*w_phi(1) 1533 w_phi(nphi) = w_phi(1) 1534 else 1535 angle_phi(1) = 0.0d0 1536 w_phi(1) = 2.0d0*pi 1537 end if 1538 return 1539 end 1540 1541 1542 1543c ********************************************* 1544c * * 1545c * nwpw_xc_density_solve2 * 1546c * * 1547c ********************************************* 1548c 1549c Calculates the atomic density lm expansions from the 1550c overlap coefficients. 1551 1552 subroutine nwpw_xc_density_solve2(ispin,ne,nprj,sw1, 1553 > n1dgrid,nbasis,lmax2, 1554 > nindx,lm_rholm, 1555 > iprj_rholm,jprj_rholm, 1556 > bi_rholm,bj_rholm, 1557 > coeff_rholm, 1558 > phi_ae,phi_ps, 1559 > rgrid, 1560 > rholm_ae,rholm_ps) 1561 implicit none 1562 integer ispin,ne(2),nprj 1563 real*8 sw1(ne(1)+ne(2),nprj) 1564 integer n1dgrid,nbasis,lmax2 1565 integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*) 1566 integer bi_rholm(*),bj_rholm(*) 1567 real*8 coeff_rholm(*) 1568 real*8 phi_ae(n1dgrid,nbasis),phi_ps(n1dgrid,nbasis) 1569 real*8 rgrid(n1dgrid) 1570 real*8 rholm_ae(n1dgrid,lmax2,ispin) 1571 real*8 rholm_ps(n1dgrid,lmax2,ispin) 1572 1573 integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2) 1574 real*8 tmp_ms,coeff,w,scal 1575 1576* **** parallel_wshared common block - used for simple OMP reductions *** 1577 real*8 wshared(100),wshared1 1578 common /parallel_wshared/ wshared,wshared1 1579 1580* **** external functions **** 1581 real*8 lattice_omega 1582 external lattice_omega 1583 1584 call nwpw_timing_start(21) 1585 n1(1) = 1 1586 n1(2) = ne(1)+1 1587 n2(1) = ne(1) 1588 n2(2) = ne(1)+ne(2) 1589 scal = 1.0d0/lattice_omega() 1590 1591* ***init to zero*** 1592 !call dcopy(ispin*lmax2*n1dgrid,0.0d0,0,rholm_ae,1) 1593 !call dcopy(ispin*lmax2*n1dgrid,0.0d0,0,rholm_ps,1) 1594 call Parallel_shared_vector_zero(.false.,ispin*lmax2*n1dgrid, 1595 > rholm_ae) 1596 call Parallel_shared_vector_zero(.true.,ispin*lmax2*n1dgrid, 1597 > rholm_ps) 1598 1599 do i=1,nindx 1600 lm = lm_rholm(i) 1601 iprj = iprj_rholm(i) 1602 jprj = jprj_rholm(i) 1603 bi = bi_rholm(i) 1604 bj = bj_rholm(i) 1605 coeff = coeff_rholm(i)*scal 1606 do ms=1,ispin 1607 if (ne(ms).gt.0) then 1608!$OMP MASTER 1609 wshared1 = 0.0d0 1610!$OMP END MASTER 1611!$OMP BARRIER 1612!$OMP DO REDUCTION(+:wshared1) 1613 do n=n1(ms),n2(ms) 1614 wshared1 = wshared1 + sw1(n,iprj)*sw1(n,jprj) 1615 end do 1616!$OMP END DO 1617 call D1dB_SumAll(wshared1) 1618 tmp_ms = wshared1*coeff 1619 call nwpw_xc_density_gen_rho(n1dgrid,tmp_ms, 1620 > rgrid,phi_ae(1,bi),phi_ae(1,bj), 1621 > rholm_ae(1,lm+1,ms)) 1622 call nwpw_xc_density_gen_rho(n1dgrid,tmp_ms, 1623 > rgrid,phi_ps(1,bi),phi_ps(1,bj), 1624 > rholm_ps(1,lm+1,ms)) 1625 end if 1626 end do 1627 end do 1628 1629 1630 call nwpw_timing_end(21) 1631 return 1632 end 1633 1634 1635 subroutine nwpw_xc_density_gen_rho(ic,alpha,r,phi1,phi2,rho) 1636 implicit none 1637 integer ic 1638 double precision alpha 1639 double precision r(ic) 1640 double precision phi1(ic) 1641 double precision phi2(ic) 1642 double precision rho(ic) 1643 1644* **** local variables **** 1645 integer i 1646 double precision tmp 1647 1648!$OMP DO 1649 do i=1,ic 1650 tmp = (phi1(i)*phi2(i))/(r(i)**2) 1651 rho(i) = rho(i) + alpha*tmp 1652 end do 1653!$OMP END DO 1654 1655 return 1656 end 1657 1658 1659c ********************************************* 1660c * * 1661c * nwpw_xc_density_prime_solve2 * 1662c * * 1663c ********************************************* 1664c 1665c Calculates the atomic density lm expansions from the 1666c overlap coefficients. 1667 subroutine nwpw_xc_density_prime_solve2(ispin,ne,nprj,sw1, 1668 > n1dgrid,nbasis,lmax2, 1669 > nindx,lm_rholm, 1670 > iprj_rholm,jprj_rholm, 1671 > bi_rholm,bj_rholm, 1672 > coeff_rholm, 1673 > phi_ae,phi_ps, 1674 > dphi_ae,dphi_ps,rgrid, 1675 > rholm_ae_prime,rholm_ps_prime) 1676 implicit none 1677 integer ispin,ne(2),nprj 1678 real*8 sw1(ne(1)+ne(2),nprj) 1679 integer n1dgrid,nbasis,lmax2 1680 integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*) 1681 integer bi_rholm(*),bj_rholm(*) 1682 real*8 coeff_rholm(*) 1683 real*8 phi_ae(n1dgrid,nbasis),phi_ps(n1dgrid,nbasis) 1684 real*8 dphi_ae(n1dgrid,nbasis),dphi_ps(n1dgrid,nbasis) 1685 real*8 rgrid(n1dgrid) 1686 real*8 rholm_ae_prime(n1dgrid,lmax2,ispin) 1687 real*8 rholm_ps_prime(n1dgrid,lmax2,ispin) 1688 1689 integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2) 1690 real*8 tmp_ms,coeff,w,scal 1691 1692* **** parallel_wshared common block - used for simple OMP reductions *** 1693 real*8 wshared(100),wshared1 1694 common /parallel_wshared/ wshared,wshared1 1695 1696 1697* **** external functions **** 1698 real*8 lattice_omega 1699 external lattice_omega 1700 1701 call nwpw_timing_start(21) 1702 n1(1) = 1 1703 n1(2) = ne(1)+1 1704 n2(1) = ne(1) 1705 n2(2) = ne(1)+ne(2) 1706 scal = 1.0d0/lattice_omega() 1707 1708* ***init to zero*** 1709 !call dcopy(n1dgrid*lmax2*ispin,0.0d0,0,rholm_ae_prime,1) 1710 !call dcopy(n1dgrid*lmax2*ispin,0.0d0,0,rholm_ps_prime,1) 1711 call Parallel_shared_vector_zero(.false.,n1dgrid*lmax2*ispin, 1712 > rholm_ae_prime) 1713 call Parallel_shared_vector_zero(.true.,n1dgrid*lmax2*ispin, 1714 > rholm_ps_prime) 1715 1716 do i=1,nindx 1717 lm = lm_rholm(i) 1718 iprj = iprj_rholm(i) 1719 jprj = jprj_rholm(i) 1720 bi = bi_rholm(i) 1721 bj = bj_rholm(i) 1722 coeff = coeff_rholm(i)*scal 1723 do ms=1,ispin 1724 if (ne(ms).gt.0) then 1725!$OMP MASTER 1726 wshared1 = 0.0d0 1727!$OMP END MASTER 1728!$OMP BARRIER 1729!$OMP DO REDUCTION(+:wshared1) 1730 do n=n1(ms),n2(ms) 1731 wshared1 = wshared1 + sw1(n,iprj)*sw1(n,jprj) 1732 end do 1733!$OMP END DO 1734 call D1dB_SumAll(wshared1) 1735 tmp_ms = wshared1*coeff 1736 call nwpw_xc_density_gen_drho(n1dgrid,tmp_ms, 1737 > rgrid, 1738 > phi_ae(1,bi),dphi_ae(1,bi), 1739 > phi_ae(1,bj),dphi_ae(1,bj), 1740 > rholm_ae_prime(1,lm+1,ms)) 1741 call nwpw_xc_density_gen_drho(n1dgrid,tmp_ms, 1742 > rgrid, 1743 > phi_ps(1,bi),dphi_ps(1,bi), 1744 > phi_ps(1,bj),dphi_ps(1,bj), 1745 > rholm_ps_prime(1,lm+1,ms)) 1746 end if 1747 end do 1748 end do 1749 1750 call nwpw_timing_end(21) 1751 return 1752 end 1753 1754 subroutine nwpw_xc_density_gen_drho(ic,alpha,r, 1755 > phi1,dphi1, 1756 > phi2,dphi2, 1757 > drho) 1758 implicit none 1759 integer ic 1760 double precision alpha 1761 double precision r(ic) 1762 double precision phi1(ic),dphi1(ic) 1763 double precision phi2(ic),dphi2(ic) 1764 double precision drho(ic) 1765 1766* **** local variables **** 1767 integer i 1768 double precision tmp 1769 1770!$OMP DO 1771 do i=1,ic 1772 tmp = (dphi1(i)*phi2(i)+phi1(i)*dphi2(i))/(r(i)**2) 1773 > - 2.0d0*(phi1(i)*phi2(i))/(r(i)**3) 1774 drho(i) = drho(i) + alpha*tmp 1775 end do 1776!$OMP END DO 1777 return 1778 end 1779 1780 1781c ********************************************* 1782c * * 1783c * nwpw_xc_sw1tosw2 * 1784c * * 1785c ********************************************* 1786c 1787 subroutine nwpw_xc_sw1tosw2(ispin,ne,nprj,sw1,sw2, 1788 > n1dgrid,nbasis,lmax2, 1789 > nindx,lm_rholm, 1790 > iprj_rholm,jprj_rholm, 1791 > bi_rholm,bj_rholm, 1792 > coeff_rholm, 1793 > matr) 1794 implicit none 1795 integer ispin,ne(2),nprj 1796 real*8 sw1(ne(1)+ne(2),nprj) 1797 real*8 sw2(ne(1)+ne(2),nprj) 1798 integer n1dgrid,nbasis,lmax2 1799 integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*) 1800 integer bi_rholm(*),bj_rholm(*) 1801 real*8 coeff_rholm(*) 1802 real*8 matr(nbasis,nbasis,lmax2,ispin) 1803 1804 integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2) 1805 real*8 tmp_ms,coeff,w,scal 1806 1807* **** external functions **** 1808 real*8 lattice_omega 1809 external lattice_omega 1810 1811 call nwpw_timing_start(21) 1812 n1(1) = 1 1813 n1(2) = ne(1)+1 1814 n2(1) = ne(1) 1815 n2(2) = ne(1)+ne(2) 1816 !scal = 1.0d0/dsqrt(lattice_omega()) 1817 1818* ***init to zero*** 1819 do i=1,nindx 1820 lm = lm_rholm(i) 1821 iprj = iprj_rholm(i) 1822 jprj = jprj_rholm(i) 1823 bi = bi_rholm(i) 1824 bj = bj_rholm(i) 1825 !coeff = coeff_rholm(i)*scal 1826 coeff = coeff_rholm(i) 1827 do ms=1,ispin 1828 if (ne(ms).gt.0) then 1829!$OMP DO 1830 do n=n1(ms),n2(ms) 1831 sw2(n,iprj) = sw2(n,iprj) 1832 > + coeff*matr(bi,bj,lm+1,ms)*sw1(n,jprj) 1833 end do 1834!$OMP END DO 1835 end if 1836 end do 1837 end do 1838 call nwpw_timing_end(21) 1839 return 1840 end 1841 1842 1843 1844c ********************************************* 1845c * * 1846c * nwpw_xc_sw1tosw2_dmatr * 1847c * * 1848c ********************************************* 1849c 1850 subroutine nwpw_xc_sw1tosw2_dmatr(ispin,ne,nprj,sw1,sw2, 1851 > n1dgrid,nbasis,lmax2, 1852 > nindx,lm_rholm, 1853 > iprj_rholm,jprj_rholm, 1854 > bi_rholm,bj_rholm, 1855 > coeff_rholm, 1856 > coeff_rholm2, 1857 > coeff_rholm3, 1858 > dmatr) 1859 implicit none 1860 integer ispin,ne(2),nprj 1861 real*8 sw1(ne(1)+ne(2),nprj) 1862 real*8 sw2(ne(1)+ne(2),nprj) 1863 integer n1dgrid,nbasis,lmax2 1864 integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*) 1865 integer bi_rholm(*),bj_rholm(*) 1866 real*8 coeff_rholm(*) 1867 real*8 coeff_rholm2(*) 1868 real*8 coeff_rholm3(*) 1869 real*8 dmatr(nbasis,nbasis,lmax2,3,ispin) 1870 1871 integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2) 1872 real*8 tmp_ms,coeff,coeff2,coeff3,w,scal 1873 1874* **** external functions **** 1875 real*8 lattice_omega 1876 external lattice_omega 1877 1878 call nwpw_timing_start(21) 1879 n1(1) = 1 1880 n1(2) = ne(1)+1 1881 n2(1) = ne(1) 1882 n2(2) = ne(1)+ne(2) 1883 !scal = 1.0d0/dsqrt(lattice_omega()) 1884 1885* ***init to zero*** 1886 do i=1,nindx 1887 lm = lm_rholm(i) 1888 iprj = iprj_rholm(i) 1889 jprj = jprj_rholm(i) 1890 bi = bi_rholm(i) 1891 bj = bj_rholm(i) 1892 !coeff = coeff_rholm(i)*scal 1893 coeff = coeff_rholm(i) 1894 coeff2 = coeff_rholm2(i) 1895 coeff3 = coeff_rholm3(i) 1896 do ms=1,ispin 1897 if (ne(ms).gt.0) then 1898!$OMP DO 1899 do n=n1(ms),n2(ms) 1900 sw2(n,iprj) = sw2(n,iprj) 1901 > + coeff *dmatr(bi,bj,lm+1,1,ms)*sw1(n,jprj) 1902 > + coeff2*dmatr(bi,bj,lm+1,2,ms)*sw1(n,jprj) 1903 > + coeff3*dmatr(bi,bj,lm+1,3,ms)*sw1(n,jprj) 1904 end do 1905!$OMP END DO 1906 end if 1907 end do 1908 end do 1909 call nwpw_timing_end(21) 1910 return 1911 end 1912 1913 1914