1* 2* $Id$ 3* 4 5 subroutine lattice_min_difference(x,y,z) 6 implicit none 7 real*8 x,y,z 8 9* **** common block **** 10 real*8 ecut,wcut,omega 11 real*8 ua(3,3),unitg(3,3) 12 common / lattice_block / ua,unitg,ecut,wcut,omega 13 14 real*8 ub(3,3) 15 common / lattice_block2 / ub 16 17* *** local variables **** 18 real*8 c1,c2,c3 19 20 c1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1) 21 c2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2) 22 c3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3) 23 c1 = c1 - DNINT(c1) 24 c2 = c2 - DNINT(c2) 25 c3 = c3 - DNINT(c3) 26 x = ua(1,1)*c1 + ua(1,2)*c2 + ua(1,3)*c3 27 y = ua(2,1)*c1 + ua(2,2)*c2 + ua(2,3)*c3 28 z = ua(3,1)*c1 + ua(3,2)*c2 + ua(3,3)*c3 29 30 return 31 end 32 33 subroutine lattice_frac_to_r1(n,f1,r1) 34 implicit none 35 integer n 36 real*8 f1(3,*),r1(3,*) 37 38* **** common block **** 39 real*8 ecut,wcut,omega 40 real*8 ua(3,3),unitg(3,3) 41 common / lattice_block / ua,unitg,ecut,wcut,omega 42 43* **** local variables *** 44 integer i 45 46 do i=1,n 47 r1(1,i) = f1(1,i)*ua(1,1) + f1(2,i)*ua(2,1) + f1(3,i)*ua(3,1) 48 r1(2,i) = f1(1,i)*ua(1,2) + f1(2,i)*ua(2,2) + f1(3,i)*ua(3,2) 49 r1(3,i) = f1(1,i)*ua(1,3) + f1(2,i)*ua(2,3) + f1(3,i)*ua(3,3) 50 end do 51 return 52 end 53 54 subroutine lattice_r1_to_frac(n,r1,f1) 55 implicit none 56 integer n 57 real*8 r1(3,*),f1(3,*) 58 59* **** common block **** 60 real*8 ub(3,3) 61 common / lattice_block2 / ub 62 63* **** local variables *** 64 integer i 65 66 do i=1,n 67 f1(1,i) = r1(1,i)*ub(1,1) + r1(2,i)*ub(2,1) + r1(3,i)*ub(3,1) 68 f1(2,i) = r1(1,i)*ub(1,2) + r1(2,i)*ub(2,2) + r1(3,i)*ub(3,2) 69 f1(3,i) = r1(1,i)*ub(1,3) + r1(2,i)*ub(2,3) + r1(3,i)*ub(3,3) 70 end do 71 return 72 end 73 74 75 subroutine lattice_center_xyz_to_ijk(x,y,z,c1,c2,c3) 76 implicit none 77 real*8 x,y,z 78 integer c1,c2,c3 79 80 integer np1,np2,np3 81 real*8 f1,f2,f3 82 83* **** common block **** 84 real*8 ub(3,3) 85 common / lattice_block2 / ub 86 87 call D3dB_nx(1,np1) 88 call D3dB_ny(1,np2) 89 call D3dB_nz(1,np3) 90 f1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1) 91 f2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2) 92 f3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3) 93 c1 = (f1-0.5d0+0.5d0/dble(np1))*np1 94 c2 = (f2-0.5d0+0.5d0/dble(np2))*np2 95 c3 = (f3-0.5d0+0.5d0/dble(np3))*np3 96 return 97 end 98 99 subroutine lattice_center0_xyz_to_ijk(x,y,z,c1,c2,c3) 100 implicit none 101 real*8 x,y,z 102 integer c1,c2,c3 103 104 integer np1,np2,np3 105 real*8 f1,f2,f3 106 107* **** common block **** 108 real*8 ub(3,3) 109 common / lattice_block2 / ub 110 111 call D3dB_nx(1,np1) 112 call D3dB_ny(1,np2) 113 call D3dB_nz(1,np3) 114 f1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1) 115 f2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2) 116 f3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3) 117 c1 = (f1+0.5d0/dble(np1))*np1 118 c2 = (f2+0.5d0/dble(np2))*np2 119 c3 = (f3+0.5d0/dble(np3))*np3 120 return 121 end 122 123 124 subroutine lattice_fragcell(n1,r1) 125 implicit none 126 real*8 rcm(3) 127 integer n1 128 real*8 r1(3,*) 129 130* **** common block **** 131 real*8 ecut,wcut,omega 132 real*8 ua(3,3),unitg(3,3) 133 common / lattice_block / ua,unitg,ecut,wcut,omega 134 135* **** common block **** 136 real*8 ub(3,3) 137 common / lattice_block2 / ub 138 139* **** local variables **** 140 integer i,j 141 real*8 fcm 142 143 do i=2,n1 144 do j=1,3 145 fcm = (r1(1,i)-r1(1,1))*ub(1,j) 146 > + (r1(2,i)-r1(2,1))*ub(2,j) 147 > + (r1(3,i)-r1(3,1))*ub(3,j) 148 do while (fcm.gt.0.5) 149 r1(1,i) = r1(1,i) - ua(1,j) 150 r1(2,i) = r1(2,i) - ua(2,j) 151 r1(3,i) = r1(3,i) - ua(3,j) 152 fcm = (r1(1,i)-r1(1,1))*ub(1,j) 153 > + (r1(2,i)-r1(2,1))*ub(2,j) 154 > + (r1(3,i)-r1(3,1))*ub(3,j) 155 end do 156 do while (fcm.le.(-0.5)) 157 r1(1,i) = r1(1,i) + ua(1,j) 158 r1(2,i) = r1(2,i) + ua(2,j) 159 r1(3,i) = r1(3,i) + ua(3,j) 160 fcm = (r1(1,i)-r1(1,1))*ub(1,j) 161 > + (r1(2,i)-r1(2,1))*ub(2,j) 162 > + (r1(3,i)-r1(3,1))*ub(3,j) 163 end do 164 end do 165 end do 166 return 167 end 168 169 170 subroutine lattice_incell1_frag(rcm,n1,r1) 171 implicit none 172 real*8 rcm(3) 173 integer n1 174 real*8 r1(3,*) 175 176* **** common block **** 177 real*8 ecut,wcut,omega 178 real*8 ua(3,3),unitg(3,3) 179 common / lattice_block / ua,unitg,ecut,wcut,omega 180 181* **** common block **** 182 real*8 ub(3,3) 183 common / lattice_block2 / ub 184 185* **** local variables **** 186 integer i 187 real*8 fcm(3) 188 189 call lattice_r1_to_frac(1,rcm,fcm) 190 191 do while (fcm(1).gt.(0.5d0)) 192 rcm(1) = rcm(1) - ua(1,1) 193 rcm(2) = rcm(2) - ua(2,1) 194 rcm(3) = rcm(3) - ua(3,1) 195 do i=1,n1 196 r1(1,i) = r1(1,i) - ua(1,1) 197 r1(2,i) = r1(2,i) - ua(2,1) 198 r1(3,i) = r1(3,i) - ua(3,1) 199 end do 200 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 201 end do 202 do while (fcm(1).le.(-0.5d0)) 203 rcm(1) = rcm(1) + ua(1,1) 204 rcm(2) = rcm(2) + ua(2,1) 205 rcm(3) = rcm(3) + ua(3,1) 206 do i=1,n1 207 r1(1,i) = r1(1,i) + ua(1,1) 208 r1(2,i) = r1(2,i) + ua(2,1) 209 r1(3,i) = r1(3,i) + ua(3,1) 210 end do 211 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 212 end do 213 214 do while (fcm(2).gt.(0.5d0)) 215 rcm(1) = rcm(1) - ua(1,2) 216 rcm(2) = rcm(2) - ua(2,2) 217 rcm(3) = rcm(3) - ua(3,2) 218 do i=1,n1 219 r1(1,i) = r1(1,i) - ua(1,2) 220 r1(2,i) = r1(2,i) - ua(2,2) 221 r1(3,i) = r1(3,i) - ua(3,2) 222 end do 223 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 224 end do 225 do while (fcm(2).le.(-0.5d0)) 226 rcm(1) = rcm(1) + ua(1,2) 227 rcm(2) = rcm(2) + ua(2,2) 228 rcm(3) = rcm(3) + ua(3,2) 229 do i=1,n1 230 r1(1,i) = r1(1,i) + ua(1,2) 231 r1(2,i) = r1(2,i) + ua(2,2) 232 r1(3,i) = r1(3,i) + ua(3,2) 233 end do 234 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 235 end do 236 237 do while (fcm(3).gt.(0.5d0)) 238 rcm(1) = rcm(1) - ua(1,3) 239 rcm(2) = rcm(2) - ua(2,3) 240 rcm(3) = rcm(3) - ua(3,3) 241 do i=1,n1 242 r1(1,i) = r1(1,i) - ua(1,3) 243 r1(2,i) = r1(2,i) - ua(2,3) 244 r1(3,i) = r1(3,i) - ua(3,3) 245 end do 246 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 247 end do 248 do while (fcm(3).le.(-0.5d0)) 249 rcm(1) = rcm(1) + ua(1,3) 250 rcm(2) = rcm(2) + ua(2,3) 251 rcm(3) = rcm(3) + ua(3,3) 252 do i=1,n1 253 r1(1,i) = r1(1,i) + ua(1,3) 254 r1(2,i) = r1(2,i) + ua(2,3) 255 r1(3,i) = r1(3,i) + ua(3,3) 256 end do 257 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 258 end do 259 return 260 end 261 262 subroutine lattice_incell2_frag(rcm,n1,r1,r2) 263 implicit none 264 real*8 rcm(3) 265 integer n1 266 real*8 r1(3,*) 267 real*8 r2(3,*) 268 269* **** common block **** 270 real*8 ecut,wcut,omega 271 real*8 ua(3,3),unitg(3,3) 272 common / lattice_block / ua,unitg,ecut,wcut,omega 273 274* **** common block **** 275 real*8 ub(3,3) 276 common / lattice_block2 / ub 277 278* **** local variables **** 279 integer i 280 real*8 fcm(3) 281 282 call lattice_r1_to_frac(1,rcm,fcm) 283 284 do while (fcm(1).gt.(0.5d0)) 285 rcm(1) = rcm(1) - ua(1,1) 286 rcm(2) = rcm(2) - ua(2,1) 287 rcm(3) = rcm(3) - ua(3,1) 288 do i=1,n1 289 r1(1,i) = r1(1,i) - ua(1,1) 290 r1(2,i) = r1(2,i) - ua(2,1) 291 r1(3,i) = r1(3,i) - ua(3,1) 292 293 r2(1,i) = r2(1,i) - ua(1,1) 294 r2(2,i) = r2(2,i) - ua(2,1) 295 r2(3,i) = r2(3,i) - ua(3,1) 296 end do 297 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 298 end do 299 do while (fcm(1).le.(-0.5d0)) 300 rcm(1) = rcm(1) + ua(1,1) 301 rcm(2) = rcm(2) + ua(2,1) 302 rcm(3) = rcm(3) + ua(3,1) 303 do i=1,n1 304 r1(1,i) = r1(1,i) + ua(1,1) 305 r1(2,i) = r1(2,i) + ua(2,1) 306 r1(3,i) = r1(3,i) + ua(3,1) 307 308 r2(1,i) = r2(1,i) + ua(1,1) 309 r2(2,i) = r2(2,i) + ua(2,1) 310 r2(3,i) = r2(3,i) + ua(3,1) 311 end do 312 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 313 end do 314 315 do while (fcm(2).gt.(0.5d0)) 316 rcm(1) = rcm(1) - ua(1,2) 317 rcm(2) = rcm(2) - ua(2,2) 318 rcm(3) = rcm(3) - ua(3,2) 319 do i=1,n1 320 r1(1,i) = r1(1,i) - ua(1,2) 321 r1(2,i) = r1(2,i) - ua(2,2) 322 r1(3,i) = r1(3,i) - ua(3,2) 323 324 r2(1,i) = r2(1,i) - ua(1,2) 325 r2(2,i) = r2(2,i) - ua(2,2) 326 r2(3,i) = r2(3,i) - ua(3,2) 327 end do 328 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 329 end do 330 do while (fcm(2).le.(-0.5d0)) 331 rcm(1) = rcm(1) + ua(1,2) 332 rcm(2) = rcm(2) + ua(2,2) 333 rcm(3) = rcm(3) + ua(3,2) 334 do i=1,n1 335 r1(1,i) = r1(1,i) + ua(1,2) 336 r1(2,i) = r1(2,i) + ua(2,2) 337 r1(3,i) = r1(3,i) + ua(3,2) 338 339 r2(1,i) = r2(1,i) + ua(1,2) 340 r2(2,i) = r2(2,i) + ua(2,2) 341 r2(3,i) = r2(3,i) + ua(3,2) 342 end do 343 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 344 end do 345 346 do while (fcm(3).gt.(0.5d0)) 347 rcm(1) = rcm(1) - ua(1,3) 348 rcm(2) = rcm(2) - ua(2,3) 349 rcm(3) = rcm(3) - ua(3,3) 350 do i=1,n1 351 r1(1,i) = r1(1,i) - ua(1,3) 352 r1(2,i) = r1(2,i) - ua(2,3) 353 r1(3,i) = r1(3,i) - ua(3,3) 354 355 r2(1,i) = r2(1,i) - ua(1,3) 356 r2(2,i) = r2(2,i) - ua(2,3) 357 r2(3,i) = r2(3,i) - ua(3,3) 358 end do 359 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 360 end do 361 do while (fcm(3).le.(-0.5d0)) 362 rcm(1) = rcm(1) + ua(1,3) 363 rcm(2) = rcm(2) + ua(2,3) 364 rcm(3) = rcm(3) + ua(3,3) 365 do i=1,n1 366 r1(1,i) = r1(1,i) + ua(1,3) 367 r1(2,i) = r1(2,i) + ua(2,3) 368 r1(3,i) = r1(3,i) + ua(3,3) 369 370 r2(1,i) = r2(1,i) + ua(1,3) 371 r2(2,i) = r2(2,i) + ua(2,3) 372 r2(3,i) = r2(3,i) + ua(3,3) 373 end do 374 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 375 end do 376 return 377 end 378 379 subroutine lattice_incell3_frag(rcm,n1,r1,r2,r3) 380 implicit none 381 real*8 rcm(3) 382 integer n1 383 real*8 r1(3,*) 384 real*8 r2(3,*) 385 real*8 r3(3,*) 386 387* **** common block **** 388 real*8 ecut,wcut,omega 389 real*8 ua(3,3),unitg(3,3) 390 common / lattice_block / ua,unitg,ecut,wcut,omega 391 392* **** common block **** 393 real*8 ub(3,3) 394 common / lattice_block2 / ub 395 396* **** local variables **** 397 integer i 398 real*8 fcm(3) 399 400 call lattice_r1_to_frac(1,rcm,fcm) 401 402 do while (fcm(1).gt.(0.5d0)) 403 rcm(1) = rcm(1) - ua(1,1) 404 rcm(2) = rcm(2) - ua(2,1) 405 rcm(3) = rcm(3) - ua(3,1) 406 do i=1,n1 407 r1(1,i) = r1(1,i) - ua(1,1) 408 r1(2,i) = r1(2,i) - ua(2,1) 409 r1(3,i) = r1(3,i) - ua(3,1) 410 411 r2(1,i) = r2(1,i) - ua(1,1) 412 r2(2,i) = r2(2,i) - ua(2,1) 413 r2(3,i) = r2(3,i) - ua(3,1) 414 415 r3(1,i) = r3(1,i) - ua(1,1) 416 r3(2,i) = r3(2,i) - ua(2,1) 417 r3(3,i) = r3(3,i) - ua(3,1) 418 end do 419 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 420 end do 421 do while (fcm(1).le.(-0.5d0)) 422 rcm(1) = rcm(1) + ua(1,1) 423 rcm(2) = rcm(2) + ua(2,1) 424 rcm(3) = rcm(3) + ua(3,1) 425 do i=1,n1 426 r1(1,i) = r1(1,i) + ua(1,1) 427 r1(2,i) = r1(2,i) + ua(2,1) 428 r1(3,i) = r1(3,i) + ua(3,1) 429 430 r2(1,i) = r2(1,i) + ua(1,1) 431 r2(2,i) = r2(2,i) + ua(2,1) 432 r2(3,i) = r2(3,i) + ua(3,1) 433 434 r3(1,i) = r3(1,i) + ua(1,1) 435 r3(2,i) = r3(2,i) + ua(2,1) 436 r3(3,i) = r3(3,i) + ua(3,1) 437 end do 438 fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1) 439 end do 440 441 do while (fcm(2).gt.(0.5d0)) 442 rcm(1) = rcm(1) - ua(1,2) 443 rcm(2) = rcm(2) - ua(2,2) 444 rcm(3) = rcm(3) - ua(3,2) 445 do i=1,n1 446 r1(1,i) = r1(1,i) - ua(1,2) 447 r1(2,i) = r1(2,i) - ua(2,2) 448 r1(3,i) = r1(3,i) - ua(3,2) 449 450 r2(1,i) = r2(1,i) - ua(1,2) 451 r2(2,i) = r2(2,i) - ua(2,2) 452 r2(3,i) = r2(3,i) - ua(3,2) 453 454 r3(1,i) = r3(1,i) - ua(1,2) 455 r3(2,i) = r3(2,i) - ua(2,2) 456 r3(3,i) = r3(3,i) - ua(3,2) 457 end do 458 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 459 end do 460 do while (fcm(2).le.(-0.5d0)) 461 rcm(1) = rcm(1) + ua(1,2) 462 rcm(2) = rcm(2) + ua(2,2) 463 rcm(3) = rcm(3) + ua(3,2) 464 do i=1,n1 465 r1(1,i) = r1(1,i) + ua(1,2) 466 r1(2,i) = r1(2,i) + ua(2,2) 467 r1(3,i) = r1(3,i) + ua(3,2) 468 469 r2(1,i) = r2(1,i) + ua(1,2) 470 r2(2,i) = r2(2,i) + ua(2,2) 471 r2(3,i) = r2(3,i) + ua(3,2) 472 473 r3(1,i) = r3(1,i) + ua(1,2) 474 r3(2,i) = r3(2,i) + ua(2,2) 475 r3(3,i) = r3(3,i) + ua(3,2) 476 end do 477 fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2) 478 end do 479 480 do while (fcm(3).gt.(0.5d0)) 481 rcm(1) = rcm(1) - ua(1,3) 482 rcm(2) = rcm(2) - ua(2,3) 483 rcm(3) = rcm(3) - ua(3,3) 484 do i=1,n1 485 r1(1,i) = r1(1,i) - ua(1,3) 486 r1(2,i) = r1(2,i) - ua(2,3) 487 r1(3,i) = r1(3,i) - ua(3,3) 488 489 r2(1,i) = r2(1,i) - ua(1,3) 490 r2(2,i) = r2(2,i) - ua(2,3) 491 r2(3,i) = r2(3,i) - ua(3,3) 492 493 r3(1,i) = r3(1,i) - ua(1,3) 494 r3(2,i) = r3(2,i) - ua(2,3) 495 r3(3,i) = r3(3,i) - ua(3,3) 496 end do 497 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 498 end do 499 do while (fcm(3).le.(-0.5d0)) 500 rcm(1) = rcm(1) + ua(1,3) 501 rcm(2) = rcm(2) + ua(2,3) 502 rcm(3) = rcm(3) + ua(3,3) 503 do i=1,n1 504 r1(1,i) = r1(1,i) + ua(1,3) 505 r1(2,i) = r1(2,i) + ua(2,3) 506 r1(3,i) = r1(3,i) + ua(3,3) 507 508 r2(1,i) = r2(1,i) + ua(1,3) 509 r2(2,i) = r2(2,i) + ua(2,3) 510 r2(3,i) = r2(3,i) + ua(3,3) 511 512 r3(1,i) = r3(1,i) + ua(1,3) 513 r3(2,i) = r3(2,i) + ua(2,3) 514 r3(3,i) = r3(3,i) + ua(3,3) 515 end do 516 fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3) 517 end do 518 return 519 end 520 521 522 real*8 function lattice_wcut() 523 implicit none 524 525* **** common block **** 526 real*8 ecut,wcut,omega 527 real*8 unita(3,3),unitg(3,3) 528 common / lattice_block / unita,unitg,ecut,wcut,omega 529 530 lattice_wcut = wcut 531 return 532 end 533 534 real*8 function lattice_ecut() 535 implicit none 536 537* **** common block **** 538 real*8 ecut,wcut,omega 539 real*8 unita(3,3),unitg(3,3) 540 common / lattice_block / unita,unitg,ecut,wcut,omega 541 542 lattice_ecut = ecut 543 return 544 end 545 546 real*8 function lattice_ggcut() 547 implicit none 548 549* **** common block **** 550 real*8 ecut,wcut,omega 551 real*8 unita(3,3),unitg(3,3) 552 common / lattice_block / unita,unitg,ecut,wcut,omega 553 554 lattice_ggcut = 2.0d0*ecut 555 return 556 end 557 558 real*8 function lattice_wggcut() 559 implicit none 560 561* **** common block **** 562 real*8 ecut,wcut,omega 563 real*8 unita(3,3),unitg(3,3) 564 common / lattice_block / unita,unitg,ecut,wcut,omega 565 566 lattice_wggcut = 2.0d0*wcut 567 return 568 end 569 570 571 572 real*8 function lattice_omega() 573 implicit none 574 575* **** common block **** 576 real*8 ecut,wcut,omega 577 real*8 unita(3,3),unitg(3,3) 578 common / lattice_block / unita,unitg,ecut,wcut,omega 579 580 lattice_omega = omega 581 return 582 end 583 584 real*8 function lattice_unita(i,j) 585 implicit none 586 integer i,j 587 588* **** common block **** 589 real*8 ecut,wcut,omega 590 real*8 unita(3,3),unitg(3,3) 591 common / lattice_block / unita,unitg,ecut,wcut,omega 592 593 lattice_unita = unita(i,j) 594 return 595 end 596 597 598 real*8 function lattice_unitg(i,j) 599 implicit none 600 integer i,j 601 602* **** common block **** 603 real*8 ecut,wcut,omega 604 real*8 unita(3,3),unitg(3,3) 605 common / lattice_block / unita,unitg,ecut,wcut,omega 606 607 lattice_unitg = unitg(i,j) 608 return 609 end 610 611* ******************************************** 612* * * 613* * lattice_unita_small * 614* * * 615* ******************************************** 616 real*8 function lattice_unita_small(i,j) 617 implicit none 618 integer i,j 619 620* **** common block **** 621 logical has_small 622 real*8 omega_small 623 real*8 unita_small(3,3),unitg_small(3,3) 624 common /lattice_small_block/ unita_small,unitg_small, 625 > omega_small,has_small 626 627 lattice_unita_small = unita_small(i,j) 628 return 629 end 630 631* ******************************************** 632* * * 633* * lattice_unita_frozen_small * 634* * * 635* ******************************************** 636 real*8 function lattice_unita_frozen_small(i,j) 637 implicit none 638 integer i,j 639 640* **** common block **** 641 real*8 omega_frozen_small 642 real*8 unita_frozen_small(3,3),unitg_frozen_small(3,3) 643 common /lattice_small_frozen_block/ unita_frozen_small, 644 > unitg_frozen_small, 645 > omega_frozen_small 646 647 lattice_unita_frozen_small = unita_frozen_small(i,j) 648 return 649 end 650 651 652* ******************************************** 653* * * 654* * lattice_unitg_small * 655* * * 656* ******************************************** 657 real*8 function lattice_unitg_small(i,j) 658 implicit none 659 integer i,j 660 661* **** common block **** 662 logical has_small 663 real*8 omega_small 664 real*8 unita_small(3,3),unitg_small(3,3) 665 common /lattice_small_block/ unita_small,unitg_small, 666 > omega_small,has_small 667 668 lattice_unitg_small = unitg_small(i,j) 669 return 670 end 671 672* ******************************************** 673* * * 674* * lattice_unitg_frozen_small * 675* * * 676* ******************************************** 677 real*8 function lattice_unitg_frozen_small(i,j) 678 implicit none 679 integer i,j 680 681* **** common block **** 682 real*8 omega_frozen_small 683 real*8 unita_frozen_small(3,3),unitg_frozen_small(3,3) 684 common /lattice_small_frozen_block/ unita_frozen_small, 685 > unitg_frozen_small, 686 > omega_frozen_small 687 688 lattice_unitg_frozen_small = unitg_frozen_small(i,j) 689 return 690 end 691 692 693 694* ******************************************** 695* * * 696* * lattice_omega_small * 697* * * 698* ******************************************** 699 real*8 function lattice_omega_small(i,j) 700 implicit none 701 integer i,j 702 703* **** common block **** 704 logical has_small 705 real*8 omega_small 706 real*8 unita_small(3,3),unitg_small(3,3) 707 common /lattice_small_block/ unita_small,unitg_small, 708 > omega_small,has_small 709 710 lattice_omega_small = omega_small 711 return 712 end 713 714* ******************************************** 715* * * 716* * lattice_omega_frozen_small * 717* * * 718* ******************************************** 719 real*8 function lattice_omega_frozen_small() 720 implicit none 721 integer i,j 722 723* **** common block **** 724 real*8 omega_frozen_small 725 real*8 unita_frozen_small(3,3),unitg_frozen_small(3,3) 726 common /lattice_small_frozen_block/ unita_frozen_small, 727 > unitg_frozen_small, 728 > omega_frozen_small 729 730 lattice_omega_frozen_small = omega_frozen_small 731 return 732 end 733 734* ******************************************** 735* * * 736* * lattice_has_small * 737* * * 738* ******************************************** 739 logical function lattice_has_small() 740 implicit none 741 742* **** common block **** 743 logical has_small 744 real*8 omega_small 745 real*8 unita_small(3,3),unitg_small(3,3) 746 common /lattice_small_block/ unita_small,unitg_small, 747 > omega_small,has_small 748 749 lattice_has_small = has_small 750 return 751 end 752 753 754 755 756 757* ******************************************** 758* * * 759* * lattice_unitg_frozen * 760* * * 761* ******************************************** 762 763* frozen lattice structure - used to determine masking arrays 764 765 real*8 function lattice_unitg_frozen(i,j) 766 implicit none 767 integer i,j 768 769* **** common blocks **** 770 logical frozen 771 real*8 ecut_frozen,wcut_frozen,omega_frozen 772 real*8 unita_frozen(3,3),unitg_frozen(3,3) 773 common / lattice_froze / unita_frozen,unitg_frozen, 774 > ecut_frozen,wcut_frozen, 775 > omega_frozen,frozen 776 777 lattice_unitg_frozen = unitg_frozen(i,j) 778 return 779 end 780 781* ******************************************** 782* * * 783* * lattice_unita_frozen * 784* * * 785* ******************************************** 786 787* frozen lattice structure - used to determine masking arrays 788 789 real*8 function lattice_unita_frozen(i,j) 790 implicit none 791 integer i,j 792 793* **** common blocks **** 794 logical frozen 795 real*8 ecut_frozen,wcut_frozen,omega_frozen 796 real*8 unita_frozen(3,3),unitg_frozen(3,3) 797 common / lattice_froze / unita_frozen,unitg_frozen, 798 > ecut_frozen,wcut_frozen, 799 > omega_frozen,frozen 800 801 lattice_unita_frozen = unita_frozen(i,j) 802 return 803 end 804 805* ******************************************** 806* * * 807* * lattice_omega_frozen * 808* * * 809* ******************************************** 810 811* frozen lattice structure - used to determine masking arrays 812 813 real*8 function lattice_omega_frozen() 814 implicit none 815 816* **** common blocks **** 817 logical frozen 818 real*8 ecut_frozen,wcut_frozen,omega_frozen 819 real*8 unita_frozen(3,3),unitg_frozen(3,3) 820 common / lattice_froze / unita_frozen,unitg_frozen, 821 > ecut_frozen,wcut_frozen, 822 > omega_frozen,frozen 823 824 lattice_omega_frozen = omega_frozen 825 return 826 end 827 828* ******************************************** 829* * * 830* * lattice_wcut_frozen * 831* * * 832* ******************************************** 833 834* frozen lattice structure - used to determine masking arrays 835 836 real*8 function lattice_wcut_frozen() 837 implicit none 838 839* **** common blocks **** 840 logical frozen 841 real*8 ecut_frozen,wcut_frozen,omega_frozen 842 real*8 unita_frozen(3,3),unitg_frozen(3,3) 843 common / lattice_froze / unita_frozen,unitg_frozen, 844 > ecut_frozen,wcut_frozen, 845 > omega_frozen,frozen 846 847 lattice_wcut_frozen = wcut_frozen 848 return 849 end 850 851* ******************************************** 852* * * 853* * lattice_ecut_frozen * 854* * * 855* ******************************************** 856 857* frozen lattice structure - used to determine masking arrays 858 859 real*8 function lattice_ecut_frozen() 860 implicit none 861 862* **** common blocks **** 863 logical frozen 864 real*8 ecut_frozen,wcut_frozen,omega_frozen 865 real*8 unita_frozen(3,3),unitg_frozen(3,3) 866 common / lattice_froze / unita_frozen,unitg_frozen, 867 > ecut_frozen,wcut_frozen, 868 > omega_frozen,frozen 869 870 lattice_ecut_frozen = ecut_frozen 871 return 872 end 873 874* ******************************************** 875* * * 876* * lattice_ggcut_frozen * 877* * * 878* ******************************************** 879 880* frozen lattice structure - used to determine masking arrays 881 882 real*8 function lattice_ggcut_frozen() 883 implicit none 884 885* **** common blocks **** 886 logical frozen 887 real*8 ecut_frozen,wcut_frozen,omega_frozen 888 real*8 unita_frozen(3,3),unitg_frozen(3,3) 889 common / lattice_froze / unita_frozen,unitg_frozen, 890 > ecut_frozen,wcut_frozen, 891 > omega_frozen,frozen 892 893 lattice_ggcut_frozen = 2.0d0*ecut_frozen 894 return 895 end 896 897* ******************************************** 898* * * 899* * lattice_wggcut_frozen * 900* * * 901* ******************************************** 902 903* frozen lattice structure - used to determine masking arrays 904 905 real*8 function lattice_wggcut_frozen() 906 implicit none 907 908* **** common blocks **** 909 logical frozen 910 real*8 ecut_frozen,wcut_frozen,omega_frozen 911 real*8 unita_frozen(3,3),unitg_frozen(3,3) 912 common / lattice_froze / unita_frozen,unitg_frozen, 913 > ecut_frozen,wcut_frozen, 914 > omega_frozen,frozen 915 916 lattice_wggcut_frozen = 2.0d0*wcut_frozen 917 return 918 end 919 920* ******************************************** 921* * * 922* * lattice_frozen * 923* * * 924* ******************************************** 925 926* frozen lattice structure - used to determine masking arrays 927 928 logical function lattice_frozen() 929 implicit none 930 931* **** common blocks **** 932 logical frozen 933 real*8 ecut_frozen,wcut_frozen,omega_frozen 934 real*8 unita_frozen(3,3),unitg_frozen(3,3) 935 common / lattice_froze / unita_frozen,unitg_frozen, 936 > ecut_frozen,wcut_frozen, 937 > omega_frozen,frozen 938 939 lattice_frozen = frozen 940 return 941 end 942 943 944 945 946* ********************************************* 947* * * 948* * lattice_init * 949* * * 950* ********************************************* 951 952 subroutine lattice_init() 953 implicit none 954 955* **** common block **** 956 real*8 ecut,wcut,omega 957 real*8 unita(3,3),unitg(3,3) 958 common / lattice_block / unita,unitg,ecut,wcut,omega 959 960 real*8 ub(3,3) 961 common / lattice_block2 / ub 962 963 logical frozen 964 real*8 ecut_frozen,wcut_frozen,omega_frozen 965 real*8 unita_frozen(3,3),unitg_frozen(3,3) 966 common / lattice_froze / unita_frozen,unitg_frozen, 967 > ecut_frozen,wcut_frozen, 968 > omega_frozen,frozen 969 970 logical has_small 971 real*8 omega_small 972 real*8 unita_small(3,3),unitg_small(3,3) 973 common /lattice_small_block/ unita_small,unitg_small, 974 > omega_small,has_small 975 976 real*8 omega_frozen_small 977 real*8 unita_frozen_small(3,3),unitg_frozen_small(3,3) 978 common /lattice_small_frozen_block/ unita_frozen_small, 979 > unitg_frozen_small, 980 > omega_frozen_small 981 982 983* **** local variables **** 984 integer nx,ny,nz 985 integer nxh,nyh,nzh 986 real*8 gx,gy,gz,gg 987 real*8 gg1,gg2,gg3 988 real*8 ecut0,wcut0 989 990* **** external functions **** 991 logical control_frozen,control_has_ngrid_small 992 integer control_ngrid,control_ngrid_small 993 real*8 control_unita,control_ecut,control_wcut 994 real*8 control_unita_frozen 995 external control_frozen,control_has_ngrid_small 996 external control_ngrid,control_ngrid_small 997 external control_unita,control_ecut,control_wcut 998 external control_unita_frozen 999 1000 ecut0 = control_ecut() 1001 wcut0 = control_wcut() 1002 1003* **** define lattice **** 1004 unita(1,1) = control_unita(1,1) 1005 unita(2,1) = control_unita(2,1) 1006 unita(3,1) = control_unita(3,1) 1007 unita(1,2) = control_unita(1,2) 1008 unita(2,2) = control_unita(2,2) 1009 unita(3,2) = control_unita(3,2) 1010 unita(1,3) = control_unita(1,3) 1011 unita(2,3) = control_unita(2,3) 1012 unita(3,3) = control_unita(3,3) 1013 call get_cube(unita,unitg,omega) 1014 1015* **** define frozen lattice - note this is only different from unita if nwpw:frozen_lattice is set **** 1016 frozen = control_frozen() 1017 unita_frozen(1,1) = control_unita_frozen(1,1) 1018 unita_frozen(2,1) = control_unita_frozen(2,1) 1019 unita_frozen(3,1) = control_unita_frozen(3,1) 1020 unita_frozen(1,2) = control_unita_frozen(1,2) 1021 unita_frozen(2,2) = control_unita_frozen(2,2) 1022 unita_frozen(3,2) = control_unita_frozen(3,2) 1023 unita_frozen(1,3) = control_unita_frozen(1,3) 1024 unita_frozen(2,3) = control_unita_frozen(2,3) 1025 unita_frozen(3,3) = control_unita_frozen(3,3) 1026 call get_cube(unita_frozen,unitg_frozen,omega_frozen) 1027 1028* **** define ub **** 1029 call dcopy(9,unitg,1,ub,1) 1030 call dscal(9,(1.0d0/(8.0d0*datan(1.0d0))),ub,1) 1031 1032 1033 1034* *** set the ecut variable using the frozen lattice*** 1035c call D3dB_nx(1,nx) 1036c call D3dB_ny(1,ny) 1037c call D3dB_nz(1,nz) 1038 nx = control_ngrid(1) 1039 ny = control_ngrid(2) 1040 nz = control_ngrid(3) 1041 nxh = nx/2 1042 nyh = ny/2 1043 nzh = nz/2 1044 1045 gx = unitg_frozen(1,1)*dble(nxh) 1046 gy = unitg_frozen(2,1)*dble(nxh) 1047 gz = unitg_frozen(3,1)*dble(nxh) 1048 gg1 = gx*gx + gy*gy + gz*gz 1049 1050 gx = unitg_frozen(1,2)*dble(nyh) 1051 gy = unitg_frozen(2,2)*dble(nyh) 1052 gz = unitg_frozen(3,2)*dble(nyh) 1053 gg2 = gx*gx + gy*gy + gz*gz 1054 1055 gx = unitg_frozen(1,3)*dble(nzh) 1056 gy = unitg_frozen(2,3)*dble(nzh) 1057 gz = unitg_frozen(3,3)*dble(nzh) 1058 gg3 = gx*gx + gy*gy + gz*gz 1059 1060 gg = gg1 1061 if (gg2.lt.gg) gg=gg2 1062 if (gg3.lt.gg) gg=gg3 1063 1064 ecut = 0.5d0*gg 1065 if (ecut0.lt.ecut) then 1066 ecut = ecut0 1067 end if 1068 1069 wcut = ecut 1070 if (wcut0.lt.wcut) then 1071 wcut = wcut0 1072 end if 1073 1074* **** set ecut,wcut for frozen lattice **** 1075 wcut_frozen = wcut 1076 ecut_frozen = ecut 1077 1078 1079* **** small cell **** 1080 has_small = control_has_ngrid_small() 1081 if (has_small) then 1082 gg1 = dble(control_ngrid_small(1))/dble(nx) 1083 gg2 = dble(control_ngrid_small(2))/dble(ny) 1084 gg3 = dble(control_ngrid_small(3))/dble(nz) 1085 unita_small(1,1) = unita(1,1)*gg1 1086 unita_small(2,1) = unita(2,1)*gg1 1087 unita_small(3,1) = unita(3,1)*gg1 1088 unita_small(1,2) = unita(1,2)*gg2 1089 unita_small(2,2) = unita(2,2)*gg2 1090 unita_small(3,2) = unita(3,2)*gg2 1091 unita_small(1,3) = unita(1,3)*gg3 1092 unita_small(2,3) = unita(2,3)*gg3 1093 unita_small(3,3) = unita(3,3)*gg3 1094 call get_cube(unita_small,unitg_small,omega_small) 1095 unita_frozen_small(1,1) = unita_frozen(1,1)*gg1 1096 unita_frozen_small(2,1) = unita_frozen(2,1)*gg1 1097 unita_frozen_small(3,1) = unita_frozen(3,1)*gg1 1098 unita_frozen_small(1,2) = unita_frozen(1,2)*gg2 1099 unita_frozen_small(2,2) = unita_frozen(2,2)*gg2 1100 unita_frozen_small(3,2) = unita_frozen(3,2)*gg2 1101 unita_frozen_small(1,3) = unita_frozen(1,3)*gg3 1102 unita_frozen_small(2,3) = unita_frozen(2,3)*gg3 1103 unita_frozen_small(3,3) = unita_frozen(3,3)*gg3 1104 call get_cube(unita_frozen_small, 1105 > unitg_frozen_small, 1106 > omega_frozen_small) 1107 end if 1108 1109 return 1110 end 1111 1112 1113* ******************************* 1114* * * 1115* * lattice_p_grid * 1116* * * 1117* ******************************* 1118* 1119* This routine computes coordinates of grid points in 1120* the unit cell 1121* 1122* Uses - 1123* Parallel_taskid --- processor number 1124* D3dB_nx --- number of grid points in direction 1 1125* D3dB_ny --- number of grid points in direction 2 1126* D3dB_nz --- number of grid points in direction 2 1127* lattice_unita -- primitive lattice vectors in real space 1128* 1129* Exit - 1130* xs --- coordinates of grid points (Rx,Ry,Rz) 1131* 1132* 1133 subroutine lattice_p_grid(xs) 1134 implicit none 1135 real*8 xs(*) 1136 1137* **** local variables **** 1138 integer nfft3d,n2ft3d 1139 integer i,j,k,p,taskid 1140 integer idx,k1,k2,k3 1141 integer np1,np2,np3 1142 integer nph1,nph2,nph3 1143 real*8 a(3,3),dk1,dk2,dk3,twopi 1144 1145 1146* **** constants **** 1147 call Parallel2d_taskid_i(taskid) 1148 call D3dB_nfft3d(1,nfft3d) 1149 n2ft3d = 2*nfft3d 1150 call D3dB_nx(1,np1) 1151 call D3dB_ny(1,np2) 1152 call D3dB_nz(1,np3) 1153 twopi = 8.0d0*datan(1.0d0) 1154 1155 nph1 = np1/2 1156 nph2 = np2/2 1157 nph3 = np3/2 1158 1159* **** elemental vectors **** 1160 call dcopy(9,0.0d0,0,a,1) 1161 a(1,1) = twopi/dble(np1) 1162 a(2,2) = twopi/dble(np2) 1163 a(3,3) = twopi/dble(np3) 1164 1165 call dcopy(6*n2ft3d,0.0d0,0,xs,1) 1166 1167* **** grid points in coordination space **** 1168 do k3 = -nph3, nph3-1 1169 do k2 = -nph2, nph2-1 1170 do k1 = -nph1, nph1-1 1171 1172 i = k1 + nph1 1173 j = k2 + nph2 1174 k = k3 + nph3 1175 1176 call D3dB_ijktoindex2p(1,i+1,j+1,k+1,idx,p) 1177 if (p .eq. taskid) then 1178 dk1=dble(k1) 1179 dk2=dble(k2) 1180 dk3=dble(k3) 1181 xs(idx) =dcos(a(1,1)*dk1+a(1,2)*dk2+a(1,3)*dk3) 1182 xs(idx+n2ft3d) =dsin(a(1,1)*dk1+a(1,2)*dk2+a(1,3)*dk3) 1183 xs(idx+2*n2ft3d)=dcos(a(2,1)*dk1+a(2,2)*dk2+a(2,3)*dk3) 1184 xs(idx+3*n2ft3d)=dsin(a(2,1)*dk1+a(2,2)*dk2+a(2,3)*dk3) 1185 xs(idx+4*n2ft3d)=dcos(a(3,1)*dk1+a(3,2)*dk2+a(3,3)*dk3) 1186 xs(idx+5*n2ft3d)=dsin(a(3,1)*dk1+a(3,2)*dk2+a(3,3)*dk3) 1187 end if 1188 end do 1189 end do 1190 end do 1191 1192 return 1193 end 1194 1195 1196 1197 1198* ******************************* 1199* * * 1200* * lattice_r_grid * 1201* * * 1202* ******************************* 1203* 1204* This routine computes coordinates of grid points in 1205* the unit cell 1206* 1207* Uses - 1208* Parallel_taskid --- processor number 1209* D3dB_nx --- number of grid points in direction 1 1210* D3dB_ny --- number of grid points in direction 2 1211* D3dB_nz --- number of grid points in direction 2 1212* lattice_unita -- primitive lattice vectors in real space 1213* 1214* Exit - 1215* r --- coordinates of grid points (Rx,Ry,Rz) 1216* 1217* 1218 subroutine lattice_r_grid(r) 1219 implicit none 1220 real*8 r(3,*) 1221 1222* **** local variables **** 1223 integer nfft3d,n2ft3d 1224 integer i,j,k,p,taskid,tid,nthreads 1225 integer index,k1,k2,k3,it 1226 integer np1,np2,np3 1227 integer nph1,nph2,nph3 1228 real*8 a(3,3),dk1,dk2,dk3 1229 1230* **** external functions **** 1231 real*8 lattice_unita 1232 external lattice_unita 1233 integer Parallel_threadid,Parallel_nthreads 1234 external Parallel_threadid,Parallel_nthreads 1235 1236 1237* **** constants **** 1238 call Parallel2d_taskid_i(taskid) 1239 tid = Parallel_threadid() 1240 nthreads = Parallel_nthreads() 1241 call D3dB_nfft3d(1,nfft3d) 1242 n2ft3d = 2*nfft3d 1243 call D3dB_nx(1,np1) 1244 call D3dB_ny(1,np2) 1245 call D3dB_nz(1,np3) 1246 1247 nph1 = np1/2 1248 nph2 = np2/2 1249 nph3 = np3/2 1250 1251* **** elemental vectors **** 1252 do i=1,3 1253 a(i,1) = lattice_unita(i,1)/dble(np1) 1254 a(i,2) = lattice_unita(i,2)/dble(np2) 1255 a(i,3) = lattice_unita(i,3)/dble(np3) 1256 end do 1257 1258 !call dcopy(3*n2ft3d,0.0d0,0,r,1) 1259 call Parallel_shared_vector_zero(.true.,2*n2ft3d,r) 1260 1261* **** grid points in coordination space **** 1262 it = 0 1263 do k3 = -nph3, nph3-1 1264 do k2 = -nph2, nph2-1 1265 do k1 = -nph1, nph1-1 1266 1267 i = k1 + nph1 1268 j = k2 + nph2 1269 k = k3 + nph3 1270 1271 !call D3dB_ktoqp(1,k+1,q,p) 1272 call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p) 1273 if ((p.eq.taskid).and.(it.eq.tid)) then 1274c index = (q-1)*(np1+2)*np2 1275c > + j *(np1+2) 1276c > + i+1 1277 dk1=dble(k1) 1278 dk2=dble(k2) 1279 dk3=dble(k3) 1280 r(1,index) = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3 1281 r(2,index) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3 1282 r(3,index) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3 1283 1284c* **** reverse y and z **** 1285c r(1,index) = a(1,1)*k1 + a(1,2)*k3 + a(1,3)*k2 1286c r(2,index) = a(2,1)*k1 + a(2,2)*k3 + a(2,3)*k2 1287c r(3,index) = a(3,1)*k1 + a(3,2)*k3 + a(3,3)*k2 1288 1289 end if 1290 it = mod(it+1,nthreads) 1291 end do 1292 end do 1293 end do 1294 1295 return 1296 end 1297 1298 1299 subroutine lattice_r_grid_sym(r) 1300 implicit none 1301 real*8 r(3,*) 1302 1303* **** local variables **** 1304 integer nfft3d,n2ft3d 1305 integer i,j,k,p,taskid 1306 integer index,k1,k2,k3 1307 integer np1,np2,np3 1308 integer nph1,nph2,nph3 1309 real*8 a(3,3),dk1,dk2,dk3 1310 1311* **** external functions **** 1312 real*8 lattice_unita 1313 external lattice_unita 1314 1315 1316* **** constants **** 1317 call Parallel2d_taskid_i(taskid) 1318 call D3dB_nfft3d(1,nfft3d) 1319 n2ft3d = 2*nfft3d 1320 call D3dB_nx(1,np1) 1321 call D3dB_ny(1,np2) 1322 call D3dB_nz(1,np3) 1323 1324 nph1 = np1/2 1325 nph2 = np2/2 1326 nph3 = np3/2 1327 1328* **** elemental vectors **** 1329 do i=1,3 1330 a(i,1) = lattice_unita(i,1)/dble(np1) 1331 a(i,2) = lattice_unita(i,2)/dble(np2) 1332 a(i,3) = lattice_unita(i,3)/dble(np3) 1333 end do 1334 1335 call dcopy(3*n2ft3d,0.0d0,0,r,1) 1336 1337* **** grid points in coordination space **** 1338 do k3 = -nph3+1, nph3-1 1339 do k2 = -nph2+1, nph2-1 1340 do k1 = -nph1+1, nph1-1 1341 1342 i = k1 + nph1 1343 j = k2 + nph2 1344 k = k3 + nph3 1345 1346 !call D3dB_ktoqp(1,k+1,q,p) 1347 call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p) 1348 if (p .eq. taskid) then 1349c index = (q-1)*(np1+2)*np2 1350c > + j *(np1+2) 1351c > + i+1 1352 dk1=dble(k1) 1353 dk2=dble(k2) 1354 dk3=dble(k3) 1355 r(1,index) = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3 1356 r(2,index) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3 1357 r(3,index) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3 1358 1359c* **** reverse y and z **** 1360c r(1,index) = a(1,1)*k1 + a(1,2)*k3 + a(1,3)*k2 1361c r(2,index) = a(2,1)*k1 + a(2,2)*k3 + a(2,3)*k2 1362c r(3,index) = a(3,1)*k1 + a(3,2)*k3 + a(3,3)*k2 1363 1364 end if 1365 end do 1366 end do 1367 end do 1368 1369 return 1370 end 1371 1372 1373 subroutine lattice_r_grid_sym0(r) 1374 implicit none 1375 real*8 r(*) 1376 1377* **** local variables **** 1378 integer nfft3d,n2ft3d 1379 integer i,j,k,p,taskid 1380 integer index,k1,k2,k3 1381 integer np1,np2,np3 1382 integer nph1,nph2,nph3 1383 real*8 a(3,3),dk1,dk2,dk3 1384 1385* **** external functions **** 1386 real*8 lattice_unita 1387 external lattice_unita 1388 1389 1390* **** constants **** 1391 call Parallel2d_taskid_i(taskid) 1392 call D3dB_nfft3d(1,nfft3d) 1393 n2ft3d = 2*nfft3d 1394 call D3dB_nx(1,np1) 1395 call D3dB_ny(1,np2) 1396 call D3dB_nz(1,np3) 1397 1398 nph1 = np1/2 1399 nph2 = np2/2 1400 nph3 = np3/2 1401 1402* **** elemental vectors **** 1403 do i=1,3 1404 a(i,1) = lattice_unita(i,1)/dble(np1) 1405 a(i,2) = lattice_unita(i,2)/dble(np2) 1406 a(i,3) = lattice_unita(i,3)/dble(np3) 1407 end do 1408 1409 call dcopy(3*n2ft3d,0.0d0,0,r,1) 1410 1411* **** grid points in coordination space **** 1412 do k3 = -nph3+1, nph3-1 1413 do k2 = -nph2+1, nph2-1 1414 do k1 = -nph1+1, nph1-1 1415 1416 i = k1 + nph1 1417 j = k2 + nph2 1418 k = k3 + nph3 1419 1420 !call D3dB_ktoqp(1,k+1,q,p) 1421 call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p) 1422 if (p .eq. taskid) then 1423c index = (q-1)*(np1+2)*np2 1424c > + j *(np1+2) 1425c > + i+1 1426 dk1=dble(k1) 1427 dk2=dble(k2) 1428 dk3=dble(k3) 1429 r(index) = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3 1430 r(index+ n2ft3d) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3 1431 r(index+2*n2ft3d) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3 1432 end if 1433 end do 1434 end do 1435 end do 1436 1437 return 1438 end 1439 1440 1441* ******************************* 1442* * * 1443* * lattice_i_grid * 1444* * * 1445* ******************************* 1446* 1447* This routine computes coordinates of grid points in 1448* the unit cell 1449* 1450* Uses - 1451* Parallel_taskid --- processor number 1452* D3dB_nx --- number of grid points in direction 1 1453* D3dB_ny --- number of grid points in direction 2 1454* D3dB_nz --- number of grid points in direction 2 1455* 1456* Entry - nb 1457* Exit - 1458* ijk --- coordinates of grid points (k1,k2,k3) 1459* 1460* 1461 subroutine lattice_i_grid(nb,ijk) 1462 implicit none 1463 integer nb 1464 integer ijk(4,*) 1465 1466* **** local variables **** 1467 integer nfft3d,n2ft3d 1468 integer i,j,k,p,taskid 1469 integer index,k1,k2,k3 1470 integer np1,np2,np3 1471 integer nph1,nph2,nph3 1472 1473* **** constants **** 1474 call Parallel2d_taskid_i(taskid) 1475 call D3dB_nfft3d(nb,nfft3d) 1476 1477 n2ft3d = 2*nfft3d 1478 call D3dB_nx(nb,np1) 1479 call D3dB_ny(nb,np2) 1480 call D3dB_nz(nb,np3) 1481 1482 nph1 = np1/2 1483 nph2 = np2/2 1484 nph3 = np3/2 1485 1486 call icopy(4*n2ft3d,0,0,ijk,1) 1487 1488* **** grid points **** 1489 do k3 = -nph3, nph3-1 1490 do k2 = -nph2, nph2-1 1491 do k1 = -nph1, nph1-1 1492 1493 i = k1 + nph1 1494 j = k2 + nph2 1495 k = k3 + nph3 1496 1497 call D3dB_ijktoindex2p(nb,i+1,j+1,k+1,index,p) 1498 if (p .eq. taskid) then 1499 ijk(1,index) = k1 1500 ijk(2,index) = k2 1501 ijk(3,index) = k3 1502 ijk(4,index) = 1 1503 end if 1504 end do 1505 end do 1506 end do 1507 1508 return 1509 end 1510 1511 1512 1513 1514 1515 1516 1517 subroutine lattice_mask_sym(r) 1518 implicit none 1519 real*8 r(*) 1520 1521* **** local variables **** 1522 integer nfft3d,n2ft3d 1523 integer i,j,k,p,taskid 1524 integer index,k1,k2,k3 1525 integer np1,np2,np3 1526 integer nph1,nph2,nph3 1527 1528 1529* **** constants **** 1530 call Parallel2d_taskid_i(taskid) 1531 call D3dB_nfft3d(1,nfft3d) 1532 n2ft3d = 2*nfft3d 1533 call D3dB_nx(1,np1) 1534 call D3dB_ny(1,np2) 1535 call D3dB_nz(1,np3) 1536 1537 nph1 = np1/2 1538 nph2 = np2/2 1539 nph3 = np3/2 1540 1541 1542 call dcopy(n2ft3d,0.0d0,0,r,1) 1543 1544* **** grid points in coordination space **** 1545 do k3 = -nph3+1, nph3-1 1546 do k2 = -nph2+1, nph2-1 1547 do k1 = -nph1+1, nph1-1 1548 1549 i = k1 + nph1 1550 j = k2 + nph2 1551 k = k3 + nph3 1552 1553 !call D3dB_ktoqp(1,k+1,q,p) 1554 call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p) 1555 if (p .eq. taskid) then 1556c index = (q-1)*(np1+2)*np2 1557c > + j *(np1+2) 1558c > + i+1 1559 1560 r(index) = 1.0d0 1561 end if 1562 end do 1563 end do 1564 end do 1565 1566 return 1567 end 1568 1569 1570 1571 1572 subroutine get_cube(unita,unitg,volume) 1573 1574****************************************************************************** 1575* * 1576* This routine computes primitive vectors both in coordination * 1577* space and in reciporocal space and the volume of primitive cell. * 1578* * 1579* Inputs: * 1580* type --- type of cube (1=SC, 2=FCC, 3=BCC, 4=linear) * 1581* unit --- lattice constants * 1582* * 1583* Outputs: * 1584* volume --- volume of primitive cell * 1585* unita --- primitive vectors in coordination space * 1586* unitg --- primitive vectors in reciprocal space * 1587* * 1588* Library: DSCAL from BLAS * 1589* * 1590* Last modification: 7/03/93 by R. Kawai * 1591* * 1592****************************************************************************** 1593 1594 implicit none 1595 1596* ------------------ 1597* argument variables 1598* ------------------ 1599 double precision unita(3,3), unitg(3,3) 1600 double precision volume 1601 1602* --------------- 1603* local variables 1604* --------------- 1605 double precision twopi 1606 1607 twopi = 8.0d0*datan(1.0d0) 1608 1609 1610* ----------------------------------------- 1611* primitive vectors in the reciprocal space 1612* ----------------------------------------- 1613 unitg(1,1) = unita(2,2)*unita(3,3) - unita(3,2)*unita(2,3) 1614 unitg(2,1) = unita(3,2)*unita(1,3) - unita(1,2)*unita(3,3) 1615 unitg(3,1) = unita(1,2)*unita(2,3) - unita(2,2)*unita(1,3) 1616 unitg(1,2) = unita(2,3)*unita(3,1) - unita(3,3)*unita(2,1) 1617 unitg(2,2) = unita(3,3)*unita(1,1) - unita(1,3)*unita(3,1) 1618 unitg(3,2) = unita(1,3)*unita(2,1) - unita(2,3)*unita(1,1) 1619 unitg(1,3) = unita(2,1)*unita(3,2) - unita(3,1)*unita(2,2) 1620 unitg(2,3) = unita(3,1)*unita(1,2) - unita(1,1)*unita(3,2) 1621 unitg(3,3) = unita(1,1)*unita(2,2) - unita(2,1)*unita(1,2) 1622* --------------------- 1623* volume of a unit cell 1624* --------------------- 1625 volume = unita(1,1)*unitg(1,1) 1626 > + unita(2,1)*unitg(2,1) 1627 > + unita(3,1)*unitg(3,1) 1628 call dscal(9,twopi/volume,unitg,1) 1629 1630 volume=dabs(volume) 1631 return 1632 end 1633 1634* ******************************* 1635* * * 1636* * lattice_abc_abg * 1637* * * 1638* ******************************* 1639* 1640* This routine computes a,b,c,alpha,beta,gamma. 1641* 1642 subroutine lattice_abc_abg(a,b,c,alpha,beta,gamma) 1643 implicit none 1644 real*8 a,b,c 1645 real*8 alpha,beta,gamma 1646 1647* *** local variables **** 1648 real*8 d2,pi 1649 1650* **** external functions **** 1651 real*8 lattice_unita 1652 external lattice_unita 1653 1654* **** determine a,b,c,alpha,beta,gmma *** 1655 pi = 4.0d0*datan(1.0d0) 1656 a = dsqrt(lattice_unita(1,1)**2 1657 > + lattice_unita(2,1)**2 1658 > + lattice_unita(3,1)**2) 1659 b = dsqrt(lattice_unita(1,2)**2 1660 > + lattice_unita(2,2)**2 1661 > + lattice_unita(3,2)**2) 1662 c = dsqrt(lattice_unita(1,3)**2 1663 > + lattice_unita(2,3)**2 1664 > + lattice_unita(3,3)**2) 1665 1666 d2 = (lattice_unita(1,2)-lattice_unita(1,3))**2 1667 > + (lattice_unita(2,2)-lattice_unita(2,3))**2 1668 > + (lattice_unita(3,2)-lattice_unita(3,3))**2 1669 alpha = (b*b + c*c - d2)/(2.0d0*b*c) 1670 alpha = dacos(alpha)*180.0d0/pi 1671 1672 d2 = (lattice_unita(1,3)-lattice_unita(1,1))**2 1673 > + (lattice_unita(2,3)-lattice_unita(2,1))**2 1674 > + (lattice_unita(3,3)-lattice_unita(3,1))**2 1675 beta = (c*c + a*a - d2)/(2.0d0*c*a) 1676 beta = dacos(beta)*180.0d0/pi 1677 1678 d2 = (lattice_unita(1,1)-lattice_unita(1,2))**2 1679 > + (lattice_unita(2,1)-lattice_unita(2,2))**2 1680 > + (lattice_unita(3,1)-lattice_unita(3,2))**2 1681 gamma = (a*a + b*b - d2)/(2.0d0*a*b) 1682 gamma = dacos(gamma)*180.0d0/pi 1683 1684 return 1685 end 1686 1687 1688 1689* ************************************* 1690* * * 1691* * lattice_small_abc_abg * 1692* * * 1693* ************************************* 1694* 1695* This routine computes a,b,c,alpha,beta,gamma. 1696* 1697 subroutine lattice_small_abc_abg(a,b,c,alpha,beta,gamma) 1698 implicit none 1699 real*8 a,b,c 1700 real*8 alpha,beta,gamma 1701 1702* *** local variables **** 1703 real*8 d2,pi 1704 1705* **** external functions **** 1706 real*8 lattice_unita_small 1707 external lattice_unita_small 1708 1709* **** determine a,b,c,alpha,beta,gmma *** 1710 pi = 4.0d0*datan(1.0d0) 1711 a = dsqrt(lattice_unita_small(1,1)**2 1712 > + lattice_unita_small(2,1)**2 1713 > + lattice_unita_small(3,1)**2) 1714 b = dsqrt(lattice_unita_small(1,2)**2 1715 > + lattice_unita_small(2,2)**2 1716 > + lattice_unita_small(3,2)**2) 1717 c = dsqrt(lattice_unita_small(1,3)**2 1718 > + lattice_unita_small(2,3)**2 1719 > + lattice_unita_small(3,3)**2) 1720 1721 d2 = (lattice_unita_small(1,2)-lattice_unita_small(1,3))**2 1722 > + (lattice_unita_small(2,2)-lattice_unita_small(2,3))**2 1723 > + (lattice_unita_small(3,2)-lattice_unita_small(3,3))**2 1724 alpha = (b*b + c*c - d2)/(2.0d0*b*c) 1725 alpha = dacos(alpha)*180.0d0/pi 1726 1727 d2 = (lattice_unita_small(1,3)-lattice_unita_small(1,1))**2 1728 > + (lattice_unita_small(2,3)-lattice_unita_small(2,1))**2 1729 > + (lattice_unita_small(3,3)-lattice_unita_small(3,1))**2 1730 beta = (c*c + a*a - d2)/(2.0d0*c*a) 1731 beta = dacos(beta)*180.0d0/pi 1732 1733 d2 = (lattice_unita_small(1,1)-lattice_unita_small(1,2))**2 1734 > + (lattice_unita_small(2,1)-lattice_unita_small(2,2))**2 1735 > + (lattice_unita_small(3,1)-lattice_unita_small(3,2))**2 1736 gamma = (a*a + b*b - d2)/(2.0d0*a*b) 1737 gamma = dacos(gamma)*180.0d0/pi 1738 1739 return 1740 end 1741 1742 1743 1744 1745