1* 2* $Id$ 3* 4 5* ********************************** 6* * * 7* * integrate_d_stress * 8* * * 9* ********************************** 10 11 subroutine integrate_d_stress(version,rlocal, 12 > nrho,drho,lmax,locp,zv, 13 > vp,wp,rho,f,cs,sn, 14 > nfft3d,lmmax, 15 > G,dvl,dvnl, 16 > semicore,rho_sc_r,rho_sc_k, 17 > ierr) 18 implicit none 19 integer version 20 double precision rlocal 21 integer nrho 22 double precision drho 23 integer lmax 24 integer locp 25 double precision zv 26 double precision vp(nrho,0:lmax) 27 double precision wp(nrho,0:lmax) 28 double precision rho(nrho) 29 double precision f(nrho) 30 double precision cs(nrho) 31 double precision sn(nrho) 32 33 integer nfft3d,lmmax 34 double precision G(nfft3d,3) 35 double precision dvl(nfft3d) 36 double precision dvnl(nfft3d,3,lmmax) 37 38 logical semicore 39 double precision rho_sc_r(nrho,2) 40 double precision rho_sc_k(nfft3d,4) 41 42 integer ierr 43 44#include "errquit.fh" 45 46 integer np,taskid,MASTER 47 integer np_i,np_j,taskid_i,taskid_j,countj 48 parameter (MASTER=0) 49 50* *** local variables **** 51 integer lcount 52 integer k1,k2,k3,i,l,pzero,zero 53 double precision pi,twopi,forpi 54 double precision p0,p1,p2,p3,p,pp 55 double precision gx,gy,gz,a,q,d,dd 56 double precision duxdGx,duxdGy,duxdGz 57 double precision duydGx,duydGy,duydGz 58 double precision duzdGx,duzdGy,duzdGz 59 double precision sumx,sumy,sumz 60 double precision T,dTdux,dTduy,dTduz 61 62* **** external functions **** 63 double precision dsum,simp 64 external dsum,simp 65 double precision tollz 66 parameter(tollz=1d-16) 67 68 if (version.ne.3) then 69 call errquit('integrate_stress - unit cell is aperiodic',0, 70 & INPUT_ERR) 71 end if 72 call Parallel2d_np_i(np_i) 73 call Parallel2d_np_j(np_j) 74 call Parallel2d_taskid_i(taskid_i) 75 call Parallel2d_taskid_j(taskid_j) 76 77 pi=4.0d0*datan(1.0d0) 78 twopi=2.0d0*pi 79 forpi=4.0d0*pi 80 81 IF(LMMAX.GT.16) THEN 82 IERR=1 83 RETURN 84 ENDIF 85 IF((nrho/2)*2.Eq.nrho) THEN 86 IERR=2 87 RETURN 88 ENDIF 89 90 P0=dsqrt(forpi) 91 P1=dsqrt(3.0d0*forpi) 92 P2=dsqrt(15.0d0*forpi) 93 P3=dsqrt(105.0d0*forpi) 94 95*:::::::::::::::::: Define non-local pseudopotential :::::::::::::::: 96 do l=0,lmax 97 if (l.ne.locp) then 98 do I=1,nrho 99 vp(i,l)=vp(i,l)-vp(i,locp) 100 end do 101 end if 102 end do 103 104*====================== Fourier transformation ====================== 105 call dcopy(nfft3d,0.0d0,0,dvl,1) 106 call dcopy(3*lmmax*nfft3d,0.0d0,0,dvnl,1) 107 call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1) 108 109* ***** find the G==0 point in the lattice ***** 110 call D3dB_ijktoindexp(1,1,1,1,zero,pzero) 111 112 countj = -1 113 DO 700 k1=1,nfft3d 114 115 countj = mod(countj+1,np_j) 116 if (countj.ne.taskid_j) go to 700 117 if ((pzero.eq.taskid_i).and.(k1.eq.zero)) go to 700 118 119 q=DSqRT(G(k1,1)**2 120 > +G(k1,2)**2 121 > +G(k1,3)**2) 122 if(abs(q).lt.tollz) goto 700 123 124 gx=G(k1,1)/q 125 gy=G(k1,2)/q 126 gz=G(k1,3)/q 127 DO i=1,nrho 128 cs(i)=DCOS(q*rho(i)) 129 sn(i)=DSIN(q*rho(i)) 130 END DO 131 132* **** calculate du_r/dG_s **** 133 duxdGx = 1.0d0/q -gx*gx/q 134 duxdGy = -gx*gy/q 135 duxdGz = -gx*gz/q 136 137 duydGx = -gy*gx/q 138 duydGy = 1.0d0/q - gy*gy/q 139 duydGz = -gy*gz/q 140 141 duzdGx = -gz*gx/q 142 duzdGy = -gz*gy/q 143 duzdGz = 1.0d0/q - gz*gz/q 144 145 lcount = lmmax+1 146 GO TO (500,400,300,200), LMAX+1 147 148 149*:::::::::::::::::::::::::::::: f-wave :::::::::::::::::::::::::::::: 150 200 CONTINUE 151 if (locp.ne.3) then 152 F(1)=0.0d0 153 do i=2,nrho 154 A=sn(i)/(q*rho(i)) 155 A=15.0d0*(A-cs(i))/(q*rho(i))**2 - 6*A + cs(i) 156 f(i)=A*wp(i,3)*vp(i,3) 157 end do 158 D=P3*SIMP(nrho,F,drho)/q 159 160 F(1)=0.0d0 161 do i=2,nrho 162 A= -60.0d0*sn(i)/(rho(i)**3 * q**5) 163 > + 60.0d0*cs(i)/(rho(i)**2 * q**4) 164 > + 27.0d0*sn(i)/(rho(i) * q**3) 165 > - 7.0d0*cs(i)/(q**2) 166 > - rho(i)*sn(i)/q 167 f(i)=A*wp(i,3)*vp(i,3) 168 end do 169 DD=P3*SIMP(nrho,F,drho) 170 171 lcount = lcount-1 172 T = gy*(3.0d0*(1.0d0-gz*gz)-4.0d0*gy*gy)/dsqrt(24.0d0) 173 dTdux = 0.0d0 174 dTduy = (3.0d0*(1.0d0-gz*gz)-12.0d0*gy*gy)/dsqrt(24.0d0) 175 dTduz = -6.0d0*gy*gz/dsqrt(24.0d0) 176 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 177 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 178 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 179 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 180 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 181 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 182 183 lcount = lcount-1 184 T =gx*gy*gz 185 dTdux = gy*gz 186 dTduy = gx*gz 187 dTduz = gx*gy 188 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 189 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 190 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 191 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 192 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 193 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 194 195 lcount = lcount-1 196 T = gy*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 197 dTdux = 0.0d0 198 dTduy =(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 199 dTduz =10.0d0*gy*gz/dsqrt(40.0d0) 200 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 201 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 202 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 203 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 204 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 205 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 206 207 lcount = lcount-1 208 T =gz*(5.0d0*gz*gz-3.0d0)/dsqrt(60.0d0) 209 dTdux = 0.0d0 210 dTduy = 0.0d0 211 dTduz =(15.0d0*gz*gz -3.0d0)/dsqrt(60.0d0) 212 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 213 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 214 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 215 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 216 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 217 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 218 219 lcount = lcount-1 220 T = gx*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 221 dTdux = (5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 222 dTduy = 0.0d0 223 dTduz = 10.0d0*gx*gz/dsqrt(40.0d0) 224 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 225 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 226 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 227 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 228 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 229 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 230 231 lcount = lcount-1 232 T =gz*(gx*gx - gy*gy)/2.0d0 233 dTdux = gx*gz 234 dTduy = -gy*gz 235 dTduz = (gx*gx-gy*gy)/2.0d0 236 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 237 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 238 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 239 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 240 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 241 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 242 243 lcount = lcount-1 244 T = gx*(4.0d0*gx*gx - 3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0) 245 dTdux = (12.0d0*gx*gx-3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0) 246 dTduy = 0.0d0 247 dTduz = 6.0d0*gx*gz/dsqrt(24.0d0) 248 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 249 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 250 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 251 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 252 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 253 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 254 end if 255 256 257 258*:::::::::::::::::::::::::::::: d-wave :::::::::::::::::::::::::::::: 259 300 CONTINUE 260 if (locp.ne.2) then 261 F(1)=0.0d0 262 DO i=2,nrho 263 A=3.0d0*(sn(i)/(q*rho(i))-cs(i))/(q*rho(i))-sn(i) 264 f(i)=A*wp(i,2)*vp(i,2) 265 END DO 266 D=P2*SIMP(nrho,F,drho)/q 267 268 F(1)=0.0d0 269 DO i=2,nrho 270 A= -9.0d0*sn(i)/(rho(i)**2 * q**4) 271 > + 9.0d0*cs(i)/(rho(i) * q**3) 272 > + 4.0d0*sn(i)/(q**2) 273 > - rho(i)*cs(i)/q 274 f(i)=A*wp(i,2)*vp(i,2) 275 END DO 276 DD=P2*SIMP(nrho,F,drho) 277 278 lcount = lcount-1 279 T = gx*gy 280 dTdux = gy 281 dTduy = gx 282 dTduz = 0.0d0 283 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 284 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 285 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 286 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 287 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 288 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 289 290 lcount = lcount-1 291 T = gy*gz 292 dTdux = 0.0d0 293 dTduy = gz 294 dTduz = gy 295 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 296 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 297 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 298 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 299 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 300 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 301 302 lcount = lcount-1 303 T = (3.0d0*gz*gz-1.0d0)/(2.0d0*dsqrt(3.0d0)) 304 dTdux = 0.0d0 305 dTduy = 0.0d0 306 dTduz = 6.0d0*gz/(2.0d0*dsqrt(3.0d0)) 307 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 308 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 309 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 310 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 311 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 312 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 313 314 lcount = lcount-1 315 T = gz*gx 316 dTdux = gz 317 dTduy = 0.0d0 318 dTduz = gx 319 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 320 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 321 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 322 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 323 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 324 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 325 326 lcount = lcount-1 327 T = (gx*gx-gy*gy)/2.0d0 328 dTdux = gx 329 dTduy = -gy 330 dTduz = 0.0d0 331 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 332 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 333 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 334 dvnl(k1,1,lcount)=DD*T*gx + D*sumx 335 dvnl(k1,2,lcount)=DD*T*gy + D*sumy 336 dvnl(k1,3,lcount)=DD*T*gz + D*sumz 337 end if 338 339*:::::::::::::::::::::::::::::: p-wave :::::::::::::::::::::::::::::: 340 400 CONTINUE 341 if (locp.ne.1) then 342 F(1)=0.0d0 343 DO i=2,nrho 344 f(i)=(sn(i)/(q*rho(i)) - cs(i)) * wp(i,1)*vp(i,1) 345 END DO 346 P=P1*SIMP(nrho,F,drho)/q 347 348 F(1)=0.0d0 349 DO i=2,nrho 350 f(i)=wp(i,1)*vp(i,1)* ( -2.0d0*sn(i)/(rho(i) * q**3) 351 > + 2.0d0*cs(i)/(q**2) 352 > + rho(i)*sn(i)/q) 353 END DO 354 PP=P1*SIMP(nrho,F,drho) 355 356 lcount = lcount-1 357 T = gy 358 dTdux = 0.0d0 359 dTduy = 1.0d0 360 dTduz = 0.0d0 361 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 362 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 363 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 364 dvnl(k1,1,lcount)= PP*T*gx + P*sumx 365 dvnl(k1,2,lcount)= PP*T*gy + P*sumy 366 dvnl(k1,3,lcount)= PP*T*gz + P*sumz 367 368 lcount = lcount-1 369 T = gz 370 dTdux = 0.0d0 371 dTduy = 0.0d0 372 dTduz = 1.0d0 373 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 374 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 375 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 376 dvnl(k1,1,lcount)= PP*T*gx + P*sumx 377 dvnl(k1,2,lcount)= PP*T*gy + P*sumy 378 dvnl(k1,3,lcount)= PP*T*gz + P*sumz 379 380 lcount = lcount-1 381 T = gx 382 dTdux = 1.0d0 383 dTduy = 0.0d0 384 dTduz = 0.0d0 385 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 386 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 387 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 388 dvnl(k1,1,lcount)= PP*T*gx + P*sumx 389 dvnl(k1,2,lcount)= PP*T*gy + P*sumy 390 dvnl(k1,3,lcount)= PP*T*gz + P*sumz 391 392 end if 393 394*:::::::::::::::::::::::::::::: s-wave ::::::::::::::::::::::::::::::: 395 500 CONTINUE 396 if (locp.ne.0) then 397 DO i=1,nrho 398 f(i)=wp(i,0)*vp(i,0) * ( -sn(i)/(q**2) 399 > + rho(i)*cs(i)/q) 400 END DO 401 P = P0*SIMP(nrho,F,drho) 402 lcount = lcount-1 403 dvnl(k1,1,lcount) = P *gx 404 dvnl(k1,2,lcount) = P *gy 405 dvnl(k1,3,lcount) = P *gz 406 end if 407 408*:::::::::::::::::::::::::::::: local ::::::::::::::::::::::::::::::: 409 600 CONTINUE 410 411 do i=1,nrho 412 f(i)=rho(i)*vp(i,locp)*(rho(i)*cs(i)-sn(i)/q) 413 end do 414 dvl(k1)= SIMP(nrho,f,drho)*forpi/q 415 > + zv*forpi/(q*q)*(2.0d0*cs(nrho)/q + rho(nrho)*sn(nrho)) 416 417 418 419*::::::::::::::::::::: semicore density ::::::::::::::::::::::::::::::: 420 if (semicore) then 421 422 do i=1,nrho 423 f(i)=rho(i)*dsqrt(rho_sc_r(i,1))*(rho(i)*cs(i)-sn(i)/q) 424 end do 425 rho_sc_k(k1,1)= SIMP(nrho,f,drho)*forpi/q 426 end if 427 428 700 CONTINUE 429 call D1dB_Vector_SumAll(nfft3d,dvl) 430 call D1dB_Vector_Sumall(3*lmmax*nfft3d,dvnl) 431 call D1dB_Vector_SumAll(nfft3d,rho_sc_k) 432 433 434 435*::::::::::::::::::::::::::::::: G=0 :::::::::::::::::::::::::::::::: 436 if (pzero.eq.taskid_i) then 437 dvl(zero)= 0.0d0 438 do l=1,lmmax 439 dvnl(zero,1,l)=0.0d0 440 dvnl(zero,2,l)=0.0d0 441 dvnl(zero,3,l)=0.0d0 442 end do 443 end if 444 445 446 IERR=0 447 RETURN 448 END 449 450 451 452