1* 2* $Id$ 3* 4 5 6 subroutine integrate_kbppv3(version,rlocal, 7 > nrho,drho,lmax,locp,zv, 8 > vp,wp,rho,f,cs,sn, 9 > nfft1,nfft2,nfft3,lmmax, 10 > G,vl,vnl, 11 > n_prj,l_prj,m_prj,b_prj,vnlnrm, 12 > semicore,rho_sc_r,rho_sc_k, 13 > ierr) 14 implicit none 15 integer version 16 double precision rlocal 17 integer nrho 18 double precision drho 19 integer lmax 20 integer locp 21 double precision zv 22 double precision vp(nrho,0:lmax) 23 double precision wp(nrho,0:lmax) 24 double precision rho(nrho) 25 double precision f(nrho) 26 double precision cs(nrho) 27 double precision sn(nrho) 28 29 integer nfft1,nfft2,nfft3,lmmax 30 double precision G(nfft1/2+1,nfft2,nfft3,3) 31 double precision vl(nfft1/2+1,nfft2,nfft3) 32 double precision vnl(nfft1/2+1,nfft2,nfft3,lmmax) 33 integer n_prj(lmmax),l_prj(lmmax),m_prj(lmmax) 34 integer b_prj(lmmax) 35 double precision vnlnrm(0:lmax) 36 37 logical semicore 38 double precision rho_sc_r(nrho,2) 39 double precision rho_sc_k(nfft1/2+1,nfft2,nfft3,4) 40 41 integer ierr 42 43 integer np,taskid,MASTER 44 parameter (MASTER=0) 45 46* *** local variables **** 47 logical fast_erf 48 integer lcount,task_count,nfft3d 49 integer k1,k2,k3,i,l 50 double precision pi,twopi,forpi 51 double precision p0,p1,p2,p3,p 52 double precision gx,gy,gz,a,q,d 53 54* **** Error function parameters **** 55 real*8 yerf,xerf 56 real*8 c1,c2,c3,c4,c5,c6 57 parameter (c1=0.07052307840d0,c2=0.04228201230d0) 58 parameter (c3=0.00927052720d0) 59 parameter (c4=0.00015201430d0,c5=0.00027656720d0) 60 parameter (c6=0.00004306380d0) 61 62* **** external functions **** 63 logical control_fast_erf 64 external control_fast_erf 65 double precision dsum,simp,util_erf 66 external dsum,simp,util_erf 67 68 call Parallel_np(np) 69 call Parallel_taskid(taskid) 70 71 fast_erf = control_fast_erf() 72 73 nfft3d = (nfft1/2+1)*nfft2*nfft3 74 pi=4.0d0*datan(1.0d0) 75 twopi=2.0d0*pi 76 forpi=4.0d0*pi 77 78 IF(LMMAX.GT.16) THEN 79 IERR=1 80 RETURN 81 ENDIF 82 IF((NRHO/2)*2.EQ.NRHO) THEN 83 IERR=2 84 RETURN 85 ENDIF 86 87 P0=DSQRT(FORPI) 88 P1=DSQRT(3.0d0*FORPI) 89 P2=DSQRT(15.0d0*FORPI) 90 P3=DSQRT(105.0d0*FORPI) 91 92*:::::::::::::::::: Define non-local pseudopotential :::::::::::::::: 93 do l=0,lmax 94 if (l.ne.locp) then 95 do I=1,nrho 96 vp(i,l)=vp(i,l)-vp(i,locp) 97 end do 98 end if 99 end do 100 101*::::::::::::::::::::: Normarization constants :::::::::::::::::::::: 102 lcount = 0 103 do l=0,lmax 104 if (l.ne.locp) then 105 do i=1,nrho 106 f(I)=vp(I,L)*wp(I,L)**2 107 end do 108 a=simp(nrho,f,drho) 109 vnlnrm(l) = (1.0d0/a) 110 else 111 vnlnrm(l) = 0.0d0 112 end if 113 end do 114 115*====================== Fourier transformation ====================== 116 call dcopy(nfft3d,0.0d0,0,vl,1) 117 call dcopy(lmmax*nfft3d,0.0d0,0,vnl,1) 118 call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1) 119 task_count = -1 120 DO 700 k3=1,nfft3 121 DO 700 k2=1,nfft2 122 DO 700 k1=1,(nfft1/2+1) 123 task_count = task_count + 1 124 if (mod(task_count,np).ne.taskid) go to 700 125 126 Q=DSQRT(G(k1,k2,k3,1)**2 127 > +G(k1,k2,k3,2)**2 128 > +G(k1,k2,k3,3)**2) 129 130 if ((k1.eq.1).and.(k2.eq.1).and.(k3.eq.1)) go to 700 131 132 133 GX=G(k1,k2,k3,1)/Q 134 GY=G(k1,k2,k3,2)/Q 135 GZ=G(k1,k2,k3,3)/Q 136 DO I=1,NRHO 137 CS(I)=DCOS(Q*RHO(I)) 138 SN(I)=DSIN(Q*RHO(I)) 139 END DO 140 141 lcount = lmmax+1 142 GO TO (500,400,300,200), LMAX+1 143 144 145*:::::::::::::::::::::::::::::: f-wave :::::::::::::::::::::::::::::: 146 200 CONTINUE 147 if (locp.ne.3) then 148 F(1)=0.0d0 149 do I=2,NRHO 150 A=SN(I)/(Q*RHO(I)) 151 A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I) 152 F(I)=A*WP(I,3)*VP(I,3) 153 end do 154 D=P3*SIMP(NRHO,F,DRHO)/Q 155 156 lcount = lcount-1 157 vnl(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY) 158 > /dsqrt(24.0d0) 159 lcount = lcount-1 160 vnl(k1,k2,k3,lcount)=D*GX*GY*GZ 161 lcount = lcount-1 162 vnl(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0) 163 > /dsqrt(40.0d0) 164 lcount = lcount-1 165 vnl(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0) 166 > /dsqrt(60.0d0) 167 lcount = lcount-1 168 vnl(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0) 169 > /dsqrt(40.0d0) 170 lcount = lcount-1 171 vnl(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY) 172 > /2.0d0 173 lcount = lcount-1 174 vnl(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ)) 175 > /dsqrt(24.0d0) 176c lcount = lcount-1 177c vnl(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ)) 178c > /dsqrt(24.0d0) 179c lcount = lcount-1 180c vnl(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY) 181c > /dsqrt(24.0d0) 182c lcount = lcount-1 183c vnl(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY) 184c > /2.0d0 185c lcount = lcount-1 186c vnl(k1,k2,k3,lcount)=D*GX*GY*GZ 187c lcount = lcount-1 188c vnl(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0) 189c > /dsqrt(40.0d0) 190c lcount = lcount-1 191c vnl(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0) 192c > /dsqrt(40.0d0) 193c lcount = lcount-1 194c vnl(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0) 195c > /dsqrt(60.0d0) 196 end if 197 198 199 200*:::::::::::::::::::::::::::::: d-wave :::::::::::::::::::::::::::::: 201 300 CONTINUE 202 if (locp.ne.2) then 203 F(1)=0.0d0 204 DO I=2,NRHO 205 A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I) 206 F(I)=A*WP(I,2)*VP(I,2) 207 END DO 208 D=P2*SIMP(NRHO,F,DRHO)/Q 209 210 lcount = lcount-1 211 vnl(k1,k2,k3,lcount)=D*GX*GY 212 lcount = lcount-1 213 vnl(k1,k2,k3,lcount)=D*GY*GZ 214 lcount = lcount-1 215 vnl(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0) 216 > /(2.0d0*dsqrt(3.0d0)) 217 lcount = lcount-1 218 vnl(k1,k2,k3,lcount)=D*GZ*GX 219 lcount = lcount-1 220 vnl(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0) 221 222c lcount = lcount-1 223c vnl(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0) 224c > /(2.0d0*dsqrt(3.0d0)) 225c lcount = lcount-1 226c vnl(k1,k2,k3,lcount)=D*GX*GY 227c lcount = lcount-1 228c vnl(k1,k2,k3,lcount)=D*GY*GZ 229c lcount = lcount-1 230c vnl(k1,k2,k3,lcount)=D*GZ*GX 231c lcount = lcount-1 232c vnl(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0) 233 end if 234 235*:::::::::::::::::::::::::::::: p-wave :::::::::::::::::::::::::::::: 236 400 CONTINUE 237 if (locp.ne.1) then 238 F(1)=0.0d0 239 DO I=2,NRHO 240 F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1)*VP(I,1) 241 END DO 242 P=P1*SIMP(NRHO,F,DRHO)/Q 243 lcount = lcount-1 244 vnl(k1,k2,k3,lcount)=P*GY 245 lcount = lcount-1 246 vnl(k1,k2,k3,lcount)=P*GZ 247 lcount = lcount-1 248 vnl(k1,k2,k3,lcount)=P*GX 249c lcount = lcount-1 250c vnl(k1,k2,k3,lcount)=P*GX 251c lcount = lcount-1 252c vnl(k1,k2,k3,lcount)=P*GY 253c lcount = lcount-1 254c vnl(k1,k2,k3,lcount)=P*GZ 255 end if 256 257*:::::::::::::::::::::::::::::: s-wave ::::::::::::::::::::::::::::::: 258 500 CONTINUE 259 if (locp.ne.0) then 260 DO I=1,NRHO 261 F(I)=SN(I)*WP(I,0)*VP(I,0) 262 END DO 263 lcount = lcount-1 264 vnl(k1,k2,k3,lcount)=P0*SIMP(NRHO,F,DRHO)/Q 265 end if 266 267*:::::::::::::::::::::::::::::: local ::::::::::::::::::::::::::::::: 268 600 CONTINUE 269 270 271 if (version.eq.3) then 272 DO I=1,NRHO 273 F(I)=RHO(I)*VP(I,locp)*SN(I) 274 END DO 275 vl(k1,k2,k3)=SIMP(NRHO,F,DRHO)*FORPI/Q-ZV*FORPI*CS(NRHO)/(Q*Q) 276 end if 277 278 if (version.eq.4) then 279 if (fast_erf) then 280 do I=1,NRHO 281 xerf=RHO(I)/rlocal 282 yerf = (1.0d0 283 > + xerf*(c1 + xerf*(c2 284 > + xerf*(c3 + xerf*(c4 285 > + xerf*(c5 + xerf*c6))))))**4 286 yerf = (1.0d0 - 1.0d0/yerf**4) 287 F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I) 288 end do 289 else 290 do I=1,NRHO 291 xerf=RHO(I)/rlocal 292 yerf = util_erf(xerf) 293 F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I) 294 end do 295 end if 296 vl(k1,k2,k3)=SIMP(NRHO,F,DRHO)*FORPI/Q 297 end if 298 299 300*::::::::::::::::::::: semicore density ::::::::::::::::::::::::::::::: 301 if (semicore) then 302 do i=1,nrho 303 f(i) = rho(i)*dsqrt(rho_sc_r(i,1))*sn(i) 304 end do 305 rho_sc_k(k1,k2,k3,1) = SIMP(nrho,f,drho)*forpi/Q 306 307 do i=1,nrho 308 f(i)=(sn(i)/(Q*rho(i))-cs(i))*rho_sc_r(i,2)*rho(i) 309 end do 310 P = SIMP(nrho,f,drho)*forpi/Q 311 rho_sc_k(k1,k2,k3,2)=P*GX 312 rho_sc_k(k1,k2,k3,3)=P*GY 313 rho_sc_k(k1,k2,k3,4)=P*GZ 314 315 end if 316 317 700 CONTINUE 318 319 call Parallel_Vector_SumAll(4*nfft3d,rho_sc_k) 320 call Parallel_Vector_SumAll(nfft3d,vl) 321 call Parallel_Vector_Sumall(lmmax*nfft3d,vnl) 322*::::::::::::::::::::::::::::::: G=0 :::::::::::::::::::::::::::::::: 323 324 if (version.eq.3) then 325 DO I=1,NRHO 326 F(I)=VP(I,locp)*RHO(I)**2 327 END DO 328 vl(1,1,1)=FORPI*SIMP(NRHO,F,DRHO)+TWOPI*ZV*RHO(NRHO)**2 329 end if 330 331 if (version.eq.4) then 332 if (fast_erf) then 333 do I=1,NRHO 334 xerf=RHO(I)/rlocal 335 yerf = (1.0d0 336 > + xerf*(c1 + xerf*(c2 337 > + xerf*(c3 + xerf*(c4 338 > + xerf*(c5 + xerf*c6))))))**4 339 yerf = (1.0d0 - 1.0d0/yerf**4) 340 F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I) 341 end do 342 else 343 do I=1,NRHO 344 xerf=RHO(I)/rlocal 345 yerf = util_erf(xerf) 346 F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I) 347 end do 348 end if 349 vl(1,1,1)=FORPI*SIMP(NRHO,F,DRHO) 350 end if 351 352* **** semicore density **** 353 if (semicore) then 354 do i=1,nrho 355 f(i) = dsqrt(rho_sc_r(i,1))*rho(i)**2 356 end do 357 rho_sc_k(1,1,1,1) = forpi*SIMP(nrho,f,drho) 358 rho_sc_k(1,1,1,2) = 0.0d0 359 rho_sc_k(1,1,1,3) = 0.0d0 360 rho_sc_k(1,1,1,4) = 0.0d0 361 end if 362 363 do l=1,lmmax 364 vnl(1,1,1,l)=0.0d0 365 end do 366* *** only j0 is non-zero at zero **** 367 if (locp.ne.0) then 368 DO I=1,NRHO 369 F(I)=RHO(I)*WP(I,0)*VP(I,0) 370 END DO 371 vnl(1,1,1,1)=P0*SIMP(NRHO,F,DRHO) 372 end if 373 374 375* ******************************** 376* **** define n_prj and l_prj **** 377* ******************************** 378 lcount = lmmax+1 379 GO TO (950,940,930,920), lmax+1 380 381 !:::::: f-wave ::::::: 382 920 CONTINUE 383 if (locp.ne.3) then 384 lcount = lcount-1 385 n_prj(lcount) = 1 386 l_prj(lcount) = 3 387 m_prj(lcount) = -3 388 b_prj(lcount) = 4 389 390 lcount = lcount-1 391 n_prj(lcount) = 1 392 l_prj(lcount) = 3 393 m_prj(lcount) = -2 394 b_prj(lcount) = 4 395 396 lcount = lcount-1 397 n_prj(lcount) = 1 398 l_prj(lcount) = 3 399 m_prj(lcount) = -1 400 b_prj(lcount) = 4 401 402 lcount = lcount-1 403 n_prj(lcount) = 1 404 l_prj(lcount) = 3 405 m_prj(lcount) = 0 406 b_prj(lcount) = 4 407 408 lcount = lcount-1 409 n_prj(lcount) = 1 410 l_prj(lcount) = 3 411 m_prj(lcount) = 1 412 b_prj(lcount) = 4 413 414 lcount = lcount-1 415 n_prj(lcount) = 1 416 l_prj(lcount) = 3 417 m_prj(lcount) = 2 418 b_prj(lcount) = 4 419 420 lcount = lcount-1 421 n_prj(lcount) = 1 422 l_prj(lcount) = 3 423 m_prj(lcount) = 3 424 b_prj(lcount) = 4 425 end if 426 427 428 !:::: d-wave :::: 429 930 CONTINUE 430 if (locp.ne.2) then 431 lcount = lcount-1 432 n_prj(lcount) = 1 433 l_prj(lcount) = 2 434 m_prj(lcount) = -2 435 b_prj(lcount) = 3 436 437 lcount = lcount-1 438 n_prj(lcount) = 1 439 l_prj(lcount) = 2 440 m_prj(lcount) = -1 441 b_prj(lcount) = 3 442 443 lcount = lcount-1 444 n_prj(lcount) = 1 445 l_prj(lcount) = 2 446 m_prj(lcount) = 0 447 b_prj(lcount) = 3 448 449 lcount = lcount-1 450 n_prj(lcount) = 1 451 l_prj(lcount) = 2 452 m_prj(lcount) = 1 453 b_prj(lcount) = 3 454 455 lcount = lcount-1 456 n_prj(lcount) = 1 457 l_prj(lcount) = 2 458 m_prj(lcount) = 2 459 b_prj(lcount) = 3 460 end if 461 462 463 !:::: p-wave :::: 464 940 CONTINUE 465 if (locp.ne.1) then 466 lcount = lcount-1 467 n_prj(lcount) = 1 468 l_prj(lcount) = 1 469 m_prj(lcount) = -1 470 b_prj(lcount) = 2 471 472 lcount = lcount-1 473 n_prj(lcount) = 1 474 l_prj(lcount) = 1 475 m_prj(lcount) = 0 476 b_prj(lcount) = 2 477 478 lcount = lcount-1 479 n_prj(lcount) = 1 480 l_prj(lcount) = 1 481 m_prj(lcount) = 1 482 b_prj(lcount) = 2 483 end if 484 485 486 !:::: s-wave :::: 487 950 CONTINUE 488 if (locp.ne.0) then 489 lcount = lcount-1 490 n_prj(lcount) = 1 491 l_prj(lcount) = 0 492 m_prj(lcount) = 0 493 b_prj(lcount) = 1 494 end if 495 496 497 IERR=0 498 RETURN 499 END 500 501 502 503