1* 2* $Id$ 3* 4 5* ****************************** 6* * * 7* * Calculate_Dipole * 8* * * 9* ****************************** 10 11 subroutine Calculate_Dipole(ispin,ne,n2ft3d,dn,dipole) 12 implicit none 13 integer ispin,ne(2) 14 integer n2ft3d 15 real*8 dn(n2ft3d,ispin) 16 real*8 dipole(3) 17 18#include "bafdecls.fh" 19#include "errquit.fh" 20#include "util.fh" 21#include "stdio.fh" 22 23 24* **** local variables **** 25 logical value,oprint 26 integer ii 27 integer nx,ny,nz 28 integer n1,n2,n3,ncut 29 real*8 GX,GY,GZ 30 real*8 qGX,qGY,qGZ 31 real*8 cqGX,cqGY,cqGZ,x,y,z,r,rmax 32 real*8 cdx1,cdy1,cdz1 33 real*8 cdx2,cdy2,cdz2 34 real*8 cdx3,cdy3,cdz3 35 real*8 tmass,tcharge,ncharge,pcharge 36 real*8 dv 37 real*8 dipole_crystal(3),dipole_molecule(3) 38 real*8 dtmp(3) 39 40 integer rgrid(2) 41 integer rgx(2),rgy(2),rgz(2) 42 43 integer taskid,MASTER 44 parameter (MASTER=0) 45 46 real*8 autoDebye 47 parameter (autoDebye=2.5416d0) 48 49* **** external functions **** 50 logical control_print 51 integer ion_katm,ion_nion,control_ncut,control_version 52 real*8 ion_amass,psp_zv,ion_rion,lattice_omega 53 real*8 lattice_unita 54 external control_print 55 external ion_katm,ion_nion,control_ncut,control_version 56 external ion_amass,psp_zv,ion_rion,lattice_omega 57 external lattice_unita 58 59 call Parallel_taskid(taskid) 60 oprint= ((taskid.eq.MASTER).and.control_print(print_medium)) 61 62 63* ***** center of mass **** 64 GX=0.0d0 65 GY=0.0d0 66 GZ=0.0d0 67 tmass=0.0d0 68 DO ii=1,ion_nion() 69 tmass=tmass+ion_amass(ii) 70 GX=GX+ion_amass(ii)*ion_rion(1,ii) 71 GY=GY+ion_amass(ii)*ion_rion(2,ii) 72 GZ=GZ+ion_amass(ii)*ion_rion(3,ii) 73 END DO 74 GX=GX/tmass 75 GY=GY/tmass 76 GZ=GZ/tmass 77 78 !*** crystal center of ionic charge *** 79 ncut = 20 80 n1 = ncut-2 81 n2 = ncut-2 82 n3 = ncut-2 83 84 x = n1*lattice_unita(1,1) 85 y = n1*lattice_unita(2,1) 86 z = n1*lattice_unita(3,1) 87 rmax = dsqrt(x*x + y*y + z*z) 88 89 x = n2*lattice_unita(1,2) 90 y = n2*lattice_unita(2,2) 91 z = n2*lattice_unita(3,2) 92 r = dsqrt(x*x + y*y + z*z) 93 if (r.lt.rmax) rmax = r 94 95 x = n3*lattice_unita(1,3) 96 y = n3*lattice_unita(2,3) 97 z = n3*lattice_unita(3,3) 98 r = dsqrt(x*x + y*y + z*z) 99 if (r.lt.rmax) rmax = r 100 101 cqGX=0.0d0 102 cqGY=0.0d0 103 cqGZ=0.0d0 104 tcharge=0.0d0 105 do ii=1,ion_nion() 106 107 do n3= -ncut, ncut 108 do n2= -ncut, ncut 109 do n1= -ncut, ncut 110 x = ion_rion(1,ii) 111 > + n1*lattice_unita(1,1) 112 > + n2*lattice_unita(1,2) 113 > + n3*lattice_unita(1,3) 114 y = ion_rion(2,ii) 115 > + n1*lattice_unita(2,1) 116 > + n2*lattice_unita(2,2) 117 > + n3*lattice_unita(2,3) 118 z = ion_rion(3,ii) 119 > + n1*lattice_unita(3,1) 120 > + n2*lattice_unita(3,2) 121 > + n3*lattice_unita(3,3) 122 123 r = dsqrt(x*x+y*y+z*z) 124 125 if (r.le.rmax) then 126 cqGX=cqGX+psp_zv(ion_katm(ii))*x 127 cqGY=cqGY+psp_zv(ion_katm(ii))*y 128 cqGZ=cqGZ+psp_zv(ion_katm(ii))*z 129 tcharge=tcharge+psp_zv(ion_katm(ii)) 130 end if 131 end do 132 end do 133 end do 134 END DO 135 cqGX=cqGX/tcharge 136 cqGY=cqGY/tcharge 137 cqGZ=cqGZ/tcharge 138 139 140 141 !*** molecular center of ionic charge *** 142 qGX=0.0d0 143 qGY=0.0d0 144 qGZ=0.0d0 145 tcharge=0.0d0 146 DO ii=1,ion_nion() 147 tcharge=tcharge+psp_zv(ion_katm(ii)) 148 qGX=qGX+psp_zv(ion_katm(ii))*ion_rion(1,ii) 149 qGY=qGY+psp_zv(ion_katm(ii))*ion_rion(2,ii) 150 qGZ=qGZ+psp_zv(ion_katm(ii))*ion_rion(3,ii) 151 END DO 152 qGX=qGX/tcharge 153 qGY=qGY/tcharge 154 qGZ=qGZ/tcharge 155 156 157 158* **** calculate the center of density **** 159 value = BA_push_get(mt_dbl,3*n2ft3d,'rgrid', rgrid(2), rgrid(1)) 160 value = value.and. 161 > BA_push_get(mt_dbl, n2ft3d,'rgx',rgx(2),rgx(1)) 162 value = value.and. 163 > BA_push_get(mt_dbl, n2ft3d,'rgy',rgy(2),rgy(1)) 164 value = value.and. 165 > BA_push_get(mt_dbl, n2ft3d,'rgz',rgz(2),rgz(1)) 166 if (.not. value) 167 > call errquit('Calculate_Dipole: out of stack memory',0, MA_ERR) 168 169 call D3dB_nx(1,nx) 170 call D3dB_ny(1,ny) 171 call D3dB_nz(1,nz) 172 dv=lattice_omega()/dble(nx*ny*nz) 173 call lattice_r_grid_sym(dbl_mb(rgrid(1))) 174 call dcopy(n2ft3d,dbl_mb(rgrid(1)+0),3,dbl_mb(rgx(1)),1) 175 call dcopy(n2ft3d,dbl_mb(rgrid(1)+1),3,dbl_mb(rgy(1)),1) 176 call dcopy(n2ft3d,dbl_mb(rgrid(1)+2),3,dbl_mb(rgz(1)),1) 177 call D3dB_r_Zero_Ends(1,dbl_mb(rgx(1))) 178 call D3dB_r_Zero_Ends(1,dbl_mb(rgy(1))) 179 call D3dB_r_Zero_Ends(1,dbl_mb(rgz(1))) 180 181 call D3dB_rr_dot(1,dbl_mb(rgx(1)),dn(1,1),cdx1) 182 call D3dB_rr_dot(1,dbl_mb(rgy(1)),dn(1,1),cdy1) 183 call D3dB_rr_dot(1,dbl_mb(rgz(1)),dn(1,1),cdz1) 184 cdx1 = cdx1*dv 185 cdy1 = cdy1*dv 186 cdz1 = cdz1*dv 187 188 189* *** check for ferromagnetic case *** 190 if (ne(ispin).ne.0) then 191 call D3dB_rr_dot(1,dbl_mb(rgx(1)),dn(1,ispin),cdx2) 192 call D3dB_rr_dot(1,dbl_mb(rgy(1)),dn(1,ispin),cdy2) 193 call D3dB_rr_dot(1,dbl_mb(rgz(1)),dn(1,ispin),cdz2) 194 cdx2 = cdx2*dv 195 cdy2 = cdy2*dv 196 cdz2 = cdz2*dv 197 else 198 cdx2 = 0.0d0 199 cdy2 = 0.0d0 200 cdz2 = 0.0d0 201 end if 202 203 cdx3=cdx1+cdx2 204 cdy3=cdy1+cdy2 205 cdz3=cdz1+cdz2 206 207 call lattice_mask_sym(dbl_mb(rgrid(1))) 208 !cdx1=cdx1/ne(1) 209 !cdy1=cdy1/ne(1) 210 !cdz1=cdz1/ne(1) 211 call D3dB_rr_dot(1,dbl_mb(rgrid(1)),dn(1,1),rmax) 212 rmax = rmax*dv 213 cdx1=cdx1/rmax 214 cdy1=cdy1/rmax 215 cdz1=cdz1/rmax 216 if (ne(ispin).ne.0) then 217 !cdx2=cdx2/ne(ispin) 218 !cdy2=cdy2/ne(ispin) 219 !cdz2=cdz2/ne(ispin) 220 call D3dB_rr_dot(1,dbl_mb(rgrid(1)),dn(1,ispin),rmax) 221 rmax = rmax*dv 222 cdx2=cdx2/rmax 223 cdy2=cdy2/rmax 224 cdz2=cdz2/rmax 225 226 end if 227 !cdx3=cdx3/dble(ne(1)+ne(ispin)) 228 !cdy3=cdy3/dble(ne(1)+ne(ispin)) 229 !cdz3=cdz3/dble(ne(1)+ne(ispin)) 230 call D3dB_rr_Sum(1,dn(1,1), 231 > dn(1,ispin), 232 > dbl_mb(rgrid(1)+n2ft3d)) 233 call D3dB_rr_dot(1,dbl_mb(rgrid(1)), 234 > dbl_mb(rgrid(1)+n2ft3d), 235 > rmax) 236 rmax = rmax*dv 237 cdx3=cdx3/rmax 238 cdy3=cdy3/rmax 239 cdz3=cdz3/rmax 240 241 if (oprint) then 242 write(LuOut,1200) 243 write(LuOut,1220) 'spin up ',CDX1,CDY1,CDZ1 244 if (ne(ispin).ne.0) 245 > write(LuOut,1220) 'spin down ',CDX2,CDY2,CDZ2 246 write(LuOut,1220) ' total',CDX3,CDY3,CDZ3 247 write(LuOut,1220) 'ionic ',qGX,qGY,qGZ 248 !write(LuOut,1220) 'crystal ',cqGX,cqGY,cqGZ 249 end if 250 dtmp(1) = CDX3 251 dtmp(2) = CDY3 252 dtmp(3) = CDZ3 253 call ecce_print1('total dipole',mt_dbl,dtmp,3) 254 dtmp(1) = CDX1 255 dtmp(2) = CDY1 256 dtmp(3) = CDZ1 257 call ecce_print1('alpha dipole',mt_dbl,dtmp,3) 258 if (ne(ispin).ne.0) then 259 dtmp(1) = CDX2 260 dtmp(2) = CDY2 261 dtmp(3) = CDZ2 262 call ecce_print1('beta dipole',mt_dbl,dtmp,3) 263 endif 264 dtmp(1) = qGX 265 dtmp(2) = qGy 266 dtmp(3) = qGz 267 call ecce_print1('nuclear dipole',mt_dbl,dtmp,3) 268 269 270 1200 FORMAT(//'== Center of Charge =='/) 271 1220 FORMAT(A10,' (',F10.4,',',F10.4,',',F10.4,' )') 272 273c* ***** calculate crystal dipole with respect to center of cell **** 274c pcharge = tcharge 275c ncharge = dble(ne(1)+ne(ispin)) 276c dipole_crystal(1) = -ncharge*cdx3 + pcharge*cqGX 277c dipole_crystal(2) = -ncharge*cdy3 + pcharge*cqGY 278c dipole_crystal(3) = -ncharge*cdz3 + pcharge*cqGZ 279c cdx1 = dsqrt( dipole_crystal(1)**2 280c > + dipole_crystal(2)**2 281c > + dipole_crystal(3)**2) 282c if (oprint) then 283c write(LuOut,1240) 284c write(LuOut,1231) dipole_crystal 285c write(LuOut,1232) cdx1,cdx1*autoDebye 286c end if 287 288* ***** calculate dipole with respect to center of mass **** 289 pcharge = tcharge 290 ncharge = dble(ne(1)+ne(ispin)) 291 dipole_molecule(1) = -ncharge*cdx3 + pcharge*qGX 292 > - GX*(pcharge-ncharge) 293 dipole_molecule(2) = -ncharge*cdy3 + pcharge*qGY 294 > - GY*(pcharge-ncharge) 295 dipole_molecule(3) = -ncharge*cdz3 + pcharge*qGZ 296 > - GZ*(pcharge-ncharge) 297 cdx1 = dsqrt( dipole_molecule(1)**2 298 > + dipole_molecule(2)**2 299 > + dipole_molecule(3)**2) 300 if (oprint) then 301 write(LuOut,1230) 302 write(LuOut,1231) dipole_molecule 303 write(LuOut,1232) cdx1,cdx1*autoDebye 304 end if 305 1230 FORMAT(//'== Molecular Dipole wrt Center of Mass =='/) 306 1231 FORMAT('mu = (',F10.4,',',F10.4,',',F10.4,' ) au') 307 1232 FORMAT('|mu| = ',F10.4,' au, ',F10.4,' Debye') 308 1240 FORMAT(//'== Crystal Dipole =='/) 309 310* **** pop stack memory **** 311 value = value.and.BA_pop_stack(rgz(2)) 312 value = value.and.BA_pop_stack(rgy(2)) 313 value = value.and.BA_pop_stack(rgx(2)) 314 value = value.and.BA_pop_stack(rgrid(2)) 315 if (.not. value) 316 > call errquit('Calculate_Dipole: cannot pop stack memory',0, 317 & MA_ERR) 318 319c if (control_version().eq.3) then 320c dipole(1) = dipole_crystal(1) 321c dipole(2) = dipole_crystal(2) 322c dipole(3) = dipole_crystal(3) 323c else 324c dipole(1) = dipole_molecule(1) 325c dipole(2) = dipole_molecule(2) 326c dipole(3) = dipole_molecule(3) 327c end if 328 dipole(1) = dipole_molecule(1) 329 dipole(2) = dipole_molecule(2) 330 dipole(3) = dipole_molecule(3) 331 332 return 333 end 334 335 336 337* ****************************** 338* * * 339* * Calculate_Molecular_Dipole * 340* * * 341* ****************************** 342 343 subroutine Calculate_Molecular_Dipole(ispin,ne,n2ft3d,dn,dipole) 344 implicit none 345 integer ispin,ne(2) 346 integer n2ft3d 347 real*8 dn(n2ft3d,ispin) 348 real*8 dipole(3) 349 350#include "bafdecls.fh" 351#include "errquit.fh" 352#include "util.fh" 353 354 355* **** local variables **** 356 logical value,oprint 357 integer ii 358 integer nx,ny,nz 359 integer n1,n2,n3,ncut 360 real*8 GX,GY,GZ 361 real*8 qGX,qGY,qGZ 362 real*8 x,y,z,r,rmax 363 real*8 cdx1,cdy1,cdz1 364 real*8 cdx2,cdy2,cdz2 365 real*8 cdx3,cdy3,cdz3 366 real*8 tmass,tcharge,ncharge,pcharge 367 real*8 dv 368 369 integer rgrid(2) 370 integer rgx(2),rgy(2),rgz(2) 371 372 real*8 autoDebye 373 parameter (autoDebye=2.5416d0) 374 375* **** external functions **** 376 integer ion_katm,ion_nion 377 real*8 ion_amass,psp_zv,ion_rion,lattice_omega 378 external ion_katm,ion_nion 379 external ion_amass,psp_zv,ion_rion,lattice_omega 380 381 382* ***** center of mass **** 383 GX=0.0d0 384 GY=0.0d0 385 GZ=0.0d0 386 tmass=0.0d0 387 DO ii=1,ion_nion() 388 tmass=tmass+ion_amass(ii) 389 GX=GX+ion_amass(ii)*ion_rion(1,ii) 390 GY=GY+ion_amass(ii)*ion_rion(2,ii) 391 GZ=GZ+ion_amass(ii)*ion_rion(3,ii) 392 END DO 393 GX=GX/tmass 394 GY=GY/tmass 395 GZ=GZ/tmass 396 397 398 !*** molecular center of ionic charge *** 399 qGX=0.0d0 400 qGY=0.0d0 401 qGZ=0.0d0 402 tcharge=0.0d0 403 DO ii=1,ion_nion() 404 tcharge=tcharge+psp_zv(ion_katm(ii)) 405 qGX=qGX+psp_zv(ion_katm(ii))*ion_rion(1,ii) 406 qGY=qGY+psp_zv(ion_katm(ii))*ion_rion(2,ii) 407 qGZ=qGZ+psp_zv(ion_katm(ii))*ion_rion(3,ii) 408 END DO 409 qGX=qGX/tcharge 410 qGY=qGY/tcharge 411 qGZ=qGZ/tcharge 412 413* **** calculate the center of density **** 414 value = BA_push_get(mt_dbl,3*n2ft3d,'rgrid',rgrid(2),rgrid(1)) 415 value = value.and. 416 > BA_push_get(mt_dbl, n2ft3d,'rgx',rgx(2),rgx(1)) 417 value = value.and. 418 > BA_push_get(mt_dbl, n2ft3d,'rgy',rgy(2),rgy(1)) 419 value = value.and. 420 > BA_push_get(mt_dbl, n2ft3d,'rgz',rgz(2),rgz(1)) 421 if (.not. value) 422 > call errquit('Calculate_Molecular Dipole:out of stack',0,MA_ERR) 423 424 call D3dB_nx(1,nx) 425 call D3dB_ny(1,ny) 426 call D3dB_nz(1,nz) 427 dv=lattice_omega()/dble(nx*ny*nz) 428 call lattice_r_grid_sym(dbl_mb(rgrid(1))) 429 call dcopy(n2ft3d,dbl_mb(rgrid(1)+0),3,dbl_mb(rgx(1)),1) 430 call dcopy(n2ft3d,dbl_mb(rgrid(1)+1),3,dbl_mb(rgy(1)),1) 431 call dcopy(n2ft3d,dbl_mb(rgrid(1)+2),3,dbl_mb(rgz(1)),1) 432 call D3dB_r_Zero_Ends(1,dbl_mb(rgx(1))) 433 call D3dB_r_Zero_Ends(1,dbl_mb(rgy(1))) 434 call D3dB_r_Zero_Ends(1,dbl_mb(rgz(1))) 435 436 call D3dB_rr_dot(1,dbl_mb(rgx(1)),dn(1,1),cdx1) 437 call D3dB_rr_dot(1,dbl_mb(rgy(1)),dn(1,1),cdy1) 438 call D3dB_rr_dot(1,dbl_mb(rgz(1)),dn(1,1),cdz1) 439 cdx1 = cdx1*dv 440 cdy1 = cdy1*dv 441 cdz1 = cdz1*dv 442 443* *** check for ferromagnetic case *** 444 if (ne(ispin).ne.0) then 445 call D3dB_rr_dot(1,dbl_mb(rgx(1)),dn(1,ispin),cdx2) 446 call D3dB_rr_dot(1,dbl_mb(rgy(1)),dn(1,ispin),cdy2) 447 call D3dB_rr_dot(1,dbl_mb(rgz(1)),dn(1,ispin),cdz2) 448 cdx2 = cdx2*dv 449 cdy2 = cdy2*dv 450 cdz2 = cdz2*dv 451 else 452 cdx2 = 0.0d0 453 cdy2 = 0.0d0 454 cdz2 = 0.0d0 455 end if 456 457 cdx3=cdx1+cdx2 458 cdy3=cdy1+cdy2 459 cdz3=cdz1+cdz2 460 461 call lattice_mask_sym(dbl_mb(rgrid(1))) 462 !cdx1=cdx1/ne(1) 463 !cdy1=cdy1/ne(1) 464 !cdz1=cdz1/ne(1) 465 call D3dB_rr_dot(1,dbl_mb(rgrid(1)),dn(1,1),rmax) 466 rmax = rmax*dv 467 cdx1=cdx1/rmax 468 cdy1=cdy1/rmax 469 cdz1=cdz1/rmax 470 if (ne(ispin).ne.0) then 471 !cdx2=cdx2/ne(ispin) 472 !cdy2=cdy2/ne(ispin) 473 !cdz2=cdz2/ne(ispin) 474 call D3dB_rr_dot(1,dbl_mb(rgrid(1)),dn(1,ispin),rmax) 475 rmax = rmax*dv 476 cdx2=cdx2/rmax 477 cdy2=cdy2/rmax 478 cdz2=cdz2/rmax 479 480 end if 481 !cdx3=cdx3/dble(ne(1)+ne(ispin)) 482 !cdy3=cdy3/dble(ne(1)+ne(ispin)) 483 !cdz3=cdz3/dble(ne(1)+ne(ispin)) 484 call D3dB_rr_Sum(1,dn(1,1), 485 > dn(1,ispin), 486 > dbl_mb(rgrid(1)+n2ft3d)) 487 call D3dB_rr_dot(1,dbl_mb(rgrid(1)), 488 > dbl_mb(rgrid(1)+n2ft3d), 489 > rmax) 490 rmax = rmax*dv 491 cdx3=cdx3/rmax 492 cdy3=cdy3/rmax 493 cdz3=cdz3/rmax 494 495 496* ***** calculate dipole with respect to center of mass **** 497 pcharge = tcharge 498 ncharge = dble(ne(1)+ne(ispin)) 499 dipole(1) = -ncharge*cdx3 + pcharge*qGX 500 > - GX*(pcharge-ncharge) 501 dipole(2) = -ncharge*cdy3 + pcharge*qGY 502 > - GY*(pcharge-ncharge) 503 dipole(3) = -ncharge*cdz3 + pcharge*qGZ 504 > - GZ*(pcharge-ncharge) 505 506 507* **** pop stack memory **** 508 value = BA_pop_stack(rgz(2)) 509 value = value.and.BA_pop_stack(rgy(2)) 510 value = value.and.BA_pop_stack(rgx(2)) 511 value = value.and.BA_pop_stack(rgrid(2)) 512 if (.not.value) 513 >call errquit('Calculate_Molecular_Dipole:cannot pop stack memory', 514 > 0,MA_ERR) 515 516 return 517 end 518 519 520* ****************************** 521* * * 522* * Calculate_Resta_Dipole * 523* * * 524* ****************************** 525 526 subroutine Calculate_Resta_Dipole(doprint, 527 > ispin,ne,neq,npack1,nfft3d,psi1, 528 > dipole) 529 implicit none 530 logical doprint 531 integer ispin,ne(2),neq(2) 532 integer npack1,nfft3d 533 complex*16 psi1(npack1,*) 534 real*8 dipole(3) 535 536#include "bafdecls.fh" 537#include "errquit.fh" 538#include "util.fh" 539#include "stdio.fh" 540 541* **** local variables **** 542 integer MASTER,taskid,tmp_len 543 parameter (MASTER=0,tmp_len=140) 544 545 real*8 autoDebye 546 parameter (autoDebye=2.5416d0) 547 548 logical value,oprint 549 integer i,j,ms,n,n2ft3d 550 integer shift,rank 551 integer info 552 integer ii 553 554 integer X(2,6),Xeig(2,6),psi_r(2),psi_r2(2) 555 real*8 bv(3,6),wrk(6,6),wts(6),bmat(3,3) 556 real*8 b(3),ixmat(3,6) 557 real*8 xx,yy,zz,tmp1(tmp_len),maxtime,pcharge 558 real*8 dx(2),dy(2),dz(2),nx,ny,nz,tx,ty,tz,alpha,scal 559 complex*16 arg,wx 560 character*2 labels(6) 561 562* **** external functions **** 563 integer ion_katm,ion_nion 564 external ion_katm,ion_nion 565 real*8 ion_rion,psp_zv 566 external ion_rion,psp_zv 567 real*8 lattice_unitg,lattice_omega 568 external lattice_unitg,lattice_omega 569 logical Dneall_w_push_get,Dneall_w_pop_stack,control_print 570 external Dneall_w_push_get,Dneall_w_pop_stack,control_print 571 572 573 call Parallel_taskid(taskid) 574 oprint= ((taskid.eq.MASTER).and.control_print(print_medium) 575 > .and.doprint) 576 577* ***** allocate X,Y,Z **** 578 n2ft3d = 2*nfft3d 579 n=ne(1) 580 if (n.lt.ne(2)) n=ne(2) 581 value = BA_push_get(mt_dbl,(neq(1)+neq(2))*n2ft3d, 582 > 'psi_r',psi_r(2),psi_r(1)) 583 value = value.and. 584 > BA_push_get(mt_dbl,(neq(1)+neq(2))*n2ft3d, 585 > 'psi_r2',psi_r2(2),psi_r2(1)) 586 do j=1,6 587 value = value.and.Dneall_w_push_get(1,X(1,j)) 588 value = value.and. 589 > BA_push_get(mt_dcpl,n,'Xeig',Xeig(2,j),Xeig(1,j)) 590 end do 591 if (.not. value) 592 > call errquit('dipole_resta:out of stack',1,0) 593 594* **** transform psi1 to realspace **** 595 do n=1,neq(1)+neq(2) 596 call Pack_c_Copy(1,psi1(1,n),dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 597 end do 598 call Grsm_gh_fftb(nfft3d,neq(1)+neq(2),dbl_mb(psi_r(1))) 599 600 601 602c *** Silvestrelli G1 *** 603 ixmat(1,1)=1.0d0 604 ixmat(2,1)=0.0d0 605 ixmat(3,1)=0.0d0 606 607c *** Silvestrelli G4 *** 608 ixmat(1,2)=1.0d0 609 ixmat(2,2)=1.0d0 610 ixmat(3,2)=0.0d0 611 612c *** Silvestrelli G5 *** 613 ixmat(1,3)=1.0d0 614 ixmat(2,3)=0.0d0 615 ixmat(3,3)=1.0d0 616 617c *** Silvestrelli G2 *** 618 ixmat(1,4)=0.0d0 619 ixmat(2,4)=1.0d0 620 ixmat(3,4)=0.0d0 621 622c *** Silvestrelli G6 *** 623 ixmat(1,5)=0.0d0 624 ixmat(2,5)=1.0d0 625 ixmat(3,5)=1.0d0 626 627c *** Silvestrelli G3 *** 628 ixmat(1,6)=0.0d0 629 ixmat(2,6)=0.0d0 630 ixmat(3,6)=1.0d0 631 632 do i=1,3 633 bmat(i,1)=lattice_unitg(1,i) 634 bmat(i,2)=lattice_unitg(2,i) 635 bmat(i,3)=lattice_unitg(3,i) 636 end do 637 638 do i=1,6 639 xx=0.0d0 640 yy=0.0d0 641 zz=0.0d0 642 do j=1,3 643 xx=xx+bmat(j,1)*ixmat(j,i) 644 yy=yy+bmat(j,2)*ixmat(j,i) 645 zz=zz+bmat(j,3)*ixmat(j,i) 646 end do 647 bv(1,i)=xx 648 bv(2,i)=yy 649 bv(3,i)=zz 650 end do 651 652 do i=1,6 653 wrk(1,i)=bv(1,i)*bv(1,i) 654 wrk(2,i)=bv(1,i)*bv(2,i) 655 wrk(3,i)=bv(1,i)*bv(3,i) 656 wrk(4,i)=bv(2,i)*bv(2,i) 657 wrk(5,i)=bv(2,i)*bv(3,i) 658 wrk(6,i)=bv(3,i)*bv(3,i) 659 wts(i)=0.0d0 660 end do 661 662* *** scal=(2*pi/L)**2 *** 663 scal = lattice_omega()**(1.0d0/3.0d0) 664 scal = 8.0*datan(1.0d0)/scal 665 scal = scal*scal 666 wts(1)=scal 667 wts(4)=scal 668 wts(6)=scal 669 call dgels('N',6,6,1,wrk,6,wts,6,tmp1,tmp_len,info) 670 if (info.ne.0) then 671 write(*,*)"Illegal argument in call to dgels" 672 call flush(6) 673 end if 674 rank=0 675 do i=1,6 676 if (dabs(wts(i)).gt.1.e-6) then 677 rank=rank+1 678 wrk(1,rank)=bv(1,i) 679 wrk(2,rank)=bv(2,i) 680 wrk(3,rank)=bv(3,i) 681 wrk(4,rank)=wts(i) 682 end if 683 end do 684 do i=1,rank 685 bv(1,i)=wrk(1,i) 686 bv(2,i)=wrk(2,i) 687 bv(3,i)=wrk(3,i) 688 wts(i)=wrk(4,i) 689 end do 690 691 nx=0.0d0 692 ny=0.0d0 693 nz=0.0d0 694 pcharge = 0.0d0 695 do i=1,ion_nion() 696 j=ion_katm(i) 697 nx=nx+psp_zv(j)*ion_rion(1,i) 698 ny=ny+psp_zv(j)*ion_rion(2,i) 699 nz=nz+psp_zv(j)*ion_rion(3,i) 700 pcharge = pcharge + psp_zv(j) 701 end do 702 703 dx(1)=0.0d0 704 dx(2)=0.0d0 705 dy(1)=0.0d0 706 dy(2)=0.0d0 707 dz(1)=0.0d0 708 dz(2)=0.0d0 709 710 !if (oprint) then 711 ! write(*,1260) 712 ! write(*,1261) rank 713 ! do i=1,rank 714 ! write(*,1262) i,bv(1,i),bv(2,i),bv(3,i),wts(i) 715 ! end do 716 !end if 717 718 do ms=1,ispin 719 do i=1,rank 720 b(1) = bv(1,i) 721 b(2) = bv(2,i) 722 b(3) = bv(3,i) 723 call silvestrelli_overlap( 724 > b,ms,ne,neq, 725 > dbl_mb(psi_r(1)), 726 > dbl_mb(psi_r2(1)), 727 > dcpl_mb(X(1,i))) 728 call Dneall_w_eigenvalues(ms,dcpl_mb(X(1,i)), 729 > dcpl_mb(Xeig(1,i))) 730 end do 731 do i=1,ne(ms) 732 shift=(i-1)*ne(ms)+(i-1) 733 xx=0.0d0 734 yy=0.0d0 735 zz=0.0d0 736 do j=1,rank 737 738 !*** really just want complex eigenvalues of X here *** 739 !arg=dcpl_mb(X(1,j)+shift) 740 arg=dcpl_mb(Xeig(1,j)+(i-1)) 741 arg= -wts(j)*datan2(dimag(arg),dble(arg)) 742 xx=xx+bv(1,j)*arg/scal 743 yy=yy+bv(2,j)*arg/scal 744 zz=zz+bv(3,j)*arg/scal 745 end do 746 dx(ms)=dx(ms)+xx 747 dy(ms)=dy(ms)+yy 748 dz(ms)=dz(ms)+zz 749 end do 750 end do !* ms * 751 752cccccccccccccccccccccccccccccccccccccccccccccc 753c Molecular dipoles from Resta's theory! 754ccccccccccccccccccccccccccccccccccccccccccccc 755 756 tx=nx-dx(1)-dx(ispin) 757 ty=ny-dy(1)-dy(ispin) 758 tz=nz-dz(1)-dz(ispin) 759 xx = dsqrt(tx*tx + ty*ty + tz*tz) 760 dipole(1) = tx 761 dipole(2) = ty 762 dipole(3) = tz 763 764 if (oprint) then 765 write(luout,1771) 766 write(luout,1772) 'spin up ', 767 > dx(1)/dble(ne(1)), 768 > dy(1)/dble(ne(1)), 769 > dz(1)/dble(ne(1)) 770 if (ne(ispin).ne.0) 771 > write(luout,1772) 'spin down ', 772 > dx(ispin)/dble(ne(ispin)), 773 > dy(ispin)/dble(ne(ispin)), 774 > dz(ispin)/dble(ne(ispin)) 775 write(luout,1772) 'electronic', 776 > (dx(1)+dx(ispin))/dble(ne(1)+ne(ispin)), 777 > (dy(1)+dy(ispin))/dble(ne(1)+ne(ispin)), 778 > (dz(1)+dz(ispin))/dble(ne(1)+ne(ispin)) 779 write(luout,1772) 'ionic ', 780 > nx/pcharge, 781 > ny/pcharge, 782 > nz/pcharge 783 write(luout,1778) 784 write(luout,1774) tx,ty,tz 785 write(luout,1775) xx,xx*autoDebye 786 end if 787 788 value = .true. 789 do j=6,1,-1 790 value = value.and.BA_pop_stack(Xeig(2,j)) 791 value = value.and.Dneall_w_pop_stack(X(1,j)) 792 end do 793 value = value.and.BA_pop_stack(psi_r2(2)) 794 value = value.and.BA_pop_stack(psi_r(2)) 795 if (.not. value) 796 > call errquit('dipole_resta:pop stack',2,0) 797 798 799 return 800*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 801 802 1771 FORMAT(//'== Center of Charge =='/) 803 1772 FORMAT(A10,' (',F10.4,',',F10.4,',',F10.4,' )') 804 1773 FORMAT(//'== Wannier Crystal Dipole =='/) 805 1774 FORMAT('mu = (',F10.4,',',F10.4,',',F10.4,' ) au') 806 1775 FORMAT('|mu| = ',F10.4,' au, ',F10.4,' Debye') 807 1776 FORMAT(/"ELECTRONIC DIPOLES") 808 1777 FORMAT("DX =",F11.5," DY= ",F11.5," DZ= ",F11.5) 809 1778 FORMAT(//'== Crystal Dipole (Resta) =='/) 810 1780 FORMAT("NUCLEAR DIPOLES") 811 1785 FORMAT("TOTAL DIPOLES") 812 813 end 814 815 816 817 818