1* 2* $Id$ 3* 4* 5* $Log: not supported by cvs2svn $ 6* Revision 1.27 2006/01/12 00:54:01 bylaska 7* Added charge component analysis wavefunctions, i.e. the following is now outputed in the results section: 8* 9* number of electrons: spin up= 8.00000 down= 6.00000 (real space) 10* plane-wave part: 7.96734 6.02652 (real space) 11* augmented part: 0.03266 -0.02652 (real space) 12* 13* ...EJB 14* 15* Revision 1.26 2005/02/23 21:40:18 edo 16* fixed x1 fpe caused by 0**0 17* 18* Revision 1.25.4.1 2005/02/23 21:39:50 edo 19* fixed x1 fpe caused by 0**0 20* 21* Revision 1.25 2003/10/28 19:50:51 edo 22* errquizzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz 23* 24* Revision 1.24 2003/10/21 02:05:17 marat 25* switched to new errquit by running global replace operation 26* see the script below (note it will not work on multiline errquit calls) 27* ********************************************************* 28* #!/bin/sh 29* 30* e=`find . -name "*F" -print` 31* 32* for f in $e 33* do 34* cp $f $f.bak 35* sed 's|\(^[ ].*call[ ]*errquit([^,]*\)\(,[^,]*\)\()\)|\1,0\2\3|' $f.bak > $f 36* #rm $f.bak 37* done 38* ********************************************************** 39* 40* Revision 1.23 2003/03/21 23:41:13 bylaska 41* 42* paw updates ...EJB 43* 44* Revision 1.22 2003/03/11 17:56:03 bylaska 45* rcut algorithm has been reinstated....EJB 46* 47* Revision 1.21 2003/03/05 23:16:32 bylaska 48* Commented out write statements and other minor fixes..... 49* self-consistent loop looks like it is working..... 50* ....EJB 51* 52* Revision 1.20 2003/03/05 20:35:04 bylaska 53* bug fix.....(m.eq.(mi-mj)) changed to (m.eq.(mj-mi)) in the paw_pot_mult loop... 54* Eigevalues now agree! 55* 56* ....EJB 57* 58* Revision 1.19 2003/03/04 00:04:04 marat 59* added printouts for atomic potenitials 60* for debug purposes 61* MV 62* 63* Revision 1.18 2003/02/27 02:19:25 bylaska 64* The electron gradient has been added. 65* The eigenvalues are not yet agreeing with the F90 code....EJB 66* 67* Revision 1.17 2003/02/26 20:34:06 marat 68* fixed bug related to calculation of 69* comp_coeff, the i and j loop were incorrectly switched 70* MV 71* 72* Revision 1.16 2003/02/25 00:21:00 bylaska 73* debug write statements commented out...EJB 74* 75* Revision 1.15 2003/02/24 20:59:58 bylaska 76* paw_mult_rcut and paw_mult_ncut functions have been added....EJB 77* 78* Revision 1.14 2003/02/24 20:52:52 bylaska 79* fixed initialization of rcut and ncut.....EJB 80* 81* Revision 1.13 2003/02/23 21:37:07 bylaska 82* routines for calculating atomic multipole energies have been added....EJB 83* 84* Revision 1.12 2003/02/23 20:53:08 bylaska 85* bug fix in find_comp_coeff ....EJB 86* 87* Revision 1.11 2003/02/22 03:10:44 bylaska 88* debugging multipole coefficients...There is currently a bug in 89* find_comp_coeff...EJB 90* 91* Revision 1.10 2003/02/21 22:37:26 bylaska 92* find_comp_coeff subroutine has been added....EJB 93* 94* Revision 1.9 2003/02/21 19:44:22 bylaska 95* Routines for computing the mult_energy_coeff have been added to paw_mult 96* ...EJB 97* 98 99!************************************************** 100! 101! Name: paw_mult_init 102! 103! Purpose: initializes paw_mult 104! 105! Created: 2/16/2003 106!************************************************** 107 subroutine paw_mult_init() 108 implicit none 109 110#include "bafdecls.fh" 111#include "paw_mult_data.fh" 112#include "paw_geom.fh" 113#include "paw_ma.fh" 114#include "paw_proj.fh" 115#include "paw_basis.fh" 116 117 !**** local variables *** 118 logical ok 119 integer i,j,k,nfft3d,npack0,nion,nkatm 120 integer lm_size,lmax,v_mult_size 121 integer ia,ii,jj 122 integer l,m 123 integer Gx(2),Gy(2),Gz(2),Ylm(2) 124 integer paw_pot_mult_size 125 real*8 scal,gg,fourpi,omega,rs,w 126 complex*16 itol,cscalf 127 128 !**** external functions **** 129 integer control_ncut,G_indx,control_version 130 real*8 control_rcut,double_factorial,lattice_omega 131 real*8 lattice_unita 132 external control_ncut,G_indx,control_version 133 external control_rcut,double_factorial,lattice_omega 134 external lattice_unita 135 136 137 fourpi = 16.0d0*datan(1.0d0) 138 omega = lattice_omega() 139 140 !*** allocate paw mult memory from heap *** 141 call D3dB_nfft3d(1,nfft3d) 142 nion = ion_nion() 143 nkatm = ion_nkatm() 144 call Pack_npack(0,npack0) 145 lmax = paw_basis_max_mult_l() 146 lm_size = (paw_basis_max_mult_l()+1)**2 147 148 !*** allocate gk_smooth,gk,and glm *** 149 ok = my_alloc(mt_dbl,npack0,"gk_smooth",gk_smooth) 150 ok = ok.and.my_alloc(mt_dbl,npack0*nkatm,"gk",gk) 151 ok = ok.and.my_alloc(mt_dcpl,npack0*lm_size,"g_lm",g_lm) 152 153 !*** allocate paw mult arrays *** 154 ok = ok.and.my_alloc(mt_int,nion,"i_v_mult",i_v_mult) !** same as i_paw_qlm ** 155 if (.not.ok) 156 > call errquit("paw_mult_init: out of heap memory",0,0) 157 158 v_mult_size = 0 159 do ii=1,nion 160 ia = ion_katm(ii) 161 162 int_mb(i_v_mult(1) + ii - 1) = v_mult_size 163 v_mult_size = v_mult_size 164 > + (paw_basis_mult_l(ia)+1)**2 165 end do 166 ok = my_alloc(mt_dcpl,v_mult_size,"v_mult",v_mult) 167 ok = ok.and. 168 > my_alloc(mt_dcpl,v_mult_size,"comp_coeff",comp_coeff) 169 if (.not.ok) 170 > call errquit("paw_mult_init: out of heap memory",0,1) 171 172 173 !*** allocate self_energy_coeff and mult_energy_coeff arrays *** 174 ok = my_alloc(mt_dbl,(lmax+1)*nkatm, 175 > 'self_energy_coeff',self_energy_coeff) 176 177 lm_size = (lmax+1)**2 178 ok = ok.and. 179 > my_alloc(mt_dcpl,nion*lm_size*nion*lm_size, 180 > 'mult_energy_coeff',mult_energy_coeff) 181 if (.not.ok) 182 > call errquit("paw_mult_init: out of heap memory",0,2) 183 184 !*** allocate multiple potential nonlocal matrix *** 185 ok = my_alloc(mt_int,nion,"i_paw_pot_mult", 186 > i_paw_pot_mult) 187 if (.not.ok) 188 > call errquit('init_paw_pot_mult:out of heap memory',0,3) 189 190 paw_pot_mult_size = 0 191 do ii=1,nion 192 int_mb(i_paw_pot_mult(1)+ii-1) = paw_pot_mult_size 193 ia = ion_katm(ii) 194 paw_pot_mult_size = paw_pot_mult_size 195 > + paw_proj_nbasis(ia)**2 196 end do 197 ok = my_alloc(mt_dcpl,paw_pot_mult_size, 198 > "paw_pot_mult",paw_pot_mult) 199 if (.not.ok) 200 > call errquit("paw_mult_init:out of heap memory",0,4) 201 202 203 204 205 206 !**** initialize sigma_smooth and ncut*** 207 sigma_smooth = control_rcut() 208 if (control_version().eq.3) then 209 ncut = control_ncut() 210 else 211 ncut = 0 212 if (sigma_smooth.le.0.0d0) sigma_smooth=1.0d0 213 end if 214 215 if (ncut.lt.0) ncut=0 216 if (sigma_smooth.le.0.0d0) then 217 rs = lattice_unita(1,1)**2 218 > + lattice_unita(2,1)**2 219 > + lattice_unita(3,1)**2 220 rs = dsqrt(rs) 221 sigma_smooth=4.0d0*rs/fourpi 222 223 rs = lattice_unita(1,2)**2 224 > + lattice_unita(2,2)**2 225 > + lattice_unita(3,2)**2 226 rs = dsqrt(rs) 227 w=4.0d0*rs/fourpi 228 if (w.lt.sigma_smooth) sigma_smooth = w 229 230 rs = lattice_unita(1,3)**2 231 > + lattice_unita(2,3)**2 232 > + lattice_unita(3,3)**2 233 rs = dsqrt(rs) 234 w=4.0d0*rs/fourpi 235 if (w.lt.sigma_smooth) sigma_smooth = w 236 end if 237 238 239 240 !**** initialize gk_smooth, gk, and g_lm **** 241 242 !**** allocate stack memory **** 243 ok = BA_push_get(mt_dbl,nfft3d,'Gx',Gx(2),Gx(1)) 244 ok = ok.and. 245 > BA_push_get(mt_dbl,nfft3d,'Gy',Gy(2),Gy(1)) 246 ok = ok.and. 247 > BA_push_get(mt_dbl,nfft3d,'Gz',Gz(2),Gz(1)) 248 ok = ok.and. 249 > BA_push_get(mt_dcpl,npack0,'Ylm',Ylm(2),Ylm(1)) 250 if (.not.ok) 251 > call errquit("paw_mult_init: out of stack memory",0,2) 252 253 call D3dB_t_Copy(1,dbl_mb(G_indx(1)),dbl_mb(Gx(1))) 254 call D3dB_t_Copy(1,dbl_mb(G_indx(2)),dbl_mb(Gy(1))) 255 call D3dB_t_Copy(1,dbl_mb(G_indx(3)),dbl_mb(Gz(1))) 256 call Pack_t_pack(0,dbl_mb(Gx(1))) 257 call Pack_t_pack(0,dbl_mb(Gy(1))) 258 call Pack_t_pack(0,dbl_mb(Gz(1))) 259 260 !**** g_lm **** 261 jj = 0 262 itol = dcmplx(1.0d0,0.0d0) 263 do l=0,paw_basis_max_mult_l() 264 k = 2*l+1 265 cscalf = itol/double_factorial(k) 266 itol = itol*dcmplx(0.0d0,-1.0d0) 267 268 do m = -l,l 269 270 !*** generate Ylm *** 271 call spher_harmonics_generate(l,m,npack0, 272 > dbl_mb(Gx(1)), 273 > dbl_mb(Gy(1)), 274 > dbl_mb(Gz(1)), 275 > dcpl_mb(Ylm(1))) 276 277 !**** (-i)**l*Ylm(k)*|k|**l/(2l+1)!! **** 278 if(l.eq.0) then 279 do k=1,npack0 280 dcpl_mb(g_lm(1)+k-1+(jj)*npack0) 281 > =cscalf 282 > *dcpl_mb(Ylm(1)+k-1) 283 enddo 284 else 285 do k=1,npack0 286 gg = dbl_mb(Gx(1)+k-1)**2 287 > + dbl_mb(Gy(1)+k-1)**2 288 > + dbl_mb(Gz(1)+k-1)**2 289 dcpl_mb(g_lm(1)+k-1+(jj)*npack0) 290 > =cscalf 291 > *dsqrt(gg)**l 292 > *dcpl_mb(Ylm(1)+k-1) 293 > 294 end do !*k* 295 endif 296 jj = jj + 1 297 298 end do !*m* 299 end do !*l* 300 301 !**** gk_smooth and gk **** 302 do k=1,npack0 303 gg = dbl_mb(Gx(1)+k-1)**2 304 > + dbl_mb(Gy(1)+k-1)**2 305 > + dbl_mb(Gz(1)+k-1)**2 306 scal = 0.25d0 * sigma_smooth**2 307 dbl_mb(gk_smooth(1)+k-1) = fourpi*dexp(-scal*gg) 308 > /omega 309 do ia=1,nkatm 310 scal = 0.25d0 * paw_basis_sigma(ia)**2 311 dbl_mb(gk(1)+k-1+(ia-1)*npack0) = fourpi*dexp(-scal*gg) 312 > /omega 313 end do !*ia* 314 end do !*k* 315 316 !**** deallocate stack memory **** 317 ok = BA_pop_stack(Ylm(2)) 318 ok = ok.and.BA_pop_stack(Gz(2)) 319 ok = ok.and.BA_pop_stack(Gy(2)) 320 ok = ok.and.BA_pop_stack(Gx(2)) 321 if (.not.ok) 322 > call errquit('paw_mult_init:error popping stack',0,2) 323 324 325 326* **** allocate rcell memory **** 327 nshl3d=(2*ncut+1)**3 328 ok = my_alloc(mt_dbl,(3*nshl3d),'rcell',rcell) 329 if (.not. ok) call errquit('out of heap memory',0,0) 330 331* **** get lattice vectors in real space **** 332 l=0 333 do k=-ncut,ncut 334 do j=-ncut,ncut 335 do i=-ncut,ncut 336 l = l+1 337 dbl_mb(rcell(1)+3*(l-1) ) 338 > = i*lattice_unita(1,1) 339 > + j*lattice_unita(1,2) 340 > + k*lattice_unita(1,3) 341 dbl_mb(rcell(1)+3*(l-1)+1) 342 > = i*lattice_unita(2,1) 343 > + j*lattice_unita(2,2) 344 > + k*lattice_unita(2,3) 345 dbl_mb(rcell(1)+3*(l-1)+2) 346 > = i*lattice_unita(3,1) 347 > + j*lattice_unita(3,2) 348 > + k*lattice_unita(3,3) 349 end do 350 end do 351 end do 352 353 354 !*** intitalize self_energy_coeff and mult_energy_coeff **** 355 call find_self_energy_coeff(lmax,nkatm, 356 > dbl_mb(self_energy_coeff(1))) 357 call paw_set_mult_energy_coeff() !*needs to be recalled when geometry changes* 358 359 360 return 361 end 362 363 364************************************************* 365! 366! Name: paw_mult_end 367! 368! Purpose: deallocates heap memory 369! 370! Created: 2/16/2003 371!************************************************** 372 subroutine paw_mult_end() 373 implicit none 374 375#include "paw_mult_data.fh" 376#include "paw_ma.fh" 377 378 !**** local varables **** 379 logical ok 380 381 ok = my_dealloc(i_v_mult) 382 ok = ok.and.my_dealloc(v_mult) 383 ok = ok.and.my_dealloc(comp_coeff) 384 ok = ok.and.my_dealloc(g_lm) 385 ok = ok.and.my_dealloc(gk) 386 ok = ok.and.my_dealloc(gk_smooth) 387 ok = ok.and.my_dealloc(self_energy_coeff) 388 ok = ok.and.my_dealloc(mult_energy_coeff) 389 ok = ok.and.my_dealloc(rcell) 390 ok = ok.and.my_dealloc(paw_pot_mult) 391 ok = ok.and.my_dealloc(i_paw_pot_mult) 392 if (.not.ok) 393 > call errquit("paw_mult_end: error freeing heap",0,0) 394 395 return 396 end 397 398 399 400 401!************************************************** 402! 403! Name: paw_mult_dn_cmp_get 404! 405! Purpose: returns dn_cmp and dn_cmp_smooth 406! 407! Created: 2/16/2003 408!************************************************** 409 subroutine paw_mult_dn_cmp_get(dn_cmp, 410 > dn_cmp_smooth) 411 implicit none 412 complex*16 dn_cmp(*) 413 complex*16 dn_cmp_smooth(*) 414 415#include "bafdecls.fh" 416#include "paw_mult_data.fh" 417#include "paw_comp_charge_data.fh" 418#include "paw_geom.fh" 419#include "paw_basis.fh" 420 421 !**** local variables *** 422 logical ok 423 integer ia,ii,jj,kk,l,m,nion,npack0,mult_l 424 integer exi(2),tmp(2),QYlm(2) 425 real*8 sum 426 427 !**** allocate stack memory **** 428 call Pack_npack(0,npack0) 429 ok = BA_push_get(mt_dcpl,npack0,'exi',exi(2),exi(1)) 430 ok = ok.and. 431 > BA_push_get(mt_dcpl,npack0,'QYlm',QYlm(2),QYlm(1)) 432 tmp(1) = exi(1) 433 tmp(2) = exi(2) 434 if (.not.ok) 435 > call errquit( 436 > 'paw_mult_dn_cmp_get: out of stack memory',0,0) 437 438 439 call dcopy(2*npack0,0.0d0,0,dn_cmp,1) 440 call dcopy(2*npack0,0.0d0,0,dn_cmp_smooth,1) 441 sum = 0.0d0 442 nion = ion_nion() 443 do ii=1,nion 444 ia = ion_katm(ii) 445 446 !**** Define QYlm **** 447 mult_l = paw_basis_mult_l(ia) 448 jj = int_mb(i_paw_qlm(1)+ii-1) 449 sum = sum + dble(dcpl_mb(paw_qlm(1)+jj)) 450 kk = 0 451 call dcopy(2*npack0,0.0d0,0,dcpl_mb(QYlm(1)),1) 452 do l=0,mult_l 453 do m=-l,l 454 455 call Pack_cc_zaxpy(0, 456 > dcpl_mb(paw_qlm(1)+jj), 457 > dcpl_mb(g_lm(1)+(kk)*npack0), 458 > dcpl_mb(QYlm(1))) 459 460 jj = jj + 1 461 kk = kk + 1 462 end do 463 end do 464 465 !**** Multiply by Structure Factor **** 466 call strfac_pack(0,ii,dcpl_mb(exi(1))) 467c call Pack_cc_Mul(0, 468c > dcpl_mb(exi(1)), 469c > dcpl_mb(QYlm(1)), 470c > dcpl_mb(QYlm(1))) 471 call Pack_cc_Mul2(0,dcpl_mb(exi(1)), 472 > dcpl_mb(QYlm(1))) 473 474 475 !**** add up ncmp_smooth^ii **** 476 call Pack_tc_Mul(0, 477 > dbl_mb(gk_smooth(1)), 478 > dcpl_mb(QYlm(1)), 479 > dcpl_mb(tmp(1))) 480c call Pack_cc_Sum(0, 481c > dcpl_mb(tmp(1)), 482c > dn_cmp_smooth, 483c > dn_cmp_smooth) 484 call Pack_cc_Sum2(0, 485 > dcpl_mb(tmp(1)), 486 > dn_cmp_smooth) 487 488 !**** add up ncmp^ii *** 489 call Pack_tc_Mul(0, 490 > dbl_mb(gk(1)+(ia-1)*npack0), 491 > dcpl_mb(QYlm(1)), 492 > dcpl_mb(tmp(1))) 493c call Pack_cc_Sum(0, 494c > dcpl_mb(tmp(1)), 495c > dn_cmp, 496c > dn_cmp) 497 call Pack_cc_Sum2(0, 498 > dcpl_mb(tmp(1)), 499 > dn_cmp) 500 501 502 end do !*ii* 503 504 505 !**** deallocate stack memory **** 506 ok = BA_pop_stack(QYlm(2)) 507 ok = ok.and.BA_pop_stack(exi(2)) 508 if (.not.ok) 509 > call errquit('paw_mult_dn_cmp_get: error popping stack',0,1) 510 511 return 512 end 513 514 515!************************************************** 516! 517! Name: paw_mult_pw_force 518! 519! Purpose: returns 520! Sum(G) (i*G*vh(G)*ncmp^a(G)+i*G*vcmp(G)*n_cmp_smooth^a(G)) 521! 522! Created: 2/16/2003 523!************************************************** 524 subroutine paw_mult_pw_force(vh,vcmp,fion) 525 implicit none 526 complex*16 vh(*) 527 complex*16 vcmp(*) 528 real*8 fion(3,*) 529 530#include "bafdecls.fh" 531#include "paw_mult_data.fh" 532#include "paw_comp_charge_data.fh" 533#include "paw_geom.fh" 534#include "paw_basis.fh" 535 536 !**** local variables *** 537 logical ok 538 integer i,ia,ii,jj,kk,l,m,nion,npack0,nfft3d,mult_l 539 integer exi(2),QYlm(2),ncmp(2),ncmp_smooth(2) 540 integer xtmp(2),Gx(2),Gy(2),Gz(2),G(3) 541 real*8 fx,fy,fz,omega 542 543 !*** external functions 544 integer G_indx 545 real*8 lattice_omega 546 external G_indx 547 external lattice_omega 548 549 550 omega = lattice_omega() 551 552 !**** allocate stack memory **** 553 call D3dB_nfft3d(1,nfft3d) 554 call Pack_npack(0,npack0) 555 ok = BA_push_get(mt_dcpl,npack0,'exi',exi(2),exi(1)) 556 ok = ok.and. 557 > BA_push_get(mt_dcpl,npack0,'QYlm',QYlm(2),QYlm(1)) 558 ok = ok.and. 559 > BA_push_get(mt_dcpl,npack0,'ncmp',ncmp(2),ncmp(1)) 560 ok = ok.and. 561 > BA_push_get(mt_dcpl,npack0, 562 > 'ncmp_smooth',ncmp_smooth(2),ncmp_smooth(1)) 563 ok = ok.and. 564 > BA_push_get(mt_dbl, npack0,'xtmp',xtmp(2),xtmp(1)) 565 ok = ok.and. 566 > BA_push_get(mt_dbl, nfft3d,'Gx',Gx(2),Gx(1)) 567 ok = ok.and. 568 > BA_push_get(mt_dbl, nfft3d,'Gy',Gy(2),Gy(1)) 569 ok = ok.and. 570 > BA_push_get(mt_dbl, nfft3d,'Gz',Gz(2),Gz(1)) 571 if (.not.ok) 572 > call errquit( 573 > 'paw_mult_pw_force: out of stack memory',0,0) 574 575 !**** define Gx,Gy and Gz in packed space **** 576 G(1) = G_indx(1) 577 G(2) = G_indx(2) 578 G(3) = G_indx(3) 579 call D3dB_t_Copy(1,dbl_mb(G(1)),dbl_mb(Gx(1))) 580 call D3dB_t_Copy(1,dbl_mb(G(2)),dbl_mb(Gy(1))) 581 call D3dB_t_Copy(1,dbl_mb(G(3)),dbl_mb(Gz(1))) 582 call Pack_t_pack(0,dbl_mb(Gx(1))) 583 call Pack_t_pack(0,dbl_mb(Gy(1))) 584 call Pack_t_pack(0,dbl_mb(Gz(1))) 585 586 587 nion = ion_nion() 588 do ii=1,nion 589 ia = ion_katm(ii) 590 591 !**** Define QYlm **** 592 mult_l = paw_basis_mult_l(ia) 593 jj = int_mb(i_paw_qlm(1)+ii-1) 594 kk = 0 595 call dcopy(2*npack0,0.0d0,0,dcpl_mb(QYlm(1)),1) 596 do l=0,mult_l 597 do m=-l,l 598 call Pack_cc_zaxpy(0, 599 > dcpl_mb(paw_qlm(1)+jj), 600 > dcpl_mb(g_lm(1)+(kk)*npack0), 601 > dcpl_mb(QYlm(1))) 602 jj = jj + 1 603 kk = kk + 1 604 end do 605 end do 606 607 !**** Multiply by Structure Factor **** 608 call strfac_pack(0,ii,dcpl_mb(exi(1))) 609c call Pack_cc_Mul(0, 610c > dcpl_mb(exi(1)), 611c > dcpl_mb(QYlm(1)), 612c > dcpl_mb(QYlm(1))) 613 call Pack_cc_Mul2(0, 614 > dcpl_mb(exi(1)), 615 > dcpl_mb(QYlm(1))) 616 617 618 !**** add up ncmp_smooth^ii **** 619 call Pack_tc_Mul(0, 620 > dbl_mb(gk_smooth(1)), 621 > dcpl_mb(QYlm(1)), 622 > dcpl_mb(ncmp_smooth(1))) 623 624 !**** add up ncmp^ii *** 625 call Pack_tc_Mul(0, 626 > dbl_mb(gk(1)+(ia-1)*npack0), 627 > dcpl_mb(QYlm(1)), 628 > dcpl_mb(ncmp(1))) 629 do i=1,npack0 630 dbl_mb(xtmp(1)+i-1) 631 > = dimag(vh(i))* dble(dcpl_mb(ncmp(1)+i-1)) 632 > - dble(vh(i))*dimag(dcpl_mb(ncmp(1)+i-1)) 633 > + dimag(vcmp(i))* dble(dcpl_mb(ncmp_smooth(1)+i-1)) 634 > - dble(vcmp(i))*dimag(dcpl_mb(ncmp_smooth(1)+i-1)) 635 end do 636 637 call Pack_tt_dot(0,dbl_mb(Gx(1)),dbl_mb(xtmp(1)),fx) 638 call Pack_tt_dot(0,dbl_mb(Gy(1)),dbl_mb(xtmp(1)),fy) 639 call Pack_tt_dot(0,dbl_mb(Gz(1)),dbl_mb(xtmp(1)),fz) 640 fion(1,ii) = fion(1,ii) + fx*omega 641 fion(2,ii) = fion(2,ii) + fy*omega 642 fion(3,ii) = fion(3,ii) + fz*omega 643 644 end do !*ii* 645 646 !**** deallocate stack memory **** 647 ok = BA_pop_stack(Gz(2)) 648 ok = ok.and.BA_pop_stack(Gy(2)) 649 ok = ok.and.BA_pop_stack(Gx(2)) 650 ok = ok.and.BA_pop_stack(xtmp(2)) 651 ok = ok.and.BA_pop_stack(ncmp_smooth(2)) 652 ok = ok.and.BA_pop_stack(ncmp(2)) 653 ok = ok.and.BA_pop_stack(QYlm(2)) 654 ok = ok.and.BA_pop_stack(exi(2)) 655 if (.not.ok) 656 > call errquit('paw_mult_pw_force: error popping stack',0,1) 657 658 return 659 end 660 661 662 663 664!************************************************** 665! 666! Name: paw_mult_coeff_set 667! 668! Purpose: 669! 670! Created: 2/16/2003 671!************************************************** 672 subroutine paw_mult_coeff_set(vh,vcmp) 673 implicit none 674 complex*16 vh(*) 675 complex*16 vcmp(*) 676 677#include "bafdecls.fh" 678#include "paw_mult_data.fh" 679#include "paw_comp_charge_data.fh" 680#include "paw_comp_charge_matrix.fh" 681#include "paw_geom.fh" 682#include "paw_basis.fh" 683#include "paw_proj.fh" 684 685 686 !**** local variables *** 687 logical ok 688 integer ia,ii,ja,jj,kk,l,m,nion,npack0,mult_l,mabs 689 integer i,il,li,mi,mult_li,ill,jll 690 integer j,jl,lj,mj,mult_lj 691 integer isgn,lmax,lmax2,indx 692 integer exi(2),tmp1(2),gls(2),gl(2),t_mult(2) 693 real*8 omega 694 complex*16 csum1,csum2 695 696 integer i_mtr,i_mtr0 697 integer nb,nb2 698 integer ilm 699 integer i_cp,i_cp0 700 integer nilm,njlm 701 integer i_qlm,i_qlm0 702 complex*16 tmp_mult_pot 703 704 !**** external functions **** 705 real*8 lattice_omega,gen_gaunt_coeff 706 external lattice_omega,gen_gaunt_coeff 707 708 omega = lattice_omega() 709 710 !**** allocate stack memory **** 711 call Pack_npack(0,npack0) 712 ok = BA_push_get(mt_dcpl,npack0,'exi',exi(2),exi(1)) 713 ok = ok.and. 714 > BA_push_get(mt_dcpl,npack0,'tmp1',tmp1(2),tmp1(1)) 715 ok = ok.and. 716 > BA_push_get(mt_dcpl,npack0,'gl',gl(2),gl(1)) 717 ok = ok.and. 718 > BA_push_get(mt_dcpl,npack0,'gls',gls(2),gls(1)) 719 ok = ok.and. 720 > BA_push_get(mt_dcpl,v_mult(3), 721 > 't_mult',t_mult(2),t_mult(1)) 722 if (.not.ok) 723 > call errquit( 724 > 'paw_mult_coeff_set: out of stack memory',0,0) 725 726 nion = ion_nion() 727 do ii=1,nion 728 call strfac_pack(0,ii,dcpl_mb(exi(1))) 729 ia = ion_katm(ii) 730 731 call Pack_tc_Mul(0, 732 > dbl_mb(gk(1)+(ia-1)*npack0), 733 > dcpl_mb(exi(1)), 734 > dcpl_mb(gl(1))) 735 call Pack_tc_Mul(0, 736 > dbl_mb(gk_smooth(1)), 737 > dcpl_mb(exi(1)), 738 > dcpl_mb(gls(1))) 739 740 mult_l = paw_basis_mult_l(ia) 741 jj = int_mb(i_v_mult(1) + ii - 1) 742 kk = 0 743 do l=0,mult_l 744 do m=-l,l 745 call Pack_cc_Mul(0, 746 > dcpl_mb(gl(1)), 747 > dcpl_mb(g_lm(1)+(kk)*npack0), 748 > dcpl_mb(tmp1(1))) 749 call Pack_cc_izdot(0, 750 > vh, 751 > dcpl_mb(tmp1(1)), 752 > csum1) 753 call Pack_cc_Mul(0, 754 > dcpl_mb(gls(1)), 755 > dcpl_mb(g_lm(1)+(kk)*npack0), 756 > dcpl_mb(tmp1(1))) 757 call Pack_cc_izdot(0, 758 > vcmp, 759 > dcpl_mb(tmp1(1)), 760 > csum2) 761 762 !**** v_mult(l,m,ii) = <g_lm^a|vh> + <tg_lm^a|vcmp> **** 763 dcpl_mb(t_mult(1)+jj) = (csum1 + csum2)*omega 764 765 jj = jj + 1 766 kk = kk + 1 767 end do !*m* 768 end do !*l* 769 end do !*ii* 770 771 772 !**** unscramble multipole coefficients **** 773 !write(*,*) 774 !write(*,*) 775 do ii=1,nion 776 ia = ion_katm(ii) 777 mult_l = paw_basis_mult_l(ia) 778 jj = int_mb(i_v_mult(1) + ii - 1) 779 do l=0,mult_l 780 do m=-l,l 781 mabs = abs(m) 782 if (mod(mabs,2).eq.0) then 783 isgn = 1 784 else 785 isgn = -1 786 end if 787 kk = jj - 2*m 788 dcpl_mb(v_mult(1)+jj) = 0.5d0*(dcpl_mb(t_mult(1)+jj) 789 > + isgn*dconjg(dcpl_mb(t_mult(1)+kk))) 790 791 !write(*,*) "v_mult 1:",ii,l,m,dcpl_mb(v_mult(1)+jj) 792 jj = jj + 1 793 end do 794 end do 795 end do 796 call D3dB_Vector_SumAll(2*v_mult(3),dcpl_mb(v_mult(1))) 797 798 799 800 lmax = paw_basis_max_mult_l() 801 lmax2 = (lmax+1)**2 802 do ii=1,nion 803 ia = ion_katm(ii) 804 mult_li = paw_basis_mult_l(ia) 805 il = int_mb(i_v_mult(1) + ii - 1) 806 ill = 0 807 do li=0,mult_li 808 do mi=-li,li 809 810 do jj=1,nion 811 ja = ion_katm(jj) 812 mult_lj = paw_basis_mult_l(ja) 813 jl = int_mb(i_paw_qlm(1)+jj-1) 814 jll = 0 815 do lj=0,mult_lj 816 do mj=-lj,lj 817 818 indx = (ii-1) 819 > + ill*nion 820 > + (jj-1)*nion*lmax2 821 > + jll*nion*lmax2*nion 822 823 dcpl_mb(v_mult(1)+il) 824 > = dcpl_mb(v_mult(1)+il) 825 > + dcpl_mb(paw_qlm(1)+jl) 826 > *dcpl_mb(mult_energy_coeff(1)+indx) 827 828 829 jl = jl+1 830 jll = jll+1 831 end do !*mj* 832 end do !*lj* 833 end do !*jj* 834 835 il = il + 1 836 ill = ill + 1 837 end do !*mi* 838 end do !*li* 839 end do !*ii* 840 841 842 843 844 845 call find_comp_coeff() 846 847 do ii=1,nion 848 ia = ion_katm(ii) 849 mult_li = paw_basis_mult_l(ia) 850 il = int_mb(i_v_mult(1)+ii-1) 851 do li=0,mult_li 852 indx = li + (ia-1)*(lmax+1) 853 do mi=-li,li 854 855 dcpl_mb(v_mult(1)+il) 856 > = dcpl_mb(v_mult(1)+il) 857 > - dconjg(dcpl_mb(paw_qlm(1)+il)) 858 > *dbl_mb(self_energy_coeff(1)+indx) 859 > + dcpl_mb(comp_coeff(1)+il) 860 861 862 il = il + 1 863 end do !*mi* 864 end do !*li* 865 end do !*ii* 866 867 868 869 870 871 872 !**** multipole potential nonlocal operator **** 873 call dcopy(2*paw_pot_mult(3), 874 > 0.0d0,0, 875 > dcpl_mb(paw_pot_mult(1)),1) 876 877 do ii=1,nion 878 ia = ion_katm(ii) 879 nb = paw_basis_nbasis(ia) 880 nb2 = nb*nb 881 mult_l = paw_basis_mult_l(ia) 882 883 i_mtr0 = int_mb(i_comp_charge_matrix(1) + ia - 1) 884 i_qlm0 = int_mb(i_v_mult(1) + ii - 1) 885 i_cp0 = int_mb(i_paw_pot_mult(1) + ii - 1) 886 887 i_qlm = i_qlm0 888 do l=0,mult_l 889 do m=-l,l 890 891 i_cp = i_cp0 892 do j=1,nb 893 lj = paw_basis_orb_l(j,ia) 894 do mj=-lj,lj 895 896 do i=1,nb 897 li = paw_basis_orb_l(i,ia) 898 do mi=-li,li 899 900 if ( (l.le.(li+lj) ) .and. 901 > (l.ge.abs(li-lj)) .and. 902 > (m.eq.(mj-mi) ) ) then 903 904 905 i_mtr = i_mtr0 + (i-1) + (j-1)*nb + l*nb2 906 907 tmp_mult_pot = dbl_mb(comp_charge_matrix(1)+i_mtr)* 908 > gen_gaunt_coeff(l,m,lj,mj,li,mi)* 909 > dcpl_mb(v_mult(1) + i_qlm) 910 911 dcpl_mb(paw_pot_mult(1)+i_cp) = 912 > dcpl_mb(paw_pot_mult(1)+i_cp)+ 913 > tmp_mult_pot 914 end if 915 916 i_cp = i_cp + 1 917 end do !*mj* 918 end do !*j* 919 920 end do !*mi* 921 end do !*i* 922 923 i_qlm = i_qlm + 1 924 end do !*m* 925 end do !*l* 926 927 end do !*ii* 928 929 930 931 !**** deallocate stack memory **** 932 ok = BA_pop_stack(t_mult(2)) 933 ok = ok.and.BA_pop_stack(gls(2)) 934 ok = ok.and.BA_pop_stack(gl(2)) 935 ok = ok.and.BA_pop_stack(tmp1(2)) 936 ok = ok.and.BA_pop_stack(exi(2)) 937 if (.not.ok) 938 > call errquit('paw_mult_coeff_set: error popping stack',0,1) 939 940 return 941 end 942 943 944! subroutine paw_pot_mult_print() 945! implicit none 946! 947!#include "bafdecls.fh" 948!#include "paw_mult_data.fh" 949!#include "paw_comp_charge_data.fh" 950!#include "paw_comp_charge_matrix.fh" 951!#include "paw_geom.fh" 952!#include "paw_basis.fh" 953!#include "paw_proj.fh" 954! 955! 956! !**** local variables *** 957! logical ok 958! integer ia,ii,ja,jj,kk,l,m,nion,npack0,mult_l,mabs 959! integer i,il,li,mi,mult_li,ill,jll 960! integer j,jl,lj,mj,mult_lj 961! integer isgn,lmax,lmax2,indx 962! integer exi(2),tmp1(2),gls(2),gl(2),t_mult(2) 963! real*8 omega 964! complex*16 csum1,csum2 965! 966! integer i_mtr,i_mtr0 967! integer nb,nb2 968! integer ilm 969! integer i_cp,i_cp0 970! integer nilm,njlm 971! integer i_qlm,i_qlm0 972! complex*16 tmp_mult_pot 973! 974! 975! 976! write(48,*) paw_pot_mult(3) 977! do ii=1,ion_nion() 978! ia = ion_katm(ii) 979! nb = paw_basis_nbasis(ia) 980! nb2 = nb*nb 981! mult_l = paw_basis_mult_l(ia) 982! i_cp0 = int_mb(i_paw_pot_mult(1) + ii - 1) 983! nilm = 0 984! do i=1,nb 985! li = paw_basis_orb_l(i,ia) 986! njlm = 0 987! do j=1,nb 988! lj = paw_basis_orb_l(j,ia) 989! do mi=-li,li 990! do mj=-lj,lj 991! i_cp = i_cp0-1+(njlm+lj+mj+1)+ 992! > (nilm+li+mi)*paw_proj_nbasis(ia) 993! 994! 995! write(48,*) i,mi,j,mj,ii,dcpl_mb(paw_pot_mult(1)+i_cp) 996! 997! 998! end do !mi 999! end do !mj 1000! 1001! njlm = njlm + 2*lj+1 1002! end do !j 1003! nilm = nilm + 2*li+1 1004! end do !i 1005! end do !ii 1006! 1007! 1008! return 1009! end 1010 1011 1012 1013!************************************************** 1014! 1015! Name: paw_mult_rcut 1016! 1017! Purpose: 1018! 1019! Created: 2/16/2003 1020!************************************************** 1021 function paw_mult_rcut() 1022 implicit none 1023 real*8 paw_mult_rcut !*RESULT* 1024 1025#include "paw_mult_data.fh" 1026 1027 paw_mult_rcut = sigma_smooth 1028 return 1029 end 1030 1031 1032!************************************************** 1033! 1034! Name: paw_mult_ncut 1035! 1036! Purpose: 1037! 1038! Created: 2/16/2003 1039!************************************************** 1040 function paw_mult_ncut() 1041 implicit none 1042 integer paw_mult_ncut !*RESULT* 1043 1044#include "paw_mult_data.fh" 1045 1046 paw_mult_ncut = ncut 1047 return 1048 end 1049 1050 1051 1052!************************************************** 1053! 1054! Name: paw_mult_vzero 1055! 1056! Purpose: 1057! 1058! Created: 2/16/2003 1059!************************************************** 1060 subroutine paw_mult_vzero(vzero) 1061 implicit none 1062 real*8 vzero 1063 1064#include "bafdecls.fh" 1065#include "paw_mult_data.fh" 1066#include "paw_comp_charge_data.fh" 1067#include "paw_geom.fh" 1068#include "paw_basis.fh" 1069 1070 !**** local variables **** 1071 integer ia,ii,jj,nion 1072 real*8 fourpi 1073 1074 !**** external functions **** 1075 real*8 lattice_omega 1076 external lattice_omega 1077 1078 fourpi = 16.0d0*datan(1.0d0) 1079 nion = ion_nion() 1080 vzero = 0.0d0 1081 do ii=1,nion 1082 ia = ion_katm(ii) 1083 jj = int_mb(i_paw_qlm(1)+ii-1) 1084 vzero = vzero 1085 > + dble(dcpl_mb(paw_qlm(1)+jj)) 1086 > *(sigma_smooth**2-paw_basis_sigma(ia)**2) 1087 end do 1088 vzero = vzero*fourpi*dsqrt(fourpi)/lattice_omega()/4.0d0 1089 1090 return 1091 end 1092 1093 1094 !************************************************* 1095 ! 1096 ! Name : find_self_energy_coeff 1097 ! 1098 ! Purpose : 1099 ! 1100 ! Created : 1101 !************************************************* 1102 subroutine find_self_energy_coeff(lmax,nkatm,coeff) 1103 implicit none 1104 integer lmax,nkatm 1105 real*8 coeff(lmax+1,nkatm) 1106 1107#include "paw_basis.fh" 1108 1109 !*** local variables *** 1110 integer ia,l,mult_l 1111 real*8 sigma_tmp,twopi 1112 1113 !*** external functions *** 1114 integer paw_double_factorial 1115 external paw_double_factorial 1116 1117 twopi = 8.0d0*datan(1.0d0) 1118 call dcopy((lmax+1)*nkatm,0.0d0,0,coeff,1) 1119 1120 do ia=1,nkatm 1121 sigma_tmp = paw_basis_sigma(ia) 1122 mult_l = paw_basis_mult_l(ia) 1123 do l=0,mult_l 1124 1125 coeff(l+1,ia) = 4.0d0*dsqrt(twopi) 1126 > /(dble((2*l+1)*paw_double_factorial(2*l+1)) 1127 > *sigma_tmp**(2*l+1) ) 1128 end do 1129 end do 1130 1131 return 1132 end 1133 1134 !************************************************* 1135 ! 1136 ! Name : find_comp_coeff 1137 ! 1138 ! Purpose : 1139 ! 1140 ! Created : 1141 !************************************************* 1142 subroutine find_comp_coeff() 1143 implicit none 1144 1145#include "bafdecls.fh" 1146#include "paw_gaunt.fh" 1147#include "paw_geom.fh" 1148#include "paw_ovlp_data.fh" 1149#include "paw_mult_data.fh" 1150#include "paw_matrix_comp_pot.fh" 1151#include "paw_proj.fh" 1152#include "paw_basis.fh" 1153 1154 1155 !*** local variables ** 1156 integer i,j,ii,ia,jj 1157 integer l,m,li,mi,lj,mj,mult_l 1158 integer indx,indx_2,mtrx_ptr,mtrx_ptr2 1159 integer nion,nbasis,basis_nbasis 1160 1161 1162 call dcopy(2*comp_coeff(3),0.0d0,0, 1163 > dcpl_mb(comp_coeff(1)),1) 1164 1165 nion = ion_nion() 1166 do ii=1,nion 1167 ia = ion_katm(ii) 1168 mult_l = paw_basis_mult_l(ia) 1169 basis_nbasis = paw_basis_nbasis(ia) 1170 nbasis = paw_proj_nbasis(ia) 1171 1172 jj = int_mb(i_v_mult(1)+ii-1) 1173 mtrx_ptr = int_mb(i_paw_comp_pot_matrix(1)+ia-1) 1174 mtrx_ptr2 = int_mb(i_paw_ovlp_w(1)+ii-1) 1175 do l=0,mult_l 1176 do m=-l,l 1177 1178 indx_2 = mtrx_ptr2 1179 do i=1,basis_nbasis 1180 li=paw_basis_orb_l(i,ia) 1181 do mi=-li,li 1182 1183 do j=1,basis_nbasis 1184 lj=paw_basis_orb_l(j,ia) 1185 do mj=-lj,lj 1186 1187 !*** check for non-zero gaunt coefficient *** 1188 if ( (l.le.(li+lj)) .and. 1189 > (l.ge.abs(li-lj)).and. 1190 > (m.eq.(mi-mj)) ) then 1191 indx = mtrx_ptr + (j-1) 1192 > + (i-1)*basis_nbasis 1193 > + l*basis_nbasis**2 1194 dcpl_mb(comp_coeff(1)+jj) 1195 > = dcpl_mb(comp_coeff(1)+jj) 1196 > - gen_gaunt_coeff(l,m,li,mi,lj,mj) 1197 > *dbl_mb(paw_comp_pot_matrix(1)+indx) 1198 > *dcpl_mb(paw_ovlp_w(1)+indx_2) 1199 end if !*non-zero gaunt* 1200 1201 indx_2 = indx_2 + 1 1202 end do !*mi* 1203 end do !*i* 1204 1205 end do !*mj* 1206 end do !*j* 1207 1208 jj = jj + 1 1209 end do !*m* 1210 end do !*l* 1211 1212 end do !*ii* 1213 1214 return 1215 end 1216 1217 1218 1219 subroutine paw_mult_pot_ptr(ptr) 1220 implicit none 1221 integer ptr 1222 1223#include "paw_mult_data.fh" 1224 1225 ptr = paw_pot_mult(1) 1226 return 1227 end 1228 1229 1230 1231 1232 1233!************************************************** 1234! 1235! Name: paw_mult_dn_cmp_smooth_spin_get 1236! 1237! Purpose: returns dn_cmp_smooth for each spin w/o core charge 1238! 1239! Created: 1/11/2006 1240!************************************************** 1241 subroutine paw_mult_dn_cmp_smooth_spin_get(ms,dn_cmp_smooth) 1242 implicit none 1243 integer ms 1244 complex*16 dn_cmp_smooth(*) 1245 1246#include "bafdecls.fh" 1247#include "paw_mult_data.fh" 1248#include "paw_comp_charge_data.fh" 1249#include "paw_geom.fh" 1250#include "paw_basis.fh" 1251 1252 !**** local variables *** 1253 logical ok 1254 integer ia,ii,jj,kk,l,m,nion,npack0,mult_l 1255 integer exi(2),tmp(2),QYlm(2) 1256 1257 !**** allocate stack memory **** 1258 call Pack_npack(0,npack0) 1259 ok = BA_push_get(mt_dcpl,npack0,'exi',exi(2),exi(1)) 1260 ok = ok.and. 1261 > BA_push_get(mt_dcpl,npack0,'QYlm',QYlm(2),QYlm(1)) 1262 tmp(1) = exi(1) 1263 tmp(2) = exi(2) 1264 if (.not.ok) 1265 > call errquit( 1266 > 'paw_mult_dn_cmp_smooth_get2: out of stack memory',0,0) 1267 1268 1269 call dcopy(2*npack0,0.0d0,0,dn_cmp_smooth,1) 1270 nion = ion_nion() 1271 do ii=1,nion 1272 ia = ion_katm(ii) 1273 1274 !**** Define QYlm **** 1275 mult_l = paw_basis_mult_l(ia) 1276 jj = int_mb(i_paw_qlm(1)+ii-1) 1277 1278 kk = 0 1279 call dcopy(2*npack0,0.0d0,0,dcpl_mb(QYlm(1)),1) 1280 do l=0,mult_l 1281 do m=-l,l 1282 call Pack_cc_zaxpy(0, 1283 > dcpl_mb(paw_qlm_spin(1,ms)+jj), 1284 > dcpl_mb(g_lm(1)+(kk)*npack0), 1285 > dcpl_mb(QYlm(1))) 1286 1287 !write(*,*) " qlm=",l,m,dcpl_mb(paw_qlm_spin(1,ms)+jj) 1288 jj = jj + 1 1289 kk = kk + 1 1290 end do 1291 end do 1292 1293 !**** Multiply by Structure Factor **** 1294 call strfac_pack(0,ii,dcpl_mb(exi(1))) 1295c call Pack_cc_Mul(0, 1296c > dcpl_mb(exi(1)), 1297c > dcpl_mb(QYlm(1)), 1298c > dcpl_mb(QYlm(1))) 1299 call Pack_cc_Mul2(0, 1300 > dcpl_mb(exi(1)), 1301 > dcpl_mb(QYlm(1))) 1302 1303 1304 !**** add up ncmp_smooth^ii **** 1305 call Pack_tc_Mul(0, 1306 > dbl_mb(gk_smooth(1)), 1307 > dcpl_mb(QYlm(1)), 1308 > dcpl_mb(tmp(1))) 1309c call Pack_cc_Sum(0, 1310c > dcpl_mb(tmp(1)), 1311c > dn_cmp_smooth, 1312c > dn_cmp_smooth) 1313 call Pack_cc_Sum2(0, 1314 > dcpl_mb(tmp(1)), 1315 > dn_cmp_smooth) 1316 1317 1318 end do !*ii* 1319 1320 1321 !**** deallocate stack memory **** 1322 ok = BA_pop_stack(QYlm(2)) 1323 ok = ok.and.BA_pop_stack(exi(2)) 1324 if (.not.ok) 1325 > call errquit( 1326 > 'paw_mult_dn_cmp_smooth_get2:error popping stack',0,1) 1327 1328 return 1329 end 1330 1331