1* 2* $Id$ 3* 4 5* ********************************************** 6* * * 7* * integrate_kbpp_band_stress_nonlocal_new * 8* * * 9* ********************************************** 10 11 subroutine integrate_kbpp_band_stress_nonlocal_new(version,kvec, 12 > nrho,drho,lmax,locp,zv, 13 > vp,wp,rho,f,cs,sn, 14 > nfft1,nfft2,nfft3,lmmax, 15 > G,dvnl, 16 > nray,G_ray,dvnl_ray, 17 > ierr) 18 implicit none 19 integer version 20 double precision kvec(3) 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 nfft1,nfft2,nfft3,lmmax 34 double precision G(nfft1,nfft2,nfft3,3) 35 double precision dvnl(nfft1,nfft2,nfft3,3,lmmax) 36 37 integer nray 38 double precision G_ray(nray) 39 double precision dvnl_ray(nray,2,0:lmax,2) 40 41 integer ierr 42 43#include "errquit.fh" 44 45* *** local variables **** 46 integer np,taskid,MASTER 47 parameter (MASTER=0) 48 49 integer lcount,task_count,nfft3d 50 integer k1,k2,k3,i,l,nx 51 double precision p,pp 52 double precision gx,gy,gz,q,d,dd 53 double precision duxdGx,duxdGy,duxdGz 54 double precision duydGx,duydGy,duydGz 55 double precision duzdGx,duzdGy,duzdGz 56 double precision sumx,sumy,sumz,dG 57 double precision T,dTdux,dTduy,dTduz 58 59 60 61* **** external functions **** 62 double precision dsum,simp,nwpw_splint 63 external dsum,simp,nwpw_splint 64 65 call Parallel_np(np) 66 call Parallel_taskid(taskid) 67 68 if(lmmax.gt.16) then 69 call errquit( 70 > 'integrate_kbpp_band_stress_nonlocal_new - lmax > f',0, 71 > INPUT_ERR) 72 end if 73 if((nrho/2)*2.eq.nrho) then 74 call errquit( 75 > 'integrate_kbpp_band_stress_nonlocal_new - psp grid not odd',0, 76 > INPUT_ERR) 77 end if 78 79 nfft3d = (nfft1)*nfft2*nfft3 80 81*====================== Fourier transformation ====================== 82 dG = G_ray(3)-G_ray(2) 83 call dcopy(3*lmmax*nfft3d,0.0d0,0,dvnl,1) 84 task_count = -1 85 do 700 k3=1,nfft3 86 do 700 k2=1,nfft2 87 do 700 k1=1,nfft1 88 task_count = task_count + 1 89 if (mod(task_count,np).ne.taskid) go to 700 90 gx=G(k1,k2,k3,1)+kvec(1) 91 gy=G(k1,k2,k3,2)+kvec(2) 92 gz=G(k1,k2,k3,3)+kvec(3) 93 94 q=dsqrt(gx**2 + gy**2 + gz**2) 95 nx = (q/dG) + 1.0d0 96 97 if (dabs(q).gt.1.0d-9) then 98 99 gx=gx/q 100 gy=gy/q 101 gz=gz/q 102 do i=1,nrho 103 cs(i)=dcos(q*rho(i)) 104 sn(i)=dsin(q*rho(i)) 105 end do 106 107* **** calculate du_r/dG_s **** 108 duxdGx = 1.0d0/q -gx*gx/q 109 duxdGy = -gx*gy/q 110 duxdGz = -gx*gz/q 111 112 duydGx = -gy*gx/q 113 duydGy = 1.0d0/q - gy*gy/q 114 duydGz = -gy*gz/q 115 116 duzdGx = -gz*gx/q 117 duzdGy = -gz*gy/q 118 duzdGz = 1.0d0/q - gz*gz/q 119 120 lcount = lmmax+1 121 GO TO (500,400,300,200), LMAX+1 122 123*:::::::::::::::::::::::::::::: f-wave :::::::::::::::::::::::::::::: 124 200 CONTINUE 125 if (locp.ne.3) then 126 D = nwpw_splint(G_ray,dvnl_ray(1,1,3,1), 127 > dvnl_ray(1,1,3,2),nray,nx,Q) 128 DD = nwpw_splint(G_ray,dvnl_ray(1,2,3,1), 129 > dvnl_ray(1,2,3,2),nray,nx,Q) 130 lcount = lcount-1 131 T = gx*(4.0d0*gx*gx - 3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0) 132 dTdux = (12.0d0*gx*gx-3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0) 133 dTduy = 0.0d0 134 dTduz = 6.0d0*gx*gz/dsqrt(24.0d0) 135 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 136 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 137 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 138 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 139 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 140 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 141 142 lcount = lcount-1 143 T = gy*(3.0d0*(1.0d0-gz*gz)-4.0d0*gy*gy)/dsqrt(24.0d0) 144 dTdux = 0.0d0 145 dTduy = (3.0d0*(1.0d0-gz*gz)-12.0d0*gy*gy)/dsqrt(24.0d0) 146 dTduz = -6.0d0*gy*gz/dsqrt(24.0d0) 147 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 148 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 149 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 150 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 151 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 152 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 153 154 lcount = lcount-1 155 T =gz*(gx*gx - gy*gy)/2.0d0 156 dTdux = gx*gz 157 dTduy = -gy*gz 158 dTduz = (gx*gx-gy*gy)/2.0d0 159 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 160 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 161 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 162 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 163 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 164 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 165 166 lcount = lcount-1 167 T =gx*gy*gz 168 dTdux = gy*gz 169 dTduy = gx*gz 170 dTduz = gx*gy 171 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 172 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 173 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 174 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 175 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 176 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 177 178 lcount = lcount-1 179 T = gx*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 180 dTdux = (5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 181 dTduy = 0.0d0 182 dTduz = 10.0d0*gx*gz/dsqrt(40.0d0) 183 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 184 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 185 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 186 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 187 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 188 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 189 190 lcount = lcount-1 191 T = gy*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 192 dTdux = 0.0d0 193 dTduy =(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0) 194 dTduz =10.0d0*gy*gz/dsqrt(40.0d0) 195 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 196 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 197 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 198 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 199 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 200 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 201 202 lcount = lcount-1 203 T =gz*(5.0d0*gz*gz-3.0d0)/dsqrt(60.0d0) 204 dTdux = 0.0d0 205 dTduy = 0.0d0 206 dTduz =(15.0d0*gz*gz -3.0d0)/dsqrt(60.0d0) 207 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 208 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 209 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 210 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 211 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 212 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 213 end if 214*:::::::::::::::::::::::::::::: d-wave :::::::::::::::::::::::::::::: 215 300 CONTINUE 216 if (locp.ne.2) then 217 D = nwpw_splint(G_ray,dvnl_ray(1,1,2,1), 218 > dvnl_ray(1,1,2,2),nray,nx,Q) 219 DD = nwpw_splint(G_ray,dvnl_ray(1,2,2,1), 220 > dvnl_ray(1,2,2,2),nray,nx,Q) 221 lcount = lcount-1 222 T = (3.0d0*gz*gz-1.0d0)/(2.0d0*dsqrt(3.0d0)) 223 dTdux = 0.0d0 224 dTduy = 0.0d0 225 dTduz = 6.0d0*gz/(2.0d0*dsqrt(3.0d0)) 226 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 227 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 228 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 229 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 230 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 231 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 232 233 lcount = lcount-1 234 T = gx*gy 235 dTdux = gy 236 dTduy = gx 237 dTduz = 0.0d0 238 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 239 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 240 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 241 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 242 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 243 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 244 245 lcount = lcount-1 246 T = gy*gz 247 dTdux = 0.0d0 248 dTduy = gz 249 dTduz = gy 250 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 251 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 252 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 253 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 254 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 255 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 256 257 lcount = lcount-1 258 T = gz*gx 259 dTdux = gz 260 dTduy = 0.0d0 261 dTduz = gx 262 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 263 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 264 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 265 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 266 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 267 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 268 269 lcount = lcount-1 270 T = (gx*gx-gy*gy)/2.0d0 271 dTdux = gx 272 dTduy = -gy 273 dTduz = 0.0d0 274 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 275 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 276 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 277 dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx 278 dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy 279 dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz 280 end if 281*:::::::::::::::::::::::::::::: p-wave :::::::::::::::::::::::::::::: 282 400 CONTINUE 283 if (locp.ne.1) then 284 P = nwpw_splint(G_ray,dvnl_ray(1,1,1,1), 285 > dvnl_ray(1,1,1,2),nray,nx,Q) 286 PP = nwpw_splint(G_ray,dvnl_ray(1,2,1,1), 287 > dvnl_ray(1,2,1,2),nray,nx,Q) 288 lcount = lcount-1 289 T = gx 290 dTdux = 1.0d0 291 dTduy = 0.0d0 292 dTduz = 0.0d0 293 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 294 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 295 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 296 dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx 297 dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy 298 dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz 299 300 301 lcount = lcount-1 302 T = gy 303 dTdux = 0.0d0 304 dTduy = 1.0d0 305 dTduz = 0.0d0 306 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 307 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 308 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 309 dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx 310 dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy 311 dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz 312 313 lcount = lcount-1 314 T = gz 315 dTdux = 0.0d0 316 dTduy = 0.0d0 317 dTduz = 1.0d0 318 sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx 319 sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy 320 sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz 321 dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx 322 dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy 323 dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz 324 end if 325*:::::::::::::::::::::::::::::: s-wave ::::::::::::::::::::::::::::::: 326 500 CONTINUE 327 if (locp.ne.0) then 328 P = nwpw_splint(G_ray,dvnl_ray(1,1,0,1), 329 > dvnl_ray(1,1,0,2),nray,nx,Q) 330 lcount = lcount-1 331 dvnl(k1,k2,k3,1,lcount) = P *gx 332 dvnl(k1,k2,k3,2,lcount) = P *gy 333 dvnl(k1,k2,k3,3,lcount) = P *gz 334 end if 335 336 337 600 CONTINUE 338*::::::::::::::::::::::::::::::: G+k=0 :::::::::::::::::::::::::::::::: 339 else 340 341 do l=1,lmmax 342 dvnl(1,1,1,1,l)=0.0d0 343 dvnl(1,1,1,2,l)=0.0d0 344 dvnl(1,1,1,3,l)=0.0d0 345 end do 346 347 end if 348 349 700 CONTINUE 350 call Parallel_Vector_Sumall(3*lmmax*nfft3d,dvnl) 351 352 ierr=0 353 RETURN 354 END 355 356 357 358