1* 2* $Id$ 3* 4 5 6* ******************************************* 7* * * 8* * Multiply_Kijl_SO * 9* * * 10* ******************************************* 11ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 12c Create the G(i,j,L)*<psi|(LS)|prj(j,L)> array 13ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 14 subroutine Multiply_Kijl_SO(nn,nprj,nmax,lmax, 15 > n_prj,l_prj,m_prj, 16 > G, 17 > zsw1,zsw2) 18 implicit none 19 integer nn 20 integer nprj,nmax,lmax,nh 21 integer n_prj(nprj) 22 integer l_prj(nprj) 23 integer m_prj(nprj) 24 real*8 G(nmax,nmax,0:lmax),xx,dl2,dm2,dlm 25 complex*16 zsw1(nn,nprj) 26 complex*16 zsw2(nn,nprj) 27 28 !**** local variables **** 29 integer a,b,na,nb,la,lb,ma,mb 30 31 nh=nn/2 32 call dcopy(2*nn*nprj,0.0d0,0,zsw2,1) 33 do b=1,nprj 34 lb = l_prj(b) 35 mb = m_prj(b) 36 37 do a=1,nprj 38 la = l_prj(a) 39 ma = m_prj(a) 40 41 if ((la.eq.lb).and.(ma.eq.mb)) then 42 na = n_prj(a) 43 nb = n_prj(b) 44 xx = 0.5d0*G(nb,na,la)*dble(ma) 45 call daxpy(nn,xx,zsw1(1,a),1,zsw2(1,b),1) 46 call daxpy(nn,(-xx),zsw1(nh+1,a),1,zsw2(nh+1,b),1) 47 end if 48 49 end do 50 end do 51 do b=1,nprj 52 lb = l_prj(b) 53 mb = m_prj(b) 54 do a=1,nprj 55 la = l_prj(a) 56 ma = m_prj(a) 57 if ((la.eq.lb).and.(ma.eq.(mb+1)).and.(mb.ne.lb)) then 58 dl2=dble(la*(la+1)) 59 dm2=dble(ma*mb) 60 dlm=dsqrt(dl2-dm2) 61 na = n_prj(a) 62 nb = n_prj(b) 63 xx = 0.5d0*G(nb,na,la)*dlm 64 call daxpy(nn,xx,zsw1(1,a),1,zsw2(nh+1,b),1) 65 end if 66 end do 67 end do 68 do b=1,nprj 69 lb = l_prj(b) 70 mb = m_prj(b) 71 do a=1,nprj 72 la = l_prj(a) 73 ma = m_prj(a) 74 if ((la.eq.lb).and.(ma.eq.(mb-1)).and.(mb.ne.(-lb))) then 75 dl2=dble(la*(la+1)) 76 dm2=dble(ma*mb) 77 dlm=dsqrt(dl2-dm2) 78 na = n_prj(a) 79 nb = n_prj(b) 80 xx = 0.5d0*G(nb,na,la)*dlm 81 call daxpy(nn,xx,zsw1(nh+1,a),1,zsw2(1,b),1) 82 end if 83 end do 84 end do 85 return 86 end 87 88 89* ******************************************* 90* * * 91* * Multiply_Kijl_SO_x * 92* * * 93* ******************************************* 94ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 95c Create the G(i,j,L)*<psi|(LS)|prj(j,L)> array 96ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 97 subroutine Multiply_Kijl_SO_x(nn,nprj,nmax,lmax, 98 > n_prj,l_prj,m_prj, 99 > G, 100 > zsw1u,zsw1d,zsw2u,zsw2d) 101 implicit none 102 integer nn 103 integer nprj,nmax,lmax 104 integer n_prj(nprj) 105 integer l_prj(nprj) 106 integer m_prj(nprj) 107 real*8 G(nmax,nmax,0:lmax),xx,dl2,dm2,dlm 108 complex*16 zsw1u(nprj),zsw1d(nprj) 109 complex*16 zsw2u(nprj),zsw2d(nprj) 110 111 !**** local variables **** 112 integer a,b,na,nb,la,lb,ma,mb 113 114 call dcopy(2*nn*nprj,0.0d0,0,zsw2u,1) 115 call dcopy(2*nn*nprj,0.0d0,0,zsw2d,1) 116 do b=1,nprj 117 lb = l_prj(b) 118 mb = m_prj(b) 119 120 do a=1,nprj 121 la = l_prj(a) 122 ma = m_prj(a) 123 124 if ((la.eq.lb).and.(ma.eq.mb)) then 125 na = n_prj(a) 126 nb = n_prj(b) 127 xx = 0.5d0*G(nb,na,la)*dble(ma) 128 call daxpy(nn,xx,zsw1u(a),1,zsw2u(b),1) 129 call daxpy(nn,(-xx),zsw1d(a),1,zsw2d(b),1) 130 end if 131 132 end do 133 end do 134 do b=1,nprj 135 lb = l_prj(b) 136 mb = m_prj(b) 137 do a=1,nprj 138 la = l_prj(a) 139 ma = m_prj(a) 140 if ((la.eq.lb).and.(ma.eq.(mb+1)).and.(mb.ne.lb)) then 141 dl2=dble(la*(la+1)) 142 dm2=dble(ma*mb) 143 dlm=dsqrt(dl2-dm2) 144 na = n_prj(a) 145 nb = n_prj(b) 146 xx = 0.5d0*G(nb,na,la)*dlm 147 call daxpy(nn,xx,zsw1u(a),1,zsw2d(b),1) 148 end if 149 end do 150 end do 151 do b=1,nprj 152 lb = l_prj(b) 153 mb = m_prj(b) 154 do a=1,nprj 155 la = l_prj(a) 156 ma = m_prj(a) 157 if ((la.eq.lb).and.(ma.eq.(mb-1)).and.(mb.ne.(-lb))) then 158 dl2=dble(la*(la+1)) 159 dm2=dble(ma*mb) 160 dlm=dsqrt(dl2-dm2) 161 na = n_prj(a) 162 nb = n_prj(b) 163 xx = 0.5d0*G(nb,na,la)*dlm 164 call daxpy(nn,xx,zsw1d(a),1,zsw2u(b),1) 165 end if 166 end do 167 end do 168 return 169 end 170 171* ******************************************* 172* * * 173* * cpsp_v_nonlocal_orb_2com * 174* * * 175* ******************************************* 176 177 subroutine cpsp_v_nonlocal_orb_2com(nb,orb1,orb2) 178 implicit none 179 integer nb 180 complex*16 orb1(*) 181 complex*16 orb2(*) 182 183#include "bafdecls.fh" 184#include "errquit.fh" 185#include "cpsp_common.fh" 186 187 188* *** local variables *** 189 complex*16 one,mone 190c parameter (one=(1.0d0,0.0d0), mone=(-1.0d0,0.0d0)) 191 192 integer nfft3d,npack1,npack,nbrill 193 integer ii,ia,l,prj_shift 194 integer shift,l_prj,nproj,shifts,ne1 195 real*8 omega,scal 196 complex*16 cxr 197 integer exi(2),zsw1u(2),zsw2u(2),zsw1d(2),zsw2d(2) 198 logical value,sd_function 199 200* **** external functions **** 201 logical is_sORd 202 integer ion_nion,ion_katm,brillioun_nbrillioun 203 integer cpsi_ne,cpsp_projector_get_ptr 204 real*8 lattice_omega 205 external is_sORd 206 external ion_nion,ion_katm,brillioun_nbrillioun 207 external lattice_omega,cpsi_ne,cpsp_projector_get_ptr 208 209 one=dcmplx(1.0d0,0.0d0) 210 mone=dcmplx(-1.0d0,0.0d0) 211 212 call nwpw_timing_start(6) 213 214 prj_shift=vso_shift 215* **** allocate local memory **** 216 call C3dB_nfft3d(1,nfft3d) 217 call Cram_max_npack(npack1) 218 nbrill = brillioun_nbrillioun() 219 220 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 221 value = value.and. 222 > BA_push_get(mt_dcpl,nprj_max, 223 > 'zsw1u',zsw1u(2),zsw1u(1)) 224 value = value.and. 225 > BA_push_get(mt_dcpl,nprj_max, 226 > 'zsw2u',zsw2u(2),zsw2u(1)) 227 value = value.and. 228 > BA_push_get(mt_dcpl,nprj_max, 229 > 'zsw1d',zsw1d(2),zsw1d(1)) 230 value = value.and. 231 > BA_push_get(mt_dcpl,nprj_max, 232 > 'zsw2d',zsw2d(2),zsw2d(1)) 233 if (.not.value) 234 > call errquit('cpsp_v_nonlocal_orb2com:pushing stack',0,MA_ERR) 235 236 ne1=cpsi_ne(1) 237 shifts = npack1*ne1 + 1 238 omega = lattice_omega() 239 scal = 1.0d0/omega 240 241 do ii=1,ion_nion() 242 ia=ion_katm(ii) 243 nproj = int_mb(nprj(1)+ia-1) 244 245 if (nproj.gt.0) then 246 247* **** structure factor **** 248 call Cram_npack(nb,npack) 249 call cstrfac_pack(nb,ii,dcpl_mb(exi(1))) 250 call cstrfac_k(ii,nb,cxr) 251c call Cram_c_ZMul(nb,cxr,dcpl_mb(exi(1)),dcpl_mb(exi(1))) 252 call zscal(npack,cxr,dcpl_mb(exi(1)),1) 253 254 do l=1,nproj 255 256c shift = vnl(1)+(l-1) *npack1 257c > +(nb-1)*npack1*vnl_stride 258c > +(ia-1)*npack1*vnl_stride*nbrill 259 shift = cpsp_projector_get_ptr( 260 > int_mb(vnl(1)+ia-1),nb,l) 261 262 l_prj = int_mb(l_projector(1)+(l-1) 263 > +(ia-1)*jmmax_max) 264 265 sd_function=.true. 266 if (mod(l_prj,2).ne.0) then 267 sd_function=.false. 268 end if 269* **** phase factor does not matter therefore **** 270* **** (-i)^l is the same as (i)^l in the **** 271* **** Rayleigh scattering formula **** 272* *** current function is s or d **** 273 if (sd_function) then 274 call Cram_rc_Mul(nb,dbl_mb(shift), 275 > dcpl_mb(exi(1)), 276 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 277* *** current function is p or f **** 278 else 279 call Cram_irc_Mul(nb,dbl_mb(shift), 280 > dcpl_mb(exi(1)), 281 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 282 end if 283 284* **** compute 1Xnproj matrix zsw1 = <orb1|prj> **** 285 call Cram_cc_izdot(nb, 286 > orb1, 287 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 288 > dcpl_mb(zsw1u(1)+(l-1))) 289 call Cram_cc_izdot(nb, 290 > orb1(shifts), 291 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 292 > dcpl_mb(zsw1d(1)+(l-1))) 293 end do !**l** 294 call C3dB_Vector_SumAll((2*nproj),dcpl_mb(zsw1u(1))) 295 call C3dB_Vector_SumAll((2*nproj),dcpl_mb(zsw1d(1))) 296 297* **** zsw2 = Gijl*zsw1 ****** 298 call Multiply_Gijl_zsw1(1, 299 > nproj, 300 > int_mb(nmax(1)+ia-1), 301 > int_mb(lmax(1)+ia-1), 302 > int_mb(n_projector(1) 303 > + (ia-1)*jmmax_max), 304 > int_mb(l_projector(1) 305 > + (ia-1)*jmmax_max), 306 > int_mb(m_projector(1) 307 > + (ia-1)*jmmax_max), 308 > dbl_mb( Gijl(1) 309 > + (ia-1)*gij_stride), 310 > dcpl_mb(zsw1u(1)), 311 > dcpl_mb(zsw2u(1))) 312 call Multiply_Gijl_zsw1(1, 313 > nproj, 314 > int_mb(nmax(1)+ia-1), 315 > int_mb(lmax(1)+ia-1), 316 > int_mb(n_projector(1) 317 > + (ia-1)*jmmax_max), 318 > int_mb(l_projector(1) 319 > + (ia-1)*jmmax_max), 320 > int_mb(m_projector(1) 321 > + (ia-1)*jmmax_max), 322 > dbl_mb(Gijl(1) 323 > + (ia-1)*gij_stride), 324 > dcpl_mb(zsw1d(1)), 325 > dcpl_mb(zsw2d(1))) 326 327* **** do Kleinman-Bylander Multiplication **** 328 call dscal(2*nproj,scal,dcpl_mb(zsw2u(1)),1) 329 call dscal(2*nproj,scal,dcpl_mb(zsw2d(1)),1) 330 call ZGEMM('N','C',npack,1,nproj, 331 > mone, 332 > dcpl_mb(prjtmp(1)), npack1, 333 > dcpl_mb(zsw2u(1)), 1, 334 > one, 335 > orb2, npack1) 336 call ZGEMM('N','C',npack,1,nproj, 337 > mone, 338 > dcpl_mb(prjtmp(1)), npack1, 339 > dcpl_mb(zsw2d(1)), 1, 340 > one, 341 > orb2(shifts), npack1) 342 343 end if !** nproj>0 ** 344 end do !** ii*** 345 346 value = BA_pop_stack(zsw2d(2)) 347 value = value.and.BA_pop_stack(zsw1d(2)) 348 value = value.and.BA_pop_stack(zsw2u(2)) 349 value = value.and.BA_pop_stack(zsw1u(2)) 350 value = value.and.BA_pop_stack(exi(2)) 351 if (.not.value) 352 > call errquit('cpsp_v_nonlocal_orb:popping stack',3,MA_ERR) 353 354 call nwpw_timing_end(6) 355 356 return 357 end 358 359 360 361* *********************************** 362* * * 363* * cpsp_v_spin_orbit * 364* * * 365* *********************************** 366 367 subroutine cpsp_v_spin_orbit(ispin,ne, 368 > psi1_tag,psi2_tag,move,fion) 369 implicit none 370 integer ispin,ne(2) 371 integer psi1_tag 372 integer psi2_tag 373c complex*16 psi1(*),psi2(*) 374 logical move 375 real*8 fion(3,*) 376 377#include "bafdecls.fh" 378#include "cpsp_common.fh" 379#include "errquit.fh" 380 381* *** local variables *** 382 complex*16 one,mone 383 integer nfft3d,G(3),npack1,npack 384 integer ii,ia,l,n,nn,nb,nbq,nbrill 385 integer shift,l_prj,nproj 386 integer psi1_shift,psi2_shift,psi_shift,nshift 387 real*8 omega,weight,scal 388 complex*16 cxr 389 integer exi(2),xtmp(2),zsw1(2),zsw2(2),sum(2) 390 integer Gx(2),Gy(2),Gz(2) 391 logical value,sd_function 392 393* **** external functions **** 394 logical is_sORd 395 integer ion_nion,ion_katm,c_G_indx,Pneb_nbrillq,Pneb_convert_nb 396 integer cpsp_projector_get_ptr,cpsi_data_get_chnk 397 real*8 lattice_omega,brillioun_weight 398 external is_sORd 399 external ion_nion,ion_katm,c_G_indx,Pneb_nbrillq,Pneb_convert_nb 400 external cpsp_projector_get_ptr,cpsi_data_get_chnk 401 external lattice_omega,brillioun_weight 402 403 404 if (.not.do_spin_orbit) return 405 406 one = dcmplx( 1.0d0,0.0d0) 407 mone = dcmplx(-1.0d0,0.0d0) 408 409 call nwpw_timing_start(6) 410 411* **** allocate local memory **** 412 nn = ne(1)+ne(2) 413 nbrill = Pneb_nbrillq() 414 call C3dB_nfft3d(1,nfft3d) 415 call Cram_max_npack(npack1) 416 nshift = 2*npack1 417 418 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 419 value = value.and. 420 > BA_push_get(mt_dcpl,nn*nprj_max, 421 > 'zsw1',zsw1(2),zsw1(1)) 422 value = value.and. 423 > BA_push_get(mt_dcpl,nn*nprj_max, 424 > 'zsw2',zsw2(2),zsw2(1)) 425 if (.not.value) call errquit('cpsp_v_nonlocal:pushing stack',0, 426 & MA_ERR) 427 428 if (move) then 429 value = BA_push_get(mt_dbl, nfft3d,'xtmp',xtmp(2),xtmp(1)) 430 value = value.and. 431 > BA_push_get(mt_dbl, nfft3d,'Gx',Gx(2),Gx(1)) 432 value = value.and. 433 > BA_push_get(mt_dbl, nfft3d,'Gy',Gy(2),Gy(1)) 434 value = value.and. 435 > BA_push_get(mt_dbl, nfft3d,'Gz',Gz(2),Gz(1)) 436 value = value.and. 437 > BA_push_get(mt_dbl,3*nn,'sum',sum(2),sum(1)) 438 if (.not.value) call errquit('cpsp_v_nonlocal:pushing stack',1, 439 & MA_ERR) 440 G(1) = c_G_indx(1) 441 G(2) = c_G_indx(2) 442 G(3) = c_G_indx(3) 443 444* **** define Gx,Gy and Gz in packed space **** 445 call C3dB_t_Copy(1,dbl_mb(G(1)),dbl_mb(Gx(1))) 446 call C3dB_t_Copy(1,dbl_mb(G(2)),dbl_mb(Gy(1))) 447 call C3dB_t_Copy(1,dbl_mb(G(3)),dbl_mb(Gz(1))) 448 end if 449 450 omega = lattice_omega() 451 scal = 1.0d0/omega 452 do 200 ii=1,ion_nion() 453 ia=ion_katm(ii) 454cccccccccc if this atom is not HGH PPOT skip it.... 455 if (int_mb(psp_type(1)+ia-1).ne.1) goto 200 456cccccccccccccccccccccccccccccccccccccccccccc 457 nproj = int_mb(nprj(1)+ia-1) 458 if (nproj.gt.0) then 459 do nbq=1,nbrill 460 nb = Pneb_convert_nb(nbq) 461 psi1_shift = cpsi_data_get_chnk(psi1_tag,nbq) 462 psi2_shift = cpsi_data_get_chnk(psi2_tag,nbq) 463 call Cram_npack(nb,npack) 464 465* **** structure factor pseudopotential **** 466 call cstrfac_pack(nb,ii,dcpl_mb(exi(1))) 467 call cstrfac_k(ii,nb,cxr) 468c call Cram_c_ZMul(nb,cxr,dcpl_mb(exi(1)),dcpl_mb(exi(1))) 469 call zscal(npack,cxr,dcpl_mb(exi(1)),1) 470 471 472* **** generate zsw1's and projectors **** 473 do 105 l=1,nproj 474 475c shift = vnlso(1)+(l-1)*npack1 476c > +(nb-1)*npack1*vso_stride 477c > +(ia-1)*npack1*vso_stride*nbrill 478 shift = cpsp_projector_get_ptr( 479 > int_mb(vnlso(1)+ia-1),nb,l) 480 481 l_prj = int_mb(l_projector(1)+(l-1) 482 > +(ia-1)*jmmax_max) 483 if (l_prj.eq.0) then 484 call dcopy(npack1*2,0.0d0,0, 485 > dcpl_mb(prjtmp(1)+(l-1)*npack1),1) 486 goto 105 487 end if 488 489 sd_function=.true. 490 if (mod(l_prj,2).ne.0) then 491 sd_function=.false. 492 end if 493 494* **** phase factor does not matter therefore **** 495* **** (-i)^l is the same as (i)^l in the **** 496* **** Rayleigh scattering formula **** 497 498c *** current function is s or d **** 499 if (sd_function) then 500 call Cram_cc_Mul(nb,dbl_mb(shift), 501 > dcpl_mb(exi(1)), 502 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 503* *** current function is p or f **** 504 else 505 call Cram_icc_Mul(nb,dbl_mb(shift), 506 > dcpl_mb(exi(1)), 507 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 508 end if 509 510* **** compute nnXnproj matrix zsw1 = <psi1|prj> **** 511 call Cram_cc_inzdot(nb,nn, 512 > dbl_mb(psi1_shift), 513 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 514 > dcpl_mb(zsw1(1)+(l-1)*nn)) 515 516 105 continue 517 call C3dB_Vector_SumAll((2*nn*nproj),dcpl_mb(zsw1(1))) 518 519 520* **** zsw2 = Kijl*zsw1 ****** 521 call Multiply_Kijl_SO(nn, 522 > nproj, 523 > int_mb(nmax(1)+ia-1), 524 > int_mb(lmax(1)+ia-1), 525 > int_mb(n_projector(1) 526 > + (ia-1)*jmmax_max), 527 > int_mb(l_projector(1) 528 > + (ia-1)*jmmax_max), 529 > int_mb(m_projector(1) 530 > + (ia-1)*jmmax_max), 531 > dbl_mb(Kijl(1)+(ia-1)*kij_stride), 532 > dcpl_mb(zsw1(1)), 533 > dcpl_mb(zsw2(1))) 534 535* **** do Kleinman-Bylander Multiplication **** 536 call dscal(2*nn*nproj,scal,dcpl_mb(zsw2(1)),1) 537 call ZGEMM('N','C',npack,nn,nproj, 538 > mone, 539 > dcpl_mb(prjtmp(1)), npack1, 540 > dcpl_mb(zsw2(1)), nn, 541 > one, 542 > dbl_mb(psi2_shift), npack1) 543 544 if (move) then 545 weight = brillioun_weight(nb) 546 if (ispin.eq.1) 547 > call dscal(2*nn*nproj,2.0d0,dcpl_mb(zsw2(1)),1) 548 549 do l=1,nproj 550 551 psi_shift = psi1_shift 552 do n=1,nn 553 554 call Cram_zccr_Multiply2(nb, 555 > dcpl_mb(zsw2(1)+(l-1)*nn+n-1), 556 > dbl_mb(psi_shift), 557 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 558 > dbl_mb(xtmp(1))) 559 psi_shift = psi_shift + nshift 560 561c do i=1,npack 562c ctmp = psi1(i+(n-1)*npack1 563c > +(nb-1)*neall*npack1) 564c > *dconjg(dcpl_mb(prjtmp(1)+(l-1)*npack1+i-1)) 565c > *dconjg(dcpl_mb(zsw2(1)+(l-1)*nn+n-1)) 566c dbl_mb(xtmp(1)+i-1) = dimag(ctmp) 567c end do 568 569* **** define Gx,Gy and Gz in packed space **** 570 call C3dB_t_Copy(1,dbl_mb(G(1)),dbl_mb(Gx(1))) 571 call C3dB_t_Copy(1,dbl_mb(G(2)),dbl_mb(Gy(1))) 572 call C3dB_t_Copy(1,dbl_mb(G(3)),dbl_mb(Gz(1))) 573 call Cram_r_pack(nb,dbl_mb(Gx(1))) 574 call Cram_r_pack(nb,dbl_mb(Gy(1))) 575 call Cram_r_pack(nb,dbl_mb(Gz(1))) 576 call Cram_rr_idot(nb,dbl_mb(Gx(1)),dbl_mb(xtmp(1)), 577 > dbl_mb(sum(1)+3*(n-1))) 578 call Cram_rr_idot(nb,dbl_mb(Gy(1)),dbl_mb(xtmp(1)), 579 > dbl_mb(sum(1)+1+3*(n-1))) 580 call Cram_rr_idot(nb,dbl_mb(Gz(1)),dbl_mb(xtmp(1)), 581 > dbl_mb(sum(1)+2+3*(n-1))) 582 end do !**n** 583 584 call C3dB_Vector_Sumall(3*(nn),dbl_mb(sum(1))) 585 586 do n=1,nn 587 fion(1,ii) = fion(1,ii) 588 > + 2.0d0*weight 589 > *dbl_mb(sum(1)+3*(n-1)) 590 fion(2,ii) = fion(2,ii) 591 > + 2.0d0*weight 592 > *dbl_mb(sum(1)+1+3*(n-1)) 593 fion(3,ii) = fion(3,ii) 594 > + 2.0d0*weight 595 > *dbl_mb(sum(1)+2+3*(n-1)) 596 end do !** nn ** 597 598 end do !** l ** 599 end if !** move ** 600 601 end do !** nb ** 602 end if !** nproj>0** 603 604200 continue !**ii** 605 606 if (move) then 607 value = BA_pop_stack(sum(2)) 608 value = value.and.BA_pop_stack(Gz(2)) 609 value = value.and.BA_pop_stack(Gy(2)) 610 value = value.and.BA_pop_stack(Gx(2)) 611 value = value.and.BA_pop_stack(xtmp(2)) 612 if (.not.value) call errquit('cpsp_v_nonlocal:popping stack',2, 613 & MA_ERR) 614 end if 615 616 value = BA_pop_stack(zsw2(2)) 617 value = value.and.BA_pop_stack(zsw1(2)) 618 value = value.and.BA_pop_stack(exi(2)) 619 if (.not.value) call errquit('cpsp_v_nonlocal:popping stack',3, 620 & MA_ERR) 621 622 call nwpw_timing_end(6) 623 624 return 625 end 626 627 628* *********************************** 629* * * 630* * cpsp_v_spin_orbit0 * 631* * * 632* *********************************** 633 subroutine cpsp_v_spin_orbit0(nbq,ispin,ne, 634 > psi1_tag,psi2_tag,move,fion) 635 implicit none 636 integer nbq 637 integer ispin,ne(2) 638 integer psi1_tag 639 integer psi2_tag 640c complex*16 psi1(*),psi2(*) 641 logical move 642 real*8 fion(3,*) 643 644#include "bafdecls.fh" 645#include "cpsp_common.fh" 646#include "errquit.fh" 647 648* *** local variables *** 649 complex*16 one,mone 650 integer nfft3d,G(3),npack1,npack 651 integer ii,ia,l,n,nn,nb,nbrill 652 integer shift,l_prj,nproj 653 integer psi1_shift,psi2_shift,psi_shift,nshift 654 real*8 omega,weight,scal 655 complex*16 cxr 656 integer exi(2),xtmp(2),zsw1(2),zsw2(2),sum(2) 657 integer Gx(2),Gy(2),Gz(2) 658 logical value,sd_function 659 660* **** external functions **** 661 logical is_sORd 662 integer ion_nion,ion_katm,c_G_indx,Pneb_nbrillq,Pneb_convert_nb 663 integer cpsp_projector_get_ptr,cpsi_data_get_chnk 664 real*8 lattice_omega,brillioun_weight 665 external is_sORd 666 external ion_nion,ion_katm,c_G_indx,Pneb_nbrillq,Pneb_convert_nb 667 external cpsp_projector_get_ptr,cpsi_data_get_chnk 668 external lattice_omega,brillioun_weight 669 670 671 if (.not.do_spin_orbit) return 672 673 one = dcmplx( 1.0d0,0.0d0) 674 mone = dcmplx(-1.0d0,0.0d0) 675 676 call nwpw_timing_start(6) 677 678* **** allocate local memory **** 679 nn = ne(1)+ne(2) 680 nbrill = Pneb_nbrillq() 681 call C3dB_nfft3d(1,nfft3d) 682 call Cram_max_npack(npack1) 683 nshift = 2*npack1 684 685 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 686 value = value.and. 687 > BA_push_get(mt_dcpl,nn*nprj_max, 688 > 'zsw1',zsw1(2),zsw1(1)) 689 value = value.and. 690 > BA_push_get(mt_dcpl,nn*nprj_max, 691 > 'zsw2',zsw2(2),zsw2(1)) 692 if (.not.value) call errquit('cpsp_v_nonlocal:pushing stack',0, 693 & MA_ERR) 694 695 if (move) then 696 value = BA_push_get(mt_dbl, nfft3d,'xtmp',xtmp(2),xtmp(1)) 697 value = value.and. 698 > BA_push_get(mt_dbl, nfft3d,'Gx',Gx(2),Gx(1)) 699 value = value.and. 700 > BA_push_get(mt_dbl, nfft3d,'Gy',Gy(2),Gy(1)) 701 value = value.and. 702 > BA_push_get(mt_dbl, nfft3d,'Gz',Gz(2),Gz(1)) 703 value = value.and. 704 > BA_push_get(mt_dbl,3*nn,'sum',sum(2),sum(1)) 705 if (.not.value) call errquit('cpsp_v_nonlocal:pushing stack',1, 706 & MA_ERR) 707 G(1) = c_G_indx(1) 708 G(2) = c_G_indx(2) 709 G(3) = c_G_indx(3) 710 711* **** define Gx,Gy and Gz in packed space **** 712 call C3dB_t_Copy(1,dbl_mb(G(1)),dbl_mb(Gx(1))) 713 call C3dB_t_Copy(1,dbl_mb(G(2)),dbl_mb(Gy(1))) 714 call C3dB_t_Copy(1,dbl_mb(G(3)),dbl_mb(Gz(1))) 715 end if 716 717 omega = lattice_omega() 718 scal = 1.0d0/omega 719 do 200 ii=1,ion_nion() 720 ia=ion_katm(ii) 721cccccccccc if this atom is not HGH PPOT skip it.... 722 if (int_mb(psp_type(1)+ia-1).ne.1) goto 200 723cccccccccccccccccccccccccccccccccccccccccccc 724 nproj = int_mb(nprj(1)+ia-1) 725 if (nproj.gt.0) then 726 !do nbq=1,nbrill 727 nb = Pneb_convert_nb(nbq) 728 psi1_shift = cpsi_data_get_chnk(psi1_tag,nbq) 729 psi2_shift = cpsi_data_get_chnk(psi2_tag,nbq) 730 call Cram_npack(nb,npack) 731 732* **** structure factor pseudopotential **** 733 call cstrfac_pack(nb,ii,dcpl_mb(exi(1))) 734 call cstrfac_k(ii,nb,cxr) 735c call Cram_c_ZMul(nb,cxr,dcpl_mb(exi(1)),dcpl_mb(exi(1))) 736 call zscal(npack,cxr,dcpl_mb(exi(1)),1) 737 738 739* **** generate zsw1's and projectors **** 740 do 105 l=1,nproj 741 742c shift = vnlso(1)+(l-1)*npack1 743c > +(nb-1)*npack1*vso_stride 744c > +(ia-1)*npack1*vso_stride*nbrill 745 shift = cpsp_projector_get_ptr( 746 > int_mb(vnlso(1)+ia-1),nb,l) 747 748 l_prj = int_mb(l_projector(1)+(l-1) 749 > +(ia-1)*jmmax_max) 750 if (l_prj.eq.0) then 751 call dcopy(npack1*2,0.0d0,0, 752 > dcpl_mb(prjtmp(1)+(l-1)*npack1),1) 753 goto 105 754 end if 755 756 sd_function=.true. 757 if (mod(l_prj,2).ne.0) then 758 sd_function=.false. 759 end if 760 761* **** phase factor does not matter therefore **** 762* **** (-i)^l is the same as (i)^l in the **** 763* **** Rayleigh scattering formula **** 764 765c *** current function is s or d **** 766 if (sd_function) then 767 call Cram_cc_Mul(nb,dbl_mb(shift), 768 > dcpl_mb(exi(1)), 769 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 770* *** current function is p or f **** 771 else 772 call Cram_icc_Mul(nb,dbl_mb(shift), 773 > dcpl_mb(exi(1)), 774 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 775 end if 776 777* **** compute nnXnproj matrix zsw1 = <psi1|prj> **** 778 call Cram_cc_inzdot(nb,nn, 779 > dbl_mb(psi1_shift), 780 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 781 > dcpl_mb(zsw1(1)+(l-1)*nn)) 782 783 105 continue 784 call C3dB_Vector_SumAll((2*nn*nproj),dcpl_mb(zsw1(1))) 785 786 787* **** zsw2 = Kijl*zsw1 ****** 788 call Multiply_Kijl_SO(nn, 789 > nproj, 790 > int_mb(nmax(1)+ia-1), 791 > int_mb(lmax(1)+ia-1), 792 > int_mb(n_projector(1) 793 > + (ia-1)*jmmax_max), 794 > int_mb(l_projector(1) 795 > + (ia-1)*jmmax_max), 796 > int_mb(m_projector(1) 797 > + (ia-1)*jmmax_max), 798 > dbl_mb(Kijl(1)+(ia-1)*kij_stride), 799 > dcpl_mb(zsw1(1)), 800 > dcpl_mb(zsw2(1))) 801 802* **** do Kleinman-Bylander Multiplication **** 803 call dscal(2*nn*nproj,scal,dcpl_mb(zsw2(1)),1) 804 call ZGEMM('N','C',npack,nn,nproj, 805 > mone, 806 > dcpl_mb(prjtmp(1)), npack1, 807 > dcpl_mb(zsw2(1)), nn, 808 > one, 809 > dbl_mb(psi2_shift), npack1) 810 811 if (move) then 812 weight = brillioun_weight(nb) 813 if (ispin.eq.1) 814 > call dscal(2*nn*nproj,2.0d0,dcpl_mb(zsw2(1)),1) 815 816 do l=1,nproj 817 818 psi_shift = psi1_shift 819 do n=1,nn 820 821 call Cram_zccr_Multiply2(nb, 822 > dcpl_mb(zsw2(1)+(l-1)*nn+n-1), 823 > dbl_mb(psi_shift), 824 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 825 > dbl_mb(xtmp(1))) 826 psi_shift = psi_shift + nshift 827 828c do i=1,npack 829c ctmp = psi1(i+(n-1)*npack1 830c > +(nb-1)*neall*npack1) 831c > *dconjg(dcpl_mb(prjtmp(1)+(l-1)*npack1+i-1)) 832c > *dconjg(dcpl_mb(zsw2(1)+(l-1)*nn+n-1)) 833c dbl_mb(xtmp(1)+i-1) = dimag(ctmp) 834c end do 835 836* **** define Gx,Gy and Gz in packed space **** 837 call C3dB_t_Copy(1,dbl_mb(G(1)),dbl_mb(Gx(1))) 838 call C3dB_t_Copy(1,dbl_mb(G(2)),dbl_mb(Gy(1))) 839 call C3dB_t_Copy(1,dbl_mb(G(3)),dbl_mb(Gz(1))) 840 call Cram_r_pack(nb,dbl_mb(Gx(1))) 841 call Cram_r_pack(nb,dbl_mb(Gy(1))) 842 call Cram_r_pack(nb,dbl_mb(Gz(1))) 843 call Cram_rr_idot(nb,dbl_mb(Gx(1)),dbl_mb(xtmp(1)), 844 > dbl_mb(sum(1)+3*(n-1))) 845 call Cram_rr_idot(nb,dbl_mb(Gy(1)),dbl_mb(xtmp(1)), 846 > dbl_mb(sum(1)+1+3*(n-1))) 847 call Cram_rr_idot(nb,dbl_mb(Gz(1)),dbl_mb(xtmp(1)), 848 > dbl_mb(sum(1)+2+3*(n-1))) 849 end do !**n** 850 851 call C3dB_Vector_Sumall(3*(nn),dbl_mb(sum(1))) 852 853 do n=1,nn 854 fion(1,ii) = fion(1,ii) 855 > + 2.0d0*weight 856 > *dbl_mb(sum(1)+3*(n-1)) 857 fion(2,ii) = fion(2,ii) 858 > + 2.0d0*weight 859 > *dbl_mb(sum(1)+1+3*(n-1)) 860 fion(3,ii) = fion(3,ii) 861 > + 2.0d0*weight 862 > *dbl_mb(sum(1)+2+3*(n-1)) 863 end do !** nn ** 864 865 end do !** l ** 866 end if !** move ** 867 868 !end do !** nbq ** 869 end if !** nproj>0** 870 871200 continue !**ii** 872 873 if (move) then 874 value = BA_pop_stack(sum(2)) 875 value = value.and.BA_pop_stack(Gz(2)) 876 value = value.and.BA_pop_stack(Gy(2)) 877 value = value.and.BA_pop_stack(Gx(2)) 878 value = value.and.BA_pop_stack(xtmp(2)) 879 if (.not.value) call errquit('cpsp_v_nonlocal:popping stack',2, 880 & MA_ERR) 881 end if 882 883 value = BA_pop_stack(zsw2(2)) 884 value = value.and.BA_pop_stack(zsw1(2)) 885 value = value.and.BA_pop_stack(exi(2)) 886 if (.not.value) call errquit('cpsp_v_nonlocal:popping stack',3, 887 & MA_ERR) 888 889 call nwpw_timing_end(6) 890 891 return 892 end 893 894 895* ******************************************* 896* * * 897* * cpsp_v_spin_orbit_orb * 898* * * 899* ******************************************* 900 901 subroutine cpsp_v_spin_orbit_orb(nb,orb1,orb2) 902 implicit none 903 integer nb 904 complex*16 orb1(*) 905 complex*16 orb2(*) 906 907#include "bafdecls.fh" 908#include "errquit.fh" 909#include "cpsp_common.fh" 910 911* *** local variables *** 912 complex*16 one,mone 913 parameter (one=(1.0d0,0.0d0), mone=(-1.0d0,0.0d0)) 914 915 integer nfft3d,npack1,npack 916 integer ii,ia,l 917 integer shift,l_prj,nproj,shifts,ne1 918 real*8 omega,scal 919 complex*16 cxr 920 integer exi(2),zsw1u(2),zsw2u(2),zsw1d(2),zsw2d(2) 921 logical value,sd_function 922 923* **** external functions **** 924 logical is_sORd 925 integer ion_nion,ion_katm 926 integer cpsi_ne,cpsp_projector_get_ptr 927 real*8 lattice_omega 928 external is_sORd 929 external ion_nion,ion_katm 930 external lattice_omega,cpsi_ne,cpsp_projector_get_ptr 931 932 call nwpw_timing_start(6) 933 934* **** allocate local memory **** 935 call C3dB_nfft3d(1,nfft3d) 936 call Cram_max_npack(npack1) 937 938 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 939 value = value.and. 940 > BA_push_get(mt_dcpl,nprj_max, 941 > 'zsw1u',zsw1u(2),zsw1u(1)) 942 value = value.and. 943 > BA_push_get(mt_dcpl,nprj_max, 944 > 'zsw2u',zsw2u(2),zsw2u(1)) 945 value = value.and. 946 > BA_push_get(mt_dcpl,nprj_max, 947 > 'zsw1d',zsw1d(2),zsw1d(1)) 948 value = value.and. 949 > BA_push_get(mt_dcpl,nprj_max, 950 > 'zsw2d',zsw2d(2),zsw2d(1)) 951 if (.not.value) 952 > call errquit('cpsp_v_spin_orbit_orb2com:pushing stack',0,MA_ERR) 953 954 ne1=cpsi_ne(1) 955 shifts = npack1*ne1 + 1 956 omega = lattice_omega() 957 scal = 1.0d0/omega 958 959 do 200 ii=1,ion_nion() 960 ia=ion_katm(ii) 961 962 if (int_mb(psp_type(1)+ia-1).ne.1) goto 200 963 964 nproj = int_mb(nprj(1)+ia-1) 965 966 if (nproj.gt.0) then 967 968* **** structure factor **** 969 call Cram_npack(nb,npack) 970 call cstrfac_pack(nb,ii,dcpl_mb(exi(1))) 971 call cstrfac_k(ii,nb,cxr) 972c call Cram_c_ZMul(nb,cxr,dcpl_mb(exi(1)),dcpl_mb(exi(1))) 973 call zscal(npack,cxr,dcpl_mb(exi(1)),1) 974 975 do l=1,nproj 976 977c shift = vnlso(1)+(l-1) *npack1 978c > +(nb-1)*npack1*vso_stride 979c > +(ia-1)*npack1*vso_stride*nbrill 980 shift = cpsp_projector_get_ptr( 981 > int_mb(vnlso(1)+ia-1),nb,l) 982 983 l_prj = int_mb(l_projector(1)+(l-1) 984 > +(ia-1)*jmmax_max) 985 986 sd_function=.true. 987 if (mod(l_prj,2).ne.0) then 988 sd_function=.false. 989 end if 990 991* **** phase factor does not matter therefore **** 992* **** (-i)^l is the same as (i)^l in the **** 993* **** Rayleigh scattering formula **** 994* *** current function is s or d **** 995 if (sd_function) then 996 call Cram_cc_Mul(nb,dbl_mb(shift), 997 > dcpl_mb(exi(1)), 998 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 999* *** current function is p or f **** 1000 else 1001 call Cram_icc_Mul(nb,dbl_mb(shift), 1002 > dcpl_mb(exi(1)), 1003 > dcpl_mb(prjtmp(1)+(l-1)*npack1)) 1004 end if 1005 1006* **** compute 1Xnproj matrix zsw1 = <psi1|prj> **** 1007 call Cram_cc_izdot(nb, 1008 > orb1, 1009 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 1010 > dcpl_mb(zsw1u(1)+(l-1))) 1011 call Cram_cc_izdot(nb, 1012 > orb1(shifts), 1013 > dcpl_mb(prjtmp(1)+(l-1)*npack1), 1014 > dcpl_mb(zsw1d(1)+(l-1))) 1015 end do !**l** 1016 call C3dB_Vector_SumAll((2*nproj),dcpl_mb(zsw1u(1))) 1017 call C3dB_Vector_SumAll((2*nproj),dcpl_mb(zsw1d(1))) 1018 1019* **** zsw2 = Gijl*zsw1 ****** 1020 call Multiply_Kijl_SO_x(1, 1021 > nproj, 1022 > int_mb(nmax(1)+ia-1), 1023 > int_mb(lmax(1)+ia-1), 1024 > int_mb(n_projector(1) 1025 > + (ia-1)*jmmax_max), 1026 > int_mb(l_projector(1) 1027 > + (ia-1)*jmmax_max), 1028 > int_mb(m_projector(1) 1029 > + (ia-1)*jmmax_max), 1030 > dbl_mb(Kijl(1) 1031 > + (ia-1)*kij_stride), 1032 > dcpl_mb(zsw1u(1)),dcpl_mb(zsw1d(1)), 1033 > dcpl_mb(zsw2u(1)),dcpl_mb(zsw2d(1))) 1034* **** do Kleinman-Bylander Multiplication **** 1035 call dscal(2*nproj,scal,dcpl_mb(zsw2u(1)),1) 1036 call dscal(2*nproj,scal,dcpl_mb(zsw2d(1)),1) 1037 call ZGEMM('N','C',npack,1,nproj, 1038 > mone, 1039 > dcpl_mb(prjtmp(1)), npack1, 1040 > dcpl_mb(zsw2u(1)), 1, 1041 > one, 1042 > orb2, npack1) 1043 call ZGEMM('N','C',npack,1,nproj, 1044 > mone, 1045 > dcpl_mb(prjtmp(1)), npack1, 1046 > dcpl_mb(zsw2d(1)), 1, 1047 > one, 1048 > orb2(shifts), npack1) 1049 1050 end if !** nproj>0 ** 1051200 continue !** ii*** 1052 1053 value = BA_pop_stack(zsw2d(2)) 1054 value = value.and.BA_pop_stack(zsw1d(2)) 1055 value = value.and.BA_pop_stack(zsw2u(2)) 1056 value = value.and.BA_pop_stack(zsw1u(2)) 1057 value = value.and.BA_pop_stack(exi(2)) 1058 if (.not.value) 1059 > call errquit('cpsp_v_spin_orbit_orb:popping stack',3,MA_ERR) 1060 1061 call nwpw_timing_end(6) 1062 1063 return 1064 end 1065 1066* *********************************************** 1067* * * 1068* * cpsp_v_nonlocal_rel * 1069* * * 1070* *********************************************** 1071 1072c Dirac Two Component Relativistic 1073c Non-Local Pseudopotential 1074c pjn --- use at your own peril 1075 1076 subroutine cpsp_v_nonlocal_rel(nb,ii,ia,nproj,npack1,nbrill,nn,ne, 1077 > exi,prtmp,zsw1,psi1,psi2, 1078 > Gij,bz_vol) 1079 implicit none 1080 integer nb,nn,npack1,nproj,ia,ii,nbrill,ne(2) 1081 complex*16 prtmp(*) 1082 complex*16 zsw1(*) 1083 complex*16 exi(*) 1084 complex*16 psi1(*) 1085 complex*16 psi2(*) 1086 real*8 bz_vol,Gij(*) 1087 1088#include "cpsp_common.fh" 1089#include "bafdecls.fh" 1090 1091* **** local variables **** 1092 complex*16 one,mone,ione 1093 parameter(one=(1.0d0,0.0d0),mone=(-1.0d0,0.0d0)) 1094 parameter(ione=(0.0d0,1.0d0)) 1095 1096 complex*16 cxr 1097 real*8 xx 1098 integer l,shifts,shiftxu,shiftxd,l_prj,npack,ne1 1099 integer zshft,vshiftu,vshiftd 1100 logical sd_function,val 1101 1102* **** external functions **** 1103 integer cpsp_projector_get_ptr 1104 external cpsp_projector_get_ptr 1105 1106cccccccccccccccccccccccccccccccccccccccccc 1107c These should be declared as parameters 1108c however there seems to be no "standard" 1109c way to declare a complex parameter that 1110c works with all compilers...gfortran,ifort, 1111c pf90,etc. 1112ccccccccccccccccccccccccccccccccccccccccc 1113c one = dcmplx( 1.0d0,0.0d0) 1114c mone = dcmplx(-1.0d0,0.0d0) 1115c ione = dcmplx( 0.0d0,1.0d0) 1116 1117 ne1=ne(1) 1118 shifts=1+npack1*ne1 1119 call Cram_npack(nb,npack) 1120 1121* **** structure factor pseudopotential **** 1122 call cstrfac_pack(nb,ii,exi) 1123 call cstrfac_k(ii,nb,cxr) 1124 call zscal(npack,cxr,exi,1) 1125 1126 1127* **** generate zsw1's and projectors **** 1128 do 105 l=1,nproj 1129 1130 vshiftu = cpsp_projector_get_ptr( 1131 > int_mb(vnl(1)+ia-1),nb,l) 1132 vshiftd = cpsp_projector_get_ptr( 1133 > int_mb(vnlso(1)+ia-1),nb,l) 1134 shiftxu=1 + (l-1)*npack1 1135 shiftxd=shiftxu + nproj*npack1 1136 zshft=1+(l-1)*ne1 1137 l_prj = int_mb(l_projector(1)+(l-1) 1138 > +(ia-1)*jmmax_max) 1139 sd_function=.true. 1140 if (mod(l_prj,2).ne.0) then 1141 sd_function=.false. 1142 end if 1143* **** phase factor does not matter therefore **** 1144* **** (-i)^l is the same as (i)^l in the **** 1145* **** Rayleigh scattering formula **** 1146 1147c *** current function is s or d **** 1148 if (sd_function) then 1149 call Cram_cc_Mul(nb,dbl_mb(vshiftu), 1150 > exi, 1151 > prtmp(shiftxu)) 1152 call Cram_cc_Mul(nb,dbl_mb(vshiftd), 1153 > exi, 1154 > prtmp(shiftxd)) 1155* *** current function is p or f **** 1156 else 1157 call Cram_icc_Mul(nb,dbl_mb(vshiftu), 1158 > exi, 1159 > prtmp(shiftxu)) 1160 1161 call Cram_icc_Mul(nb,dbl_mb(vshiftd), 1162 > exi, 1163 > prtmp(shiftxd)) 1164 end if 1165 1166* **** compute nnXnproj matrix zsw1 = <psi1|prj> **** 1167 call Cram_cc_inzdot(nb,ne1, 1168 > psi1, 1169 > prtmp(shiftxu), 1170 > zsw1(zshft)) 1171 call Cram_cc_inzdotAdd(nb,ne1, 1172 > psi1(shifts), 1173 > prtmp(shiftxd), 1174 > zsw1(zshft)) 1175 xx=Gij(l)*bz_vol 1176 call dscal(2*ne1,xx,zsw1(zshft),1) 1177 105 continue 1178 1179 1180 call C3dB_Vector_SumAll((2*ne1*nproj),zsw1) 1181 1182 call ZGEMM('N','C',npack,ne1,nproj, 1183 > mone, 1184 > prtmp, npack1, 1185 > zsw1, ne1, 1186 > one, 1187 > psi2(1), npack1) 1188 call ZGEMM('N','C',npack,ne1,nproj, 1189 > mone, 1190 > prtmp(1+npack1*nproj), npack1, 1191 > zsw1, ne1, 1192 > one, 1193 > psi2(shifts), npack1) 1194 1195 return 1196 end 1197 1198 1199ccccccccccc 1200* *********************************** 1201* * * 1202* * cpsp_v_nonlocal_orb * 1203* * * 1204* *********************************** 1205 subroutine cpsp_v_nonlocal_rel_orb(nb,orb1,orb2, 1206 > zsw1,Gij,exi,nfft3d,ia,ii,ne1,npack1,nproj) 1207 implicit none 1208 integer nb 1209 real*8 Gij(*) 1210 complex*16 orb1(*) 1211 complex*16 orb2(*) 1212 1213#include "bafdecls.fh" 1214#include "errquit.fh" 1215#include "cpsp_common.fh" 1216 1217 1218* *** local variables *** 1219 complex*16 one,mone,ione 1220 integer nfft3d,npack1,npack 1221 integer ii,ia,l,ne1 1222 integer l_prj,nproj 1223 integer shiftu,shiftd,shifts,shiftux,shiftdx 1224 real*8 omega,scal 1225 complex*16 cxr 1226 complex*16 exi(*) 1227 complex*16 zsw1(*) 1228 logical sd_function 1229 1230* **** external functions **** 1231 logical is_sORd,cpsi_spin_orbit 1232 integer ion_nion,ion_katm 1233 integer cpsp_projector_get_ptr 1234 real*8 lattice_omega 1235 external is_sORd,cpsi_spin_orbit 1236 external ion_nion,ion_katm 1237 external cpsp_projector_get_ptr 1238 external lattice_omega 1239 1240 one=dcmplx(1.0d0,0.0d0) 1241 mone=dcmplx(-1.0d0,0.d0) 1242 ione=dcmplx(0.0d0,1.d0) 1243 1244 omega = lattice_omega() 1245 scal = 1.0d0/omega 1246 shifts= ne1*npack1+1 1247* **** structure factor **** 1248 1249 call Cram_npack(nb,npack) 1250 call cstrfac_pack(nb,ii,exi) 1251 call cstrfac_k(ii,nb,cxr) 1252c call Cram_c_ZMul(nb,cxr,exi,exi) 1253 call zscal(npack,cxr,exi,1) 1254 1255 do l=1,nproj 1256 shiftu = cpsp_projector_get_ptr( 1257 > int_mb(vnl(1)+ia-1),nb,l) 1258 shiftd = cpsp_projector_get_ptr( 1259 > int_mb(vnlso(1)+ia-1),nb,l) 1260 l_prj = int_mb(l_projector(1)+(l-1) 1261 > +(ia-1)*jmmax_max) 1262 shiftux=(l-1)*npack1 1263 shiftdx=shiftux+npack1*nproj 1264 sd_function=.true. 1265 if (mod(l_prj,2).ne.0) then 1266 sd_function=.false. 1267 end if 1268 1269* **** phase factor does not matter therefore **** 1270* **** (-i)^l is the same as (i)^l in the **** 1271* **** Rayleigh scattering formula **** 1272* *** current function is s or d **** 1273 if (sd_function) then 1274 call Cram_cc_Mul(nb,dbl_mb(shiftu),exi, 1275 > dcpl_mb(prjtmp(1)+shiftux)) 1276 call Cram_cc_Mul(nb,dbl_mb(shiftd),exi, 1277 > dcpl_mb(prjtmp(1)+shiftdx)) 1278 1279* *** current function is p or f **** 1280 else 1281 call Cram_icc_Mul(nb,dbl_mb(shiftu),exi, 1282 > dcpl_mb(prjtmp(1)+shiftux)) 1283 call Cram_icc_Mul(nb,dbl_mb(shiftd),exi, 1284 > dcpl_mb(prjtmp(1)+shiftdx)) 1285 end if 1286 1287* **** compute 1Xnproj matrix zsw1 = <psi1|prj> **** 1288 call Cram_cc_izdot(nb, 1289 > orb1, 1290 > dcpl_mb(prjtmp(1)+shiftux), 1291 > zsw1(1+(l-1))) 1292 call Cram_cc_izdotAdd(nb, 1293 > orb1(shifts), 1294 > dcpl_mb(prjtmp(1)+shiftdx), 1295 > zsw1(1+(l-1))) 1296 zsw1(1+(l-1))=zsw1(1+(l-1))*Gij(l)*scal 1297 end do !**l** 1298 1299 call C3dB_Vector_SumAll((2*nproj),zsw1) 1300 1301* **** do Kleinman-Bylander Multiplication **** 1302 call ZGEMM('N','C',npack,1,nproj, 1303 > mone, 1304 > dcpl_mb(prjtmp(1)), npack1, 1305 > zsw1, 1, 1306 > one, 1307 > orb2, npack1) 1308 1309 call ZGEMM('N','C',npack,1,nproj, 1310 > mone, 1311 > dcpl_mb(prjtmp(1)+nproj*npack1), npack1, 1312 > zsw1, 1, 1313 > one, 1314 > orb2(shifts), npack1) 1315 1316 return 1317 end 1318ccccccccccccccccc 1319 subroutine cpsp_f_nonlocal_rel(nb,ii,ia,nproj,npack1,ne1, 1320 > exi,zsw1,psi,prtmp,xtmp,sum1,Gij,Gx,Gy,Gz, 1321 > fx,fy,fz,weight,scal) 1322 implicit none 1323 integer ia,ii,nproj,ne1,npack1,nb,npack 1324 complex*16 exi(*),zsw1(*),psi(*),prtmp(*),cxr 1325 real*8 xtmp(*),sum1(*),fx,fy,fz,xx 1326 real*8 weight,scal,Gij(*),Gx(*),Gy(*),Gz(*) 1327#include "cpsp_common.fh" 1328#include "bafdecls.fh" 1329ccccccccccc locals 1330 integer shiftxu,shiftxd,vshiftu,vshiftd,zshft 1331 integer l_prj,l,n,pshftu,pshftd,sshft,shifts 1332 logical sd_function 1333ccccccccccc external 1334 integer cpsp_projector_get_ptr 1335 external cpsp_projector_get_ptr 1336cccccccccccccccccccccccccccccccccccccccccccc 1337* **** structure factor **** 1338 call Cram_npack(nb,npack) 1339 call cstrfac_pack(nb,ii,exi) 1340 call cstrfac_k(ii,nb,cxr) 1341c call Cram_c_ZMul(nb,cxr,exi,exi) 1342 call zscal(npack,cxr,exi,1) 1343 1344 shifts=1+npack1*ne1 1345 do l=1,nproj 1346 vshiftu = cpsp_projector_get_ptr( 1347 > int_mb(vnl(1)+ia-1),nb,l) 1348 vshiftd = cpsp_projector_get_ptr( 1349 > int_mb(vnlso(1)+ia-1),nb,l) 1350 shiftxu=1 + (l-1)*npack1 1351 shiftxd=shiftxu + nproj*npack1 1352 1353 l_prj = int_mb(l_projector(1)+(l-1) 1354 > +(ia-1)*jmmax_max) 1355 1356 sd_function=.true. 1357 if (mod(l_prj,2).ne.0) then 1358 sd_function=.false. 1359 end if 1360 1361* **** phase factor does not matter therefore **** 1362* **** (-i)^l is the same as (i)^l in the **** 1363* **** Rayleigh scattering formula **** 1364* *** current function is s or d **** 1365 if (sd_function) then 1366 call Cram_cc_Mul(nb,dbl_mb(vshiftu),exi,prtmp(shiftxu)) 1367 call Cram_cc_Mul(nb,dbl_mb(vshiftd),exi,prtmp(shiftxd)) 1368 else 1369 call Cram_icc_Mul(nb,dbl_mb(vshiftu),exi,prtmp(shiftxu)) 1370 call Cram_icc_Mul(nb,dbl_mb(vshiftd),exi,prtmp(shiftxd)) 1371 end if 1372 1373 zshft=1+(l-1)*ne1 1374 1375 call Cram_cc_inzdot(nb,ne1,psi,prtmp(shiftxu), 1376 > zsw1(zshft)) 1377 call Cram_cc_inzdotAdd(nb,ne1,psi(shifts),prtmp(shiftxd), 1378 > zsw1(zshft)) 1379 1380 xx=Gij(l)*scal 1381 call dscal(ne1*2,xx,zsw1(zshft),1) 1382 1383 end do 1384 call C3db_Vector_SumAll(2*ne1*nproj,zsw1) 1385 do l=1,nproj 1386 shiftxu=1+(l-1)*npack1 1387 shiftxd=shiftxu+nproj*npack1 1388 do n=1,ne1 1389 zshft=(l-1)*ne1+n 1390 pshftu=1+(n-1)*npack1 1391 pshftd=pshftu+ne1*npack1 1392 call Cram_zccr_Multiply2(nb, 1393 > zsw1(zshft), 1394 > psi(pshftu), 1395 > prtmp(shiftxu), 1396 > xtmp) 1397 call Cram_zccr_Multiply2Add(nb, 1398 > zsw1(zshft), 1399 > psi(pshftd), 1400 > prtmp(shiftxd), 1401 > xtmp) 1402 1403 sshft=1+3*(n-1) 1404 call Cram_rr_idot(nb,Gx,xtmp,sum1(sshft)) 1405 call Cram_rr_idot(nb,Gy,xtmp,sum1(sshft+1)) 1406 call Cram_rr_idot(nb,Gz,xtmp,sum1(sshft+2)) 1407 end do 1408 call C3db_Vector_Sumall(3*ne1,sum1) 1409 do n=1,ne1 1410 sshft=1+3*(n-1) 1411 fx=fx+2.0d0*weight*sum1(sshft) 1412 fy=fy+2.0d0*weight*sum1(sshft+1) 1413 fz=fz+2.0d0*weight*sum1(sshft+2) 1414 end do 1415 end do 1416 return 1417 end 1418ccccccccccccccccccccccccc 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428