1* 2* $Id$ 3* 4 5 6************************ f_orb orbitals Part ************************ 7 8************************ virtural orbital Part ************************ 9 10 11************************ KS orbital Part ************************ 12 13* *********************************** 14* * * 15* * paw_psi_KS_update * 16* * * 17* *********************************** 18 19* This routine (approximately) diagonalizes the KS matrix. 20* 21 subroutine paw_psi_KS_update(paw_psi_number,precondition,maxerror) 22 implicit none 23 integer paw_psi_number 24 logical precondition 25 real*8 maxerror 26 27#include "bafdecls.fh" 28#include "paw_psi.fh" 29 30* **** local variables **** 31 logical done 32 integer i,j,neall,maxit_orb,maxit_orbs 33 real*8 error,error_out,tim1,tim2,tim,sum 34 35* **** external functions **** 36 integer control_ks_maxit_orb,control_ks_maxit_orbs 37 external control_ks_maxit_orb,control_ks_maxit_orbs 38 39 tim = 0.0d0 40 neall = ne(1)+ne(2) 41 maxit_orb = control_ks_maxit_orb() !*** should be read from rtdb *** 42 maxit_orbs = control_ks_maxit_orbs() !*** should be read from rtdb *** 43 j = 0 44 2 j = j+1 45 error = 0.0d0 46 !do i=neall,1,-1 47 do i=1,neall 48 call current_second(tim1) 49 50 !*** orthogonalize to lower orbitals **** 51! call paw_psi_project_out_f_orb1( 52! > i, 53! > dcpl_mb(psi1(1)+(i-1)*npack1)) 54 55 !*** normalize **** 56 call Pack_cc_dot(1, 57 > dcpl_mb(psi1(1) +(i-1)*npack1), 58 > dcpl_mb(psi1(1) +(i-1)*npack1), 59 > sum) 60 sum = 1.0d0/dsqrt(sum) 61c call Pack_c_SMul(1,sum, 62c > dcpl_mb(psi1(1) +(i-1)*npack1), 63c > dcpl_mb(psi1(1) +(i-1)*npack1)) 64 call Pack_c_SMul1(1,sum,dcpl_mb(psi1(1) +(i-1)*npack1)) 65 66 67 68 call paw_psi_KS_update_orb(paw_psi_number,precondition, 69 > maxit_orb,maxerror, 70 > 0.1d0,i,error_out) 71 call current_second(tim2) 72 tim = tim + (tim2-tim1) 73 error = error+error_out 74 end do 75 error = error/dble(neall) 76 77 done = ((j.gt.maxit_orbs).or.(error.lt.maxerror)) 78 if (.not.done) go to 2 79c write(*,*) "TIME ALL=",tim 80 81 return 82 end 83 84 85* *********************************** 86* * * 87* * paw_psi_KS_update_orb * 88* * * 89* *********************************** 90 91* This routine performs a KS update on orbital i 92* 93 subroutine paw_psi_KS_update_orb(paw_psi_number, 94 > precondition,maxiteration, 95 > maxerror,perror,i, 96 > error_out) 97 implicit none 98#include "errquit.fh" 99 integer paw_psi_number 100 logical precondition 101 integer maxiteration 102 real*8 maxerror,perror 103 integer i 104 real*8 error_out 105 106#include "bafdecls.fh" 107#include "paw_psi.fh" 108 109* **** local variables **** 110 integer MASTER,taskid 111 parameter (MASTER=0) 112 113 logical value,done,oneloop 114 integer it 115 real*8 e0,eold,percent_error,error0,de0,lmbda_r0,lmbda_r1 116 real*8 theta,sigma 117 integer r1(2),t0(2),t(2),g(2) 118 integer paw_psi_ptr 119 120 if (paw_psi_number.eq.1) then 121 paw_psi_ptr=psi1(1) 122 else 123 paw_psi_ptr=psi2(1) 124 end if 125 126 call Parallel_taskid(taskid) 127 128 value = BA_push_get(mt_dcpl,npack1,'t0',t0(2),t0(1)) 129 value = value.and. 130 > BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 131 value = value.and. 132 > BA_push_get(mt_dcpl,npack1,'g',g(2),g(1)) 133 value = value.and. 134 > BA_push_get(mt_dcpl,npack1,'t',t(2),t(1)) 135 if (.not. value) call errquit( 136 > 'paw_psi_KS_update_orb: out of stack memory',0, MA_ERR) 137 138 done = .false. 139 error0 = 0.0d0 140 e0 = 0.0d0 141 theta = -3.14159d0/600.0d0 142 it = 0 143 2 continue 144 145 it = it + 1 146 eold = e0 147 148* *** calculate residual (steepest descent) direction for a single band *** 149 call paw_psi_get_gradient_orb(paw_psi_number,i,dcpl_mb(g(1))) 150 call Pack_cc_dot(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 151 > dcpl_mb(g(1)), 152 > e0) 153 e0 = -e0 154 155 156 done = ((it.gt.maxiteration) 157 > .or. 158 > (dabs(e0-eold).lt.maxerror)) 159 160 if (done) go to 4 161 162* **** preconditioning **** 163 if (precondition) then 164 call ke_Precondition(npack1,1, 165 > dcpl_mb(g(1)), 166 > dcpl_mb(g(1))) 167 168 end if 169 170c call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 171c call Pack_cc_daxpy(1,(e0), 172c > dcpl_mb(paw_psi_ptr+(i-1)*npack1), 173c > dcpl_mb(r1(1))) 174 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 175 call paw_psi_project_out_orb(paw_psi_number,i,dcpl_mb(r1(1))) 176 177 178 179 180* *** determine conjuagate direction *** 181 call Pack_cc_dot(1,dcpl_mb(r1(1)), 182 > dcpl_mb(r1(1)), 183 > lmbda_r1) 184 call Pack_c_Copy(1,dcpl_mb(r1(1)),dcpl_mb(t(1))) 185 186 if (it.gt.1) then 187 call Pack_cc_daxpy(1,(lmbda_r1/lmbda_r0), 188 > dcpl_mb(t0(1)), 189 > dcpl_mb(t(1))) 190 end if 191 lmbda_r0 = lmbda_r1 192 oneloop = .true. 193 3 call Pack_c_Copy(1,dcpl_mb(t(1)),dcpl_mb(t0(1))) 194 195 196 197c!* **** project out psi components from t **** 198c! call paw_psi_project_out_orb(paw_psi_number,i,dcpl_mb(t(1))) 199c! call Pack_cc_dot(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 200c! > dcpl_mb(t(1)), 201c! > de0) 202c! de0 = -de0 203c! call Pack_cc_daxpy(1,(de0), 204c! > dcpl_mb(paw_psi_ptr+(i-1)*npack1), 205c! > dcpl_mb(t(1))) 206 207 208* *** normalize search direction, t **** 209 call Pack_cc_dot(1,dcpl_mb(t(1)), 210 > dcpl_mb(t(1)), 211 > sigma) 212 sigma = dsqrt(sigma) 213 de0 = 1.0d0/sigma 214c call Pack_c_SMul(1,de0,dcpl_mb(t(1)),dcpl_mb(t(1))) 215 call Pack_c_SMul1(1,de0,dcpl_mb(t(1))) 216 217 218* **** compute de0 = <t|g> **** 219 call Pack_cc_dot(1,dcpl_mb(t(1)), 220 > dcpl_mb(g(1)), 221 > de0) 222 223* *** bad direction *** 224 if ((de0.lt.0.0d0).and.oneloop) then 225 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(t(1))) 226 oneloop = .false. 227 go to 3 228 end if 229 230 de0 = -2.0d0*de0 231 call paw_psi_linesearch_update2(paw_psi_number,i, 232 > theta,e0,de0, 233 > dcpl_mb(t(1)), 234 > sigma, 235 > dcpl_mb(t0(1))) 236 237 go to 2 238 239 240* **** release stack memory **** 241 4 value = BA_pop_stack(t(2)) 242 value = value.and.BA_pop_stack(g(2)) 243 value = value.and.BA_pop_stack(r1(2)) 244 value = value.and.BA_pop_stack(t0(2)) 245 if (.not. value) call errquit( 246 > 'paw_psi_KS_update_orb: popping stack memory',1, MA_ERR) 247 248c write(*,*) "iterations=",it," eig=",e0," error=",error_out, 249c > theta 250 error_out = dabs(e0-eold) 251 return 252 end 253 254 255 256* *********************************** 257* * * 258* * paw_psi_linesearch_update * 259* * * 260* *********************************** 261 262* This routine performs a linesearch on orbital i, in the direction t. 263* This routine is needed for a KS minimizer. 264* e0 = <orb|g> 265* de0 = 2*<t|g> 266* 267 subroutine paw_psi_linesearch_update(paw_psi_number,i, 268 > theta,e0,de0,t) 269 implicit none 270 integer paw_psi_number 271 integer i 272 real*8 theta 273 real*8 e0,de0 274 complex*16 t(*) !search direction 275 276#include "bafdecls.fh" 277#include "paw_psi.fh" 278#include "errquit.fh" 279 280* **** local variables **** 281 logical value 282 integer orb(2),g(2),paw_psi_ptr 283 real*8 x,y,pi,e1,de1 284 real*8 theta2,e2,de2 285 286 if (paw_psi_number.eq.1) then 287 paw_psi_ptr=psi1(1) 288 else 289 paw_psi_ptr=psi2(1) 290 end if 291 292 pi = 4.0d0*datan(1.0d0) 293 294* **** allocate stack memory **** 295 value = BA_push_get(mt_dcpl,npack1,'orb', 296 > orb(2),orb(1)) 297 value = value.and. 298 > BA_push_get(mt_dcpl,npack1,'g', 299 > g(2),g(1)) 300 if (.not. value) call errquit( 301 > 'paw_psi_linesearch_update: out of stack memory',0, MA_ERR) 302 303 304 call Pack_c_Copy(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 305 > dcpl_mb(orb(1))) 306 307* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 308 !theta = pi/300.0d0 309 x = cos(theta) 310 y = sin(theta) 311 call Pack_c_SMul(1,x, 312 > dcpl_mb(orb(1)), 313 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 314 call Pack_cc_daxpy(1,y, 315 > t, 316 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 317 318* *** determine theta *** 319 call paw_psi_get_gradient_orb(paw_psi_number,i,dcpl_mb(g(1))) 320 321 call Pack_cc_dot(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 322 > dcpl_mb(g(1)), 323 > e1) 324 e1 = -e1 325 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 326 > /(1.0d0-cos(2*theta)) 327 theta = 0.5d0*datan(0.5d0*de0/x) 328 329c call Pack_cc_dot(1,t, 330c > dcpl_mb(g(1)), 331c > de1) 332c de1 = -2.0d0*de1 333c theta = -de1*(pi/300.0d0)/(de1-de0) 334 335 !write(*,*) "i,theta,e1:",i,theta,e1 336 337 338* **** orb2 = orb*cos(theta) + t*sin(theta) **** 339 x = cos(theta) 340 y = sin(theta) 341 call Pack_c_SMul(1,x, 342 > dcpl_mb(orb(1)), 343 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 344 call Pack_cc_daxpy(1,y, 345 > t, 346 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 347 348* **** update orb2_r and H*orb2 **** 349 !call paw_electron_run_orb(i,dcpl_mb(paw_psi_ptr)) 350c call paw_psi_get_gradient_orb(paw_psi_number,i,dcpl_mb(g(1))) 351c call Pack_cc_dot(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 352c > dcpl_mb(g(1)), 353c > e2) 354c e2 = -e2 355c call Pack_cc_dot(1,t, 356c > dcpl_mb(g(1)), 357c > de2) 358c de2 = -2.0d0*de2 359 360c write(*,*) "i,theta,es:",i,theta,e0,e1,e2 361c write(*,*) 362 363* **** release stack memory **** 364 value = BA_pop_stack(g(2)) 365 value = value.and.BA_pop_stack(orb(2)) 366 if (.not. value) call errquit( 367 > 'paw_psi_linesearch_update: popping stack memory',1, MA_ERR) 368 369 return 370 end 371 372* *********************************** 373* * * 374* * paw_psi_linesearch_update2 * 375* * * 376* *********************************** 377 378* This routine performs a linesearch on orbital i, in the direction t. 379* This routine is needed for a KS minimizer. 380* e0 = <orb|g> 381* de0 = 2*<t|g> 382* 383 subroutine paw_psi_linesearch_update2(paw_psi_number, 384 > i,theta,e0,de0,t, 385 > sigma,tau_t) 386 implicit none 387#include "errquit.fh" 388 integer paw_psi_number 389 integer i 390 real*8 theta 391 real*8 e0,de0 392 complex*16 t(*) !search direction 393 394 real*8 sigma 395 complex*16 tau_t(*) !parallel transported search direction 396 397#include "bafdecls.fh" 398#include "paw_psi.fh" 399 400* **** local variables **** 401 logical value 402 integer orb(2),g(2),paw_psi_ptr 403 real*8 x,y,pi,e1,de1 404 real*8 theta2,e2,de2 405 406 if (paw_psi_number.eq.1) then 407 paw_psi_ptr=psi1(1) 408 else 409 paw_psi_ptr=psi2(1) 410 end if 411 412 pi = 4.0d0*datan(1.0d0) 413 414* **** allocate stack memory **** 415 value = BA_push_get(mt_dcpl,npack1,'orb', 416 > orb(2),orb(1)) 417 value = value.and. 418 > BA_push_get(mt_dcpl,npack1,'g', 419 > g(2),g(1)) 420 if (.not. value) call errquit( 421 > 'paw_psi_linesearch_update: out of stack memory',0,MA_ERR) 422 423 424 call Pack_c_Copy(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 425 > dcpl_mb(orb(1))) 426 427* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 428 !theta = pi/300.0d0 429 x = cos(theta) 430 y = sin(theta) 431 call Pack_c_SMul(1,x, 432 > dcpl_mb(orb(1)), 433 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 434 call Pack_cc_daxpy(1,y, 435 > t, 436 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 437 438* *** determine theta *** 439 call paw_psi_get_gradient_orb(paw_psi_number,i,dcpl_mb(g(1))) 440 441 call Pack_cc_dot(1,dcpl_mb(paw_psi_ptr+(i-1)*npack1), 442 > dcpl_mb(g(1)), 443 > e1) 444 e1 = -e1 445 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 446 > /(1.0d0-cos(2*theta)) 447 theta = 0.5d0*datan(0.5d0*de0/x) 448 449 x = cos(theta) 450 y = sin(theta) 451 452* **** tau_t = (-orb*sin(theta) + t*cos(theta))*sigma **** 453 call Pack_c_SMul(1,(-y), 454 > dcpl_mb(orb(1)), 455 > tau_t) 456 call Pack_cc_daxpy(1,x, 457 > t, 458 > tau_t) 459c call Pack_c_SMul(1,sigma, 460c > tau_t, 461c > tau_t) 462 call Pack_c_SMul1(1,sigma,tau_t) 463 464* **** orb2 = orb*cos(theta) + t*sin(theta) **** 465 call Pack_c_SMul(1,x, 466 > dcpl_mb(orb(1)), 467 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 468 call Pack_cc_daxpy(1,y, 469 > t, 470 > dcpl_mb(paw_psi_ptr+(i-1)*npack1)) 471 472 473* **** release stack memory **** 474 value = BA_pop_stack(g(2)) 475 value = value.and.BA_pop_stack(orb(2)) 476 if (.not. value) call errquit( 477 > 'paw_psi_linesearch_update: popping stack memory',1,MA_ERR) 478 479 return 480 end 481 482 483 484* *************************** 485* * * 486* * paw_psi_set_orb * 487* * * 488* *************************** 489 490* This routine copies an orbital, orb, into the ith psi of psi1. 491* This routine is needed for a KS minimizer. 492* 493 subroutine paw_psi_set_orb(paw_psi_number,i,orb) 494 implicit none 495 integer paw_psi_number 496 integer i 497 complex*16 orb(*) 498 499#include "bafdecls.fh" 500#include "paw_psi.fh" 501 502* **** local variables **** 503 integer index,paw_psi_ptr 504 505 if (paw_psi_number.eq.1) then 506 paw_psi_ptr=psi1(1) 507 else 508 paw_psi_ptr=psi2(1) 509 end if 510 511 index = (i-1)*npack1 512 513 call zcopy(npack1, 514 > orb, 1, 515 > dcpl_mb(paw_psi_ptr+index),1) 516 return 517 end 518 519 520* *************************** 521* * * 522* * paw_psi_get_orb * 523* * * 524* *************************** 525 526* This routine copies the ith psi of psi1 into an orbital, orb. 527* This routine is needed for a KS minimizer. 528* 529 subroutine paw_psi_get_orb(paw_psi_number,i,orb) 530 implicit none 531 integer paw_psi_number 532 integer i 533 complex*16 orb(*) 534 535#include "bafdecls.fh" 536#include "paw_psi.fh" 537 538* **** local variables **** 539 integer index,paw_psi_ptr 540 541 542 if (paw_psi_number.eq.1) then 543 paw_psi_ptr=psi1(1) 544 else 545 paw_psi_ptr=psi2(1) 546 end if 547 548 index = (i-1)*npack1 549 550 call zcopy(npack1, 551 > dcpl_mb(paw_psi_ptr+index), 1, 552 > orb, 1) 553 return 554 end 555 556* *********************************** 557* * * 558* * paw_psi_get_gradient_orb * 559* * * 560* *********************************** 561 562* This routine returns the Hpsi(i). 563* This routine is needed for a KS minimizer. 564* 565 subroutine paw_psi_get_gradient_orb(paw_psi_number,i,Horb) 566 implicit none 567 integer paw_psi_number 568 integer i 569 complex*16 Horb(*) 570 571#include "bafdecls.fh" 572#include "paw_psi.fh" 573 574* **** local variables **** 575 integer paw_psi_ptr 576 577 if (paw_psi_number.eq.1) then 578 paw_psi_ptr=psi1(1) 579 else 580 paw_psi_ptr=psi2(1) 581 end if 582 583 call paw_electron_run_orb(i,dcpl_mb(paw_psi_ptr)) 584 call paw_electron_get_gradient_orb(i,Horb) 585 586 return 587 end 588 589 590* ******************************************* 591* * * 592* * paw_psi_project_out_orb * 593* * * 594* ******************************************* 595* 596* This routine projects out non-orthogonal components of Horb. 597* This routine is needed for a KS minimizer. 598* 599 subroutine paw_psi_project_out_orb(paw_psi_number,i,Horb) 600 implicit none 601#include "errquit.fh" 602 integer paw_psi_number 603 integer i 604 complex*16 Horb(*) 605 606#include "bafdecls.fh" 607#include "paw_psi.fh" 608 609* **** local variables **** 610 logical ok 611 integer ii,n,paw_psi_ptr,np 612 integer x(2) 613 real*8 sum 614 615 call Parallel_np(np) 616 617* **** allocate stack memory **** 618 ok = BA_push_get(mt_dbl,ne(1),'x',x(2),x(1)) 619 if (.not.ok) 620 > call errquit('paw_psi_project_out_orb: out of stack memory',0, 621 & MA_ERR) 622 623 624 if (paw_psi_number.eq.1) then 625 paw_psi_ptr=psi1(1) 626 else 627 paw_psi_ptr=psi2(1) 628 end if 629 630* **** spin up orbital **** 631 if (i.le.ne(1)) then 632 633 ii = i 634! do n=1,(ii) 635! call Pack_cc_dot(1, 636! > dcpl_mb(paw_psi_ptr +(n-1)*npack1), 637! > Horb, 638! > sum) 639! call daxpy(2*npack1, 640! > (-sum), 641! > dcpl_mb(paw_psi_ptr+(n-1)*npack1),1, 642! > Horb,1) 643! end do 644 call Pack_cc_ndot(1,ii, 645 > dcpl_mb(paw_psi_ptr), 646 > Horb, 647 > dbl_mb(x(1))) 648 do n=1,(ii) 649 call daxpy(2*npack1, 650 > (-dbl_mb(x(1)+n-1)), 651 > dcpl_mb(paw_psi_ptr+(n-1)*npack1),1, 652 > Horb,1) 653 end do 654 655 656 657* **** spin down orbital **** 658 else 659 660 661 ii = i - ne(1) 662 do n=(ne(1)+1),(ne(1)+ii) 663 call Pack_cc_dot(1, 664 > dcpl_mb(paw_psi_ptr +(n-1)*npack1), 665 > Horb, 666 > sum) 667 call daxpy(2*npack1, 668 > (-sum), 669 > dcpl_mb(paw_psi_ptr+(n-1)*npack1),1, 670 > Horb,1) 671 end do 672 673 674 end if 675 676* **** release stack memory **** 677 ok = BA_pop_stack(x(2)) 678 if (.not. ok) 679 > call errquit('paw_psi_project_out_orb: poping stack memory',0, 680 & MA_ERR) 681 682 return 683 end 684 685 686 687 688 689 690* *************************** 691* * * 692* * paw_psi_set_density * 693* * * 694* *************************** 695 696* This routine sets the densities and potentials in psi and electron. 697* This routine is needed for a band by band minimizer. 698* 699 subroutine paw_psi_set_density(paw_psi_number,rho) 700 implicit none 701 integer paw_psi_number 702 real*8 rho(*) 703 704 705#include "bafdecls.fh" 706#include "paw_psi.fh" 707 708 709* **** local variables **** 710 integer dng_ptr,rho_ptr 711 712 if (paw_psi_number.eq.1) then 713 rho_ptr = rho1(1) 714 dng_ptr = dng1(1) 715 else 716 rho_ptr = rho2(1) 717 dng_ptr = dng2(1) 718 end if 719 720 721 call dcopy(4*nfft3d, 722 > rho, 1, 723 > dbl_mb(rho_ptr),1) 724 725 call paw_electron_gen_dng(dbl_mb(rho_ptr), 726 > dcpl_mb(dng_ptr)) 727 call paw_electron_gen_scf_potentials(dbl_mb(rho_ptr), 728 > dcpl_mb(dng_ptr)) 729 call paw_electron_gen_vall() 730 return 731 end 732 733 734* *************************** 735* * * 736* * paw_psi_get_density * 737* * * 738* *************************** 739 740* This routine gets the densities in psi. 741* This routine is needed for a band by band minimizer. 742* 743 subroutine paw_psi_get_density(paw_psi_number,rho) 744 implicit none 745 integer paw_psi_number 746 real*8 rho(*) 747 748 749#include "bafdecls.fh" 750#include "paw_psi.fh" 751 752* **** local variables **** 753 integer rho_ptr 754 755 if (paw_psi_number.eq.1) then 756 rho_ptr = rho1(1) 757 else 758 rho_ptr = rho2(1) 759 end if 760 761 call dcopy(4*nfft3d, 762 > dbl_mb(rho_ptr),1, 763 > rho,1) 764 return 765 end 766 767 768* ************************************** 769* * * 770* * paw_psi_gen_density_potentials * 771* * * 772* ************************************** 773 774* This routine sets the densities and potentials in psi and electron. 775* This routine is needed for a band by band minimizer. 776* 777 subroutine paw_psi_gen_density_potentials(paw_psi_number) 778 implicit none 779 integer paw_psi_number 780 781 782#include "bafdecls.fh" 783#include "paw_psi.fh" 784 785* **** local variables **** 786 integer paw_psi_ptr,rho_ptr,dng_ptr 787 788 if (paw_psi_number.eq.1) then 789 paw_psi_ptr = psi1(1) 790 rho_ptr = rho1(1) 791 dng_ptr = dng1(1) 792 else 793 paw_psi_ptr = psi2(1) 794 rho_ptr = rho2(1) 795 dng_ptr = dng2(1) 796 end if 797 798 799 call paw_electron_gen_densities(dcpl_mb(paw_psi_ptr), 800 > dbl_mb(rho_ptr), 801 > dcpl_mb(dng_ptr)) 802 call paw_electron_gen_scf_potentials(dbl_mb(rho_ptr), 803 > dcpl_mb(dng_ptr)) 804 call paw_electron_gen_vall() 805 return 806 end 807 808 809************************ Grasmman orbitals Part ************************ 810 811* *************************** 812* * * 813* * paw_psi_1to2 * 814* * * 815* *************************** 816 subroutine paw_psi_1to2() 817 implicit none 818 819#include "bafdecls.fh" 820#include "paw_psi.fh" 821 822 call zcopy(npack1*(ne(1)+ne(2)), 823 > dcpl_mb(psi1(1)),1, 824 > dcpl_mb(psi2(1)),1) 825 826 return 827 end 828 829 830* *************************** 831* * * 832* * paw_psi_2to1 * 833* * * 834* *************************** 835 subroutine paw_psi_2to1() 836 implicit none 837 838#include "bafdecls.fh" 839#include "paw_psi.fh" 840 841 842 call zcopy(npack1*(ne(1)+ne(2)), 843 > dcpl_mb(psi2(1)),1, 844 > dcpl_mb(psi1(1)),1) 845 846c call OrthoCheck(ispin,ne,dcpl_mb(psi1(1))) 847 return 848 end 849 850 851* *************************** 852* * * 853* * paw_psi_1get_psi * 854* * * 855* *************************** 856 subroutine paw_psi_1get_psi(rpsi) 857 implicit none 858 complex*16 rpsi(*) 859 860#include "bafdecls.fh" 861#include "paw_psi.fh" 862 863 call zcopy(npack1*(ne(1)+ne(2)), 864 > dcpl_mb(psi1(1)),1, 865 > rpsi,1) 866 867 return 868 end 869 870 871* *************************** 872* * * 873* * paw_psi_2get_psi * 874* * * 875* *************************** 876 subroutine paw_psi_2get_psi(rpsi) 877 implicit none 878 complex*16 rpsi(*) 879 880#include "bafdecls.fh" 881#include "paw_psi.fh" 882 883 call zcopy(npack1*(ne(1)+ne(2)), 884 > dcpl_mb(psi2(1)),1, 885 > rpsi,1) 886 887 return 888 end 889 890* *************************** 891* * * 892* * paw_psi_check * 893* * * 894* *************************** 895 subroutine paw_psi_check() 896 implicit none 897 898#include "bafdecls.fh" 899#include "paw_psi.fh" 900 901 902 call OrthoCheck(ispin,ne,dcpl_mb(psi1(1))) 903 return 904 end 905 906 907 908* *************************** 909* * * 910* * paw_rho_2to1 * 911* * * 912* *************************** 913 subroutine paw_rho_2to1() 914 implicit none 915 916#include "bafdecls.fh" 917#include "paw_psi.fh" 918 919 call dcopy(4*nfft3d, 920 > dbl_mb(rho2(1)),1, 921 > dbl_mb(rho1(1)),1) 922 923 return 924 end 925 926* *************************** 927* * * 928* * paw_rho_1to2 * 929* * * 930* *************************** 931 subroutine paw_rho_1to2() 932 implicit none 933 934#include "bafdecls.fh" 935#include "paw_psi.fh" 936 937 938 call dcopy(4*nfft3d, 939 > dbl_mb(rho1(1)),1, 940 > dbl_mb(rho2(1)),1) 941 942 return 943 end 944 945* *************************** 946* * * 947* * paw_dng_2to1 * 948* * * 949* *************************** 950 subroutine paw_dng_2to1() 951 implicit none 952 953#include "bafdecls.fh" 954#include "paw_psi.fh" 955 956 call zcopy(npack0, 957 > dcpl_mb(dng2(1)),1, 958 > dcpl_mb(dng1(1)),1) 959 960 return 961 end 962 963 964 965 966* *********************************** 967* * * 968* * paw_psi_1toelectron * 969* * * 970* *********************************** 971 subroutine paw_psi_1toelectron() 972 implicit none 973 974#include "bafdecls.fh" 975#include "paw_psi.fh" 976 977 978 call paw_electron_run(dcpl_mb(psi1(1)), 979 > dbl_mb(rho1(1)), 980 > dcpl_mb(dng1(1))) 981 982 return 983 end 984 985* *********************************** 986* * * 987* * paw_psi_1energy * 988* * * 989* *********************************** 990 real*8 function paw_psi_1energy() 991 implicit none 992 993#include "bafdecls.fh" 994#include "paw_psi.fh" 995 996* **** external functions **** 997 real*8 paw_electron_energy 998 external paw_electron_energy 999 1000 call paw_electron_run(dcpl_mb(psi1(1)), 1001 > dbl_mb(rho1(1)), 1002 > dcpl_mb(dng1(1))) 1003 paw_psi_1energy = paw_electron_energy(dcpl_mb(psi1(1)), 1004 > dbl_mb(rho1(1)), 1005 > dcpl_mb(dng1(1))) 1006 return 1007 end 1008 1009* *********************************** 1010* * * 1011* * paw_psi_1_noupdate_energy * 1012* * * 1013* *********************************** 1014 real*8 function paw_psi_1_noupdate_energy() 1015 implicit none 1016 1017#include "bafdecls.fh" 1018#include "paw_psi.fh" 1019 1020* **** external functions **** 1021 real*8 paw_electron_energy 1022 external paw_electron_energy 1023 1024 1025 !call electron_gen_Hpaw_psi_k(dcpl_mb(psi1(1))) 1026 paw_psi_1_noupdate_energy = paw_electron_energy(dcpl_mb(psi1(1)), 1027 > dbl_mb(rho1(1)), 1028 > dcpl_mb(dng1(1)) ) 1029 return 1030 end 1031 1032 1033* *********************************** 1034* * * 1035* * paw_psi_2energy * 1036* * * 1037* *********************************** 1038 real*8 function paw_psi_2energy() 1039 implicit none 1040 1041#include "bafdecls.fh" 1042#include "paw_psi.fh" 1043 1044 1045* **** external functions **** 1046 real*8 paw_electron_energy 1047 external paw_electron_energy 1048 1049 call paw_electron_run(dcpl_mb(psi2(1)), 1050 > dbl_mb(rho2(1)), 1051 > dcpl_mb(dng2(1))) 1052 paw_psi_2energy = paw_electron_energy(dcpl_mb(psi2(1)), 1053 > dbl_mb(rho2(1)), 1054 > dcpl_mb(dng2(1))) 1055 1056 return 1057 end 1058 1059 1060 1061* *********************************** 1062* * * 1063* * paw_psi_1eorbit * 1064* * * 1065* *********************************** 1066 real*8 function paw_psi_1eorbit() 1067 implicit none 1068 1069#include "bafdecls.fh" 1070#include "paw_psi.fh" 1071 1072* **** external functions **** 1073 real*8 paw_electron_eorbit 1074 external paw_electron_eorbit 1075 1076 paw_psi_1eorbit = paw_electron_eorbit(dcpl_mb(psi1(1))) 1077 1078 return 1079 end 1080 1081 1082* *********************************** 1083* * * 1084* * paw_psi_1ke * 1085* * * 1086* *********************************** 1087 real*8 function paw_psi_1ke() 1088 implicit none 1089 1090#include "bafdecls.fh" 1091#include "paw_psi.fh" 1092 1093* **** local variables *** 1094 real*8 ave,occ(1) 1095 1096* **** external functions **** 1097 real*8 paw_electron_eorbit 1098 external paw_electron_eorbit 1099 1100 call ke_ave(ispin,ne,dcpl_mb(psi1(1)),ave,.false.,occ) 1101 paw_psi_1ke = ave 1102 1103 return 1104 end 1105 1106 1107* *********************************** 1108* * * 1109* * paw_psi_1vl * 1110* * * 1111* *********************************** 1112 real*8 function paw_psi_1vl() 1113 implicit none 1114 1115#include "bafdecls.fh" 1116#include "paw_psi.fh" 1117 1118 1119* **** external functions **** 1120 real*8 paw_electron_psi_vl_ave 1121 external paw_electron_psi_vl_ave 1122 1123 paw_psi_1vl 1124 > = paw_electron_psi_vl_ave(dcpl_mb(psi1(1)),dbl_mb(rho1(1))) 1125 1126 return 1127 end 1128 1129 1130 1131 1132 1133* *********************************** 1134* * * 1135* * paw_rho_1exc * 1136* * * 1137* *********************************** 1138 real*8 function paw_rho_1exc() 1139 implicit none 1140 1141#include "bafdecls.fh" 1142#include "paw_psi.fh" 1143 1144 1145* **** external functions **** 1146 real*8 paw_electron_exc 1147 external paw_electron_exc 1148 1149 paw_rho_1exc = paw_electron_exc(dbl_mb(rho1(1))) 1150 return 1151 end 1152 1153* *********************************** 1154* * * 1155* * paw_rho_1pxc * 1156* * * 1157* *********************************** 1158 real*8 function paw_rho_1pxc() 1159 implicit none 1160 1161#include "bafdecls.fh" 1162#include "paw_psi.fh" 1163 1164* **** external functions **** 1165 real*8 paw_electron_pxc 1166 external paw_electron_pxc 1167 1168 paw_rho_1pxc = paw_electron_pxc(dbl_mb(rho1(1))) 1169 return 1170 end 1171 1172 1173* *********************************** 1174* * * 1175* * paw_dng_1ehartree * 1176* * * 1177* *********************************** 1178 real*8 function paw_dng_1ehartree() 1179 implicit none 1180 1181#include "bafdecls.fh" 1182#include "paw_psi.fh" 1183 1184* **** external functions **** 1185 real*8 paw_electron_ehartree 1186 external paw_electron_ehartree 1187 1188 paw_dng_1ehartree = paw_electron_ehartree(dcpl_mb(dng1(1))) 1189 return 1190 end 1191 1192 1193 1194* *********************************** 1195* * * 1196* * paw_psi_2toelectron * 1197* * * 1198* *********************************** 1199 subroutine paw_psi_2toelectron() 1200 implicit none 1201 1202#include "bafdecls.fh" 1203#include "paw_psi.fh" 1204 1205 1206 call paw_electron_run(dcpl_mb(psi2(1)), 1207 > dbl_mb(rho2(1)), 1208 > dcpl_mb(dng2(1))) 1209 return 1210 end 1211 1212 1213* *********************************** 1214* * * 1215* * paw_psi_1check_Tangent * 1216* * * 1217* *********************************** 1218* 1219* This routine checks the accuracy of the tangent vector. 1220* MM = Yt*S*H = Yt*S*(I-Y*Yt*S)*G = Yt*S*G - Yt*S*Y*Yt*S*G = Yt*S*G - Yt*S*G == 0 1221 1222* Updated - 5-18-2002 1223* 1224 subroutine paw_psi_1check_Tangent(H) 1225 implicit none 1226 complex*16 H(*) 1227 1228#include "errquit.fh" 1229#include "bafdecls.fh" 1230#include "paw_psi.fh" 1231 1232* **** local variables **** 1233 logical value 1234 integer ms,n,indx,i,j 1235 integer MM(2) 1236 real*8 sum 1237 1238 do ms=1,ispin 1239 n = ne(ms) 1240 if (n.eq.0) go to 101 !*** ferromagnetic check *** 1241 value = BA_push_get(mt_dbl,n*n,'MM',MM(2),MM(1)) 1242 if (.not. value) 1243 > call errquit('out of stack memory in paw_psi_1check_Tangent',0, 1244 & MA_ERR) 1245 1246 indx = (ms-1)*ne(1)*npack1 1247 1248* **** calculate MM = Yt*S*H **** 1249 call paw_overlap_matrix_gen(n,n, 1250 > dcpl_mb(psi1(1)+indx), 1251 > H(1+indx), 1252 > dbl_mb(MM(1))) 1253 1254* **** write out MM matrix **** 1255 sum = 0.0d0 1256 do j=1,n 1257 do i=1,n 1258 sum = sum + dbl_mb(MM(1)+(i-1)+(j-1)*n) 1259 end do 1260 end do 1261 write(*,*) "paw_psi_1check_Tangent:",sum 1262 1263 1264 1265 value = BA_pop_stack(MM(2)) 1266 if (.not. value) 1267 > call errquit( 1268 > 'error popping stack memory in paw_psi_1check_Tangent',0, 1269 > MA_ERR) 1270 1271 101 continue 1272 end do 1273 1274 return 1275 end 1276 1277 1278 1279* *********************************** 1280* * * 1281* * paw_psi_2check_Tangent * 1282* * * 1283* *********************************** 1284* 1285* This routine checks the accuracy of the tangent vector. 1286* MM = Yt*S*H = Yt*S*(I-Y*Yt*S)*G = Yt*S*G - Yt*S*Y*Yt*S*G = Yt*S*G - Yt*S*G == 0 1287 1288* Updated - 5-18-2002 1289* 1290 subroutine paw_psi_2check_Tangent(H) 1291 implicit none 1292 complex*16 H(*) 1293 1294#include "bafdecls.fh" 1295#include "paw_psi.fh" 1296#include "errquit.fh" 1297 1298* **** local variables **** 1299 logical value 1300 integer ms,n,indx,i,j 1301 integer MM(2) 1302 real*8 sum 1303 1304 do ms=1,ispin 1305 n = ne(ms) 1306 if (n.eq.0) go to 101 !*** ferromagnetic check *** 1307 value = BA_push_get(mt_dbl,n*n,'MM',MM(2),MM(1)) 1308 if (.not. value) 1309 > call errquit('out of stack memory in paw_psi_1check_Tangent',0, 1310 > MA_ERR) 1311 1312 indx = (ms-1)*ne(1)*npack1 1313 1314* **** calculate MM = Yt*S*H **** 1315 call paw_overlap_matrix_gen(n,n, 1316 > dcpl_mb(psi2(1)+indx), 1317 > H(1+indx), 1318 > dbl_mb(MM(1))) 1319 1320* **** write out MM matrix **** 1321 sum = 0.0d0 1322 do j=1,n 1323 do i=1,n 1324 sum = sum + dbl_mb(MM(1)+(i-1)+(j-1)*n) 1325 end do 1326 end do 1327 write(*,*) "paw_psi_2check_Tangent:",sum 1328 1329 1330 1331 value = BA_pop_stack(MM(2)) 1332 if (.not. value) 1333 > call errquit( 1334 > 'error popping stack memory in paw_psi_2check_Tangent',0, 1335 & MA_ERR) 1336 1337 101 continue 1338 end do 1339 1340 return 1341 end 1342 1343 1344 1345* *********************************** 1346* * * 1347* * paw_psi_1get_Tgradient * 1348* * * 1349* *********************************** 1350 1351* THpsi = Hpsi - Y*Y^t*S*Hpsi ! used by Grassman minimizers 1352* THpsi = Hpsi - S*Y*Y^t*Hpsi ! used by Grassman minimizers 1353* THpsi = S^(-1)*Hpsi - Y*Y^t*Hpsi ! used by Grassman minimizers 1354* Y(t) = Y*V*Cos(Sigma*t)*Vt + U*Sin(Sigma*t)*Vt 1355* 1356 subroutine paw_psi_1get_Tgradient(THpsi,Eout) 1357 implicit none 1358 complex*16 THpsi(*) 1359 real*8 Eout 1360 1361#include "bafdecls.fh" 1362#include "paw_psi.fh" 1363#include "errquit.fh" 1364 1365 1366* **** local variables **** 1367 logical value 1368 integer tmp1(2) 1369 1370* **** external functions **** 1371 real*8 paw_electron_energy 1372 external paw_electron_energy 1373 1374 1375 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1376 if (.not. value) 1377 > call errquit('out of stack memory in paw_psi_1get_Tradient',0, 1378 & MA_ERR) 1379 1380 1381 call paw_electron_run(dcpl_mb(psi1(1)), 1382 > dbl_mb(rho1(1)), 1383 > dcpl_mb(dng1(1))) 1384 Eout = paw_electron_energy(dcpl_mb(psi1(1)), 1385 > dbl_mb(rho1(1)), 1386 > dcpl_mb(dng1(1))) 1387 call paw_electron_gen_hml_S(dcpl_mb(psi1(1)), 1388 > dbl_mb(tmp1(1))) 1389 call paw_electron_get_Tgradient(dcpl_mb(psi1(1)), 1390 > dbl_mb(tmp1(1)), 1391 > THpsi) 1392 1393 value = BA_pop_stack(tmp1(2)) 1394 if (.not. value) 1395 > call errquit('paw_psi_1get_Tgradient:error popping stack',1, 1396 > MA_ERR) 1397 1398 return 1399 end 1400 1401* *********************************** 1402* * * 1403* * paw_psi_1get_residual * 1404* * * 1405* *********************************** 1406 1407* 1408* R = Hpsi - S*Y*Y^t*Hpsi ! used by Grassman minimizers 1409* 1410 subroutine paw_psi_1get_residual(R,Spsi,Eout) 1411 implicit none 1412 complex*16 R(*) 1413 complex*16 Spsi(*) 1414 real*8 Eout 1415 1416#include "bafdecls.fh" 1417#include "paw_psi.fh" 1418#include "errquit.fh" 1419 1420 1421* **** local variables **** 1422 logical value 1423 integer tmp1(2) 1424 1425* **** external functions **** 1426 real*8 paw_electron_energy 1427 external paw_electron_energy 1428 1429 1430 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1431 if (.not. value) 1432 > call errquit('out of stack memory in paw_psi_1gen_residual',0, 1433 & MA_ERR) 1434 1435 1436 call paw_electron_run(dcpl_mb(psi1(1)), 1437 > dbl_mb(rho1(1)), 1438 > dcpl_mb(dng1(1))) 1439 1440 Eout = paw_electron_energy(dcpl_mb(psi1(1)), 1441 > dbl_mb(rho1(1)), 1442 > dcpl_mb(dng1(1))) 1443 1444 call paw_electron_gen_hml(dcpl_mb(psi1(1)), ! tmp = <psi|Hpsi> 1445 > dbl_mb(tmp1(1))) 1446 1447 call paw_ovlp_S((ne(1)+ne(2)),dcpl_mb(psi1(1)),Spsi) ! Spsi = S*psi1 1448 1449 call paw_electron_get_Tgradient(Spsi,dbl_mb(tmp1(1)),R) ! R = Hpsi - Spsi*tmp 1450 1451 value = BA_pop_stack(tmp1(2)) 1452 if (.not. value) 1453 > call errquit('paw_psi_1get_residual:error popping stack',1, 1454 > MA_ERR) 1455 return 1456 end 1457 1458 1459* *********************************** 1460* * * 1461* * paw_psi_2get_residual * 1462* * * 1463* *********************************** 1464 1465* 1466* R = Hpsi - S*Y*Y^t*Hpsi ! used by Grassman minimizers 1467* 1468 subroutine paw_psi_2get_residual(option,R,Spsi,Eout) 1469 implicit none 1470 integer option 1471 complex*16 R(*) 1472 complex*16 Spsi(*) 1473 real*8 Eout 1474 1475#include "bafdecls.fh" 1476#include "paw_psi.fh" 1477#include "errquit.fh" 1478 1479 1480* **** local variables **** 1481 logical value 1482 integer tmp1(2) 1483 1484* **** external functions **** 1485 real*8 paw_electron_energy 1486 external paw_electron_energy 1487 1488 1489 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1490 if (.not. value) 1491 > call errquit('out of stack memory in paw_psi_2get_residual',0, 1492 & MA_ERR) 1493 1494 1495 if (option.le.1) then 1496 call paw_electron_run(dcpl_mb(psi2(1)), 1497 > dbl_mb(rho2(1)), 1498 > dcpl_mb(dng2(1))) 1499 end if 1500 Eout = paw_electron_energy(dcpl_mb(psi2(1)), 1501 > dbl_mb(rho2(1)), 1502 > dcpl_mb(dng2(1))) 1503 1504 call paw_electron_gen_hml(dcpl_mb(psi2(1)), ! tmp = <psi|Hpsi> 1505 > dbl_mb(tmp1(1))) 1506 1507 call paw_ovlp_S((ne(1)+ne(2)),dcpl_mb(psi2(1)),Spsi) ! Spsi = S*psi1 1508 1509 call paw_electron_get_Tgradient(Spsi,dbl_mb(tmp1(1)),R) ! R = Hpsi - Spsi*tmp 1510 1511 value = BA_pop_stack(tmp1(2)) 1512 if (.not. value) 1513 > call errquit('paw_psi_2get_residual:error popping stack',1, 1514 > MA_ERR) 1515 return 1516 end 1517 1518 1519 1520 1521* *********************************** 1522* * * 1523* * paw_psi_1get_Gradient * 1524* * * 1525* *********************************** 1526 1527* THpsi = Hpsi ! used by Projected Grassman minimizers 1528* 1529 subroutine paw_psi_1get_Gradient(THpsi,Eout) 1530 implicit none 1531 complex*16 THpsi(*) 1532 real*8 Eout 1533 1534#include "bafdecls.fh" 1535#include "paw_psi.fh" 1536 1537 1538* **** local variables **** 1539 1540* **** external functions **** 1541 real*8 paw_electron_energy 1542 external paw_electron_energy 1543 1544 1545 call paw_electron_run(dcpl_mb(psi1(1)), 1546 > dbl_mb(rho1(1)), 1547 > dcpl_mb(dng1(1))) 1548 1549 Eout = paw_electron_energy(dcpl_mb(psi1(1)), 1550 > dbl_mb(rho1(1)), 1551 > dcpl_mb(dng1(1))) 1552 1553 call paw_electron_get_Gradient(THpsi) 1554 1555 return 1556 end 1557 1558 1559* *********************************** 1560* * * 1561* * paw_psi_1gen_Tangent * 1562* * * 1563* *********************************** 1564 1565* T = T - Y*Y^t*S*T ! used by Grassman minimizers 1566* 1567 subroutine paw_psi_1gen_Tangent(T) 1568 implicit none 1569 complex*16 T(*) 1570 1571#include "bafdecls.fh" 1572#include "paw_psi.fh" 1573#include "errquit.fh" 1574 1575* **** local variables **** 1576 logical value 1577 integer tmp1(2),ii,jj 1578 1579 1580 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1581 if (.not. value) 1582 >call errquit('paw_psi_1gen_Tangent:out of stack memory',0,MA_ERR) 1583 1584 call paw_elecpsitang(dcpl_mb(psi1(1)), 1585 > T, 1586 > dbl_mb(tmp1(1))) 1587 call paw_electron_gen_Tangent(dcpl_mb(psi1(1)), 1588 > dbl_mb(tmp1(1)), 1589 > T) 1590 1591 value = BA_pop_stack(tmp1(2)) 1592 if (.not. value) 1593 > call errquit( 1594 > 'error popping stack memory in paw_psi_1gen_Tangent',0,MA_ERR) 1595 1596 return 1597 end 1598 1599 1600 1601 1602 1603* *********************************** 1604* * * 1605* * paw_psi_2get_Tgradient * 1606* * * 1607* *********************************** 1608 subroutine paw_psi_2get_Tgradient(option,THpsi,Eout) 1609 implicit none 1610 integer option 1611 complex*16 THpsi(*) 1612 real*8 Eout 1613 1614#include "errquit.fh" 1615#include "bafdecls.fh" 1616#include "paw_psi.fh" 1617 1618 1619* *** local variables **** 1620 logical value 1621 integer tmp1(2) 1622 1623* **** external functions **** 1624 real*8 paw_electron_energy 1625 external paw_electron_energy 1626 1627 1628 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1629 if (.not. value) 1630 > call errquit('out of stack memory in paw_psi_1get_Tradient',0, 1631 & MA_ERR) 1632 1633 if (option.le.1) then 1634 call paw_electron_run(dcpl_mb(psi2(1)), 1635 > dbl_mb(rho2(1)), 1636 > dcpl_mb(dng2(1))) 1637 end if 1638 1639 Eout = paw_electron_energy(dcpl_mb(psi2(1)), 1640 > dbl_mb(rho2(1)), 1641 > dcpl_mb(dng2(1))) 1642 call paw_electron_gen_hml_S(dcpl_mb(psi2(1)), 1643 > dbl_mb(tmp1(1))) 1644 call paw_electron_get_Tgradient(dcpl_mb(psi2(1)), 1645 > dbl_mb(tmp1(1)), 1646 > THpsi) 1647 1648 value = BA_pop_stack(tmp1(2)) 1649 if (.not. value) 1650 > call errquit( 1651 > 'paw_psi_2get_Tgradient:error popping stack',1,MA_ERR) 1652 1653 return 1654 end 1655 1656 1657* *********************************** 1658* * * 1659* * paw_psi_2get_Gradient * 1660* * * 1661* *********************************** 1662 subroutine paw_psi_2get_Gradient(option,THpsi,Eout) 1663 implicit none 1664 integer option 1665 complex*16 THpsi(*) 1666 real*8 Eout 1667 1668#include "bafdecls.fh" 1669#include "paw_psi.fh" 1670 1671 1672* *** local variables **** 1673 1674* **** external functions **** 1675 real*8 paw_electron_energy 1676 external paw_electron_energy 1677 1678 1679 if (option.le.1) then 1680 call paw_electron_run(dcpl_mb(psi2(1)), 1681 > dbl_mb(rho2(1)), 1682 > dcpl_mb(dng2(1))) 1683 end if 1684 1685 Eout = paw_electron_energy(dcpl_mb(psi2(1)), 1686 > dbl_mb(rho2(1)), 1687 > dcpl_mb(dng2(1))) 1688 1689 call paw_electron_get_Gradient(THpsi) 1690 1691 return 1692 end 1693 1694 1695* *********************************** 1696* * * 1697* * paw_psi_2gen_Tangent * 1698* * * 1699* *********************************** 1700 1701* T = T - Y*Y^t*S*T ! used by Grassman minimizers 1702* 1703 subroutine paw_psi_2gen_Tangent(T) 1704 implicit none 1705 complex*16 T(*) 1706 1707#include "bafdecls.fh" 1708#include "paw_psi.fh" 1709#include "errquit.fh" 1710 1711* **** local variables **** 1712 logical value 1713 integer tmp1(2),ii,jj 1714 1715 1716 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1717 if (.not. value) 1718 >call errquit('paw_psi_2gen_Tangent:out of stack memory',0,MA_ERR) 1719 1720 call paw_elecpsitang(dcpl_mb(psi2(1)), 1721 > T, 1722 > dbl_mb(tmp1(1))) 1723 call paw_electron_gen_Tangent(dcpl_mb(psi2(1)), 1724 > dbl_mb(tmp1(1)), 1725 > T) 1726 1727 value = BA_pop_stack(tmp1(2)) 1728 if (.not. value) 1729 > call errquit( 1730 > 'error popping stack memory in paw_psi_2gen_Tangent',0,MA_ERR) 1731 1732 return 1733 end 1734 1735 1736 1737 1738* *********************************** 1739* * * 1740* * paw_psi_1get_TSgradient * 1741* * * 1742* *********************************** 1743 1744* THpsi = Hpsi - Y*Hpsi^t*Y ! used by Stiefel minimizers 1745* 1746 subroutine paw_psi_1get_TSgradient(THpsi,Eout) 1747 implicit none 1748 complex*16 THpsi(*) 1749 real*8 Eout 1750 1751#include "errquit.fh" 1752#include "bafdecls.fh" 1753#include "paw_psi.fh" 1754 1755* **** local variables **** 1756 logical value 1757 integer tmp1(2) 1758 1759* **** external functions **** 1760 real*8 paw_electron_energy 1761 external paw_electron_energy 1762 1763 1764 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1765 if (.not. value) 1766 > call errquit('paw_psi_1get_TSradient:pushing stack',0, MA_ERR) 1767 1768 1769 call paw_electron_run(dcpl_mb(psi1(1)), 1770 > dbl_mb(rho1(1)), 1771 > dcpl_mb(dng1(1))) 1772 1773 Eout = paw_electron_energy(dcpl_mb(psi1(1)), 1774 > dbl_mb(rho1(1)), 1775 > dcpl_mb(dng1(1))) 1776 1777 1778 call paw_electron_gen_hmlt(dcpl_mb(psi1(1)), 1779 > dbl_mb(tmp1(1))) 1780 call paw_electron_get_Tgradient(dcpl_mb(psi1(1)), 1781 > dbl_mb(tmp1(1)), 1782 > THpsi) 1783 1784 1785 value = BA_pop_stack(tmp1(2)) 1786 if (.not. value) 1787 > call errquit('paw_psi_1get_TSgradient:popping stack',1, MA_ERR) 1788 1789 return 1790 end 1791 1792 1793* *********************************** 1794* * * 1795* * paw_psi_2get_TSgradient * 1796* * * 1797* *********************************** 1798 1799* THpsi = Hpsi - Y*Hpsi^t*Y ! used by Stiefel minimizers 1800* 1801 subroutine paw_psi_2get_TSgradient(option,THpsi,Eout) 1802 implicit none 1803 integer option 1804 complex*16 THpsi(*) 1805 real*8 Eout 1806 1807#include "errquit.fh" 1808#include "bafdecls.fh" 1809#include "paw_psi.fh" 1810 1811* *** local variables **** 1812 logical value 1813 integer tmp1(2) 1814 1815* **** external functions **** 1816 real*8 paw_electron_energy 1817 external paw_electron_energy 1818 1819 1820 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 1821 if (.not. value) 1822 > call errquit('paw_psi_2get_TSgradient:pushing stack',0, MA_ERR) 1823 1824 if (option.le.1) then 1825 call paw_electron_run(dcpl_mb(psi2(1)), 1826 > dbl_mb(rho2(1)), 1827 > dcpl_mb(dng2(1))) 1828 end if 1829 1830 Eout = paw_electron_energy(dcpl_mb(psi2(1)), 1831 > dbl_mb(rho2(1)), 1832 > dcpl_mb(dng2(1))) 1833 1834 call paw_electron_gen_hmlt(dcpl_mb(psi2(1)), 1835 > dbl_mb(tmp1(1))) 1836 call paw_electron_get_Tgradient(dcpl_mb(psi2(1)), 1837 > dbl_mb(tmp1(1)), 1838 > THpsi) 1839 1840 value = BA_pop_stack(tmp1(2)) 1841 if (.not. value) 1842 > call errquit('paw_psi_2get_TSgradient:popping stack',1, MA_ERR) 1843 1844 return 1845 end 1846 1847 1848 1849 1850* *********************************** 1851* * * 1852* * paw_psi_1get_TMgradient * 1853* * * 1854* *********************************** 1855 subroutine paw_psi_1get_TMgradient(THpsi,Eout) 1856 implicit none 1857 complex*16 THpsi(*) 1858 real*8 Eout 1859 1860#include "bafdecls.fh" 1861#include "paw_psi.fh" 1862 1863 1864* **** external functions **** 1865 real*8 paw_electron_energy 1866 external paw_electron_energy 1867 1868 1869 call paw_electron_run(dcpl_mb(psi1(1)), 1870 > dbl_mb(rho1(1)), 1871 > dcpl_mb(dng1(1))) 1872 1873 Eout = paw_electron_energy(dcpl_mb(psi1(1)), 1874 > dbl_mb(rho1(1)), 1875 > dcpl_mb(dng1(1))) 1876 1877 call paw_electron_get_TMgradient(dcpl_mb(psi1(1)), 1878 > THpsi) 1879 1880 return 1881 end 1882 1883 1884 1885* *********************************** 1886* * * 1887* * paw_psi_2get_TMgradient * 1888* * * 1889* *********************************** 1890 subroutine paw_psi_2get_TMgradient(option,THpsi,Eout) 1891 implicit none 1892 integer option 1893 complex*16 THpsi(*) 1894 real*8 Eout 1895 1896#include "bafdecls.fh" 1897#include "paw_psi.fh" 1898 1899 1900* **** external functions **** 1901 real*8 paw_electron_energy 1902 external paw_electron_energy 1903 1904 if (option.le.1) then 1905 call paw_electron_run(dcpl_mb(psi2(1)), 1906 > dbl_mb(rho2(1)), 1907 > dcpl_mb(dng2(1))) 1908 end if 1909 1910 Eout = paw_electron_energy(dcpl_mb(psi2(1)), 1911 > dbl_mb(rho2(1)), 1912 > dcpl_mb(dng2(1))) 1913 1914 call paw_electron_get_TMgradient(dcpl_mb(psi2(1)), 1915 > THpsi) 1916 1917 return 1918 end 1919 1920* *********************************** 1921* * * 1922* * paw_psi_1ke_Precondition * 1923* * * 1924* *********************************** 1925 subroutine paw_psi_1ke_Precondition(Hpsi) 1926 implicit none 1927 complex*16 Hpsi(*) 1928 1929#include "bafdecls.fh" 1930#include "paw_psi.fh" 1931 1932* **** local variables **** 1933 integer neall 1934 1935 neall = ne(1)+ne(2) 1936 call ke_Precondition(npack1,neall, 1937 > dcpl_mb(psi1(1)), 1938 > Hpsi) 1939 return 1940 end 1941 1942 1943 1944* *********************************** 1945* * * 1946* * paw_psi_1geodesic_transport * 1947* * * 1948* *********************************** 1949 subroutine paw_psi_1geodesic_transport(t,H0) 1950 implicit none 1951 real*8 t 1952 complex*16 H0(*) 1953 1954#include "bafdecls.fh" 1955#include "paw_psi.fh" 1956 1957 1958 call paw_geodesic_transport(t,dcpl_mb(psi1(1)),H0) 1959 1960 return 1961 end 1962 1963 1964* *********************************** 1965* * * 1966* * paw_psi_1geodesic_Gtransport * 1967* * * 1968* *********************************** 1969 subroutine paw_psi_1geodesic_Gtransport(t,G0) 1970 implicit none 1971 real*8 t 1972 complex*16 G0(*) 1973 1974#include "bafdecls.fh" 1975#include "paw_psi.fh" 1976 1977 call paw_geodesic_Gtransport(t,dcpl_mb(psi1(1)),G0) 1978 1979 return 1980 end 1981 1982 1983 1984* *********************************** 1985* * * 1986* * paw_psi_geodesic_energy * 1987* * * 1988* *********************************** 1989 real*8 function paw_psi_geodesic_energy(t) 1990 implicit none 1991 real*8 t 1992 1993#include "bafdecls.fh" 1994#include "paw_psi.fh" 1995 1996 1997* **** local variables **** 1998 real*8 e_new 1999 2000* **** external functions **** 2001 real*8 paw_electron_energy,paw_psi_CheckOrtho 2002 external paw_electron_energy,paw_psi_CheckOrtho 2003 2004 2005 call paw_geodesic_get(t,dcpl_mb(psi1(1)), 2006 > dcpl_mb(psi2(1))) 2007 call paw_electron_run(dcpl_mb(psi2(1)), 2008 > dbl_mb(rho2(1)), 2009 > dcpl_mb(dng2(1))) 2010 e_new = paw_electron_energy(dcpl_mb(psi2(1)), 2011 > dbl_mb(rho2(1)), 2012 > dcpl_mb(dng2(1))) 2013 2014c write(*,*) "paw_psi_geodesic_energy:",t,e_new 2015 paw_psi_geodesic_energy = e_new 2016 return 2017 end 2018 2019* *********************************** 2020* * * 2021* * paw_psi_geodesic_denergy * 2022* * * 2023* *********************************** 2024 real*8 function paw_psi_geodesic_denergy(t) 2025 implicit none 2026 real*8 t 2027 2028#include "bafdecls.fh" 2029#include "paw_psi.fh" 2030 2031* **** external functions **** 2032 real*8 paw_electron_eorbit 2033 external paw_electron_eorbit 2034 2035 2036 call paw_geodesic_transport(t,dcpl_mb(psi1(1)), 2037 > dcpl_mb(psi2(1))) 2038 paw_psi_geodesic_denergy 2039 > = 2.0d0*paw_electron_eorbit(dcpl_mb(psi2(1))) 2040 2041 return 2042 end 2043 2044* *********************************** 2045* * * 2046* * paw_psi_geodesic_final * 2047* * * 2048* *********************************** 2049 subroutine paw_psi_geodesic_final(t) 2050 implicit none 2051 real*8 t 2052 2053#include "bafdecls.fh" 2054#include "paw_psi.fh" 2055 2056 integer taskid,MASTER 2057 parameter (MASTER=0) 2058c real*8 sum1,sum2 2059 2060 call Parallel_taskid(taskid) 2061 2062 call paw_geodesic_get(t,dcpl_mb(psi1(1)), 2063 > dcpl_mb(psi2(1))) 2064 return 2065 end 2066 2067 2068 2069 2070 2071* *********************************** 2072* * * 2073* * paw_psito2_sd_update * 2074* * * 2075* *********************************** 2076 subroutine paw_psi1to2_sd_update(dte) 2077 implicit none 2078 real*8 dte 2079 2080#include "bafdecls.fh" 2081#include "paw_psi.fh" 2082#include "frac_occ.fh" 2083#include "errquit.fh" 2084 2085 2086* **** local variables **** 2087 logical value 2088 integer nemax,ierr 2089 integer lmd(2),tmp_L(2) 2090 2091 2092 call paw_electron_run(dcpl_mb(psi1(1)), 2093 > dbl_mb(rho1(1)), 2094 > dcpl_mb(dng1(1))) 2095 2096* **** do a steepest descent step **** 2097 call paw_electron_sd_update(dcpl_mb(psi1(1)), 2098 > dcpl_mb(psi2(1)), 2099 > dte) 2100 2101* **** lagrange multiplier corrections **** 2102 nemax = ne(1)+ne(2) 2103 2104* **** allocate MA local variables **** 2105 value = BA_push_get(mt_dbl,(8*nemax*nemax), 2106 > 'tmp_L',tmp_L(2),tmp_L(1)) 2107 value = value.and. 2108 > BA_push_get(mt_dbl,(2*nemax*nemax), 2109 > 'lmd',lmd(2),lmd(1)) 2110 2111 call paw_ovlp_S(nemax,dcpl_mb(psi1(1)),dcpl_mb(psi1(1))) 2112 call paw_psi_lmbda(ispin,ne,nemax,npack1, 2113 > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dte, 2114 > dbl_mb(lmd(1)), 2115 > dbl_mb(tmp_L(1)),ierr) 2116 2117 2118 value = value.and.BA_pop_stack(lmd(2)) 2119 value = value.and.BA_pop_stack(tmp_L(2)) 2120 if (.not. value) 2121 > call errquit( 2122 > 'psi1to2_sd_update: stack failure', 0, MA_ERR) 2123 return 2124 end 2125 2126 2127* *********************************** 2128* * * 2129* * paw_psi_1force * 2130* * * 2131* *********************************** 2132 subroutine paw_psi_1force(fion) 2133 implicit none 2134 real*8 fion(3,*) 2135 2136#include "bafdecls.fh" 2137#include "paw_psi.fh" 2138 2139 call paw_electron_run(dcpl_mb(psi1(1)), 2140 > dbl_mb(rho1(1)), 2141 > dcpl_mb(dng1(1))) 2142 2143 call paw_electron_force(dcpl_mb(psi1(1)),dcpl_mb(dng1(1)),fion) 2144 return 2145 end 2146 2147 2148 2149* *********************************** 2150* * * 2151* * paw_psi_1Shml * 2152* * * 2153* *********************************** 2154 subroutine paw_psi_1Shml(S0,S0hml) 2155 implicit none 2156 complex*16 S0(*) 2157 complex*16 S0hml(*) 2158 2159#include "bafdecls.fh" 2160#include "paw_psi.fh" 2161 2162 integer ms,n,shift1,shift2 2163 2164 call paw_electron_gen_hml(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 2165 do ms=1,ispin 2166 n = ne(ms) 2167 if (n.le.0) go to 30 2168 shift1 = 1 + (ms-1)*ne(1)*npack1 2169 shift2 = (ms-1)*ne(1)*ne(1) 2170 call dgemm('N','N',2*npack1,n,n, 2171 > (1.0d0), 2172 > S0(shift1), 2*npack1, 2173 > dbl_mb(hml(1)+shift2), n, 2174 > (0.0d0), 2175 > S0hml(shift1), 2*npack1) 2176 30 continue 2177 end do 2178 return 2179 end 2180 2181 2182 2183* *********************************** 2184* * * 2185* * paw_psi_1gen_hml * 2186* * * 2187* *********************************** 2188 subroutine paw_psi_1gen_hml() 2189 implicit none 2190 2191#include "bafdecls.fh" 2192#include "paw_psi.fh" 2193 2194 2195 call paw_electron_gen_hml(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 2196 2197 return 2198 end 2199 2200 2201 2202 2203* *********************************** 2204* * * 2205* * paw_psi_1gen_hml_g * 2206* * * 2207* *********************************** 2208 subroutine paw_psi_1gen_hml_g() 2209 implicit none 2210 2211#include "bafdecls.fh" 2212#include "paw_psi.fh" 2213 2214 2215 call paw_electron_gen_hml_g(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 2216 2217 return 2218 end 2219 2220 2221* *********************************** 2222* * * 2223* * paw_psi_2gen_hml * 2224* * * 2225* *********************************** 2226 subroutine paw_psi_2gen_hml() 2227 implicit none 2228 2229#include "bafdecls.fh" 2230#include "paw_psi.fh" 2231 2232 2233 call paw_electron_gen_hml(dcpl_mb(psi2(1)),dbl_mb(hml(1))) 2234 2235 return 2236 end 2237 2238* *********************************** 2239* * * 2240* * paw_psi_eigenvalue * 2241* * * 2242* *********************************** 2243 real*8 function paw_psi_eigenvalue(ms,i) 2244 implicit none 2245 integer ms 2246 integer i 2247 2248#include "bafdecls.fh" 2249#include "paw_psi.fh" 2250#include "frac_occ.fh" 2251 2252 real*8 sum 2253 2254 sum = dbl_mb(eig(1)+(i-1)+(ms-1)*ne(1)) 2255 paw_psi_eigenvalue = sum 2256 2257 return 2258 end 2259 2260* *********************************** 2261* * * 2262* * paw_psi_virtual * 2263* * * 2264* *********************************** 2265 real*8 function paw_psi_virtual(ms,i) 2266 implicit none 2267 integer ms 2268 integer i 2269 2270#include "bafdecls.fh" 2271#include "paw_psi.fh" 2272 2273 paw_psi_virtual=dbl_mb(eig_excited(1)+(i-1)+(ms-1)*ne_excited(1)) 2274 2275 return 2276 end 2277 2278* *********************************** 2279* * * 2280* * paw_psi_hml * 2281* * * 2282* *********************************** 2283 real*8 function paw_psi_hml(ms,i,j) 2284 implicit none 2285 integer ms 2286 integer i,j 2287 2288#include "bafdecls.fh" 2289#include "paw_psi.fh" 2290 2291 paw_psi_hml = dbl_mb(hml(1)-1 + i 2292 > + (j-1)*ne(ms) 2293 > + (ms-1)*ne(1)*ne(1)) 2294 2295 return 2296 end 2297 2298 2299* *********************************** 2300* * * 2301* * paw_psi_spin_density * 2302* * * 2303* *********************************** 2304 subroutine paw_psi_spin_density(en1,en2) 2305 implicit none 2306 real*8 en1(2),en2(2) 2307 2308#include "bafdecls.fh" 2309#include "errquit.fh" 2310#include "paw_psi.fh" 2311 2312 2313* **** local variables **** 2314 integer ms,nx,ny,nz,n2ft3d,tmp1(2) 2315 real*8 scale,sumall 2316 2317* **** external functions **** 2318 real*8 lattice_omega 2319 external lattice_omega 2320 2321 call D3dB_nfft3d(1,n2ft3d) 2322 n2ft3d = 2*n2ft3d 2323 call D3dB_nx(1,nx) 2324 call D3dB_ny(1,ny) 2325 call D3dB_nz(1,nz) 2326 scale = lattice_omega()/dble(nx*ny*nz) 2327 2328* **** check total number of electrons **** 2329 do ms =1,ispin 2330 call D3dB_r_dsum(1,dbl_mb(rho1(1)+(ms-1)*n2ft3d),sumall) 2331 en1(ms) = sumall*scale 2332 end do 2333 2334 if (.not.BA_push_get(mt_dbl,n2ft3d,'tmp1',tmp1(2),tmp1(1))) 2335 > call errquit( 2336 > 'paw_psi_spin_density: out of stack memory',0,MA_ERR) 2337 2338 call paw_comp_charge_spin_update() 2339 2340 do ms =1,ispin 2341 call paw_mult_dn_cmp_smooth_spin_get(ms,dbl_mb(tmp1(1))) 2342 call Pack_c_unpack(0,dbl_mb(tmp1(1))) 2343 call D3dB_cr_fft3b(1,dbl_mb(tmp1(1))) 2344 call D3dB_r_Zero_Ends(1,dbl_mb(tmp1(1))) 2345 call D3dB_r_dsum(1,dbl_mb(tmp1(1)),sumall) 2346 en2(ms) = sumall*scale 2347 2348 end do 2349 if (.not.BA_pop_stack(tmp1(2))) 2350 > call errquit( 2351 > 'paw_psi_spin_density: popping stack memory',0,MA_ERR) 2352 2353 2354 2355 return 2356 end 2357 2358* *********************************** 2359* * * 2360* * paw_psi_spin2 * 2361* * * 2362* *********************************** 2363 subroutine paw_psi_spin2(Sab) 2364 implicit none 2365 real*8 Sab 2366 2367#include "bafdecls.fh" 2368#include "paw_psi.fh" 2369 2370 call Calculate_psi_spin2(ispin,ne,npack1,dcpl_mb(psi1(1)),Sab) 2371 2372 return 2373 end 2374 2375 2376 2377 2378 2379* *********************************** 2380* * * 2381* * paw_psi_1rotate2 * 2382* * * 2383* *********************************** 2384 subroutine paw_psi_1rotate2() 2385 implicit none 2386 2387#include "bafdecls.fh" 2388#include "paw_psi.fh" 2389 2390* ***** local variables ***** 2391 integer ms,index,i,j,shift1,shift2 2392 2393 2394 do ms=1,ispin 2395 if (ne(ms).le.0) go to 30 2396 shift1 = (ms-1)*ne(1) 2397 shift2 = (ms-1)*ne(1)*ne(1) 2398 2399 call dgemm('N','N',2*npack1,ne(ms),ne(ms), 2400 > (1.0d0), 2401 > dcpl_mb(psi1(1)+shift1*npack1),2*npack1, 2402 > dbl_mb(hml(1)+shift2),ne(ms), 2403 > (0.0d0), 2404 > dcpl_mb(psi2(1)+shift1*npack1),2*npack1) 2405 2406 30 continue 2407 end do 2408 2409 return 2410 end 2411 2412* *********************************** 2413* * * 2414* * paw_psi_2rotate1 * 2415* * * 2416* *********************************** 2417 subroutine paw_psi_2rotate1() 2418 implicit none 2419 2420#include "bafdecls.fh" 2421#include "paw_psi.fh" 2422 2423* ***** local variables ***** 2424 integer ms,index,i,j,shift1,shift2 2425 2426 do ms=1,ispin 2427 if (ne(ms).le.0) go to 30 2428 shift1 = (ms-1)*ne(1) 2429 shift2 = (ms-1)*ne(1)*ne(1) 2430 2431 call dgemm('N','N',2*npack1,ne(ms),ne(ms), 2432 > (1.0d0), 2433 > dcpl_mb(psi2(1)+shift1*npack1),2*npack1, 2434 > dbl_mb(hml(1)+shift2),ne(ms), 2435 > (0.0d0), 2436 > dcpl_mb(psi1(1)+shift1*npack1),2*npack1) 2437 2438 30 continue 2439 end do 2440 2441 return 2442 end 2443 2444 2445* *********************************** 2446* * * 2447* * paw_psi_diagonalize_hml * 2448* * * 2449* *********************************** 2450 subroutine paw_psi_diagonalize_hml() 2451 implicit none 2452#include "errquit.fh" 2453 2454#include "bafdecls.fh" 2455#include "paw_psi.fh" 2456 2457 2458* ***** local variables **** 2459 logical value 2460 integer ms,shift1,shift2,ierr,i,j,indx 2461 integer tmp1(2) 2462 2463 2464 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 2465 if (.not. value) 2466 > call errquit( 2467 > 'out of stack memory in paw_psi_diagonalize_hml',0,MA_ERR) 2468 2469 2470* ***** diagonalize the hamiltonian matrix ***** 2471 call dcopy((ne(1)+ne(2)),0.0d0,0,dbl_mb(eig(1)),1) 2472 do ms=1,ispin 2473 shift1 = (ms-1)*ne(1) 2474 shift2 = (ms-1)*ne(1)*ne(1) 2475 if (ne(ms).le.0) go to 30 2476 2477 call dsyev('V','U',ne(ms), 2478 > dbl_mb(hml(1)+shift2),ne(ms), 2479 > dbl_mb(eig(1)+shift1), 2480 > dbl_mb(tmp1(1)),2*ne(1)*ne(1), 2481 > ierr) 2482 2483 call eigsrt(dbl_mb(eig(1)+shift1), 2484 > dbl_mb(hml(1)+shift2), 2485 > ne(ms),ne(ms)) 2486 2487 30 continue 2488 end do 2489 2490 2491 value = BA_pop_stack(tmp1(2)) 2492 if (.not. value) 2493 > call errquit( 2494 > 'error popping stack in paw_psi_diagonalize_hml',0,MA_ERR) 2495 2496 return 2497 end 2498 2499* ******************************************* 2500* * * 2501* * paw_psi_diagonalize_hml_assending * 2502* * * 2503* ******************************************* 2504 subroutine paw_psidiaghmlasse() 2505 implicit none 2506#include "errquit.fh" 2507 2508#include "bafdecls.fh" 2509#include "paw_psi.fh" 2510 2511 2512* ***** local variables **** 2513 logical value 2514 integer ms,shift1,shift2,ierr 2515 integer tmp1(2) 2516 2517 2518 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 2519 if (.not. value) 2520 > call errquit( 2521 > 'out of stack memory in paw_psi_diagonalize_hml_assending',0, 2522 > MA_ERR) 2523 2524 2525* ***** diagonalize the hamiltonian matrix ***** 2526 call dcopy((ne(1)+ne(2)),0.0d0,0,dbl_mb(eig(1)),1) 2527 do ms=1,ispin 2528 shift1 = (ms-1)*ne(1) 2529 shift2 = (ms-1)*ne(1)*ne(1) 2530 if (ne(ms).le.0) go to 30 2531 2532 call dsyev('V','U',ne(ms), 2533 > dbl_mb(hml(1)+shift2),ne(ms), 2534 > dbl_mb(eig(1)+shift1), 2535 > dbl_mb(tmp1(1)),2*ne(1)*ne(1), 2536 > ierr) 2537 2538 30 continue 2539 end do 2540 2541 2542 value = BA_pop_stack(tmp1(2)) 2543 if (.not. value) 2544 > call errquit( 2545 > 'error popping stack in paw_psi_diagonalize_hml_assending',0, 2546 > MA_ERR) 2547 2548 return 2549 end 2550 2551 2552 2553* *************************** 2554* * * 2555* * paw_psi_error * 2556* * * 2557* *************************** 2558 real*8 function paw_psi_error() 2559 implicit none 2560#include "errquit.fh" 2561 2562#include "bafdecls.fh" 2563#include "paw_psi.fh" 2564 2565* ***** local variables **** 2566 logical value 2567 integer k,n 2568 real*8 error,sum,size 2569 integer tmp1(2) 2570 2571 value = BA_push_get(mt_dcpl,(npack1),'tmp1',tmp1(2),tmp1(1)) 2572 if (.not. value) 2573 > call errquit('out of stack memory in paw_psi_error',0, MA_ERR) 2574 2575 2576 error = 0.0d0 2577 size = dble(ne(1)+ne(2)) 2578 do n=1, (ne(1)+ne(2)) 2579 do k=1,npack1 2580 dcpl_mb(tmp1(1)+k-1) = dcpl_mb(psi2(1)+k-1+(n-1)*npack1) 2581 > - dcpl_mb(psi1(1)+k-1+(n-1)*npack1) 2582 end do 2583c call D3dB_cc_dot(1,dcpl_mb(tmp1(1)),dcpl_mb(tmp1(1)),sum) 2584 call Pack_cc_dot(1,dcpl_mb(tmp1(1)),dcpl_mb(tmp1(1)),sum) 2585 2586 error = error + sum 2587 end do 2588 error = dsqrt(error)/size 2589 2590 value = BA_pop_stack(tmp1(2)) 2591 if (.not. value) 2592 > call errquit( 2593 > 'error popping stack memory in paw_psi_error',0,MA_ERR) 2594 2595 2596 paw_psi_error = error 2597 return 2598 end 2599 2600* *************************** 2601* * * 2602* * paw_rho_error * 2603* * * 2604* *************************** 2605 real*8 function paw_rho_error() 2606 implicit none 2607 2608#include "errquit.fh" 2609#include "bafdecls.fh" 2610#include "paw_psi.fh" 2611 2612* ***** local variables **** 2613 logical value 2614 integer k,nx,ny,nz 2615 real*8 error,scale 2616 integer tmp1(2) 2617 2618* ***** external functions ***** 2619 real*8 lattice_omega 2620 external lattice_omega 2621 2622 value = BA_push_get(mt_dbl,(2*nfft3d),'tmp1',tmp1(2),tmp1(1)) 2623 if (.not. value) 2624 > call errquit('out of stack memory in rho_error',0, MA_ERR) 2625 2626 2627 call D3dB_nx(1,nx) 2628 call D3dB_ny(1,ny) 2629 call D3dB_nz(1,nz) 2630 scale = lattice_omega() 2631 2632 scale = (scale)/dble(nx*ny*nz) 2633* scale = (scale)/dble(nx*ny*nz) 2634* scale = (scale*scale) 2635 2636 do k=1,(2*nfft3d) 2637 dbl_mb(tmp1(1)+k-1) = (dbl_mb(rho2(1)+k-1) 2638 > -dbl_mb(rho1(1)+k-1)) 2639 dbl_mb(tmp1(1)+k-1) = dbl_mb(tmp1(1)+k-1) 2640 > + (dbl_mb(rho2(1)+k-1+(ispin-1)*(2*nfft3d)) 2641 > -dbl_mb(rho1(1)+k-1+(ispin-1)*(2*nfft3d))) 2642 end do 2643 call D3dB_rr_dot(1,dbl_mb(tmp1(1)),dbl_mb(tmp1(1)),error) 2644 error = error*scale 2645* error = dsqrt(error) 2646 2647 value = BA_pop_stack(tmp1(2)) 2648 if (.not. value) 2649 > call errquit('error popping stack memory in rho_error',0, MA_ERR) 2650 2651 2652 paw_rho_error = error 2653 return 2654 end 2655 2656 2657* *************************** 2658* * * 2659* * paw_rho_dipole * 2660* * * 2661* *************************** 2662* 2663* Uses - Calculate_dipole (pspw/lib/psi/dipole.f) 2664* 2665 subroutine paw_rho_dipole(dipole) 2666 implicit none 2667 real*8 dipole(3) 2668 2669#include "bafdecls.fh" 2670#include "paw_psi.fh" 2671 2672 call Calculate_Dipole(ispin,ne,2*nfft3d,dbl_mb(rho1(1)),dipole) 2673 return 2674 end 2675 2676 2677* *************************** 2678* * * 2679* * paw_psi_ispin * 2680* * * 2681* *************************** 2682 integer function paw_psi_ispin() 2683 implicit none 2684 2685#include "bafdecls.fh" 2686#include "paw_psi.fh" 2687 2688 paw_psi_ispin = ispin 2689 return 2690 end 2691 2692 2693* *************************** 2694* * * 2695* * paw_psi_ne * 2696* * * 2697* *************************** 2698 integer function paw_psi_ne(ms) 2699 implicit none 2700 integer ms 2701 2702#include "bafdecls.fh" 2703#include "paw_psi.fh" 2704 2705 paw_psi_ne = ne(ms) 2706 return 2707 end 2708 2709 2710 2711* *************************** 2712* * * 2713* * paw_psi_initialize * 2714* * * 2715* *************************** 2716 2717 logical function paw_psi_initialize() 2718 implicit none 2719 2720 2721#include "bafdecls.fh" 2722#include "errquit.fh" 2723#include "btdb.fh" 2724#include "paw_psi.fh" 2725 2726* **** local variables **** 2727 integer MASTER,taskid 2728 parameter (MASTER=0) 2729 logical value,psi_nogrid 2730 integer nemax,ms,idum 2731 real*8 sum1,sum2,dum(1) 2732 integer hversion,hnfft(3),hispin,hne(2) 2733 real*8 hunita(3,3) 2734 integer rtdb,ind 2735 2736 integer control_rtdb,control_ngrid 2737 external control_rtdb,control_ngrid 2738 character*50 filename 2739 character*50 control_input_psi 2740 external control_input_psi 2741 logical wvfnc_expander,Dneall_m_allocate 2742 external wvfnc_expander,Dneall_m_allocate 2743 real*8 paw_psi_CheckOrtho 2744 external paw_psi_CheckOrtho 2745 2746 2747 call Parallel_taskid(taskid) 2748 2749* ***** get ispin, and ne, and nfft3d **** 2750 call psi_get_ne(ispin,ne) 2751 call Dneall_neq(neq) 2752 call D3dB_nfft3d(1,nfft3d) 2753 call Pack_npack(1,npack1) 2754 call Pack_npack(0,npack0) 2755 nemax = ne(1)+ne(2) 2756 2757* **** allocate memory **** 2758 value = BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 2759 > 'psi2',psi2(2),psi2(1)) 2760 value = value.and. 2761 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 2762 > 'psi1',psi1(2),psi1(1)) 2763 value = value.and. 2764 > BA_alloc_get(mt_dbl,4*nfft3d, 2765 > 'rho1',rho1(2),rho1(1)) 2766 value = value.and. 2767 > BA_alloc_get(mt_dbl,4*nfft3d, 2768 > 'rho2',rho2(2),rho2(1)) 2769 value = value.and. 2770 > BA_alloc_get(mt_dcpl,npack0, 2771 > 'dng1',dng1(2),dng1(1)) 2772 value = value.and. 2773 > BA_alloc_get(mt_dcpl,npack0, 2774 > 'dng2',dng2(2),dng2(1)) 2775c value = value.and. 2776c > BA_alloc_get(mt_dbl,(2*nemax*nemax),'hml',hml(2),hml(1)) 2777 value = value.and.Dneall_m_allocate(0,hml) 2778 2779 value = value.and. 2780 > BA_alloc_get(mt_dbl,(2*nemax),'eig',eig(2),eig(1)) 2781 2782 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 2783 2784* ***** read initial wavefunctions into psi1 **** 2785 rtdb = control_rtdb() 2786 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 2787 > mt_log,1,psi_nogrid)) 2788 > psi_nogrid = .true. 2789 2790 if (psi_nogrid) then 2791 2792 call psi_get_header(hversion,hnfft,hunita,hispin,hne) 2793 2794 if ( (hnfft(1).ne.control_ngrid(1)) .or. 2795 > (hnfft(2).ne.control_ngrid(2)) .or. 2796 > (hnfft(3).ne.control_ngrid(3)) ) then 2797 2798 hnfft(1) = control_ngrid(1) 2799 hnfft(2) = control_ngrid(2) 2800 hnfft(3) = control_ngrid(3) 2801 2802 call ga_sync() 2803 value = btdb_parallel(.false.) 2804 call ga_sync() 2805 if (taskid.eq.MASTER) then 2806 2807 filename = control_input_psi() 2808 2809 ind = index(filename,' ') - 1 2810 if (.not. btdb_cput(rtdb,'xpndr:old_wavefunction_filename', 2811 > 1,filename(1:ind))) 2812 > call errquit( 2813 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 2814 2815 if (.not. btdb_cput(rtdb,'xpndr:new_wavefunction_filename', 2816 > 1,filename(1:ind))) 2817 > call errquit( 2818 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 2819 2820 if (.not. btdb_put(rtdb,'xpndr:ngrid',mt_int,3,hnfft)) 2821 > call errquit( 2822 > 'wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 2823 2824 write(*,*) 2825 write(*,*) "Grid is being converted:" 2826 write(*,*) "------------------------" 2827 write(*,*) 2828 write(*,*) "To turn off automatic grid conversion:" 2829 write(*,*) 2830 write(*,*) "set nwpw:psi_nogrid .false." 2831 write(*,*) 2832 value = wvfnc_expander(rtdb) 2833 2834 end if 2835 call ga_sync() 2836 value = btdb_parallel(.true.) 2837 2838 end if 2839 2840 end if 2841 2842 call psi_read(ispin,ne,dcpl_mb(psi1(1)),idum,dum) 2843 2844c call psi_history_read(ispin,ne, 2845c > dcpl_mb(psi1(1)), 2846c > dcpl_mb(psi2(1))) 2847 2848 call paw_ovlp_coeff_set(dcpl_mb(psi1(1))) 2849 call paw_ovlp_weights_set() 2850 2851 !**** Ortho Check **** 2852 do ms=1,ispin 2853 sum1=paw_psi_CheckOrtho(npack1,ne(ms), 2854 > dcpl_mb(psi1(1)+(ms-1)*ne(1)*npack1)) 2855 2856 if (sum1.gt.1.0d-10) then 2857 call paw_psi_MakeOrtho(npack1,ne(ms), 2858 > dcpl_mb(psi1(1)+(ms-1)*ne(1)*npack1)) 2859 sum2=paw_psi_CheckOrtho(npack1,ne(ms), 2860 > dcpl_mb(psi1(1)+(ms-1)*ne(1)*npack1)) 2861 if (taskid.eq.MASTER) then 2862 if (ms.eq.1) then 2863 write(*,24) sum1,sum2 2864 end if 2865 if (ms.eq.2) then 2866 write(*,25) sum1,sum2 2867 end if 2868 24 format('Gram-Schmidt performed on up spin : (old error=', 2869 > E10.3,' new error=',E10.3,')' ) 2870 25 format('Gram-Schmidt performed on down spin: (old error=', 2871 > E10.3,' new error=',E10.3,')' ) 2872 2873 end if 2874 end if 2875 2876 end do 2877 2878 2879 paw_psi_initialize = value 2880 return 2881 end 2882 2883* *************************** 2884* * * 2885* * paw_psi_tmp_write * 2886* * * 2887* *************************** 2888 subroutine paw_psi_tmp_write() 2889 implicit none 2890 2891#include "bafdecls.fh" 2892#include "paw_psi.fh" 2893 2894 real*8 dum(1) 2895 2896* ***** write psi1 wavefunctions **** 2897 call psi_write(ispin,ne,dcpl_mb(psi1(1)),-1,dum) 2898 2899 return 2900 end 2901 2902 2903 2904* *************************** 2905* * * 2906* * paw_psi_finalize * 2907* * * 2908* *************************** 2909 2910 logical function paw_psi_finalize(wpsi) 2911 implicit none 2912 logical wpsi 2913 2914#include "bafdecls.fh" 2915#include "errquit.fh" 2916#include "paw_psi.fh" 2917 2918 2919* **** local variables **** 2920 logical value 2921 real*8 dum 2922 2923* ***** write psi1 wavefunctions **** 2924 if (wpsi) then 2925 call psi_write(ispin,ne,dcpl_mb(psi1(1)),-1,dum) 2926c call psi_history_write(ispin,ne,dcpl_mb(psi1(1))) 2927 end if 2928 2929 value = BA_free_heap(eig(2)) 2930 value = value.and.BA_free_heap(hml(2)) 2931 value = value.and.BA_free_heap(dng2(2)) 2932 value = value.and.BA_free_heap(dng1(2)) 2933 value = value.and.BA_free_heap(rho2(2)) 2934 value = value.and.BA_free_heap(rho1(2)) 2935 value = value.and.BA_free_heap(psi2(2)) 2936 value = value.and.BA_free_heap(psi1(2)) 2937 if (.not. value) 2938 > call errquit('paw_psi_finalize: error freeing heap',0, MA_ERR) 2939 2940 paw_psi_finalize = value 2941 return 2942 end 2943 2944 2945* *************************** 2946* * * 2947* * paw_psi_ne_excited * 2948* * * 2949* *************************** 2950 integer function paw_psi_ne_excited(ms) 2951 implicit none 2952 integer ms 2953 2954#include "bafdecls.fh" 2955#include "paw_psi.fh" 2956 2957 paw_psi_ne_excited = ne_excited(ms) 2958 return 2959 end 2960 2961* *************************** 2962* * * 2963* * paw_epsi_initialize * 2964* * * 2965* *************************** 2966 2967 logical function paw_epsi_initialize() 2968 implicit none 2969#include "errquit.fh" 2970 2971#include "bafdecls.fh" 2972#include "btdb.fh" 2973#include "paw_psi.fh" 2974 2975* **** local variables **** 2976 integer MASTER,taskid 2977 parameter (MASTER=0) 2978 logical value,psi_nogrid 2979 integer nemax,ispin0 2980 real*8 sum1,sum2 2981 2982 integer hversion,hnfft(3),hispin,hne(2) 2983 real*8 hunita(3,3) 2984 integer rtdb,ind 2985 integer control_rtdb,control_ngrid 2986 external control_rtdb,control_ngrid 2987 character*50 filename 2988 character*50 control_input_epsi 2989 external control_input_epsi 2990 logical wvfnc_expander 2991 external wvfnc_expander 2992 2993 2994* ***** get ispin, and ne, and nfft3d **** 2995 call psi_get_ne_excited(ispin0,ne_excited) 2996 nemax = ne_excited(1)+ne_excited(2) 2997 2998* **** allocate memory **** 2999 value = BA_alloc_get(mt_dcpl,npack1*(nemax), 3000 > 'psi2_excited',psi2_excited(2),psi2_excited(1)) 3001 value = value.and. 3002 > BA_alloc_get(mt_dcpl,npack1*(nemax), 3003 > 'psi1_excited',psi1_excited(2),psi1_excited(1)) 3004 value = value.and. 3005 > BA_alloc_get(mt_dbl,(2*nemax),'eig_excited', 3006 > eig_excited(2),eig_excited(1)) 3007 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 3008 3009 call dcopy(2*npack1*nemax,0.0d0,0,dcpl_mb(psi2_excited(1)),1) 3010 call dcopy(2*npack1*nemax,0.0d0,0,dcpl_mb(psi1_excited(1)),1) 3011 call dcopy(2*nemax,0.0d0,0,dbl_mb(eig_excited(1)),1) 3012 3013 3014* ***** read initial wavefunctions into psi1 **** 3015 rtdb = control_rtdb() 3016 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 3017 > mt_log,1,psi_nogrid)) 3018 > psi_nogrid = .true. 3019 3020 if (psi_nogrid) then 3021 3022 3023 filename = control_input_epsi() 3024 call psi_get_header_filename(filename, 3025 > hversion,hnfft,hunita,hispin,hne) 3026 3027 if ( (hnfft(1).ne.control_ngrid(1)) .or. 3028 > (hnfft(2).ne.control_ngrid(2)) .or. 3029 > (hnfft(3).ne.control_ngrid(3)) ) then 3030 3031 hnfft(1) = control_ngrid(1) 3032 hnfft(2) = control_ngrid(2) 3033 hnfft(3) = control_ngrid(3) 3034 call Parallel_taskid(taskid) 3035 value = btdb_parallel(.false.) 3036 if (taskid.eq.MASTER) then 3037 3038 3039 ind = index(filename,' ') - 1 3040 if (.not. btdb_cput(rtdb,'xpndr:old_wavefunction_filename', 3041 > 1,filename(1:ind))) 3042 > call errquit( 3043 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 3044 3045 if (.not. btdb_cput(rtdb,'xpndr:new_wavefunction_filename', 3046 > 1,filename(1:ind))) 3047 > call errquit( 3048 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 3049 3050 if (.not. btdb_put(rtdb,'xpndr:ngrid',mt_int,3,hnfft)) 3051 > call errquit( 3052 > 'wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 3053 3054 write(*,*) 3055 write(*,*) "Grid is being converted:" 3056 write(*,*) "------------------------" 3057 write(*,*) 3058 write(*,*) "To turn off automatic grid conversion:" 3059 write(*,*) 3060 write(*,*) "set nwpw:psi_nogrid .false." 3061 write(*,*) 3062 value = wvfnc_expander(rtdb) 3063 3064 end if 3065 value = btdb_parallel(.true.) 3066 3067 end if 3068 3069 end if 3070 3071* ***** read initial wavefunctions into psi1 **** 3072 call epsi_read(ispin0,ne_excited,dcpl_mb(psi1_excited(1))) 3073 3074 3075 3076 3077 3078 paw_epsi_initialize = value 3079 return 3080 end 3081 3082 3083 3084* *************************** 3085* * * 3086* * paw_epsi_finalize * 3087* * * 3088* *************************** 3089 3090 logical function paw_epsi_finalize(writepsi) 3091 implicit none 3092#include "errquit.fh" 3093 logical writepsi 3094 3095#include "bafdecls.fh" 3096#include "paw_psi.fh" 3097 3098* **** local variables **** 3099 logical value 3100 3101* ***** write psi1 wavefunctions **** 3102 if (writepsi) 3103 > call epsi_write(ispin,ne_excited,dcpl_mb(psi1_excited(1))) 3104 3105 value = BA_free_heap(eig_excited(2)) 3106 value = value.and.BA_free_heap(psi2_excited(2)) 3107 value = value.and.BA_free_heap(psi1_excited(2)) 3108 if (.not. value) 3109 > call errquit('epaw_psi_finalize: error freeing heap',0, MA_ERR) 3110 3111 paw_epsi_finalize = value 3112 return 3113 end 3114 3115 3116 3117