1* 2* $Id$ 3* 4 5 6* $Log: not supported by cvs2svn $ 7* Revision 1.98 2009/03/19 20:42:17 bylaska 8* ...EJB 9* 10* Revision 1.97 2009/02/24 21:30:17 bert 11* In psi_finalize, occ1 and occ2 were deallocated in wrong order. 12* 13* Revision 1.96 2009/02/07 03:50:56 bylaska 14* Bassi Vectorization Fix...EJB 15* 16* Revision 1.95 2008/12/18 21:15:51 bylaska 17* ...updates for calculating spin contaminatio....EJB 18* 19* Revision 1.94 2008/11/17 17:25:45 bylaska 20* fractional occupation updates....EJB 21* 22* Revision 1.93 2008/10/22 23:56:43 bylaska 23* added NWCHEM_NWPW_LIBRARY to nwchemrc. fixed bug in paw...EJB 24* 25* Revision 1.92 2008/09/30 19:53:35 bylaska 26* Added Baden's exchange algorithm...EJB 27* 28* Revision 1.91 2008/09/17 00:55:36 bylaska 29* ...EJB 30* 31* Revision 1.90 2008/09/15 20:25:33 bylaska 32* ...fractional bug fix..EJB 33* 34* Revision 1.89 2008/09/11 21:26:51 bylaska 35* ...EJB 36* 37* Revision 1.88 2008/06/21 19:37:16 bylaska 38* 39* initalization error fixed with psi_Tgradient...EJB 40* 41* Revision 1.87 2008/06/02 15:20:04 bylaska 42* ..io fixes...EJB 43* 44* Revision 1.86 2008/05/13 02:10:36 bylaska 45* ...EJB 46* 47* Revision 1.85 2008/04/21 19:34:27 bylaska 48* queue fft added to cpsi_H, bug fixes in DMatrix_dgemm1_rot (MPI_Allgather routine replaced with a MPI_AllReduce for stability), np_orbital keyword replaced with np_dimensions keyword (needed for Parallel3d routines) 49* ...EJB 50* 51* Revision 1.84 2007/11/17 00:25:26 d3p708 52* ....PNJ 53* 54* Revision 1.83 2007/11/17 00:19:09 d3p708 55* pjn 56* bugfix 57* 58* Revision 1.82 2007/11/16 22:32:53 d3p708 59* pjn 60* stuff for berry phase dipole 61* 62* Revision 1.81 2007/10/01 23:02:46 bylaska 63* PAW changes..EJB 64* 65* Revision 1.80 2007/09/29 00:33:37 bylaska 66* ...EJB 67* 68* Revision 1.79 2007/09/24 16:58:14 bylaska 69* ...preliminary PAW modifications... 70* - basis file format changed 71* - .vpp formatting routines added to pspw 72* 73* - zdotc routines currently modified to tzdotc. 74* ...EJB 75* 76* Revision 1.78 2007/09/13 20:38:36 bylaska 77* occupation template added to band...EJB 78* 79* Revision 1.77 2007/03/27 02:02:49 bylaska 80* more qmmm_updates....EJB 81* 82* Revision 1.76 2007/03/22 20:46:22 bylaska 83* New implementation of QM/MM. 84* ....EJB 85* 86* Revision 1.75 2007/02/23 01:24:32 bylaska 87* ...EJB 88* 89* Revision 1.74 2007/02/10 03:56:54 bylaska 90* ...bug fix... 91* ..EJB 92* 93* Revision 1.73 2007/02/10 03:40:18 bylaska 94* replaced calls to Grsm_g_MakeOrtho with Dneall_f_ortho 95* ...EJB 96* 97* Revision 1.72 2007/01/02 18:36:52 bylaska 98* HGH pseudopotentials added to band. 99* ...EJB 100* 101* Revision 1.71 2006/10/13 01:43:58 bylaska 102* tcgmsg code for 2d grid distribution. Also cleaned up Dmatrix_ calls so 103* that they are Dneall_ calls instead. 104* ....EJB 105* 106* Revision 1.70 2006/10/07 00:10:07 bylaska 107* Initial implementation of 2d processor grid parallelization in pspw. Currently works with: 108* 109* task pspw steepest_descent 110* task pspw energy (only minimizer 1, minimizer 2?, other minimizers not yet implemented) 111* 112* Currently only works with USE_MPIF option, straight tcgmsg only partially implemented. Car-Parrinello, HFX, SIC, and various analysis codes are also not yet ported. 113* 114* 115* The number of processors along the orbital dimension, np_orbital, is entered as follows, e.g.: 116* 117* nwpw 118* np_orbital 2 119* end 120* 121* The number of processors along the grid dimension, np_grid, is currently defined using np_orbital as 122* 123* np_grid = np/np_orbital 124* 125* where np is the total number of processors. 126* 127* ...EJB 128* 129* Revision 1.69 2006/09/20 19:18:49 bylaska 130* Adding Dmatrix 131* ...EJB 132* 133* Revision 1.68 2006/08/13 01:03:28 bylaska 134* Checking in code not include in 5.0 release. 135* A chain algorithm was added to Nose-Hoover thermostats. 136* Preliminary implementation of a processor group decomposition added to pspw, i.e. parallel decomposition is over fft grid and electrons. 137* ...EJB 138* 139* Revision 1.67 2006/01/26 18:29:36 bylaska 140* bug fix for gga checking...EJB 141* 142* Revision 1.66 2006/01/06 22:52:28 bylaska 143* parallel io bug fix...EJB 144* 145* Revision 1.65 2006/01/06 21:48:43 bylaska 146* io changes for inversion symmetry....EJB 147* 148* Revision 1.64 2005/12/29 03:06:09 marat 149* qmmm interface stuff 150* 151* Revision 1.63 2005/12/22 01:35:07 bylaska 152* revPBE added and gga logic restructured....EJB 153* 154* Revision 1.62 2005/10/05 21:21:30 bylaska 155* psi_iptr_write added...EJB 156* 157* Revision 1.61 2005/05/24 17:36:27 bylaska 158* Stresses added to SIC and HFX. 159* ....BLYP functional updates 160* ....ECCE hacks 161* ...EJB 162* 163* Revision 1.60 2005/02/09 02:39:10 bylaska 164* .............................EJB 165* 166* Revision 1.59 2005/01/17 20:51:33 edo 167* fixed a couple of FPEs 168* 169* Revision 1.58.2.1 2005/01/17 20:51:06 edo 170* fixed a couple of FPEs 171* 172* Revision 1.58 2004/12/21 16:58:35 bylaska 173* various io fixes for dos...EJB 174* 175* Revision 1.57 2004/12/06 20:03:25 bylaska 176* RMM-DIIS diagonalizer added to pspw. 177* nwpw_list updated to handle multiple lists. 178* ....EJB 179* 180* Revision 1.56 2004/11/29 16:05:21 bylaska 181* Finite difference stresses added to PSPW, BAND, and PAW modules. 182* - This is currently the default for BAND and PAW 183* Fixed the analytic unrestricted gga stress term in PSPW. 184* Fixed unrestricted optimization for minimizers 1 and 2 in the PSPW and PAW modules. 185* Partial implementation of analytic stress in BAND module. 186* - kinetic, ewald, and coulomb stresses have been implemented 187* 188* ....EJB 189* 190* Revision 1.55 2004/11/08 01:32:49 bylaska 191* 192* Unrestricted and closed-shell restricted Hartree-Fock has been implemented into pspw. 193* - works with minimizers 1,2,3,4,6,and 7. 194* - band-by-band minimizer (minimizer 5) not yet implemented 195* - free-space coulomb and screened-coulomb kernels implemented. 196* - the free-space coulomb kernel has been tested. 197* - the screened-coulomb kernel needs to be debugged/tested? 198* - restricted open-shell HF has not yet been implemented 199* 200* ....EJB 201* 202* Revision 1.54 2004/03/08 22:51:27 bylaska 203* Fractional occupation working in pspw with minimizer 4, steepest descent, and Car-Parrinello. 204* 205* Bug fix in velocity initialization in liquid and solid-state Car-Parrinello simulations...incell3 instead of incell2 was used in newton step. 206* 207* Added restart capabilites to thermostat masses...Qe and Qr and eke0 are now propagated to 208* restart Car-Parrinello simulations. 209* 210* SIC input modifications. 211* 212* Wannier orbital output modifications. 213* 214* ....EJB 215* 216* Revision 1.53 2004/03/02 00:10:22 bylaska 217* ....EJB 218* 219* Revision 1.52 2004/03/01 05:14:33 bylaska 220* Mulliken and DOS fixes. 221* Added Mulliken projections based on atomic orbitals 222* Added projected density of states (based on Mulliken projections) 223* ...EJB 224* 225* Revision 1.51 2004/02/06 01:57:22 bylaska 226* 227* Option added to write out temporary psi for Kiril. 228* Tempory psi written if 229* set nwpw:psi_tmp .true. 230* 231* 232* ....EJB 233* 234* Revision 1.50 2004/01/28 00:08:47 bylaska 235* Bug fixes...EJB 236* 237* Revision 1.49 2003/12/13 21:07:41 bylaska 238* Kohn-Sham scf minimizer added to pspw. Mixing includes 239* - simple mixing 240* - Anderson mixing 241* - Johnson-Pulay mixing 242* - Kerker preconditioner 243* 244* ....EJB 245* 246* Revision 1.48 2003/12/02 23:55:34 bylaska 247* density of state generation added...EJB 248* 249* Revision 1.47 2003/12/02 19:17:10 bylaska 250* HGH pseudpotential added. 251* TM, Hamman, HGH, pspw_default, and paw_default pseudopotential libraries have been added. 252* KS minimizer updates. 253* ...EJB 254* 255 256 257************************ f_orb orbitals Part ************************ 258 259* *********************************** 260* * * 261* * psi_minimize_f_orb * 262* * * 263* *********************************** 264 subroutine psi_minimize_f_orb() 265 implicit none 266 267#include "bafdecls.fh" 268#include "psi.fh" 269 270 !*** local variables *** 271 integer maxit_orb 272 integer ii,l 273 real*8 sum,maxerror,error_out,e0 274 275 !*** external functions *** 276 real*8 control_tole 277 external control_tole 278 279 !call psi_gen_density_potentials(1) 280 maxit_orb=120 281 maxerror = control_tole() 282 283 do ii=1,(ne(1)+ne(2)) 284 285 !*** orthogonalize to lower orbitals **** 286 call psi_project_out_f_orb1( 287 > ii, 288 > dcpl_mb(psi1(1)+(ii-1)*npack1)) 289 290 !*** normalize **** 291 call Pack_cc_dot(1, 292 > dcpl_mb(psi1(1) +(ii-1)*npack1), 293 > dcpl_mb(psi1(1) +(ii-1)*npack1), 294 > sum) 295 sum = 1.0d0/dsqrt(sum) 296c call Pack_c_SMul(1,sum, 297c > dcpl_mb(psi1(1) +(ii-1)*npack1), 298c > dcpl_mb(psi1(1) +(ii-1)*npack1)) 299 call Pack_c_SMul1(1,sum,dcpl_mb(psi1(1) +(ii-1)*npack1)) 300 301 !*** minimize orbital **** 302 l = 0 303 2 call psi_KS_update_f_orb(maxit_orb, 304 > maxerror, 305 > 0.001d0,ii,error_out,e0) 306 !write(*,*) "e0:",ii,l,e0,error_out 307 l = l+1 308 if ((error_out.gt.maxerror).and.(l.le.4)) go to 2 309 310 dbl_mb(eig(1)+ii-1) = e0 311 312 end do 313 call psi_sort_f_orb() 314 315 316 return 317 end 318 319 subroutine psi_sort_f_orb() 320 implicit none 321#include "errquit.fh" 322 323#include "bafdecls.fh" 324#include "psi.fh" 325 326 logical value 327 integer i,j,ii,jj,ms 328 integer r1(2) 329 real*8 ei,ej 330 331 value = BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 332 if (.not. value) call errquit( 333 > 'psi_sort_f_orb: out of stack memory',0,MA_ERR) 334 335 do ms=1,ispin 336 337 !*** Bubble sort *** 338 do ii=1,ne(ms) 339 do jj=ii+1,ne(ms) 340 i = ii + (ms-1)*ne(1) 341 j = jj + (ms-1)*ne(1) 342 ei = dbl_mb(eig(1)+i-1) 343 ej = dbl_mb(eig(1)+j-1) 344 345 !*** swap *** 346 if (ej.lt.ei) then 347 dbl_mb(eig(1)+i-1) = ej 348 dbl_mb(eig(1)+j-1) = ei 349 call Pack_c_Copy(1,dcpl_mb(psi1(1)+(i-1)*npack1), 350 > dcpl_mb(r1(1))) 351 call Pack_c_Copy(1,dcpl_mb(psi1(1)+(j-1)*npack1), 352 > dcpl_mb(psi1(1)+(i-1)*npack1)) 353 call Pack_c_Copy(1,dcpl_mb(r1(1)), 354 > dcpl_mb(psi1(1)+(j-1)*npack1)) 355 end if 356 357 end do 358 end do 359 360 end do 361 362 value = BA_pop_stack(r1(2)) 363 if (.not. value) call errquit( 364 > 'psi_sort_f_orb: popping stack memory',1, MA_ERR) 365 return 366 end 367 368* *********************************** 369* * * 370* * psi_KS_update_f_orb * 371* * * 372* *********************************** 373 374* This routine performs a KS update on orbital ii 375* 376 subroutine psi_KS_update_f_orb(maxiteration, 377 > maxerror,perror,ii, 378 > error_out,e0) 379 implicit none 380#include "errquit.fh" 381 integer maxiteration 382 real*8 maxerror,perror 383 integer ii 384 real*8 error_out 385 real*8 e0 386 387#include "bafdecls.fh" 388#include "psi.fh" 389 390* **** local variables **** 391 integer MASTER,taskid 392 parameter (MASTER=0) 393 394 logical value,done,oneloop 395 integer it 396 real*8 eold,percent_error,error0,de0,lmbda_r0,lmbda_r1 397 real*8 theta 398 integer r1(2),t0(2),t(2),g(2) 399 integer psi_ptr 400 401 psi_ptr=psi1(1) 402 403 call Parallel_taskid(taskid) 404 405 value = BA_push_get(mt_dcpl,npack1,'t0',t0(2),t0(1)) 406 value = value.and. 407 > BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 408 value = value.and. 409 > BA_push_get(mt_dcpl,npack1,'g',g(2),g(1)) 410 value = value.and. 411 > BA_push_get(mt_dcpl,npack1,'t',t(2),t(1)) 412 if (.not. value) call errquit( 413 > 'psi_KS_update_f_orb: out of stack memory',0, MA_ERR) 414 415 done = .false. 416 error0 = 0.0d0 417 e0 = 0.0d0 418 theta = -3.14159d0/600.0d0 419 lmbda_r0 = 1.0d0 420 it = 0 421 2 continue 422 423 it = it + 1 424 eold = e0 425 426* *** calculate residual (steepest descent) direction for a single band *** 427 call psi_get_gradient_f_orb(ii,dcpl_mb(g(1))) 428 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 429 > dcpl_mb(g(1)), 430 > e0) 431 432 e0 = -e0 433 percent_error=0d0 434 if(error0.ne.0d0) 435 A percent_error = dabs(e0-eold)/error0 436 437 done = ((it.gt.maxiteration) 438 > .or. 439 > (dabs(e0-eold).lt.maxerror)) 440 441 if (done) go to 4 442 443 444 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 445 call Pack_cc_daxpy(1,(e0), 446 > dcpl_mb(psi_ptr+(ii-1)*npack1), 447 > dcpl_mb(r1(1))) 448 449 450 !*** determine conjuagate direction *** 451 call Pack_cc_dot(1,dcpl_mb(r1(1)), 452 > dcpl_mb(r1(1)), 453 > lmbda_r1) 454 call Pack_c_Copy(1,dcpl_mb(r1(1)),dcpl_mb(t(1))) 455 456 if (it.gt.1) then 457 call Pack_cc_daxpy(1,(lmbda_r1/lmbda_r0), 458 > dcpl_mb(t0(1)), 459 > dcpl_mb(t(1))) 460 end if 461 lmbda_r0 = lmbda_r1 462 oneloop = .true. 463 3 call Pack_c_Copy(1,dcpl_mb(t(1)),dcpl_mb(t0(1))) 464 465 466* *** normalize search direction, t **** 467 call psi_project_out_f_orb(ii,dcpl_mb(t(1))) 468 call Pack_cc_dot(1,dcpl_mb(t(1)), 469 > dcpl_mb(t(1)), 470 > de0) 471 de0 = 1.0d0/dsqrt(de0) 472c call Pack_c_SMul(1,de0,dcpl_mb(t(1)),dcpl_mb(t(1))) 473 call Pack_c_SMul1(1,de0,dcpl_mb(t(1))) 474 call Pack_cc_dot(1,dcpl_mb(t(1)), 475 > dcpl_mb(g(1)), 476 > de0) 477 478* *** bad direction *** 479 if ((de0.lt.0.0d0).and.oneloop) then 480 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(t(1))) 481 oneloop = .false. 482 go to 3 483 end if 484 485 de0 = -2.0d0*de0 486 call psi_linesearch_f_orb(ii, 487 > theta,e0,de0,dcpl_mb(t(1))) 488 489 go to 2 490 491 492* **** release stack memory **** 493 4 value = BA_pop_stack(t(2)) 494 value = value.and.BA_pop_stack(g(2)) 495 value = value.and.BA_pop_stack(r1(2)) 496 value = value.and.BA_pop_stack(t0(2)) 497 if (.not. value) call errquit( 498 > 'psi_KS_update_f_orb: popping stack memory',1, MA_ERR) 499 500 error_out = dabs(e0-eold) 501 e0 = -e0 502 return 503 end 504 505* *********************************** 506* * * 507* * psi_linesearch_f_orb 508* * * 509* *********************************** 510 511* This routine performs a linesearch on orbital ii, in the direction t. 512* This routine is needed for a KS minimizer. 513* e0 = <orb|g> 514* de0 = 2*<t|g> 515* 516 subroutine psi_linesearch_f_orb(ii,theta,e0,de0,t) 517 implicit none 518#include "errquit.fh" 519 integer ii 520 real*8 theta 521 real*8 e0,de0 522 complex*16 t(*) !search direction 523 524#include "bafdecls.fh" 525#include "psi.fh" 526 527* **** local variables **** 528 logical value 529 integer orb(2),g(2),psi_ptr 530 real*8 x,y,pi,dtheta_min,e1 531 532 psi_ptr=psi1(1) 533 534 pi = 4.0d0*datan(1.0d0) 535 !dtheta = pi/300.0d0 536 dtheta_min = 0.01*theta 537 538* **** allocate stack memory **** 539 value = BA_push_get(mt_dcpl,npack1,'orb', 540 > orb(2),orb(1)) 541 value = value.and. 542 > BA_push_get(mt_dcpl,npack1,'g', 543 > g(2),g(1)) 544 if (.not. value) call errquit( 545 > 'psi_linesearch_f_orb: out of stack memory',0, MA_ERR) 546 547 548 call Pack_c_Copy(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 549 > dcpl_mb(orb(1))) 550 551* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 552 10 x = cos(theta) 553 y = sin(theta) 554 call Pack_c_SMul(1,x, 555 > dcpl_mb(orb(1)), 556 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 557 call Pack_cc_daxpy(1,y, 558 > t, 559 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 560 561* *** determine theta *** 562 call psi_get_gradient_f_orb(ii,dcpl_mb(g(1))) 563 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 564 > dcpl_mb(g(1)), 565 > e1) 566 e1 = -e1 567 568 569 570 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 571 > /(1.0d0-cos(2*theta)) 572 theta = 0.5d0*datan(0.5d0*de0/x) 573 574 575 576* **** orb2 = orb*cos(theta) + t*sin(theta) **** 577 x = cos(theta) 578 y = sin(theta) 579 call Pack_c_SMul(1,x, 580 > dcpl_mb(orb(1)), 581 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 582 call Pack_cc_daxpy(1,y, 583 > t, 584 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 585 586 587* **** release stack memory **** 588 value = BA_pop_stack(g(2)) 589 value = value.and.BA_pop_stack(orb(2)) 590 if (.not. value) call errquit( 591 > 'psi_linesearch_f_orb: popping stack memory',1,MA_ERR) 592 593 return 594 end 595 596 597* *********************************** 598* * * 599* * psi_get_gradient_f_orb * 600* * * 601* *********************************** 602 603* This routine returns the Hpsi(i). 604* This routine is needed for a KS minimizer. 605* 606 subroutine psi_get_gradient_f_orb(ii,Horb) 607 implicit none 608 integer ii 609 complex*16 Horb(*) 610 611#include "bafdecls.fh" 612#include "psi.fh" 613 614* **** local variables **** 615 integer psi_ptr,ms 616 617 psi_ptr=psi1(1)+(ii-1)*npack1 618 619 if (ii.le.neq(1)) then 620 ms = 1 621 else 622 ms = 2 623 end if 624 625 call electron_get_gradient_virtual(ms,dcpl_mb(psi_ptr),Horb) 626 627 return 628 end 629 630* ******************************************* 631* * * 632* * psi_project_out_f_orb1 * 633* * * 634* ******************************************* 635* 636* This routine projects out non-orthogonal components of Horb. 637* This routine is needed for a KS minimizer. 638* 639 subroutine psi_project_out_f_orb1(ii,Horb) 640 integer ii 641 complex*16 Horb(*) 642 643#include "bafdecls.fh" 644#include "psi.fh" 645 646 integer ms,jj,kk,shift,shifte 647 real*8 sum 648 649* **** spin up orbital **** 650 if (ii.le.neq(1)) then 651 652 shift = 0 653 shifte = 0 654 ms = 1 655 kk = ii 656* **** spin down orbital **** 657 else 658 shift = neq(1)*npack1 659 shifte = neq(1)*npack1 660 ms = 2 661 kk = ii-neq(1) 662 end if 663 664 665 !**** project out orbitals **** 666 do jj=1,(kk-1) 667 call Pack_cc_dot(1, 668 > dcpl_mb(psi1(1) +(jj-1)*npack1+shifte), 669 > Horb, 670 > sum) 671 672 call Pack_cc_daxpy(1,(-sum), 673 > dcpl_mb(psi1(1) +(jj-1)*npack1+shifte), 674 > Horb) 675 end do 676 677 678 return 679 end 680 681 682 683 684* ******************************************* 685* * * 686* * psi_project_out_f_orb * 687* * * 688* ******************************************* 689* 690* This routine projects out non-orthogonal components of Horb. 691* This routine is needed for a KS minimizer. 692* 693 subroutine psi_project_out_f_orb(ii,Horb) 694 integer ii 695 complex*16 Horb(*) 696 697#include "bafdecls.fh" 698#include "psi.fh" 699 700 integer ms,jj,kk,shift,shifte 701 real*8 sum 702 703* **** spin up orbital **** 704 if (ii.le.neq(1)) then 705 706 shift = 0 707 shifte = 0 708 ms = 1 709 kk = ii 710* **** spin down orbital **** 711 else 712 shift = neq(1)*npack1 713 shifte = neq(1)*npack1 714 ms = 2 715 kk = ii-neq(1) 716 end if 717 718 719 !**** project out orbitals **** 720 do jj=1,(kk) 721 call Pack_cc_dot(1, 722 > dcpl_mb(psi1(1) +(jj-1)*npack1+shifte), 723 > Horb, 724 > sum) 725 726 call Pack_cc_daxpy(1,(-sum), 727 > dcpl_mb(psi1(1) +(jj-1)*npack1+shifte), 728 > Horb) 729 end do 730 731 732 return 733 end 734 735 736 737************************ gradient virtural orbital Part ************************ 738* *********************************** 739* * * 740* * psi_gen_hml_virtual * 741* * * 742* *********************************** 743 subroutine psi_gen_hml_virtual(assending) 744 implicit none 745 logical assending 746 747#include "bafdecls.fh" 748#include "psi.fh" 749#include "errquit.fh" 750 751 !*** local variables *** 752 integer ii,jj,jjstart,jjend,ms,mshift,indx,indxt 753 integer Horb(2) 754 real*8 tsum 755 756* **** allocate temporary memory from stack **** 757 if (.not.BA_push_get(mt_dcpl,npack1,'Horb',Horb(2),Horb(1))) 758 > call errquit('psi_gen_hml_virtual:pusing stack',0,MA_ERR) 759 760 do ii=1,(ne_excited(1)+ne_excited(2)) 761 762 !*** orthogonalize to lower orbitals **** 763 call psi_project_out_virtual1( 764 > ii, 765 > dcpl_mb(psi1_excited(1)+(ii-1)*npack1)) 766 767 !*** normalize **** 768 call Pack_cc_dot(1, 769 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 770 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 771 > tsum) 772 tsum = 1.0d0/dsqrt(tsum) 773 call Pack_c_SMul1(1,tsum, 774 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1)) 775 776 if (ii.le.ne_excited(1)) then 777 ms = 1 778 jjstart = 1 779 jjend = ii 780 mshift = 0 781 else 782 ms = 2 783 jjstart = ne_excited(1)+1 784 jjend = ii 785 mshift = ne_excited(1)*ne_excited(1) 786 end if 787c call Pack_cc_dot(1, 788c > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 789c > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 790c > tsum) 791 !write(*,*) "ii,norm=",ii,tsum 792 793 call electron_get_gradient_virtual(ms, 794 > dcpl_mb(psi1_excited(1)+(ii-1)*npack1), 795 > dcpl_mb(Horb(1))) 796 797 do jj=jjstart,jjend 798c call Pack_cc_dot(1, 799c > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 800c > dcpl_mb(psi1_excited(1) +(jj-1)*npack1), 801c > tsum) 802 !write(*,*) "ii,jj,norm=",ii,jj,tsum 803 call Pack_cc_dot(1, 804 > dcpl_mb(psi1_excited(1) +(jj-1)*npack1), 805 > dcpl_mb(Horb(1)), 806 > tsum) 807 !write(*,*) "ii,jj,e=",ii,jj,tsum 808 indx = (ii-jjstart)+(jj-jjstart)*ne_excited(ms)+mshift 809 dbl_mb(hml_excited(1)+indx) = tsum 810 if (ii.ne.jj) then 811 indxt = (jj-jjstart)+(ii-jjstart)*ne_excited(ms)+mshift 812 dbl_mb(hml_excited(1)+indxt) = tsum 813 end if 814 end do 815 end do 816 817* **** deallocate temporary memory from stack **** 818 if (.not.BA_pop_stack(Horb(2))) 819 > call errquit('psi_gen_hml_virtual:popping stack',0,MA_ERR) 820 821c call Dnexall_m_diagonalize(0,ispin,ne_excited, 822c > dbl_mb(hml_excited(1)), 823c > dbl_mb(eig_excited(1)),.true.) 824 call Dnexall_m_diagonalize(0,ispin,ne_excited, 825 > dbl_mb(hml_excited(1)), 826 > dbl_mb(eig_excited(1)),assending) 827 call Dnexall_fmf_Multiply(0,ispin,ne_excited, 828 > dcpl_mb(psi1_excited(1)),npack1, 829 > dbl_mb(hml_excited(1)),1.0d0, 830 > dcpl_mb(psi2_excited(1)),0.0) 831 indx = psi1_excited(1) 832 psi1_excited(1) = psi2_excited(1) 833 psi2_excited(1) = indx 834c call dcopy(2*npack1, 835c > dcpl_mb(psi1_excited(1)+3*npack1),1, 836c > dcpl_mb(psi1(1)),1) 837 838 839 return 840 end 841 842 843 844 845************************ virtural orbital Part ************************ 846* *********************************** 847* * * 848* * psi_minimize_virtual * 849* * * 850* *********************************** 851 subroutine psi_minimize_virtual() 852 implicit none 853 854#include "bafdecls.fh" 855#include "psi.fh" 856 857 !*** local variables *** 858 integer maxit_orb 859 integer ii,l 860 real*8 sum,maxerror,error_out,e0 861 862 !*** external functions *** 863 real*8 control_tole 864 external control_tole 865 866 !call psi_gen_density_potentials(1) 867 maxit_orb=120 868 maxerror = control_tole() 869 870 do ii=1,(ne_excited(1)+ne_excited(2)) 871 872 !*** orthogonalize to lower orbitals **** 873 call psi_project_out_virtual1( 874 > ii, 875 > dcpl_mb(psi1_excited(1)+(ii-1)*npack1)) 876 877 !*** normalize **** 878 call Pack_cc_dot(1, 879 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 880 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1), 881 > sum) 882 sum = 1.0d0/dsqrt(sum) 883 call Pack_c_SMul1(1,sum, 884 > dcpl_mb(psi1_excited(1) +(ii-1)*npack1)) 885 886 !*** minimize orbital **** 887 l = 0 888 2 call psi_KS_update_virtual(maxit_orb, 889 > maxerror, 890 > 0.001d0,ii,error_out,e0) 891 l = l+1 892 if ((error_out.gt.maxerror).and.(l.le.4)) go to 2 893 894 dbl_mb(eig_excited(1)+ii-1) = e0 895 896 end do 897 call psi_sort_virtual() 898 899 900 return 901 end 902 903 subroutine psi_sort_virtual() 904 implicit none 905 906#include "bafdecls.fh" 907#include "psi.fh" 908#include "errquit.fh" 909 910 logical value 911 integer i,j,ii,jj,ms 912 integer r1(2) 913 real*8 ei,ej 914 915 value = BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 916 if (.not. value) call errquit( 917 > 'psi_sort_virtual: out of stack memory',0, MA_ERR) 918 919 do ms=1,ispin 920 921 !*** Bubble sort *** 922 do ii=1,ne_excited(ms) 923 do jj=ii+1,ne_excited(ms) 924 i = ii + (ms-1)*ne_excited(1) 925 j = jj + (ms-1)*ne_excited(1) 926 ei = dbl_mb(eig_excited(1)+i-1) 927 ej = dbl_mb(eig_excited(1)+j-1) 928 929 !*** swap *** 930 if (ej.lt.ei) then 931 dbl_mb(eig_excited(1)+i-1) = ej 932 dbl_mb(eig_excited(1)+j-1) = ei 933 call Pack_c_Copy(1,dcpl_mb(psi1_excited(1)+(i-1)*npack1), 934 > dcpl_mb(r1(1))) 935 call Pack_c_Copy(1,dcpl_mb(psi1_excited(1)+(j-1)*npack1), 936 > dcpl_mb(psi1_excited(1)+(i-1)*npack1)) 937 call Pack_c_Copy(1,dcpl_mb(r1(1)), 938 > dcpl_mb(psi1_excited(1)+(j-1)*npack1)) 939 end if 940 941 end do 942 end do 943 944 end do 945 946 value = BA_pop_stack(r1(2)) 947 if (.not. value) call errquit( 948 > 'psi_sort_virtual: popping stack memory',1, MA_ERR) 949 return 950 end 951 952* *********************************** 953* * * 954* * psi_KS_update_virtual * 955* * * 956* *********************************** 957 958* This routine performs a KS update on virtual ii 959* 960 subroutine psi_KS_update_virtual(maxiteration, 961 > maxerror,perror,ii, 962 > error_out,e0) 963 implicit none 964 integer maxiteration 965 real*8 maxerror,perror 966 integer ii 967 real*8 error_out 968 real*8 e0 969 970#include "bafdecls.fh" 971#include "psi.fh" 972#include "errquit.fh" 973#include "util.fh" 974#include "stdio.fh" 975 976* **** local variables **** 977 integer MASTER,taskid 978 parameter (MASTER=0) 979 980 logical value,done,oneloop,precondition,oprint 981 integer it,pit 982 real*8 eold,percent_error,error0,de0,lmbda_r0,lmbda_r1 983 real*8 theta,ep,sp 984 integer r1(2),t0(2),t(2),g(2) 985 integer psi_ptr 986 987 logical control_print 988 externalcontrol_print 989 real*8 control_Ep,control_Sp 990 external control_Ep,control_Sp 991 992 psi_ptr=psi1_excited(1) 993 994 call Parallel_taskid(taskid) 995 oprint= ((taskid.eq.MASTER).and.control_print(print_medium)) 996 997 998 value = BA_push_get(mt_dcpl,npack1,'t0',t0(2),t0(1)) 999 value = value.and. 1000 > BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 1001 value = value.and. 1002 > BA_push_get(mt_dcpl,npack1,'g',g(2),g(1)) 1003 value = value.and. 1004 > BA_push_get(mt_dcpl,npack1,'t',t(2),t(1)) 1005 if (.not. value) call errquit( 1006 > 'psi_KS_update_virtual: out of stack memory',0, MA_ERR) 1007 1008 ep = control_Ep() 1009 sp = control_Sp() 1010 precondition = .true. 1011 done = .false. 1012 error0 = 0.0d0 1013 e0 = 0.0d0 1014 theta = -3.14159d0/600.0d0 1015 lmbda_r0 = 1.0d0 1016 it = 0 1017 pit = 0 1018 2 continue 1019 1020 it = it + 1 1021 eold = e0 1022 1023* *** calculate residual (steepest descent) direction for a single band *** 1024 call psi_get_gradient_virtual(ii,dcpl_mb(g(1))) 1025 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 1026 > dcpl_mb(g(1)), 1027 > e0) 1028 1029 e0 = -e0 1030 1031 !if (it.eq.1) then 1032 ! percent_error = 1.0d0 1033 !else if (it.eq.2) then 1034 ! error0 = dabs(e0-eold) 1035 ! percent_error = 1.0d0 1036 !else 1037 percent_error=0.0d0 1038 if(error0.ne.0.0d0) 1039 A percent_error = dabs(e0-eold)/error0 1040 !end if 1041 1042 precondition = (dabs(e0-eold).gt.(sp*maxerror)) 1043 1044 done = ((it.gt.maxiteration) 1045 > .or. 1046 > (dabs(e0-eold).lt.maxerror)) 1047 1048 if (done) go to 4 1049 1050 1051 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 1052 call Pack_cc_daxpy(1,(e0), 1053 > dcpl_mb(psi_ptr+(ii-1)*npack1), 1054 > dcpl_mb(r1(1))) 1055 1056* **** preconditioning **** 1057 if (precondition) then 1058 pit = pit + 1 1059 call ke_Precondition(npack1,1,dcpl_mb(r1(1)),dcpl_mb(r1(1))) 1060 end if 1061 1062 !*** determine conjuagate direction *** 1063 call Pack_cc_dot(1,dcpl_mb(r1(1)), 1064 > dcpl_mb(r1(1)), 1065 > lmbda_r1) 1066 call Pack_c_Copy(1,dcpl_mb(r1(1)),dcpl_mb(t(1))) 1067 1068 if (it.gt.1) then 1069 call Pack_cc_daxpy(1,(lmbda_r1/lmbda_r0), 1070 > dcpl_mb(t0(1)), 1071 > dcpl_mb(t(1))) 1072 end if 1073 lmbda_r0 = lmbda_r1 1074 oneloop = .true. 1075 3 call Pack_c_Copy(1,dcpl_mb(t(1)),dcpl_mb(t0(1))) 1076 1077 1078* *** normalize search direction, t **** 1079 call psi_project_out_virtual(ii,dcpl_mb(t(1))) 1080 call Pack_cc_dot(1,dcpl_mb(t(1)), 1081 > dcpl_mb(t(1)), 1082 > de0) 1083 de0 = 1.0d0/dsqrt(de0) 1084c call Pack_c_SMul(1,de0,dcpl_mb(t(1)),dcpl_mb(t(1))) 1085 call Pack_c_SMul1(1,de0,dcpl_mb(t(1))) 1086 call Pack_cc_dot(1,dcpl_mb(t(1)), 1087 > dcpl_mb(g(1)), 1088 > de0) 1089 1090* *** bad direction *** 1091 if ((de0.lt.0.0d0).and.oneloop) then 1092 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(t(1))) 1093 oneloop = .false. 1094 go to 3 1095 end if 1096 1097 de0 = -2.0d0*de0 1098 call psi_linesearch_virtual(ii, 1099 > theta,e0,de0,dcpl_mb(t(1))) 1100 1101 go to 2 1102 1103 1104* **** release stack memory **** 1105 4 value = BA_pop_stack(t(2)) 1106 value = value.and.BA_pop_stack(g(2)) 1107 value = value.and.BA_pop_stack(r1(2)) 1108 value = value.and.BA_pop_stack(t0(2)) 1109 if (.not. value) call errquit( 1110 > 'psi_KS_update_virtual: popping stack memory',1, MA_ERR) 1111 1112 if (oprint) then 1113 write(luout,921) ii,-e0,dabs(e0-eold),it,pit,ep,sp 1114 921 format(5x,"orbital",I4," current e=",E10.3, 1115 > " (error=",E9.3,")", 1116 > " iterations",I4,"(",I4, 1117 > " preconditioned, Ep,Sp=",F5.1,F7.1,")") 1118 end if 1119 1120 error_out = dabs(e0-eold) 1121 e0 = -e0 1122 return 1123 end 1124 1125* *********************************** 1126* * * 1127* * psi_linesearch_virtual 1128* * * 1129* *********************************** 1130 1131* This routine performs a linesearch on orbital ii, in the direction t. 1132* This routine is needed for a KS minimizer. 1133* e0 = <orb|g> 1134* de0 = 2*<t|g> 1135* 1136 subroutine psi_linesearch_virtual(ii,theta,e0,de0,t) 1137 implicit none 1138 integer ii 1139 real*8 theta 1140 real*8 e0,de0 1141 complex*16 t(*) !search direction 1142#include "errquit.fh" 1143#include "bafdecls.fh" 1144#include "psi.fh" 1145 1146* **** local variables **** 1147 logical value 1148 integer orb(2),g(2),psi_ptr 1149 real*8 x,y,pi,dtheta_min,e1 1150 1151 psi_ptr=psi1_excited(1) 1152 1153 pi = 4.0d0*datan(1.0d0) 1154 !dtheta = pi/300.0d0 1155 dtheta_min = 0.01*theta 1156 1157* **** allocate stack memory **** 1158 value = BA_push_get(mt_dcpl,npack1,'orb', 1159 > orb(2),orb(1)) 1160 value = value.and. 1161 > BA_push_get(mt_dcpl,npack1,'g', 1162 > g(2),g(1)) 1163 if (.not. value) call errquit( 1164 > 'psi_linesearch_virtual: out of stack memory',0, MA_ERR) 1165 1166 1167 call Pack_c_Copy(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 1168 > dcpl_mb(orb(1))) 1169 1170* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 1171 10 x = cos(theta) 1172 y = sin(theta) 1173 call Pack_c_SMul(1,x, 1174 > dcpl_mb(orb(1)), 1175 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 1176 call Pack_cc_daxpy(1,y, 1177 > t, 1178 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 1179 1180* *** determine theta *** 1181 call psi_get_gradient_virtual(ii,dcpl_mb(g(1))) 1182 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(ii-1)*npack1), 1183 > dcpl_mb(g(1)), 1184 > e1) 1185 e1 = -e1 1186 1187 1188 !if (((-e1).gt.(-e0)).and.(theta.gt.dtheta_min)) then 1189 ! theta = 0.5d0*theta 1190 ! go to 10 1191 !end if 1192 1193 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 1194 > /(1.0d0-cos(2*theta)) 1195 theta = 0.5d0*datan(0.5d0*de0/x) 1196 1197 1198 1199* **** orb2 = orb*cos(theta) + t*sin(theta) **** 1200 x = cos(theta) 1201 y = sin(theta) 1202 call Pack_c_SMul(1,x, 1203 > dcpl_mb(orb(1)), 1204 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 1205 call Pack_cc_daxpy(1,y, 1206 > t, 1207 > dcpl_mb(psi_ptr+(ii-1)*npack1)) 1208 1209 1210* **** release stack memory **** 1211 value = BA_pop_stack(g(2)) 1212 value = value.and.BA_pop_stack(orb(2)) 1213 if (.not. value) call errquit( 1214 > 'psi_linesearch_virtual: popping stack memory',1, MA_ERR) 1215 1216 return 1217 end 1218 1219 1220* *********************************** 1221* * * 1222* * psi_get_gradient_virtual * 1223* * * 1224* *********************************** 1225 1226* This routine returns the Hpsi(i). 1227* This routine is needed for a KS minimizer. 1228* 1229 subroutine psi_get_gradient_virtual(ii,Horb) 1230 implicit none 1231 integer ii 1232 complex*16 Horb(*) 1233 1234#include "bafdecls.fh" 1235#include "psi.fh" 1236 1237* **** local variables **** 1238 integer psi_ptr,ms 1239 1240 psi_ptr=psi1_excited(1)+(ii-1)*npack1 1241 1242 if (ii.le.ne_excited(1)) then 1243 ms = 1 1244 else 1245 ms = 2 1246 end if 1247 1248 call electron_get_gradient_virtual(ms,dcpl_mb(psi_ptr),Horb) 1249 1250 return 1251 end 1252 1253* ******************************************* 1254* * * 1255* * psi_project_out_virtual1 * 1256* * * 1257* ******************************************* 1258* 1259* This routine projects out non-orthogonal components of Horb. 1260* This routine is needed for a KS minimizer. 1261* 1262 subroutine psi_project_out_virtual1(ii,Horb) 1263 integer ii 1264 complex*16 Horb(*) 1265 1266#include "bafdecls.fh" 1267#include "errquit.fh" 1268#include "psi.fh" 1269 1270 integer ms,i,jj,kk,shift,shifte 1271 integer etmp(2) 1272 real*8 sum 1273 1274* **** spin up orbital **** 1275 if (ii.le.ne_excited(1)) then 1276 1277 shift = 0 1278 shifte = 0 1279 ms = 1 1280 kk = ii 1281* **** spin down orbital **** 1282 else 1283 shift = neq(1)*npack1 1284 shifte = ne_excited(1)*npack1 1285 ms = 2 1286 kk = ii-ne_excited(1) 1287 end if 1288 1289 !**** project out filled orbitals **** 1290 if (neq(ms).eq.ne(ms)) then 1291 1292 do i=1,ne(ms) 1293 call Pack_cc_dot(1, 1294 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1295 > Horb, 1296 > sum) 1297 call Pack_cc_daxpy(1,(-sum), 1298 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1299 > Horb) 1300 end do 1301 1302 else 1303 1304 if (.not.BA_push_get(mt_dcpl,npack1,'etmp',etmp(2),etmp(1))) 1305 > call errquit('psi_project_out_virtual1: out of stack',0,MA_ERR) 1306 1307 !call dcopy(2*npack1,0.0d0,0,dcpl_mb(etmp(1)),1) 1308 call Parallel_shared_vector_zero(.true., 1309 > 2*npack1,dcpl_mb(etmp(1))) 1310 do i=1,neq(ms) 1311 call Pack_cc_dot(1, 1312 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1313 > Horb, 1314 > sum) 1315 call Pack_cc_daxpy(1,(-sum), 1316 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1317 > dcpl_mb(etmp(1))) 1318 end do 1319 call D1dB_Vector_SumAll(2*npack1,dcpl_mb(etmp(1))) 1320 call daxpy_omp(2*npack1,1.0d0,dcpl_mb(etmp(1)),1,Horb,1) 1321 1322 if (.not.BA_pop_stack(etmp(2))) 1323 > call errquit('psi_project_out_virtual1:popping stack',0,MA_ERR) 1324 1325 end if 1326 1327 !**** project out virtual orbitals **** 1328 do jj=1,(kk-1) 1329 call Pack_cc_dot(1, 1330 > dcpl_mb(psi1_excited(1) +(jj-1)*npack1+shifte), 1331 > Horb, 1332 > sum) 1333 1334 call Pack_cc_daxpy(1,(-sum), 1335 > dcpl_mb(psi1_excited(1) +(jj-1)*npack1+shifte), 1336 > Horb) 1337 end do 1338 1339 1340 return 1341 end 1342 1343 1344 1345 1346* ******************************************* 1347* * * 1348* * psi_project_out_virtual * 1349* * * 1350* ******************************************* 1351* 1352* This routine projects out non-orthogonal components of Horb. 1353* This routine is needed for a KS minimizer. 1354* 1355 subroutine psi_project_out_virtual(ii,Horb) 1356 integer ii 1357 complex*16 Horb(*) 1358 1359#include "bafdecls.fh" 1360#include "errquit.fh" 1361#include "psi.fh" 1362 1363 integer ms,i,jj,kk,shift,shifte 1364 integer etmp(2) 1365 real*8 sum 1366 1367* **** spin up orbital **** 1368 if (ii.le.ne_excited(1)) then 1369 1370 shift = 0 1371 shifte = 0 1372 ms = 1 1373 kk = ii 1374* **** spin down orbital **** 1375 else 1376 shift = neq(1)*npack1 1377 shifte = ne_excited(1)*npack1 1378 ms = 2 1379 kk = ii-ne_excited(1) 1380 end if 1381 1382 !**** project out filled orbitals **** 1383 if (neq(ms).eq.ne(ms)) then 1384 1385 do i=1,ne(ms) 1386 call Pack_cc_dot(1, 1387 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1388 > Horb, 1389 > sum) 1390 1391 call Pack_cc_daxpy(1,(-sum), 1392 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1393 > Horb) 1394 end do 1395 1396 else 1397 if (.not.BA_push_get(mt_dcpl,npack1,'etmp',etmp(2),etmp(1))) 1398 > call errquit('psi_project_out_virtual1: out of stack',0,MA_ERR) 1399 1400 !call dcopy(2*npack1,0.0d0,0,dcpl_mb(etmp(1)),1) 1401 call Parallel_shared_vector_zero(.true., 1402 > 2*npack1,dcpl_mb(etmp(1))) 1403 do i=1,neq(ms) 1404 call Pack_cc_dot(1, 1405 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1406 > Horb, 1407 > sum) 1408 call Pack_cc_daxpy(1,(-sum), 1409 > dcpl_mb(psi1(1) +(i-1)*npack1+shift), 1410 > dcpl_mb(etmp(1))) 1411 end do 1412 call D1dB_Vector_SumAll(2*npack1,dcpl_mb(etmp(1))) 1413 call daxpy(2*npack1,1.0d0,dcpl_mb(etmp(1)),1,Horb,1) 1414 1415 if (.not.BA_pop_stack(etmp(2))) 1416 > call errquit('psi_project_out_virtual:popping stack',0,MA_ERR) 1417 1418 end if 1419 1420 !**** project out virtual orbitals **** 1421 do jj=1,(kk) 1422 call Pack_cc_dot(1, 1423 > dcpl_mb(psi1_excited(1) +(jj-1)*npack1+shifte), 1424 > Horb, 1425 > sum) 1426 1427 call Pack_cc_daxpy(1,(-sum), 1428 > dcpl_mb(psi1_excited(1) +(jj-1)*npack1+shifte), 1429 > Horb) 1430 end do 1431 1432 1433 return 1434 end 1435 1436 1437 1438************************ KS orbital Part ************************ 1439 1440* *********************************** 1441* * * 1442* * psi_KS_update * 1443* * * 1444* *********************************** 1445 1446* This routine (approximately) diagonalizes the KS matrix. 1447* 1448 subroutine psi_KS_update(psi_number, 1449 > ks_algorithm, 1450 > precondition, 1451 > maxerror) 1452 implicit none 1453 integer psi_number 1454 integer ks_algorithm 1455 logical precondition 1456 real*8 maxerror 1457 1458#include "bafdecls.fh" 1459#include "psi.fh" 1460 1461* **** local variables **** 1462 logical done 1463 integer i,j,neall,maxit_orb,maxit_orbs 1464 real*8 error,error_out,tim1,tim2,tim,sum 1465 1466* **** external functions **** 1467 integer control_ks_maxit_orb,control_ks_maxit_orbs 1468 external control_ks_maxit_orb,control_ks_maxit_orbs 1469 1470 tim = 0.0d0 1471 neall = neq(1)+neq(2) 1472 maxit_orb = control_ks_maxit_orb() !*** should be read from rtdb *** 1473 maxit_orbs = control_ks_maxit_orbs() !*** should be read from rtdb *** 1474 j = 0 1475 2 j = j+1 1476 error = 0.0d0 1477 !do i=neall,1,-1 1478 do i=1,neall 1479 1480 !*** orthogonalize to lower orbitals **** 1481 call psi_project_out_f_orb1( 1482 > i, 1483 > dcpl_mb(psi1(1)+(i-1)*npack1)) 1484 1485 !*** normalize **** 1486 call Pack_cc_dot(1, 1487 > dcpl_mb(psi1(1) +(i-1)*npack1), 1488 > dcpl_mb(psi1(1) +(i-1)*npack1), 1489 > sum) 1490 sum = 1.0d0/dsqrt(sum) 1491c call Pack_c_SMul(1,sum, 1492c > dcpl_mb(psi1(1) +(i-1)*npack1), 1493c > dcpl_mb(psi1(1) +(i-1)*npack1)) 1494 call Pack_c_SMul1(1,sum, 1495 > dcpl_mb(psi1(1) +(i-1)*npack1)) 1496 1497 1498 1499 if (ks_algorithm.eq.1) then 1500 call psi_KS_update_orb2(psi_number,precondition,maxit_orb, 1501 > maxerror, 1502 > 0.1d0,i,error_out) 1503 else 1504 call psi_KS_update_orb(psi_number,precondition,maxit_orb, 1505 > maxerror, 1506 > 0.1d0,i,error_out) 1507 end if 1508 1509 error = error+error_out 1510 end do 1511 error = error/dble(neall) 1512 1513 done = ((j.gt.maxit_orbs).or.(error.lt.maxerror)) 1514 if (.not.done) go to 2 1515 1516 return 1517 end 1518 1519 1520* *********************************** 1521* * * 1522* * psi_KS_update_orb * 1523* * * 1524* *********************************** 1525 1526* This routine performs a KS update on orbital i 1527* 1528 subroutine psi_KS_update_orb(psi_number, 1529 > precondition,maxiteration, 1530 > maxerror,perror,i, 1531 > error_out) 1532 implicit none 1533 integer psi_number 1534 logical precondition 1535 integer maxiteration 1536 real*8 maxerror,perror 1537 integer i 1538 real*8 error_out 1539 1540#include "bafdecls.fh" 1541#include "psi.fh" 1542#include "errquit.fh" 1543 1544* **** local variables **** 1545 integer MASTER,taskid 1546 parameter (MASTER=0) 1547 1548 logical value,done,oneloop 1549 integer it 1550 real*8 e0,eold,error0,de0,lmbda_r0,lmbda_r1 1551 real*8 theta,sigma 1552 integer r1(2),t0(2),t(2),g(2) 1553 integer psi_ptr 1554 1555 if (psi_number.eq.1) then 1556 psi_ptr=psi1(1) 1557 else 1558 psi_ptr=psi2(1) 1559 end if 1560 1561 call Parallel_taskid(taskid) 1562 1563 value = BA_push_get(mt_dcpl,npack1,'t0',t0(2),t0(1)) 1564 value = value.and. 1565 > BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 1566 value = value.and. 1567 > BA_push_get(mt_dcpl,npack1,'g',g(2),g(1)) 1568 value = value.and. 1569 > BA_push_get(mt_dcpl,npack1,'t',t(2),t(1)) 1570 if (.not. value) call errquit( 1571 > 'psi_KS_update_orb: out of stack memory',0, MA_ERR) 1572 1573 done = .false. 1574 error0 = 0.0d0 1575 e0 = 0.0d0 1576 theta = -3.14159d0/600.0d0 1577 lmbda_r0 = 1.0d0 1578 it = 0 1579 2 continue 1580 1581 it = it + 1 1582 eold = e0 1583 1584* *** calculate residual (steepest descent) direction for a single band *** 1585 call psi_get_gradient_orb(psi_number,i,dcpl_mb(g(1))) 1586 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(i-1)*npack1), 1587 > dcpl_mb(g(1)), 1588 > e0) 1589 !dbl_mb(eig(1)+i-1) = e0 1590 e0 = -e0 1591 1592 done = ((it.gt.maxiteration) 1593 > .or. 1594 > (dabs(e0-eold).lt.maxerror)) 1595 1596 if (done) go to 4 1597 1598* **** preconditioning **** 1599 if (precondition) then 1600 call ke_Precondition(npack1,1, 1601 > dcpl_mb(g(1)), 1602 > dcpl_mb(g(1))) 1603 1604 end if 1605 1606c call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 1607c call Pack_cc_daxpy(1,(e0), 1608c > dcpl_mb(psi_ptr+(i-1)*npack1), 1609c > dcpl_mb(r1(1))) 1610 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(r1(1))) 1611 call psi_project_out_orb(psi_number,i,dcpl_mb(r1(1))) 1612 1613 1614 1615 1616* *** determine conjuagate direction *** 1617 call Pack_cc_dot(1,dcpl_mb(r1(1)), 1618 > dcpl_mb(r1(1)), 1619 > lmbda_r1) 1620 call Pack_c_Copy(1,dcpl_mb(r1(1)),dcpl_mb(t(1))) 1621 1622 if (it.gt.1) then 1623 call Pack_cc_daxpy(1,(lmbda_r1/lmbda_r0), 1624 > dcpl_mb(t0(1)), 1625 > dcpl_mb(t(1))) 1626 end if 1627 lmbda_r0 = lmbda_r1 1628 oneloop = .true. 1629 3 call Pack_c_Copy(1,dcpl_mb(t(1)),dcpl_mb(t0(1))) 1630 1631 1632 1633* **** project out psi components from t - may not be needed! **** 1634 !call psi_project_out_orb(psi_number,i,dcpl_mb(t(1))) 1635c! call psi_project_out_orb(psi_number,i,dcpl_mb(t(1))) 1636c! call Pack_cc_dot(1,dcpl_mb(psi_ptr+(i-1)*npack1), 1637c! > dcpl_mb(t(1)), 1638c! > de0) 1639c! de0 = -de0 1640c! call Pack_cc_daxpy(1,(de0), 1641c! > dcpl_mb(psi_ptr+(i-1)*npack1), 1642c! > dcpl_mb(t(1))) 1643 1644 1645* *** normalize search direction, t **** 1646 call Pack_cc_dot(1,dcpl_mb(t(1)), 1647 > dcpl_mb(t(1)), 1648 > sigma) 1649 sigma = dsqrt(sigma) 1650 de0 = 1.0d0/sigma 1651c call Pack_c_SMul(1,de0,dcpl_mb(t(1)),dcpl_mb(t(1))) 1652 call Pack_c_SMul1(1,de0,dcpl_mb(t(1))) 1653 1654 1655* **** compute de0 = <t|g> **** 1656 call Pack_cc_dot(1,dcpl_mb(t(1)), 1657 > dcpl_mb(g(1)), 1658 > de0) 1659 1660* *** bad direction *** 1661 if ((de0.lt.0.0d0).and.oneloop) then 1662 call Pack_c_Copy(1,dcpl_mb(g(1)),dcpl_mb(t(1))) 1663 oneloop = .false. 1664 go to 3 1665 end if 1666 1667 de0 = -2.0d0*de0 1668 call psi_linesearch_update2(psi_number,i, 1669 > theta,e0,de0, 1670 > dcpl_mb(t(1)), 1671 > sigma, 1672 > dcpl_mb(t0(1))) 1673 1674 go to 2 1675 1676 1677* **** release stack memory **** 1678 4 value = BA_pop_stack(t(2)) 1679 value = value.and.BA_pop_stack(g(2)) 1680 value = value.and.BA_pop_stack(r1(2)) 1681 value = value.and.BA_pop_stack(t0(2)) 1682 if (.not. value) call errquit( 1683 > 'psi_KS_update_orb: popping stack memory',1, MA_ERR) 1684 1685c write(*,*) "iterations=",it," eig=",e0," error=",error_out, 1686c > theta 1687 error_out = dabs(e0-eold) 1688 return 1689 end 1690 1691 1692 1693 1694* *********************************** 1695* * * 1696* * psi_KS_update_orb2 * 1697* * * 1698* *********************************** 1699 1700* This routine performs a RMM-DIIS KS update on orbital i 1701* 1702 subroutine psi_KS_update_orb2(psi_number, 1703 > precondition,maxiteration, 1704 > maxerror,perror,i, 1705 > error_out) 1706 implicit none 1707 integer psi_number 1708 logical precondition 1709 integer maxiteration 1710 real*8 maxerror,perror 1711 integer i 1712 real*8 error_out 1713 1714#include "bafdecls.fh" 1715#include "errquit.fh" 1716#include "psi.fh" 1717 1718* **** local variables **** 1719 integer MASTER,taskid 1720 parameter (MASTER=0) 1721 1722 logical value,done 1723 integer it 1724 real*8 sigma,e0,eold,error0 1725 real*8 lambda 1726 integer r1(2),g1(2) 1727 integer psi_ptr 1728 1729 if (psi_number.eq.1) then 1730 psi_ptr=psi1(1) 1731 else 1732 psi_ptr=psi2(1) 1733 end if 1734 1735 call Parallel_taskid(taskid) 1736 1737* **** allocate memory **** 1738 value = BA_push_get(mt_dcpl,npack1,'r1',r1(2),r1(1)) 1739 value = value.and. BA_push_get(mt_dcpl,npack1,'g1',g1(2),g1(1)) 1740 if (.not. value) call errquit( 1741 > 'psi_KS_update_orb2: out of stack memory',0, MA_ERR) 1742 1743* **** set lambda *** 1744 lambda = 0.1d0 1745 1746 1747* *** calculate residual (steepest descent) direction for a single band *** 1748 call psi_get_gradient_orb(psi_number,i,dcpl_mb(g1(1))) 1749 call Pack_cc_dot(1, 1750 > dcpl_mb(psi_ptr+(i-1)*npack1), 1751 > dcpl_mb(g1(1)), 1752 > e0) 1753 call Pack_cc_dot(1, 1754 > dcpl_mb(psi_ptr+(i-1)*npack1), 1755 > dcpl_mb(psi_ptr+(i-1)*npack1), 1756 > sigma) 1757 e0 = e0/sigma 1758 call Pack_c_SMul(1,e0, 1759 > dcpl_mb(psi_ptr+(i-1)*npack1), 1760 > dcpl_mb(r1(1))) 1761 call Pack_cc_daxpy(1,(-1.0d0), 1762 > dcpl_mb(g1(1)), 1763 > dcpl_mb(r1(1))) 1764 1765 !write(*,*) "i=",i,"it=",0, " eig=",e0,sigma 1766 1767* ***** rmmdiis start ***** 1768 call pspw_rmmdiis_start(lambda, 1769 > dcpl_mb(r1(1)), 1770 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1771 1772 1773 done = .false. 1774 error0 = 0.0d0 1775 it = 0 1776 2 continue 1777 1778 it = it + 1 1779 eold = e0 1780 1781* *** calculate residual (steepest descent) direction for a single band *** 1782 call psi_get_gradient_orb(psi_number,i,dcpl_mb(g1(1))) 1783 call Pack_cc_dot(1, 1784 > dcpl_mb(psi_ptr+(i-1)*npack1), 1785 > dcpl_mb(g1(1)), 1786 > e0) 1787 call Pack_cc_dot(1, 1788 > dcpl_mb(psi_ptr+(i-1)*npack1), 1789 > dcpl_mb(psi_ptr+(i-1)*npack1), 1790 > sigma) 1791 e0 = e0/sigma 1792 call Pack_c_SMul(1,(e0), 1793 > dcpl_mb(psi_ptr+(i-1)*npack1), 1794 > dcpl_mb(r1(1))) 1795 call Pack_cc_daxpy(1,(-1.0d0), 1796 > dcpl_mb(g1(1)), 1797 > dcpl_mb(r1(1))) 1798 !e0 = -e0 1799 !write(*,*) "i=",i,"it=",it, " eig=",e0,sigma,dabs(e0-eold) 1800 1801 done = ((it.gt.maxiteration) 1802 > .or. 1803 > (dabs(e0-eold).lt.maxerror)) 1804 1805 if (done) go to 4 1806 1807* ***** rmmdiis update ***** 1808 call pspw_rmmdiis(lambda, 1809 > dcpl_mb(r1(1)), 1810 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1811 1812 go to 2 1813 1814* ***** normalize psi **** 1815 4 call Pack_cc_dot(1, 1816 > dcpl_mb(psi_ptr+(i-1)*npack1), 1817 > dcpl_mb(psi_ptr+(i-1)*npack1), 1818 > sigma) 1819c call Pack_c_SMul(1,(1.0d0/dsqrt(sigma)), 1820c > dcpl_mb(psi_ptr+(i-1)*npack1), 1821c > dcpl_mb(psi_ptr+(i-1)*npack1)) 1822 call Pack_c_SMul1(1,(1.0d0/dsqrt(sigma)), 1823 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1824 1825 1826* **** release stack memory **** 1827 value = BA_pop_stack(g1(2)) 1828 value = value.and.BA_pop_stack(r1(2)) 1829 if (.not. value) call errquit( 1830 > 'psi_KS_update_orb2: popping stack memory',1, MA_ERR) 1831 error_out = dabs(e0-eold) 1832 1833c write(*,*) "i=",i,"iterations=",it," eig=",e0, 1834c > " error=",error_out, 1835c > lambda 1836 return 1837 end 1838 1839 1840 1841 1842 1843 1844 1845 1846* *********************************** 1847* * * 1848* * psi_linesearch_update * 1849* * * 1850* *********************************** 1851 1852* This routine performs a linesearch on orbital i, in the direction t. 1853* This routine is needed for a KS minimizer. 1854* e0 = <orb|g> 1855* de0 = 2*<t|g> 1856* 1857 subroutine psi_linesearch_update(psi_number,i,theta,e0,de0,t) 1858 implicit none 1859#include "errquit.fh" 1860 integer psi_number 1861 integer i 1862 real*8 theta 1863 real*8 e0,de0 1864 complex*16 t(*) !search direction 1865 1866#include "bafdecls.fh" 1867#include "psi.fh" 1868 1869* **** local variables **** 1870 logical value 1871 integer orb(2),g(2),psi_ptr 1872 real*8 x,y,pi,e1 1873 1874 if (psi_number.eq.1) then 1875 psi_ptr=psi1(1) 1876 else 1877 psi_ptr=psi2(1) 1878 end if 1879 1880 pi = 4.0d0*datan(1.0d0) 1881 1882* **** allocate stack memory **** 1883 value = BA_push_get(mt_dcpl,npack1,'orb', 1884 > orb(2),orb(1)) 1885 value = value.and. 1886 > BA_push_get(mt_dcpl,npack1,'g', 1887 > g(2),g(1)) 1888 if (.not. value) call errquit( 1889 > 'psi_linesearch_update: out of stack memory',0, MA_ERR) 1890 1891 1892 call Pack_c_Copy(1,dcpl_mb(psi_ptr+(i-1)*npack1), 1893 > dcpl_mb(orb(1))) 1894 1895* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 1896 !theta = pi/300.0d0 1897 x = cos(theta) 1898 y = sin(theta) 1899 call Pack_c_SMul(1,x, 1900 > dcpl_mb(orb(1)), 1901 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1902 call Pack_cc_daxpy(1,y, 1903 > t, 1904 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1905 1906* *** determine theta *** 1907 call psi_get_gradient_orb(psi_number,i,dcpl_mb(g(1))) 1908 1909 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(i-1)*npack1), 1910 > dcpl_mb(g(1)), 1911 > e1) 1912 e1 = -e1 1913 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 1914 > /(1.0d0-cos(2*theta)) 1915 theta = 0.5d0*datan(0.5d0*de0/x) 1916 1917c call Pack_cc_dot(1,t, 1918c > dcpl_mb(g(1)), 1919c > de1) 1920c de1 = -2.0d0*de1 1921c theta = -de1*(pi/300.0d0)/(de1-de0) 1922 1923 !write(*,*) "i,theta,e1:",i,theta,e1 1924 1925 1926* **** orb2 = orb*cos(theta) + t*sin(theta) **** 1927 x = cos(theta) 1928 y = sin(theta) 1929 call Pack_c_SMul(1,x, 1930 > dcpl_mb(orb(1)), 1931 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1932 call Pack_cc_daxpy(1,y, 1933 > t, 1934 > dcpl_mb(psi_ptr+(i-1)*npack1)) 1935 1936* **** update orb2_r and H*orb2 **** 1937 !call electron_run_orb(i,dcpl_mb(psi_ptr)) 1938c call psi_get_gradient_orb(psi_number,i,dcpl_mb(g(1))) 1939c call Pack_cc_dot(1,dcpl_mb(psi_ptr+(i-1)*npack1), 1940c > dcpl_mb(g(1)), 1941c > e2) 1942c e2 = -e2 1943c call Pack_cc_dot(1,t, 1944c > dcpl_mb(g(1)), 1945c > de2) 1946c de2 = -2.0d0*de2 1947 1948c write(*,*) "i,theta,es:",i,theta,e0,e1,e2 1949c write(*,*) 1950 1951* **** release stack memory **** 1952 value = BA_pop_stack(g(2)) 1953 value = value.and.BA_pop_stack(orb(2)) 1954 if (.not. value) call errquit( 1955 > 'psi_linesearch_update: popping stack memory',1, MA_ERR) 1956 1957 return 1958 end 1959 1960* *********************************** 1961* * * 1962* * psi_linesearch_update2 * 1963* * * 1964* *********************************** 1965 1966* This routine performs a linesearch on orbital i, in the direction t. 1967* This routine is needed for a KS minimizer. 1968* e0 = <orb|g> 1969* de0 = 2*<t|g> 1970* 1971 subroutine psi_linesearch_update2(psi_number,i,theta,e0,de0,t, 1972 > sigma,tau_t) 1973 implicit none 1974#include "errquit.fh" 1975 integer psi_number 1976 integer i 1977 real*8 theta 1978 real*8 e0,de0 1979 complex*16 t(*) !search direction 1980 1981 real*8 sigma 1982 complex*16 tau_t(*) !parallel transported search direction 1983 1984#include "bafdecls.fh" 1985#include "psi.fh" 1986 1987* **** local variables **** 1988 logical value 1989 integer orb(2),g(2),psi_ptr 1990 real*8 x,y,pi,e1 1991 1992 if (psi_number.eq.1) then 1993 psi_ptr=psi1(1) 1994 else 1995 psi_ptr=psi2(1) 1996 end if 1997 1998 pi = 4.0d0*datan(1.0d0) 1999 2000* **** allocate stack memory **** 2001 value = BA_push_get(mt_dcpl,npack1,'orb', 2002 > orb(2),orb(1)) 2003 value = value.and. 2004 > BA_push_get(mt_dcpl,npack1,'g', 2005 > g(2),g(1)) 2006 if (.not. value) call errquit( 2007 > 'psi_linesearch_update: out of stack memory',0, MA_ERR) 2008 2009 2010 call Pack_c_Copy(1,dcpl_mb(psi_ptr+(i-1)*npack1), 2011 > dcpl_mb(orb(1))) 2012 2013* **** orb2 = orb*cos(pi/300) + t*sin(pi/300) **** 2014 !theta = pi/300.0d0 2015 x = cos(theta) 2016 y = sin(theta) 2017 call Pack_c_SMul(1,x, 2018 > dcpl_mb(orb(1)), 2019 > dcpl_mb(psi_ptr+(i-1)*npack1)) 2020 call Pack_cc_daxpy(1,y, 2021 > t, 2022 > dcpl_mb(psi_ptr+(i-1)*npack1)) 2023 2024* *** determine theta *** 2025 call psi_get_gradient_orb(psi_number,i,dcpl_mb(g(1))) 2026 2027 call Pack_cc_dot(1,dcpl_mb(psi_ptr+(i-1)*npack1), 2028 > dcpl_mb(g(1)), 2029 > e1) 2030 e1 = -e1 2031 x = (e0 - e1 + 0.5d0*de0*sin(2*theta)) 2032 > /(1.0d0-cos(2*theta)) 2033 theta = 0.5d0*datan(0.5d0*de0/x) 2034 2035 x = cos(theta) 2036 y = sin(theta) 2037 2038* **** tau_t = (-orb*sin(theta) + t*cos(theta))*sigma **** 2039 call Pack_c_SMul(1,(-y), 2040 > dcpl_mb(orb(1)), 2041 > tau_t) 2042 call Pack_cc_daxpy(1,x, 2043 > t, 2044 > tau_t) 2045c call Pack_c_SMul(1,sigma, 2046c > tau_t, 2047c > tau_t) 2048 call Pack_c_SMul1(1,sigma,tau_t) 2049 2050* **** orb2 = orb*cos(theta) + t*sin(theta) **** 2051 call Pack_c_SMul(1,x, 2052 > dcpl_mb(orb(1)), 2053 > dcpl_mb(psi_ptr+(i-1)*npack1)) 2054 call Pack_cc_daxpy(1,y, 2055 > t, 2056 > dcpl_mb(psi_ptr+(i-1)*npack1)) 2057 2058 2059* **** release stack memory **** 2060 value = BA_pop_stack(g(2)) 2061 value = value.and.BA_pop_stack(orb(2)) 2062 if (.not. value) call errquit( 2063 > 'psi_linesearch_update: popping stack memory',1, MA_ERR) 2064 2065 return 2066 end 2067 2068 2069 2070* *************************** 2071* * * 2072* * psi_set_orb * 2073* * * 2074* *************************** 2075 2076* This routine copies an orbital, orb, into the ith psi of psi1. 2077* This routine is needed for a KS minimizer. 2078* 2079 subroutine psi_set_orb(psi_number,i,orb) 2080 implicit none 2081 integer psi_number 2082 integer i 2083 complex*16 orb(*) 2084 2085#include "bafdecls.fh" 2086#include "psi.fh" 2087 2088* **** local variables **** 2089 integer index,psi_ptr 2090 2091 if (psi_number.eq.1) then 2092 psi_ptr=psi1(1) 2093 else 2094 psi_ptr=psi2(1) 2095 end if 2096 2097 index = (i-1)*npack1 2098 2099 call zcopy(npack1, 2100 > orb, 1, 2101 > dcpl_mb(psi_ptr+index),1) 2102 return 2103 end 2104 2105 2106* *************************** 2107* * * 2108* * psi_get_orb * 2109* * * 2110* *************************** 2111 2112* This routine copies the ith psi of psi1 into an orbital, orb. 2113* This routine is needed for a KS minimizer. 2114* 2115 subroutine psi_get_orb(psi_number,i,orb) 2116 implicit none 2117 integer psi_number 2118 integer i 2119 complex*16 orb(*) 2120 2121#include "bafdecls.fh" 2122#include "psi.fh" 2123 2124* **** local variables **** 2125 integer index,psi_ptr 2126 2127 2128 if (psi_number.eq.1) then 2129 psi_ptr=psi1(1) 2130 else 2131 psi_ptr=psi2(1) 2132 end if 2133 2134 index = (i-1)*npack1 2135 2136 call zcopy(npack1, 2137 > dcpl_mb(psi_ptr+index), 1, 2138 > orb, 1) 2139 return 2140 end 2141 2142* *********************************** 2143* * * 2144* * psi_get_gradient_orb * 2145* * * 2146* *********************************** 2147 2148* This routine returns the Hpsi(i). 2149* This routine is needed for a KS minimizer. 2150* 2151 subroutine psi_get_gradient_orb(psi_number,i,Horb) 2152 implicit none 2153 integer psi_number 2154 integer i 2155 complex*16 Horb(*) 2156 2157#include "bafdecls.fh" 2158#include "psi.fh" 2159 2160* **** local variables **** 2161 integer psi_ptr 2162 2163 if (psi_number.eq.1) then 2164 psi_ptr=psi1(1) 2165 else 2166 psi_ptr=psi2(1) 2167 end if 2168 2169 call electron_run_orb(i,dcpl_mb(psi_ptr)) 2170 call electron_get_gradient_orb(i,Horb) 2171 2172 return 2173 end 2174 2175 2176* ******************************************* 2177* * * 2178* * psi_project_out_orb * 2179* * * 2180* ******************************************* 2181* 2182* This routine projects out non-orthogonal components of Horb. 2183* This routine is needed for a KS minimizer. 2184* 2185 subroutine psi_project_out_orb(psi_number,i,Horb) 2186 implicit none 2187#include "errquit.fh" 2188 integer psi_number 2189 integer i 2190 complex*16 Horb(*) 2191 2192#include "bafdecls.fh" 2193#include "psi.fh" 2194 2195* **** local variables **** 2196 logical ok 2197 integer ii,n,psi_ptr,np 2198 integer x(2) 2199 real*8 sum 2200 2201 call Parallel_np(np) 2202 2203* **** allocate stack memory **** 2204 ok = BA_push_get(mt_dbl,ne(1),'x',x(2),x(1)) 2205 if (.not.ok) 2206 > call errquit('psi_project_out_orb: out of stack memory',0, 2207 & MA_ERR) 2208 2209 2210 if (psi_number.eq.1) then 2211 psi_ptr=psi1(1) 2212 else 2213 psi_ptr=psi2(1) 2214 end if 2215 2216* **** spin up orbital **** 2217 if (i.le.ne(1)) then 2218 2219 ii = i 2220! do n=1,(ii) 2221! call Pack_cc_dot(1, 2222! > dcpl_mb(psi_ptr +(n-1)*npack1), 2223! > Horb, 2224! > sum) 2225! call daxpy(2*npack1, 2226! > (-sum), 2227! > dcpl_mb(psi_ptr+(n-1)*npack1),1, 2228! > Horb,1) 2229! end do 2230 call Pack_cc_ndot(1,ii, 2231 > dcpl_mb(psi_ptr), 2232 > Horb, 2233 > dbl_mb(x(1))) 2234 do n=1,(ii) 2235 call daxpy(2*npack1, 2236 > (-dbl_mb(x(1)+n-1)), 2237 > dcpl_mb(psi_ptr+(n-1)*npack1),1, 2238 > Horb,1) 2239 end do 2240 2241 2242 2243* **** spin down orbital **** 2244 else 2245 2246 2247 ii = i - ne(1) 2248 do n=(ne(1)+1),(ne(1)+ii) 2249 call Pack_cc_dot(1, 2250 > dcpl_mb(psi_ptr +(n-1)*npack1), 2251 > Horb, 2252 > sum) 2253 call daxpy(2*npack1, 2254 > (-sum), 2255 > dcpl_mb(psi_ptr+(n-1)*npack1),1, 2256 > Horb,1) 2257 end do 2258 2259 2260 end if 2261 2262* **** release stack memory **** 2263 ok = BA_pop_stack(x(2)) 2264 if (.not. ok) 2265 > call errquit('psi_project_out_orb: poping stack memory',0, 2266 & MA_ERR) 2267 2268 return 2269 end 2270 2271 2272 2273 2274 2275 2276* *************************** 2277* * * 2278* * psi_set_density * 2279* * * 2280* *************************** 2281 2282* This routine sets the densities and potentials in psi and electron. 2283* This routine is needed for a band by band minimizer. 2284* 2285 subroutine psi_set_density(psi_number,rho) 2286 implicit none 2287 integer psi_number 2288 real*8 rho(*) 2289 2290 2291#include "bafdecls.fh" 2292#include "psi.fh" 2293 2294* ***** rhoall common block **** 2295 integer rho1_all(2) 2296 integer rho2_all(2) 2297 common / rhoall_block / rho1_all,rho2_all 2298 2299* **** local variables **** 2300 integer rho_ptr,dng_ptr,rho_all_ptr 2301 2302 if (psi_number.eq.1) then 2303 rho_ptr = rho1(1) 2304 dng_ptr = dng1(1) 2305 rho_all_ptr = rho1_all(1) 2306 else 2307 rho_ptr = rho2(1) 2308 dng_ptr = dng2(1) 2309 rho_all_ptr = rho2_all(1) 2310 end if 2311 2312c call dcopy(4*nfft3d, 2313c > rho, 1, 2314c > dbl_mb(rho_ptr),1) 2315 2316 call Parallel_shared_vector_copy(.true.,4*nfft3d, 2317 > rho,dbl_mb(rho_ptr)) 2318 2319 call electron_gen_dng_dnall(dbl_mb(rho_ptr), 2320 > dcpl_mb(dng_ptr), 2321 > dbl_mb(rho_all_ptr)) 2322 call electron_gen_scf_potentials(dbl_mb(rho_ptr), 2323 > dcpl_mb(dng_ptr), 2324 > dbl_mb(rho_all_ptr)) 2325 call electron_gen_vall() 2326 return 2327 end 2328 2329 2330* *************************** 2331* * * 2332* * psi_get_density * 2333* * * 2334* *************************** 2335 2336* This routine gets the densities in psi. 2337* This routine is needed for a band by band minimizer. 2338* 2339 subroutine psi_get_density(psi_number,rho) 2340 implicit none 2341 integer psi_number 2342 real*8 rho(*) 2343 2344 2345#include "bafdecls.fh" 2346#include "psi.fh" 2347 2348* ***** rhoall common block **** 2349 integer rho1_all(2) 2350 integer rho2_all(2) 2351 common / rhoall_block / rho1_all,rho2_all 2352 2353* **** local variables **** 2354 integer rho_ptr 2355 2356 if (psi_number.eq.1) then 2357 rho_ptr = rho1(1) 2358 else 2359 rho_ptr = rho2(1) 2360 end if 2361 2362c call dcopy(4*nfft3d, 2363c > dbl_mb(rho_ptr),1, 2364c > rho,1) 2365 call Parallel_shared_vector_copy(.true.,4*nfft3d, 2366 > dbl_mb(rho_ptr),rho) 2367 return 2368 end 2369 2370* *************************** 2371* * * 2372* * psi_write_density * 2373* * * 2374* *************************** 2375 2376* This routine writes the densities in psi to disk. 2377* This routine is needed for a band by band minimizer. 2378* 2379 subroutine psi_write_density(psi_number) 2380 implicit none 2381 integer psi_number 2382 2383#include "bafdecls.fh" 2384#include "psi.fh" 2385 2386* ***** rhoall common block **** 2387 integer rho1_all(2) 2388 integer rho2_all(2) 2389 common / rhoall_block / rho1_all,rho2_all 2390 2391* **** local variables **** 2392 integer rho_ptr 2393 2394 if (psi_number.eq.1) then 2395 rho_ptr = rho1(1) 2396 else 2397 rho_ptr = rho2(1) 2398 end if 2399 2400 call rho_write(ispin,dbl_mb(rho_ptr)) 2401 return 2402 end 2403 2404* *************************** 2405* * * 2406* * psi_try_read_density * 2407* * * 2408* *************************** 2409 2410* This routine reads the densities from disk to psi. 2411* This routine is needed for a band by band minimizer. 2412* 2413 logical function psi_try_read_density(psi_number) 2414 implicit none 2415 integer psi_number 2416 2417#include "bafdecls.fh" 2418#include "psi.fh" 2419 2420* ***** rhoall common block **** 2421 integer rho1_all(2) 2422 integer rho2_all(2) 2423 common / rhoall_block / rho1_all,rho2_all 2424 2425* **** local variables **** 2426 logical ok 2427 integer rho_ptr 2428 2429* **** external functions **** 2430 logical rho_check_header 2431 external rho_check_header 2432 2433 if (psi_number.eq.1) then 2434 rho_ptr = rho1(1) 2435 else 2436 rho_ptr = rho2(1) 2437 end if 2438 2439 if (rho_check_header(ispin,.false.)) then 2440 call rho_read(ispin,dbl_mb(rho_ptr)) 2441 ok = .true. 2442 else 2443 ok = .false. 2444 end if 2445 2446 psi_try_read_density = ok 2447 return 2448 end 2449 2450 2451 2452* ************************************** 2453* * * 2454* * psi_gen_density_potentials * 2455* * * 2456* ************************************** 2457 2458* This routine sets the densities and potentials in psi and electron. 2459* This routine is needed for a band by band minimizer. 2460* 2461 subroutine psi_gen_density_potentials(psi_number0) 2462 implicit none 2463 integer psi_number0 2464 2465 2466#include "bafdecls.fh" 2467#include "psi.fh" 2468 2469* ***** rhoall common block **** 2470 integer rho1_all(2) 2471 integer rho2_all(2) 2472 common / rhoall_block / rho1_all,rho2_all 2473 2474* **** local variables **** 2475 logical zeroscf 2476 integer psi_ptr,rho_ptr,dng_ptr,rho_all_ptr,occ_ptr 2477 integer psi_number 2478 2479* *** hacky for now *** 2480 psi_number = psi_number0 2481 zeroscf = .false. 2482 if (psi_number.gt.2) then 2483 psi_number = psi_number - 2 2484 zeroscf = .true. 2485 end if 2486 !write(*,*) "psi_number,zeroscf=",psi_number,zeroscf 2487 2488 if (psi_number.eq.1) then 2489 psi_ptr = psi1(1) 2490 rho_ptr = rho1(1) 2491 dng_ptr = dng1(1) 2492 rho_all_ptr = rho1_all(1) 2493 occ_ptr = occ1(1) 2494 else 2495 psi_ptr = psi2(1) 2496 rho_ptr = rho2(1) 2497 dng_ptr = dng2(1) 2498 rho_all_ptr = rho2_all(1) 2499 occ_ptr = occ2(1) 2500 end if 2501 2502 2503 if (zeroscf) then 2504 call electron_gen_vall0() 2505 else 2506 call electron_gen_densities(dcpl_mb(psi_ptr), 2507 > dbl_mb(rho_ptr), 2508 > dcpl_mb(dng_ptr), 2509 > dbl_mb(rho_all_ptr), 2510 > occupation_on,dbl_mb(occ_ptr)) 2511 call electron_gen_scf_potentials(dbl_mb(rho_ptr), 2512 > dcpl_mb(dng_ptr), 2513 > dbl_mb(rho_all_ptr)) 2514 !if (zeroscf) then 2515 ! call electron_zero_scf_potentials() 2516 !end if 2517 call electron_gen_vall() 2518 end if 2519 return 2520 end 2521 2522 2523************************ Grasmman orbitals Part ************************ 2524 2525* *************************** 2526* * * 2527* * psi_1to2 * 2528* * * 2529* *************************** 2530 subroutine psi_1to2() 2531 implicit none 2532 2533#include "bafdecls.fh" 2534#include "psi.fh" 2535 2536c call zcopy(npack1*(neq(1)+neq(2)), 2537c > dcpl_mb(psi1(1)),1, 2538c > dcpl_mb(psi2(1)),1) 2539 call Parallel_shared_vector_copy(.true.,2*npack1*(neq(1)+neq(2)), 2540 > dcpl_mb(psi1(1)), 2541 > dcpl_mb(psi2(1))) 2542 2543 return 2544 end 2545 2546 2547* *************************** 2548* * * 2549* * psi_2to1 * 2550* * * 2551* *************************** 2552 subroutine psi_2to1() 2553 implicit none 2554 2555#include "bafdecls.fh" 2556#include "psi.fh" 2557 2558 2559c call zcopy(npack1*(neq(1)+neq(2)), 2560c > dcpl_mb(psi2(1)),1, 2561c > dcpl_mb(psi1(1)),1) 2562 call Parallel_shared_vector_copy(.true.,2*npack1*(neq(1)+neq(2)), 2563 > dcpl_mb(psi2(1)), 2564 > dcpl_mb(psi1(1))) 2565 2566c call OrthoCheck(ispin,ne,dcpl_mb(psi1(1))) 2567 return 2568 end 2569 2570 2571* *************************** 2572* * * 2573* * epsi_2to1 * 2574* * * 2575* *************************** 2576 subroutine epsi_2to1() 2577 implicit none 2578 2579#include "bafdecls.fh" 2580#include "psi.fh" 2581 2582 call zcopy(npack1*(ne_excited(1)+ne_excited(2)), 2583 > dcpl_mb(psi2_excited(1)),1, 2584 > dcpl_mb(psi1_excited(1)),1) 2585 return 2586 end 2587 2588 2589* *************************** 2590* * * 2591* * epsi_1to2 * 2592* * * 2593* *************************** 2594 subroutine epsi_1to2() 2595 implicit none 2596 2597#include "bafdecls.fh" 2598#include "psi.fh" 2599 2600 call zcopy(npack1*(ne_excited(1)+ne_excited(2)), 2601 > dcpl_mb(psi1_excited(1)),1, 2602 > dcpl_mb(psi2_excited(1)),1) 2603 return 2604 end 2605 2606 2607 2608* *************************** 2609* * * 2610* * psi_1get_psi * 2611* * * 2612* *************************** 2613 subroutine psi_1get_psi(rpsi) 2614 implicit none 2615 complex*16 rpsi(*) 2616 2617#include "bafdecls.fh" 2618#include "psi.fh" 2619 2620 call zcopy(npack1*(neq(1)+neq(2)), 2621 > dcpl_mb(psi1(1)),1, 2622 > rpsi,1) 2623 2624 return 2625 end 2626 2627 2628* *************************** 2629* * * 2630* * psi_2get_psi * 2631* * * 2632* *************************** 2633 subroutine psi_2get_psi(rpsi) 2634 implicit none 2635 complex*16 rpsi(*) 2636 2637#include "bafdecls.fh" 2638#include "psi.fh" 2639 2640 call zcopy(npack1*(neq(1)+neq(2)), 2641 > dcpl_mb(psi2(1)),1, 2642 > rpsi,1) 2643 2644 return 2645 end 2646 2647* *************************** 2648* * * 2649* * psi_check * 2650* * * 2651* *************************** 2652 subroutine psi_check() 2653 implicit none 2654 2655#include "bafdecls.fh" 2656#include "psi.fh" 2657 2658 2659 call OrthoCheck(ispin,ne,dcpl_mb(psi1(1))) 2660 return 2661 end 2662 2663 2664 2665* *************************** 2666* * * 2667* * rho_2to1 * 2668* * * 2669* *************************** 2670 subroutine rho_2to1() 2671 implicit none 2672 2673#include "bafdecls.fh" 2674#include "psi.fh" 2675 2676* ***** rhoall common block **** 2677 integer rho1_all(2) 2678 integer rho2_all(2) 2679 common / rhoall_block / rho1_all,rho2_all 2680 2681c call dcopy(4*nfft3d, 2682c > dbl_mb(rho2(1)),1, 2683c > dbl_mb(rho1(1)),1) 2684 call Parallel_shared_vector_copy(.true.,4*nfft3d, 2685 > dbl_mb(rho2(1)), 2686 > dbl_mb(rho1(1))) 2687 2688c call dcopy(4*nfft3d, 2689c > dbl_mb(rho2_all(1)),1, 2690c > dbl_mb(rho1_all(1)),1) 2691 call Parallel_shared_vector_copy(.true.,4*nfft3d, 2692 > dbl_mb(rho2_all(1)), 2693 > dbl_mb(rho1_all(1))) 2694 2695 return 2696 end 2697 2698* *************************** 2699* * * 2700* * rho_1to2 * 2701* * * 2702* *************************** 2703 subroutine rho_1to2() 2704 implicit none 2705 2706#include "bafdecls.fh" 2707#include "psi.fh" 2708 2709* ***** rhoall common block **** 2710 integer rho1_all(2) 2711 integer rho2_all(2) 2712 common / rhoall_block / rho1_all,rho2_all 2713 2714 call dcopy(4*nfft3d, 2715 > dbl_mb(rho1(1)),1, 2716 > dbl_mb(rho2(1)),1) 2717 2718 call dcopy(4*nfft3d, 2719 > dbl_mb(rho1_all(1)),1, 2720 > dbl_mb(rho2_all(1)),1) 2721 2722 return 2723 end 2724 2725* *************************** 2726* * * 2727* * dng_2to1 * 2728* * * 2729* *************************** 2730 subroutine dng_2to1() 2731 implicit none 2732 2733#include "bafdecls.fh" 2734#include "psi.fh" 2735 2736 call zcopy(npack0, 2737 > dcpl_mb(dng2(1)),1, 2738 > dcpl_mb(dng1(1)),1) 2739 2740 return 2741 end 2742 2743* *************************** 2744* * * 2745* * dng_1to2 * 2746* * * 2747* *************************** 2748 subroutine dng_1to2() 2749 implicit none 2750 2751#include "bafdecls.fh" 2752#include "psi.fh" 2753 2754 call zcopy(npack0, 2755 > dcpl_mb(dng1(1)),1, 2756 > dcpl_mb(dng2(1)),1) 2757 2758 return 2759 end 2760 2761 2762 2763* *********************************** 2764* * * 2765* * psi_1add_oep_to_vall * 2766* * * 2767* *********************************** 2768 subroutine psi_1add_oep_to_vall() 2769 implicit none 2770 2771#include "bafdecls.fh" 2772#include "psi.fh" 2773 2774 2775 call electron_add_oep_to_vall(dbl_mb(rho1(1))) 2776 2777 return 2778 end 2779 2780 2781* *********************************** 2782* * * 2783* * psi_1toelectron * 2784* * * 2785* *********************************** 2786 subroutine psi_1toelectron() 2787 implicit none 2788 2789#include "bafdecls.fh" 2790#include "psi.fh" 2791 2792* ***** rhoall common block **** 2793 integer rho1_all(2) 2794 integer rho2_all(2) 2795 common / rhoall_block / rho1_all,rho2_all 2796 2797 call electron_run(dcpl_mb(psi1(1)), 2798 > dbl_mb(rho1(1)), 2799 > dcpl_mb(dng1(1)), 2800 > dbl_mb(rho1_all(1)), 2801 > occupation_on,dbl_mb(occ1(1))) 2802 2803 return 2804 end 2805 2806* *********************************** 2807* * * 2808* * psi_1genrho * 2809* * * 2810* *********************************** 2811 subroutine psi_1genrho() 2812 implicit none 2813 2814#include "bafdecls.fh" 2815#include "psi.fh" 2816 2817* ***** rhoall common block **** 2818 integer rho1_all(2) 2819 integer rho2_all(2) 2820 common / rhoall_block / rho1_all,rho2_all 2821 2822 call electron_genrho(dcpl_mb(psi1(1)), 2823 > dbl_mb(rho1(1)), 2824 > occupation_on,dbl_mb(occ1(1))) 2825 2826 return 2827 end 2828 2829 2830 2831* *********************************** 2832* * * 2833* * psi_1energy * 2834* * * 2835* *********************************** 2836 real*8 function psi_1energy() 2837 implicit none 2838 2839#include "bafdecls.fh" 2840#include "psi.fh" 2841 2842* ***** rhoall common block **** 2843 integer rho1_all(2) 2844 integer rho2_all(2) 2845 common / rhoall_block / rho1_all,rho2_all 2846 2847* **** external functions **** 2848 real*8 electron_energy 2849 external electron_energy 2850 2851 call electron_run(dcpl_mb(psi1(1)), 2852 > dbl_mb(rho1(1)), 2853 > dcpl_mb(dng1(1)), 2854 > dbl_mb(rho1_all(1)), 2855 > occupation_on,dbl_mb(occ1(1))) 2856 psi_1energy = electron_energy(dcpl_mb(psi1(1)), 2857 > dbl_mb(rho1(1)), 2858 > dcpl_mb(dng1(1)), 2859 > dbl_mb(rho1_all(1)), 2860 > occupation_on,dbl_mb(occ1(1))) 2861 2862 return 2863 end 2864 2865* *********************************** 2866* * * 2867* * psi_1_noupdate_energy * 2868* * * 2869* *********************************** 2870 real*8 function psi_1_noupdate_energy() 2871 implicit none 2872 2873#include "bafdecls.fh" 2874#include "psi.fh" 2875 2876* ***** rhoall common block **** 2877 integer rho1_all(2) 2878 integer rho2_all(2) 2879 common / rhoall_block / rho1_all,rho2_all 2880 2881* **** external functions **** 2882 real*8 electron_energy 2883 external electron_energy 2884 2885 2886 !call electron_gen_Hpsi_k(dcpl_mb(psi1(1))) 2887 psi_1_noupdate_energy = electron_energy(dcpl_mb(psi1(1)), 2888 > dbl_mb(rho1(1)), 2889 > dcpl_mb(dng1(1)), 2890 > dbl_mb(rho1_all(1)), 2891 > occupation_on,dbl_mb(occ1(1))) 2892 2893 return 2894 end 2895 2896 2897* *********************************** 2898* * * 2899* * psi_2energy * 2900* * * 2901* *********************************** 2902 real*8 function psi_2energy() 2903 implicit none 2904 2905#include "bafdecls.fh" 2906#include "psi.fh" 2907 2908* ***** rhoall common block **** 2909 integer rho1_all(2) 2910 integer rho2_all(2) 2911 common / rhoall_block / rho1_all,rho2_all 2912 2913* **** external functions **** 2914 real*8 electron_energy 2915 external electron_energy 2916 2917 call electron_run(dcpl_mb(psi2(1)), 2918 > dbl_mb(rho2(1)), 2919 > dcpl_mb(dng2(1)), 2920 > dbl_mb(rho2_all(1)), 2921 > occupation_on,dbl_mb(occ2(1))) 2922 psi_2energy = electron_energy(dcpl_mb(psi2(1)), 2923 > dbl_mb(rho2(1)), 2924 > dcpl_mb(dng2(1)), 2925 > dbl_mb(rho2_all(1)), 2926 > occupation_on,dbl_mb(occ2(1))) 2927 2928 return 2929 end 2930 2931 2932 2933* *********************************** 2934* * * 2935* * psi_1eorbit * 2936* * * 2937* *********************************** 2938 real*8 function psi_1eorbit() 2939 implicit none 2940 2941#include "bafdecls.fh" 2942#include "psi.fh" 2943 2944* **** external functions **** 2945 real*8 electron_eorbit 2946 external electron_eorbit 2947 2948 psi_1eorbit = electron_eorbit(dcpl_mb(psi1(1)), 2949 > occupation_on,dbl_mb(occ1(1))) 2950 2951 return 2952 end 2953 2954 2955* *********************************** 2956* * * 2957* * psi_1ke * 2958* * * 2959* *********************************** 2960 real*8 function psi_1ke() 2961 implicit none 2962 2963#include "bafdecls.fh" 2964#include "psi.fh" 2965 2966* **** local variables *** 2967 real*8 ave 2968 2969 call ke_ave(ispin,neq,dcpl_mb(psi1(1)),ave, 2970 > occupation_on,dbl_mb(occ1(1))) 2971 2972 psi_1ke = ave 2973 return 2974 end 2975 2976 2977 2978 2979 2980* *********************************** 2981* * * 2982* * psi_1ke_atom * 2983* * * 2984* *********************************** 2985 real*8 function psi_1ke_atom() 2986 implicit none 2987 2988#include "bafdecls.fh" 2989#include "psi.fh" 2990 2991* **** local variables *** 2992 real*8 ave 2993 2994* **** external functions **** 2995 real*8 psp_kinetic_atom 2996 external psp_kinetic_atom 2997 2998 ave = psp_kinetic_atom(ispin,neq,dcpl_mb(psi1(1))) 2999 3000 psi_1ke_atom = ave 3001 return 3002 end 3003 3004* *********************************** 3005* * * 3006* * psi_1valence_core_atom * 3007* * * 3008* *********************************** 3009 real*8 function psi_1valence_core_atom() 3010 implicit none 3011 3012#include "bafdecls.fh" 3013#include "psi.fh" 3014 3015* **** local variables *** 3016 real*8 ave 3017 3018* **** external functions **** 3019 real*8 psp_valence_core_atom 3020 external psp_valence_core_atom 3021 3022 ave = psp_valence_core_atom(ispin,neq,dcpl_mb(psi1(1))) 3023 3024 psi_1valence_core_atom = ave 3025 return 3026 end 3027 3028 3029 3030* *********************************** 3031* * * 3032* * psi_1vloc_atom * 3033* * * 3034* *********************************** 3035 real*8 function psi_1vloc_atom() 3036 implicit none 3037 3038#include "bafdecls.fh" 3039#include "psi.fh" 3040 3041* **** local variables *** 3042 real*8 ave 3043 3044* **** external functions **** 3045 real*8 psp_vloc_atom 3046 external psp_vloc_atom 3047 3048 ave = psp_vloc_atom(ispin,neq,dcpl_mb(psi1(1))) 3049 3050 psi_1vloc_atom = ave 3051 return 3052 end 3053 3054 3055 3056* *********************************** 3057* * * 3058* * psi_1ncmp_vloc * 3059* * * 3060* *********************************** 3061 real*8 function psi_1ncmp_vloc() 3062 implicit none 3063 3064#include "bafdecls.fh" 3065#include "psi.fh" 3066 3067* **** local variables *** 3068 real*8 ave 3069 3070* **** external functions **** 3071 real*8 psp_ncmp_vloc 3072 external psp_ncmp_vloc 3073 3074 ave = psp_ncmp_vloc(ispin) 3075 3076 psi_1ncmp_vloc = ave 3077 return 3078 end 3079 3080 3081 3082* *********************************** 3083* * * 3084* * psi_1hartree_atom * 3085* * * 3086* *********************************** 3087 real*8 function psi_1hartree_atom() 3088 implicit none 3089 3090#include "bafdecls.fh" 3091#include "psi.fh" 3092 3093* **** local variables *** 3094 real*8 ave 3095 3096* **** external functions **** 3097 real*8 psp_hartree_atom 3098 external psp_hartree_atom 3099 3100 ave = psp_hartree_atom(ispin,neq,dcpl_mb(psi1(1))) 3101 3102 psi_1hartree_atom = ave 3103 return 3104 end 3105 3106 3107* *********************************** 3108* * * 3109* * psi_1hartree_cmp_cmp * 3110* * * 3111* *********************************** 3112 real*8 function psi_1hartree_cmp_cmp() 3113 implicit none 3114 3115#include "psi.fh" 3116 3117* **** external functions **** 3118 real*8 psp_hartree_cmp_cmp 3119 external psp_hartree_cmp_cmp 3120 3121 psi_1hartree_cmp_cmp = psp_hartree_cmp_cmp(ispin) 3122 return 3123 end 3124 3125* *********************************** 3126* * * 3127* * psi_1hartree_cmp_pw * 3128* * * 3129* *********************************** 3130 real*8 function psi_1hartree_cmp_pw() 3131 implicit none 3132 3133#include "bafdecls.fh" 3134#include "psi.fh" 3135 3136* **** external functions **** 3137 real*8 psp_hartree_cmp_pw 3138 external psp_hartree_cmp_pw 3139 3140 psi_1hartree_cmp_pw = psp_hartree_cmp_pw(ispin,dcpl_mb(dng1(1)), 3141 > dbl_mb(rho1(1))) 3142 return 3143 end 3144 3145 3146c* *********************************** 3147c* * * 3148c* * dng_1vlpaw_pw * 3149c* * * 3150c* *********************************** 3151c real*8 function dng_1vlpaw_pw() 3152c implicit none 3153c 3154c#include "bafdecls.fh" 3155c#include "psi.fh" 3156c 3157c* **** external functions **** 3158c real*8 electron_dng_vlpaw_ave 3159c external electron_dng_vlpaw_ave 3160c 3161c dng_1vlpaw_pw = electron_dng_vlpaw_ave(dcpl_mb(dng1(1))) 3162c 3163c return 3164c end 3165 3166 3167 3168* *********************************** 3169* * * 3170* * psi_1xc_atom * 3171* * * 3172* *********************************** 3173 subroutine psi_1xc_atom(exc,pxc) 3174 implicit none 3175 real*8 exc,pxc 3176 3177#include "bafdecls.fh" 3178#include "psi.fh" 3179 3180 call psp_xc_atom(ispin,neq,dcpl_mb(psi1(1)),exc,pxc) 3181 !call D1dB_SumAll(exc) 3182 !call D1dB_SumAll(pxc) 3183 return 3184 end 3185 3186 3187* *********************************** 3188* * * 3189* * psi_1qlm_atom * 3190* * * 3191* *********************************** 3192 subroutine psi_1qlm_atom() 3193 implicit none 3194 3195#include "bafdecls.fh" 3196#include "psi.fh" 3197 3198 call psp_qlm_atom(ispin,neq,dcpl_mb(psi1(1))) 3199 return 3200 end 3201 3202 3203 3204 3205* *********************************** 3206* * * 3207* * psi_1vl * 3208* * * 3209* *********************************** 3210 real*8 function psi_1vl() 3211 implicit none 3212 3213#include "bafdecls.fh" 3214#include "psi.fh" 3215 3216 3217* **** external functions **** 3218 real*8 electron_psi_vl_ave 3219 external electron_psi_vl_ave 3220 3221 psi_1vl = electron_psi_vl_ave(dcpl_mb(psi1(1)),dbl_mb(rho1(1))) 3222 3223 return 3224 end 3225 3226 3227* *********************************** 3228* * * 3229* * dng_1vl_mm * 3230* * * 3231* *********************************** 3232 real*8 function dng_1vl_mm() 3233 implicit none 3234 3235#include "bafdecls.fh" 3236#include "errquit.fh" 3237#include "psi.fh" 3238 3239* **** local variables **** 3240 integer nx,ny,nz,n2ft3d,vltmp(2),r_grid(2) 3241 real*8 elocal,esum,dv 3242 3243* **** external functions **** 3244 integer control_version 3245 external control_version 3246 real*8 lattice_omega 3247 external lattice_omega 3248 3249 call D3dB_n2ft3d(1,n2ft3d) 3250 call D3dB_nx(1,nx) 3251 call D3dB_ny(1,ny) 3252 call D3dB_nz(1,nz) 3253 dv = lattice_omega()/dble(nx*ny*nz) 3254 3255 if (.not.BA_push_get(mt_dbl,(n2ft3d),'vltmp',vltmp(2),vltmp(1))) 3256 > call errquit('dng_1_vl_mm: out of stack memory',0,MA_ERR) 3257 3258* **** average Kohn-Sham v_local energy **** 3259 call v_local_mm(dbl_mb(vltmp(1))) 3260 call Pack_cc_dot(0,dcpl_mb(dng1(1)),dbl_mb(vltmp(1)),elocal) 3261 3262* *** add in long range part **** 3263 if (control_version().eq.4) then 3264 if (.not.BA_push_get(mt_dbl,(3*n2ft3d),'r_grid', 3265 > r_grid(2),r_grid(1))) 3266 > call errquit('dng_1_vl_mm: out of stack memory',1,MA_ERR) 3267 call lattice_r_grid(dbl_mb(r_grid(1))) 3268 3269 call v_lr_local_mm(dbl_mb(r_grid(1)),dbl_mb(vltmp(1))) 3270 call D3dB_rr_dot(1,dbl_mb(rho1(1)),dbl_mb(vltmp(1)),esum) 3271 elocal = elocal + esum*dv 3272 if (.not.BA_pop_stack(r_grid(2))) 3273 > call errquit('dng_1_vl_mm: popping stack',0,MA_ERR) 3274 end if 3275 3276 if (.not.BA_pop_stack(vltmp(2))) 3277 > call errquit('dng_1_vl_mm: popping stack',1,MA_ERR) 3278 3279 3280 dng_1vl_mm = elocal 3281 return 3282 end 3283 3284 3285 3286 3287 3288* *********************************** 3289* * * 3290* * psi_1vnl * 3291* * * 3292* *********************************** 3293 real*8 function psi_1vnl() 3294 implicit none 3295 3296#include "bafdecls.fh" 3297#include "psi.fh" 3298 3299 3300* **** external functions **** 3301 real*8 electron_psi_vnl_ave 3302 external electron_psi_vnl_ave 3303 3304 psi_1vnl = electron_psi_vnl_ave(dcpl_mb(psi1(1)), 3305 > occupation_on,dbl_mb(occ1(1))) 3306 3307 return 3308 end 3309 3310* ******************************* 3311* * * 3312* * psi_1v_field * 3313* * * 3314* ******************************* 3315 real*8 function psi_1v_field() 3316 implicit none 3317 3318#include "bafdecls.fh" 3319#include "psi.fh" 3320 3321 3322 3323* **** external functions **** 3324 real*8 electron_psi_v_field_ave 3325 external electron_psi_v_field_ave 3326 3327 psi_1v_field = electron_psi_v_field_ave(dcpl_mb(psi1(1)), 3328 > dbl_mb(rho1(1))) 3329 3330 return 3331 end 3332 3333 3334* *********************************** 3335* * * 3336* * rho_1Fcharge * 3337* * * 3338* *********************************** 3339 subroutine rho_1Fcharge(Fcharge) 3340 implicit none 3341 real*8 Fcharge(*) 3342 3343#include "bafdecls.fh" 3344#include "psi.fh" 3345#include "errquit.fh" 3346 3347* **** local variables **** 3348 logical value 3349 integer n2ft3d,nx,ny,nz 3350 integer r_grid(2),rho(2) 3351 real*8 dv 3352 3353* **** external functions **** 3354 real*8 lattice_omega 3355 external lattice_omega 3356 3357* **** Initializationsr **** 3358 call D3dB_n2ft3d(1,n2ft3d) 3359 call D3dB_nx(1,nx) 3360 call D3dB_ny(1,ny) 3361 call D3dB_nz(1,nz) 3362 dv = lattice_omega()/dble(nx*ny*nz) 3363 3364 3365* **** Push memory **** 3366 value = BA_push_get(mt_dbl,(3*n2ft3d),'r_grid', 3367 > r_grid(2),r_grid(1)) 3368 value = value.and. 3369 > BA_push_get(mt_dbl,(3*n2ft3d),'rho', 3370 > rho(2),rho(1)) 3371 if (.not. value) call errquit( 3372 > 'rho_1Fcharge: out of stack memory',0, MA_ERR) 3373 call dcopy(3*n2ft3d,0.0d0,dbl_mb(rho(1)),1) 3374 3375 3376* **** Get r_grid and rho **** 3377 call lattice_r_grid(dbl_mb(r_grid(1))) 3378 call D3dB_rr_Sum(1,dbl_mb(rho1(1)), 3379 > dbl_mb(rho1(1)+(ispin-1)*n2ft3d), 3380 > dbl_mb(rho(1))) 3381 3382* **** Now calculate Fcharge **** 3383 call pspw_charge_rho_Fcharge(n2ft3d,dbl_mb(r_grid(1)), 3384 > dbl_mb(rho(1)), 3385 > dv,Fcharge) 3386 3387* **** Pop memory **** 3388 value = BA_pop_stack(rho(2)) 3389 value = value.and.BA_pop_stack(r_grid(2)) 3390 if (.not. value) call errquit( 3391 > 'rho_1Fcharge: error popping stack memory',0, MA_ERR) 3392 3393 return 3394 end 3395 3396 3397 3398* *********************************** 3399* * * 3400* * rho_1exc * 3401* * * 3402* *********************************** 3403 real*8 function rho_1exc() 3404 implicit none 3405 3406#include "bafdecls.fh" 3407#include "psi.fh" 3408 3409* ***** rhoall common block **** 3410 integer rho1_all(2) 3411 integer rho2_all(2) 3412 common / rhoall_block / rho1_all,rho2_all 3413 3414* **** external functions **** 3415 real*8 electron_exc 3416 external electron_exc 3417 3418 rho_1exc = electron_exc(dbl_mb(rho1_all(1))) 3419 return 3420 end 3421 3422* *********************************** 3423* * * 3424* * rho_1pxc * 3425* * * 3426* *********************************** 3427 real*8 function rho_1pxc() 3428 implicit none 3429 3430#include "bafdecls.fh" 3431#include "psi.fh" 3432 3433* **** external functions **** 3434 real*8 electron_pxc 3435 external electron_pxc 3436 3437 rho_1pxc = electron_pxc(dbl_mb(rho1(1))) 3438 return 3439 end 3440 3441 3442* *********************************** 3443* * * 3444* * psi_1meta_gga_pxc * 3445* * * 3446* *********************************** 3447 real*8 function psi_1meta_gga_pxc() 3448 implicit none 3449 3450#include "bafdecls.fh" 3451#include "psi.fh" 3452 3453* **** external functions **** 3454 real*8 nwpw_meta_gga_pxc 3455 external nwpw_meta_gga_pxc 3456 3457 psi_1meta_gga_pxc = nwpw_meta_gga_pxc(ispin,neq,dcpl_mb(psi1(1))) 3458 return 3459 end 3460 3461 3462 3463* *********************************** 3464* * * 3465* * dng_1ehartree * 3466* * * 3467* *********************************** 3468 real*8 function dng_1ehartree() 3469 implicit none 3470 3471#include "bafdecls.fh" 3472#include "psi.fh" 3473 3474* **** external functions **** 3475 integer control_version 3476 real*8 electron_ehartree,electron_ehartree2 3477 external control_version 3478 external electron_ehartree,electron_ehartree2 3479 3480* **** local variables ***** 3481 real*8 eh 3482 3483 eh = 0.0d0 3484 if (control_version().eq.3) 3485 > eh = electron_ehartree(dcpl_mb(dng1(1))) 3486 3487 if (control_version().eq.4) 3488 > eh = electron_ehartree2(dbl_mb(rho1(1))) 3489 3490 dng_1ehartree = eh 3491 return 3492 end 3493 3494 3495 3496* *********************************** 3497* * * 3498* * psi_2toelectron * 3499* * * 3500* *********************************** 3501 subroutine psi_2toelectron() 3502 implicit none 3503 3504#include "bafdecls.fh" 3505#include "psi.fh" 3506 3507* ***** rhoall common block **** 3508 integer rho1_all(2) 3509 integer rho2_all(2) 3510 common / rhoall_block / rho1_all,rho2_all 3511 3512 call electron_run(dcpl_mb(psi2(1)), 3513 > dbl_mb(rho2(1)), 3514 > dcpl_mb(dng2(1)), 3515 > dbl_mb(rho2_all(1)), 3516 > occupation_on,dbl_mb(occ2(1))) 3517 3518 return 3519 end 3520 3521 3522* *********************************** 3523* * * 3524* * psi_1check_Tangent * 3525* * * 3526* *********************************** 3527* 3528* This routine checks the accuracy of the tangent vector. 3529* MM = Yt*H = Yt*(I-Y*Yt)*G = Yt*G - Yt*Y*Yt*G = Yt*G - Yt*G == 0 3530 3531* Updated - 5-18-2002 3532* 3533 subroutine psi_1check_Tangent(H) 3534 implicit none 3535 complex*16 H(*) 3536 3537#include "errquit.fh" 3538#include "bafdecls.fh" 3539#include "psi.fh" 3540 3541* **** local variables **** 3542 logical value 3543 integer ms,n,indx,i,j 3544 integer MM(2) 3545 real*8 sum 3546 3547 do ms=1,ispin 3548 n = ne(ms) 3549 if (n.eq.0) go to 101 !*** ferromagnetic check *** 3550 value = BA_push_get(mt_dbl,n*n,'MM',MM(2),MM(1)) 3551 if (.not. value) 3552 > call errquit('out of stack memory in psi_1check_Tangent',0, 3553 & MA_ERR) 3554 3555 indx = (ms-1)*ne(1)*npack1 3556 3557* **** calculate MM = Yt*H **** 3558 call Grsm_ggm_dot(npack1,n, 3559 > dcpl_mb(psi1(1)+indx), 3560 > H(1+indx), 3561 > dbl_mb(MM(1))) 3562 3563* **** write out MM matrix **** 3564 sum = 0.0d0 3565 do j=1,n 3566 do i=1,n 3567 sum = sum + dbl_mb(MM(1)+(i-1)+(j-1)*n) 3568 end do 3569 end do 3570 write(*,*) "psi_1check_Tangent:",sum 3571 3572 3573 3574 value = BA_pop_stack(MM(2)) 3575 if (.not. value) 3576 > call errquit( 3577 > 'error popping stack memory in psi_1check_Tangent',0, 3578 > MA_ERR) 3579 3580 101 continue 3581 end do 3582 3583 return 3584 end 3585 3586 3587 3588* *********************************** 3589* * * 3590* * psi_2check_Tangent * 3591* * * 3592* *********************************** 3593* 3594* This routine checks the accuracy of the tangent vector. 3595* MM = Yt*H = Yt*(I-Y*Yt)*G = Yt*G - Yt*Y*Yt*G = Yt*G - Yt*G == 0 3596 3597* Updated - 5-18-2002 3598* 3599 subroutine psi_2check_Tangent(H) 3600 implicit none 3601 complex*16 H(*) 3602 3603#include "bafdecls.fh" 3604#include "psi.fh" 3605#include "errquit.fh" 3606 3607* **** local variables **** 3608 logical value 3609 integer ms,n,indx,i,j 3610 integer MM(2) 3611 real*8 sum 3612 3613 do ms=1,ispin 3614 n = ne(ms) 3615 if (n.eq.0) go to 101 !*** ferromagnetic check *** 3616 value = BA_push_get(mt_dbl,n*n,'MM',MM(2),MM(1)) 3617 if (.not. value) 3618 > call errquit('out of stack memory in psi_1check_Tangent',0, 3619 > MA_ERR) 3620 3621 indx = (ms-1)*ne(1)*npack1 3622 3623* **** calculate MM = Yt*H **** 3624 call Grsm_ggm_dot(npack1,n, 3625 > dcpl_mb(psi2(1)+indx), 3626 > H(1+indx), 3627 > dbl_mb(MM(1))) 3628 3629* **** write out MM matrix **** 3630 sum = 0.0d0 3631 do j=1,n 3632 do i=1,n 3633 sum = sum + dbl_mb(MM(1)+(i-1)+(j-1)*n) 3634 end do 3635 end do 3636 write(*,*) "psi_2check_Tangent:",sum 3637 3638 3639 3640 value = BA_pop_stack(MM(2)) 3641 if (.not. value) 3642 > call errquit( 3643 > 'error popping stack memory in psi_2check_Tangent',0, 3644 & MA_ERR) 3645 3646 101 continue 3647 end do 3648 3649 return 3650 end 3651 3652 3653 3654* *********************************** 3655* * * 3656* * psi_1get_Tgradient * 3657* * * 3658* *********************************** 3659 3660* THpsi = Hpsi - Y*Y^t*Hpsi ! used by Grassman minimizers 3661* 3662 subroutine psi_1get_Tgradient(THpsi,Eout) 3663 implicit none 3664 complex*16 THpsi(*) 3665 real*8 Eout 3666 3667#include "bafdecls.fh" 3668#include "errquit.fh" 3669#include "psi.fh" 3670 3671* ***** rhoall common block **** 3672 integer rho1_all(2) 3673 integer rho2_all(2) 3674 common / rhoall_block / rho1_all,rho2_all 3675 3676* **** local variables **** 3677 integer tmp1(2),i,n 3678 logical value 3679 3680* **** external functions **** 3681 logical Dneall_m_push_get,Dneall_m_pop_stack 3682 external Dneall_m_push_get,Dneall_m_pop_stack 3683 3684 real*8 electron_energy 3685 external electron_energy 3686 3687 if (.not.Dneall_m_push_get(0,tmp1)) 3688 > call errquit('out of stack memory in psi_1get_Tradient',0, 3689 > MA_ERR) 3690 3691 call electron_run(dcpl_mb(psi1(1)), 3692 > dbl_mb(rho1(1)), 3693 > dcpl_mb(dng1(1)), 3694 > dbl_mb(rho1_all(1)), 3695 > occupation_on,dbl_mb(occ1(1))) 3696 Eout = electron_energy(dcpl_mb(psi1(1)), 3697 > dbl_mb(rho1(1)), 3698 > dcpl_mb(dng1(1)), 3699 > dbl_mb(rho1_all(1)), 3700 > occupation_on,dbl_mb(occ1(1))) 3701 call electron_gen_hml(dcpl_mb(psi1(1)), 3702 > dbl_mb(tmp1(1))) 3703 3704* **** get the tangent THpsi = Hpsi - psi1*tmp1 , TY = HY - Y*Y'HY **** 3705 call electron_get_Tgradient(dcpl_mb(psi1(1)), 3706 > dbl_mb(tmp1(1)), 3707 > THpsi) 3708 3709 if (.not.Dneall_m_pop_stack(tmp1)) 3710 > call errquit('psi_1get_Tgradient:error popping stack',1, 3711 > MA_ERR) 3712 3713 return 3714 end 3715 3716 3717* *********************************** 3718* * * 3719* * psi_1get_STgradient * 3720* * * 3721* *********************************** 3722 3723* THpsi = Hpsi - Y*Y^t*Hpsi ! used by Grassman minimizers 3724* 3725 subroutine psi_1get_STgradient(Rpsi,THpsi,Eout) 3726 implicit none 3727 complex*16 Rpsi(*) 3728 complex*16 THpsi(*) 3729 real*8 Eout 3730 3731#include "bafdecls.fh" 3732#include "errquit.fh" 3733#include "psi.fh" 3734 3735* ***** rhoall common block **** 3736 integer rho1_all(2) 3737 integer rho2_all(2) 3738 common / rhoall_block / rho1_all,rho2_all 3739 3740* **** local variables **** 3741 integer tmp1(2),i,n 3742 logical value 3743 3744* **** external functions **** 3745 logical Dneall_m_push_get,Dneall_m_pop_stack 3746 external Dneall_m_push_get,Dneall_m_pop_stack 3747 3748 real*8 electron_energy 3749 external electron_energy 3750 3751 if (.not.Dneall_m_push_get(0,tmp1)) 3752 > call errquit('out of stack memory in psi_1get_STradient',0, 3753 > MA_ERR) 3754 3755 call electron_run(dcpl_mb(psi1(1)), 3756 > dbl_mb(rho1(1)), 3757 > dcpl_mb(dng1(1)), 3758 > dbl_mb(rho1_all(1)), 3759 > occupation_on,dbl_mb(occ1(1))) 3760 Eout = electron_energy(dcpl_mb(psi1(1)), 3761 > dbl_mb(rho1(1)), 3762 > dcpl_mb(dng1(1)), 3763 > dbl_mb(rho1_all(1)), 3764 > occupation_on,dbl_mb(occ1(1))) 3765 call electron_gen_hml(dcpl_mb(psi1(1)), 3766 > dbl_mb(tmp1(1))) 3767 3768* **** get the tangent Rpsi = Hpsi - Spsi1*tmp1 , TY = HY - SY*Y'HY **** 3769 call psp_overlap_S(ispin,neq,dcpl_mb(psi1(1)),dcpl_mb(spsi1(1))) 3770 call electron_get_Tgradient(dcpl_mb(spsi1(1)), 3771 > dbl_mb(tmp1(1)), 3772 > Rpsi) 3773 3774* **** THpsi = Rpsi - psi1 * <spsi1|Rpsi> **** 3775 call Grsm_gg_Copy(npack1,neq(1)+neq(2),Rpsi,THpsi) 3776 call Dneall_ffm_sym_Multiply(0,dcpl_mb(spsi1(1)),Rpsi,npack1, 3777 > dbl_mb(tmp1(1))) 3778 call Dneall_fmf_Multiply(0,dcpl_mb(psi1(1)),npack1, 3779 > dbl_mb(tmp1(1)), 3780 > -1.0d0,THpsi,1.0d0) 3781 3782 if (.not.Dneall_m_pop_stack(tmp1)) 3783 > call errquit('psi_1get_STgradient:error popping stack',1, 3784 > MA_ERR) 3785 3786 return 3787 end 3788 3789 3790* *********************************** 3791* * * 3792* * psi_2get_STgradient * 3793* * * 3794* *********************************** 3795 3796* THpsi = Hpsi - Y*Y^t*Hpsi ! used by Grassman minimizers 3797* 3798 subroutine psi_2get_STgradient(option,Rpsi,THpsi,Eout) 3799 implicit none 3800 integer option 3801 complex*16 Rpsi(*) 3802 complex*16 THpsi(*) 3803 real*8 Eout 3804 3805#include "bafdecls.fh" 3806#include "errquit.fh" 3807#include "psi.fh" 3808 3809* ***** rhoall common block **** 3810 integer rho1_all(2) 3811 integer rho2_all(2) 3812 common / rhoall_block / rho1_all,rho2_all 3813 3814* **** local variables **** 3815 integer tmp1(2),i,n 3816 logical value 3817 3818* **** external functions **** 3819 logical Dneall_m_push_get,Dneall_m_pop_stack 3820 external Dneall_m_push_get,Dneall_m_pop_stack 3821 3822 real*8 electron_energy 3823 external electron_energy 3824 3825 if (.not.Dneall_m_push_get(0,tmp1)) 3826 > call errquit('out of stack memory in psi_2get_STradient',0, 3827 > MA_ERR) 3828 3829 if (option.le.1) then 3830 call electron_run(dcpl_mb(psi2(1)), 3831 > dbl_mb(rho2(1)), 3832 > dcpl_mb(dng2(1)), 3833 > dbl_mb(rho2_all(1)), 3834 > occupation_on,dbl_mb(occ2(1))) 3835 end if 3836 Eout = electron_energy(dcpl_mb(psi2(1)), 3837 > dbl_mb(rho2(1)), 3838 > dcpl_mb(dng2(1)), 3839 > dbl_mb(rho2_all(1)), 3840 > occupation_on,dbl_mb(occ2(1))) 3841 call electron_gen_hml(dcpl_mb(psi2(1)), 3842 > dbl_mb(tmp1(1))) 3843 3844* **** get the tangent Rpsi = Hpsi - Spsi2*tmp1 , TY = HY - SY*Y'HY **** 3845 call psp_overlap_S(ispin,neq,dcpl_mb(psi2(1)),dcpl_mb(spsi1(1))) 3846 call electron_get_Tgradient(dcpl_mb(spsi1(1)), 3847 > dbl_mb(tmp1(1)), 3848 > Rpsi) 3849 3850* **** THpsi = Rpsi - psi2 * <spsi2|Rpsi> **** 3851 call Grsm_gg_Copy(npack1,neq(1)+neq(2),Rpsi,THpsi) 3852 call Dneall_ffm_sym_Multiply(0,dcpl_mb(spsi1(1)),Rpsi,npack1, 3853 > dbl_mb(tmp1(1))) 3854 call Dneall_fmf_Multiply(0,dcpl_mb(psi2(1)),npack1, 3855 > dbl_mb(tmp1(1)), 3856 > -1.0d0,THpsi,1.0d0) 3857 3858 if (.not.Dneall_m_pop_stack(tmp1)) 3859 > call errquit('psi_2get_STgradient:error popping stack',1, 3860 > MA_ERR) 3861 3862 return 3863 end 3864 3865 3866 3867 3868 3869* *********************************** 3870* * * 3871* * psi_1get_Gradient * 3872* * * 3873* *********************************** 3874 3875* THpsi = Hpsi ! used by Projected Grassman minimizers 3876* 3877 subroutine psi_1get_Gradient(THpsi,Eout) 3878 implicit none 3879 complex*16 THpsi(*) 3880 real*8 Eout 3881 3882#include "bafdecls.fh" 3883#include "psi.fh" 3884 3885* ***** rhoall common block **** 3886 integer rho1_all(2) 3887 integer rho2_all(2) 3888 common / rhoall_block / rho1_all,rho2_all 3889 3890* **** local variables **** 3891 3892* **** external functions **** 3893 real*8 electron_energy 3894 external electron_energy 3895 3896 3897 call electron_run(dcpl_mb(psi1(1)), 3898 > dbl_mb(rho1(1)), 3899 > dcpl_mb(dng1(1)), 3900 > dbl_mb(rho1_all(1)), 3901 > occupation_on,dbl_mb(occ1(1))) 3902 3903 Eout = electron_energy(dcpl_mb(psi1(1)), 3904 > dbl_mb(rho1(1)), 3905 > dcpl_mb(dng1(1)), 3906 > dbl_mb(rho1_all(1)), 3907 > occupation_on,dbl_mb(occ1(1))) 3908 3909 call electron_get_Gradient(THpsi) 3910 3911 return 3912 end 3913 3914 3915* *********************************** 3916* * * 3917* * psi_1gen_Tangent * 3918* * * 3919* *********************************** 3920 3921* THpsi = Hpsi - Y*Y^t*Hpsi ! used by Grassman minimizers 3922* 3923 subroutine psi_1gen_Tangent(THpsi) 3924 implicit none 3925 complex*16 THpsi(*) 3926 3927#include "bafdecls.fh" 3928#include "psi.fh" 3929#include "errquit.fh" 3930 3931* **** local variables **** 3932 integer tmp1(2) 3933 3934* **** external functions **** 3935 logical Dneall_m_push_get,Dneall_m_pop_stack 3936 external Dneall_m_push_get,Dneall_m_pop_stack 3937 3938 if (.not. Dneall_m_push_get(0,tmp1)) 3939 > call errquit('psi_1gen_Tangent:out of stack memory',0, MA_ERR) 3940 3941 call electron_gen_psiTangenthml(dcpl_mb(psi1(1)), 3942 > THpsi, 3943 > dbl_mb(tmp1(1))) 3944 call electron_gen_Tangent(dcpl_mb(psi1(1)), 3945 > dbl_mb(tmp1(1)), 3946 > THpsi) 3947 3948 if (.not. Dneall_m_pop_stack(tmp1)) 3949 > call errquit('error popping stack memory in psi_1get_Tradient',0, 3950 & MA_ERR) 3951 3952 return 3953 end 3954 3955 3956 3957 3958 3959* *********************************** 3960* * * 3961* * psi_2get_Tgradient * 3962* * * 3963* *********************************** 3964 subroutine psi_2get_Tgradient(option,THpsi,Eout) 3965 implicit none 3966 integer option 3967 complex*16 THpsi(*) 3968 real*8 Eout 3969 3970#include "errquit.fh" 3971#include "bafdecls.fh" 3972#include "psi.fh" 3973 3974* ***** rhoall common block **** 3975 integer rho1_all(2) 3976 integer rho2_all(2) 3977 common / rhoall_block / rho1_all,rho2_all 3978 3979* *** local variables **** 3980 integer tmp1(2) 3981 3982 3983* **** external functions **** 3984 logical Dneall_m_push_get,Dneall_m_pop_stack 3985 external Dneall_m_push_get,Dneall_m_pop_stack 3986 3987 real*8 electron_energy 3988 external electron_energy 3989 3990 3991!$OMP BARRIER 3992 3993 if (.not.Dneall_m_push_get(0,tmp1)) 3994 > call errquit('out of stack memory in psi_2get_Tradient',0, 3995 > MA_ERR) 3996 3997 if (option.le.1) then 3998 call electron_run(dcpl_mb(psi2(1)), 3999 > dbl_mb(rho2(1)), 4000 > dcpl_mb(dng2(1)), 4001 > dbl_mb(rho2_all(1)), 4002 > occupation_on,dbl_mb(occ2(1))) 4003 end if 4004 4005 Eout = electron_energy(dcpl_mb(psi2(1)), 4006 > dbl_mb(rho2(1)), 4007 > dcpl_mb(dng2(1)), 4008 > dbl_mb(rho2_all(1)), 4009 > occupation_on,dbl_mb(occ2(1))) 4010!$OMP BARRIER 4011 4012 call electron_gen_hml(dcpl_mb(psi2(1)), 4013 > dbl_mb(tmp1(1))) 4014 call electron_get_Tgradient(dcpl_mb(psi2(1)), 4015 > dbl_mb(tmp1(1)), 4016 > THpsi) 4017 4018 if (.not. Dneall_m_pop_stack(tmp1)) 4019 >call errquit('psi_2get_Tgradient:error popping stack',1,MA_ERR) 4020 4021 return 4022 end 4023 4024 4025* *********************************** 4026* * * 4027* * psi_2get_Gradient * 4028* * * 4029* *********************************** 4030 subroutine psi_2get_Gradient(option,THpsi,Eout) 4031 implicit none 4032 integer option 4033 complex*16 THpsi(*) 4034 real*8 Eout 4035 4036#include "bafdecls.fh" 4037#include "psi.fh" 4038 4039* ***** rhoall common block **** 4040 integer rho1_all(2) 4041 integer rho2_all(2) 4042 common / rhoall_block / rho1_all,rho2_all 4043 4044* *** local variables **** 4045 4046* **** external functions **** 4047 real*8 electron_energy 4048 external electron_energy 4049 4050 4051 if (option.le.1) then 4052 call electron_run(dcpl_mb(psi2(1)), 4053 > dbl_mb(rho2(1)), 4054 > dcpl_mb(dng2(1)), 4055 > dbl_mb(rho2_all(1)), 4056 > occupation_on,dbl_mb(occ2(1))) 4057 end if 4058 4059 Eout = electron_energy(dcpl_mb(psi2(1)), 4060 > dbl_mb(rho2(1)), 4061 > dcpl_mb(dng2(1)), 4062 > dbl_mb(rho2_all(1)), 4063 > occupation_on,dbl_mb(occ2(1))) 4064 4065 call electron_get_Gradient(THpsi) 4066 4067 return 4068 end 4069 4070* *********************************** 4071* * * 4072* * psi_2gen_Tangent * 4073* * * 4074* *********************************** 4075 subroutine psi_2gen_Tangent(THpsi) 4076 implicit none 4077 complex*16 THpsi(*) 4078 4079#include "bafdecls.fh" 4080#include "psi.fh" 4081#include "errquit.fh" 4082 4083* *** local variables **** 4084 integer tmp1(2) 4085 4086* **** external functions **** 4087 logical Dneall_m_push_get,Dneall_m_pop_stack 4088 external Dneall_m_push_get,Dneall_m_pop_stack 4089 4090 4091 if (.not. Dneall_m_push_get(0,tmp1)) 4092 > call errquit('psi_2gen_Tangent: out of stack memory',0,MA_ERR) 4093 4094 4095 call electron_gen_psiTangenthml(dcpl_mb(psi2(1)), 4096 > THpsi, 4097 > dbl_mb(tmp1(1))) 4098 call electron_gen_Tangent(dcpl_mb(psi2(1)), 4099 > dbl_mb(tmp1(1)), 4100 > THpsi) 4101 4102 if (.not. Dneall_m_pop_stack(tmp1)) 4103 > call errquit('error popping stack memory in psi_1get_Tradient',0, 4104 & MA_ERR) 4105 4106 return 4107 end 4108 4109 4110 4111 4112* *********************************** 4113* * * 4114* * psi_1get_TSgradient * 4115* * * 4116* *********************************** 4117 4118* THpsi = Hpsi - Y*Hpsi^t*Y ! used by Stiefel minimizers 4119* 4120 subroutine psi_1get_TSgradient(THpsi,Eout) 4121 implicit none 4122 complex*16 THpsi(*) 4123 real*8 Eout 4124 4125#include "errquit.fh" 4126#include "bafdecls.fh" 4127#include "psi.fh" 4128 4129* ***** rhoall common block **** 4130 integer rho1_all(2) 4131 integer rho2_all(2) 4132 common / rhoall_block / rho1_all,rho2_all 4133 4134* **** local variables **** 4135 integer tmp1(2) 4136 4137* **** external functions **** 4138 logical Dneall_m_push_get,Dneall_m_pop_stack 4139 external Dneall_m_push_get,Dneall_m_pop_stack 4140 4141 real*8 electron_energy 4142 external electron_energy 4143 4144 4145 if (.not. Dneall_m_push_get(0,tmp1)) 4146 > call errquit('psi_1get_TSradient:pushing stack',0, MA_ERR) 4147 4148 4149 call electron_run(dcpl_mb(psi1(1)), 4150 > dbl_mb(rho1(1)), 4151 > dcpl_mb(dng1(1)), 4152 > dbl_mb(rho1_all(1)), 4153 > occupation_on,dbl_mb(occ1(1))) 4154 4155 Eout = electron_energy(dcpl_mb(psi1(1)), 4156 > dbl_mb(rho1(1)), 4157 > dcpl_mb(dng1(1)), 4158 > dbl_mb(rho1_all(1)), 4159 > occupation_on,dbl_mb(occ1(1))) 4160 4161 4162 call electron_gen_hmlt(dcpl_mb(psi1(1)), 4163 > dbl_mb(tmp1(1))) 4164 call electron_get_Tgradient(dcpl_mb(psi1(1)), 4165 > dbl_mb(tmp1(1)), 4166 > THpsi) 4167 4168 4169 if (.not. Dneall_m_pop_stack(tmp1)) 4170 > call errquit('psi_1get_TSgradient:popping stack',1, MA_ERR) 4171 4172 return 4173 end 4174 4175 4176* *********************************** 4177* * * 4178* * psi_2get_TSgradient * 4179* * * 4180* *********************************** 4181 4182* THpsi = Hpsi - Y*Hpsi^t*Y ! used by Stiefel minimizers 4183* 4184 subroutine psi_2get_TSgradient(option,THpsi,Eout) 4185 implicit none 4186 integer option 4187 complex*16 THpsi(*) 4188 real*8 Eout 4189 4190#include "errquit.fh" 4191#include "bafdecls.fh" 4192#include "psi.fh" 4193 4194* ***** rhoall common block **** 4195 integer rho1_all(2) 4196 integer rho2_all(2) 4197 common / rhoall_block / rho1_all,rho2_all 4198 4199* *** local variables **** 4200 integer tmp1(2) 4201 4202* **** external functions **** 4203 logical Dneall_m_push_get,Dneall_m_pop_stack 4204 external Dneall_m_push_get,Dneall_m_pop_stack 4205 4206 real*8 electron_energy 4207 external electron_energy 4208 4209 4210 if (.not. Dneall_m_push_get(0,tmp1)) 4211 > call errquit('psi_2get_TSgradient:pushing stack',0, MA_ERR) 4212 4213 if (option.le.1) then 4214 call electron_run(dcpl_mb(psi2(1)), 4215 > dbl_mb(rho2(1)), 4216 > dcpl_mb(dng2(1)), 4217 > dbl_mb(rho2_all(1)), 4218 > occupation_on,dbl_mb(occ2(1))) 4219 end if 4220 4221 Eout = electron_energy(dcpl_mb(psi2(1)), 4222 > dbl_mb(rho2(1)), 4223 > dcpl_mb(dng2(1)), 4224 > dbl_mb(rho2_all(1)), 4225 > occupation_on,dbl_mb(occ2(1))) 4226 4227 call electron_gen_hmlt(dcpl_mb(psi2(1)), 4228 > dbl_mb(tmp1(1))) 4229 call electron_get_Tgradient(dcpl_mb(psi2(1)), 4230 > dbl_mb(tmp1(1)), 4231 > THpsi) 4232 4233 if (.not. Dneall_m_pop_stack(tmp1)) 4234 > call errquit('psi_2get_TSgradient:popping stack',1, MA_ERR) 4235 4236 return 4237 end 4238 4239 4240 4241 4242* *********************************** 4243* * * 4244* * psi_1get_TMgradient * 4245* * * 4246* *********************************** 4247 subroutine psi_1get_TMgradient(THpsi,Eout) 4248 implicit none 4249 complex*16 THpsi(*) 4250 real*8 Eout 4251 4252#include "bafdecls.fh" 4253#include "psi.fh" 4254 4255* ***** rhoall common block **** 4256 integer rho1_all(2) 4257 integer rho2_all(2) 4258 common / rhoall_block / rho1_all,rho2_all 4259 4260* **** external functions **** 4261 real*8 electron_energy 4262 external electron_energy 4263 4264 4265 call electron_run(dcpl_mb(psi1(1)), 4266 > dbl_mb(rho1(1)), 4267 > dcpl_mb(dng1(1)), 4268 > dbl_mb(rho1_all(1)), 4269 > occupation_on,dbl_mb(occ1(1))) 4270 4271 Eout = electron_energy(dcpl_mb(psi1(1)), 4272 > dbl_mb(rho1(1)), 4273 > dcpl_mb(dng1(1)), 4274 > dbl_mb(rho1_all(1)), 4275 > occupation_on,dbl_mb(occ1(1))) 4276 4277 call electron_get_TMgradient(dcpl_mb(psi1(1)), 4278 > THpsi) 4279 4280 return 4281 end 4282 4283 4284 4285* *********************************** 4286* * * 4287* * psi_2get_TMgradient * 4288* * * 4289* *********************************** 4290 subroutine psi_2get_TMgradient(option,THpsi,Eout) 4291 implicit none 4292 integer option 4293 complex*16 THpsi(*) 4294 real*8 Eout 4295 4296#include "bafdecls.fh" 4297#include "psi.fh" 4298 4299* ***** rhoall common block **** 4300 integer rho1_all(2) 4301 integer rho2_all(2) 4302 common / rhoall_block / rho1_all,rho2_all 4303 4304* **** external functions **** 4305 real*8 electron_energy 4306 external electron_energy 4307 4308 if (option.le.1) then 4309 call electron_run(dcpl_mb(psi2(1)), 4310 > dbl_mb(rho2(1)), 4311 > dcpl_mb(dng2(1)), 4312 > dbl_mb(rho2_all(1)), 4313 > occupation_on,dbl_mb(occ2(1))) 4314 end if 4315 4316 Eout = electron_energy(dcpl_mb(psi2(1)), 4317 > dbl_mb(rho2(1)), 4318 > dcpl_mb(dng2(1)), 4319 > dbl_mb(rho2_all(1)), 4320 > occupation_on,dbl_mb(occ2(1))) 4321 4322 call electron_get_TMgradient(dcpl_mb(psi2(1)), 4323 > THpsi) 4324 4325 return 4326 end 4327 4328* *********************************** 4329* * * 4330* * psi_1ke_Precondition * 4331* * * 4332* *********************************** 4333 subroutine psi_1ke_Precondition(Hpsi) 4334 implicit none 4335 complex*16 Hpsi(*) 4336 4337#include "bafdecls.fh" 4338#include "psi.fh" 4339 4340* **** local variables **** 4341 integer neall 4342 4343 neall = neq(1)+neq(2) 4344 call ke_Precondition(npack1,neall, 4345 > dcpl_mb(psi1(1)), 4346 > Hpsi) 4347 return 4348 end 4349 4350 4351 4352* *********************************** 4353* * * 4354* * psi_1geodesic_transport * 4355* * * 4356* *********************************** 4357 subroutine psi_1geodesic_transport(t,H0) 4358 implicit none 4359 real*8 t 4360 complex*16 H0(*) 4361 4362#include "bafdecls.fh" 4363#include "psi.fh" 4364 4365 4366 call geodesic_transport(t,dcpl_mb(psi1(1)),H0) 4367 4368 return 4369 end 4370 4371 4372* *********************************** 4373* * * 4374* * psi_1geodesic_Gtransport * 4375* * * 4376* *********************************** 4377 subroutine psi_1geodesic_Gtransport(t,G0) 4378 implicit none 4379 real*8 t 4380 complex*16 G0(*) 4381 4382#include "bafdecls.fh" 4383#include "psi.fh" 4384 4385 call geodesic_Gtransport(t,dcpl_mb(psi1(1)),G0) 4386 4387 return 4388 end 4389 4390 4391 4392* *********************************** 4393* * * 4394* * psi_geodesic_energy * 4395* * * 4396* *********************************** 4397 real*8 function psi_geodesic_energy(t) 4398 implicit none 4399 real*8 t 4400 4401#include "bafdecls.fh" 4402#include "psi.fh" 4403 4404* ***** rhoall common block **** 4405 integer rho1_all(2) 4406 integer rho2_all(2) 4407 common / rhoall_block / rho1_all,rho2_all 4408 4409 real*8 e_new 4410* **** external functions **** 4411 real*8 electron_energy 4412 external electron_energy 4413 4414 call geodesic_get(t,dcpl_mb(psi1(1)), 4415 > dcpl_mb(psi2(1))) 4416 call electron_run(dcpl_mb(psi2(1)), 4417 > dbl_mb(rho2(1)), 4418 > dcpl_mb(dng2(1)), 4419 > dbl_mb(rho2_all(1)), 4420 > occupation_on,dbl_mb(occ2(1))) 4421 e_new = electron_energy(dcpl_mb(psi2(1)), 4422 > dbl_mb(rho2(1)), 4423 > dcpl_mb(dng2(1)), 4424 > dbl_mb(rho2_all(1)), 4425 > occupation_on,dbl_mb(occ2(1))) 4426 4427 psi_geodesic_energy = e_new 4428 return 4429 end 4430 4431* *********************************** 4432* * * 4433* * psi_geodesic_denergy * 4434* * * 4435* *********************************** 4436 real*8 function psi_geodesic_denergy(t) 4437 implicit none 4438 real*8 t 4439 4440#include "bafdecls.fh" 4441#include "psi.fh" 4442 4443* **** external functions **** 4444 real*8 electron_eorbit_noocc 4445 external electron_eorbit_noocc 4446 4447 4448 call geodesic_transport(t,dcpl_mb(psi1(1)), 4449 > dcpl_mb(psi2(1))) 4450 psi_geodesic_denergy 4451 > = 2.0d0*electron_eorbit_noocc(dcpl_mb(psi2(1))) 4452 4453 return 4454 end 4455 4456* *********************************** 4457* * * 4458* * psi_geodesic_final * 4459* * * 4460* *********************************** 4461 subroutine psi_geodesic_final(t) 4462 implicit none 4463 real*8 t 4464 4465#include "bafdecls.fh" 4466#include "psi.fh" 4467 4468 call geodesic_get(t,dcpl_mb(psi1(1)), 4469 > dcpl_mb(psi2(1))) 4470 return 4471 end 4472 4473 4474 4475* *********************************** 4476* * * 4477* * psi_1geodesic2_start * 4478* * * 4479* *********************************** 4480 subroutine psi_1geodesic2_start(H0,max_sigma,dE0) 4481 implicit none 4482 complex*16 H0(*) 4483 real*8 max_sigma 4484 real*8 dE0 4485 4486#include "bafdecls.fh" 4487#include "psi.fh" 4488 4489 call geodesic2_start(dcpl_mb(psi1(1)),H0,max_sigma,dE0) 4490 4491 return 4492 end 4493 4494* *********************************** 4495* * * 4496* * psi_1geodesic2_transport * 4497* * * 4498* *********************************** 4499 subroutine psi_1geodesic2_transport(t,Hnew) 4500 implicit none 4501 real*8 t 4502 complex*16 Hnew(*) 4503 4504#include "bafdecls.fh" 4505#include "psi.fh" 4506 4507 call geodesic2_transport(t,dcpl_mb(psi1(1)),Hnew) 4508 4509 return 4510 end 4511 4512 4513 4514* *********************************** 4515* * * 4516* * psi_geodesic2_energy * 4517* * * 4518* *********************************** 4519 real*8 function psi_geodesic2_energy(t) 4520 implicit none 4521 real*8 t 4522 4523#include "bafdecls.fh" 4524#include "psi.fh" 4525 4526* ***** rhoall common block **** 4527 integer rho1_all(2) 4528 integer rho2_all(2) 4529 common / rhoall_block / rho1_all,rho2_all 4530 4531 real*8 e_new 4532 4533* **** external functions **** 4534 real*8 electron_energy 4535 external electron_energy 4536 4537 4538 call geodesic2_get(t,dcpl_mb(psi1(1)), 4539 > dcpl_mb(psi2(1))) 4540 4541 if (occupation_on) 4542 > call dcopy((ne(1)+ne(2)),dbl_mb(occ1(1)),1,dbl_mb(occ2(1)),1) 4543 4544* **** check Orthogonality of psi2 **** !debug 4545* call OrthoCheck_geo(ispin,ne,dcpl_mb(psi2(1))) !debug 4546 4547 4548 call electron_run(dcpl_mb(psi2(1)), 4549 > dbl_mb(rho2(1)), 4550 > dcpl_mb(dng2(1)), 4551 > dbl_mb(rho2_all(1)), 4552 > occupation_on,dbl_mb(occ2(1))) 4553 e_new = electron_energy(dcpl_mb(psi2(1)), 4554 > dbl_mb(rho2(1)), 4555 > dcpl_mb(dng2(1)), 4556 > dbl_mb(rho2_all(1)), 4557 > occupation_on,dbl_mb(occ2(1))) 4558 4559 psi_geodesic2_energy = e_new 4560 return 4561 end 4562 4563* *********************************** 4564* * * 4565* * psi_geodesic2_denergy * 4566* * * 4567* *********************************** 4568 real*8 function psi_geodesic2_denergy(t) 4569 implicit none 4570 real*8 t 4571 4572#include "bafdecls.fh" 4573#include "psi.fh" 4574 4575* **** external functions **** 4576 real*8 electron_eorbit 4577 external electron_eorbit 4578 4579 4580 call geodesic2_transport(t,dcpl_mb(psi1(1)), 4581 > dcpl_mb(psi2(1))) 4582 if (occupation_on) 4583 > call dcopy((ne(1)+ne(2)),dbl_mb(occ1(1)),1,dbl_mb(occ2(1)),1) 4584 4585 psi_geodesic2_denergy = 2.0d0*electron_eorbit(dcpl_mb(psi2(1)), 4586 > occupation_on,dbl_mb(occ2(1))) 4587 4588 return 4589 end 4590 4591* *********************************** 4592* * * 4593* * psi_geodesic2_final * 4594* * * 4595* *********************************** 4596 subroutine psi_geodesic2_final(t) 4597 implicit none 4598 real*8 t 4599 4600#include "bafdecls.fh" 4601#include "psi.fh" 4602 4603 integer taskid,MASTER 4604 parameter (MASTER=0) 4605c real*8 sum1,sum2 4606 4607 call Parallel_taskid(taskid) 4608 4609 call geodesic2_get(t,dcpl_mb(psi1(1)), 4610 > dcpl_mb(psi2(1))) 4611 if (occupation_on) 4612 > call dcopy((ne(1)+ne(2)),dbl_mb(occ1(1)),1,dbl_mb(occ2(1)),1) 4613 return 4614 end 4615 4616 4617 4618* *********************************** 4619* * * 4620* * psito2_sd_update * 4621* * * 4622* *********************************** 4623 subroutine psi1to2_sd_update(dte) 4624 implicit none 4625 real*8 dte 4626 4627#include "bafdecls.fh" 4628#include "psi.fh" 4629#include "errquit.fh" 4630 4631 4632* ***** rhoall common block **** 4633 integer rho1_all(2) 4634 integer rho2_all(2) 4635 common / rhoall_block / rho1_all,rho2_all 4636 4637* **** local variables **** 4638 logical value 4639 integer nemaxq,ierr 4640 integer lmd(2),tmp_L(2) 4641 4642* **** external functions **** 4643 logical pspw_SIC,Dneall_m_push_get,Dneall_m_push_get_block 4644 logical Dneall_m_pop_stack 4645 external pspw_SIC,Dneall_m_push_get,Dneall_m_push_get_block 4646 external Dneall_m_pop_stack 4647 4648 call electron_run(dcpl_mb(psi1(1)), 4649 > dbl_mb(rho1(1)), 4650 > dcpl_mb(dng1(1)), 4651 > dbl_mb(rho1_all(1)), 4652 > occupation_on,dbl_mb(occ1(1))) 4653 4654* **** do a steepest descent step **** 4655 call electron_sd_update(dcpl_mb(psi1(1)), 4656 > dcpl_mb(psi2(1)), 4657 > dte) 4658 4659* **** lagrange multiplier corrections **** 4660 nemaxq = neq(1)+neq(2) 4661 4662* **** allocate MA local variables **** 4663 value = Dneall_m_push_get_block(1,8,tmp_L) 4664 value = value.and.Dneall_m_push_get(0,lmd) 4665 4666c if (occupation_on) then 4667c call psi_lmbda2(ispin,ne,nemaxq,npack1, 4668c > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)), 4669c > dte,dbl_mb(occ1(1)), 4670c > dbl_mb(lmd(1)), 4671c > dbl_mb(tmp_L(1)),ierr) 4672 4673 if (pawexist) then 4674 call psp_overlap_S(ispin,neq, 4675 > dcpl_mb(psi1(1)), 4676 > dcpl_mb(spsi1(1))) 4677 call psi_lmbda_paw(ispin,neq,nemaxq,npack1, 4678 > dcpl_mb(spsi1(1)), 4679 > dcpl_mb(psi2(1)), 4680 > dte, 4681 > dbl_mb(lmd(1)), 4682 > dbl_mb(tmp_L(1)),ierr) 4683 else if (occupation_on) then 4684 call psi_lmbda2(ispin,ne,nemaxq,npack1, 4685 > dcpl_mb(psi1(1)), 4686 > dcpl_mb(psi2(1)), 4687 > dte,dbl_mb(occ1(1)), 4688 > dbl_mb(lmd(1)), 4689 > dbl_mb(tmp_L(1)),ierr) 4690 4691 else if (pspw_SIC()) then 4692 call psi_lmbda_sic(ispin,ne,nemaxq,npack1, 4693 > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dte, 4694 > dbl_mb(lmd(1)), 4695 > dbl_mb(tmp_L(1)),ierr) 4696 else 4697 call psi_lmbda(ispin,neq,nemaxq,npack1, 4698 > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dte, 4699 > dbl_mb(lmd(1)), 4700 > dbl_mb(tmp_L(1)),ierr) 4701 end if 4702 4703 4704 value = value.and.Dneall_m_pop_stack(lmd) 4705 value = value.and.Dneall_m_pop_stack(tmp_L) 4706 if (.not. value) 4707 > call errquit( 4708 > 'psi1to2_sd_update:stack failure', 0, MA_ERR) 4709 return 4710 end 4711 4712 4713* *********************************** 4714* * * 4715* * psi_1force * 4716* * * 4717* *********************************** 4718 subroutine psi_1force(fion) 4719 implicit none 4720 real*8 fion(3,*) 4721 4722#include "bafdecls.fh" 4723#include "errquit.fh" 4724#include "psi.fh" 4725 4726* **** local variables **** 4727 logical value 4728 integer r_grid(2),tmp(2) 4729 4730* **** external functions **** 4731 integer control_version 4732 external control_version 4733 logical psp_U_psputerm 4734 external psp_U_psputerm 4735 logical Dneall_m_push_get, Dneall_m_pop_stack 4736 external Dneall_m_push_get, Dneall_m_pop_stack 4737 4738c call electron_gen_psi_r(dcpl_mb(psi1(1))) 4739c call electron_gen_densities(dcpl_mb(psi1(1)), 4740c > dbl_mb(rho1(1)), 4741c > dcpl_mb(dng1(1))) 4742 4743 call f_vlocal(dcpl_mb(dng1(1)),fion) 4744 4745 if (control_version().eq.4) then 4746 value = BA_push_get(mt_dbl,(2*nfft3d),'tmp', 4747 > tmp(2),tmp(1)) 4748 value = value.and. 4749 > BA_push_get(mt_dbl,(6*nfft3d),'r_grid', 4750 > r_grid(2),r_grid(1)) 4751 if (.not. value) 4752 > call errquit('psi_1force:out of stack memory',0, MA_ERR) 4753 !call dcopy(2*nfft3d,0.0d0,0,dbl_mb(tmp(1)),1) 4754 call Parallel_shared_vector_zero(.true., 4755 > 2*nfft3d,dbl_mb(tmp(1))) 4756 4757 call D3dB_rr_Sum(1,dbl_mb(rho1(1)), 4758 > dbl_mb(rho1(1)+(ispin-1)*2*nfft3d), 4759 > dbl_mb(tmp(1))) 4760 call lattice_r_grid(dbl_mb(r_grid(1))) 4761 call grad_v_lr_local(dbl_mb(r_grid(1)), 4762 > dbl_mb(tmp(1)), 4763 > fion) 4764 4765 value = BA_pop_stack(r_grid(2)) 4766 value = value.and.BA_pop_stack(tmp(2)) 4767 if (.not.value) 4768 > call errquit('psi_1force:error popping stack memory',0, 4769 & MA_ERR) 4770 end if 4771 4772 call f_vnonlocal(ispin, 4773 > neq, 4774 > dcpl_mb(psi1(1)), 4775 > fion, 4776 > occupation_on,dbl_mb(occ1(1))) 4777 4778 if (pawexist) then 4779 value = Dneall_m_push_get(0,tmp) 4780 if (.not.value) 4781 > call errquit('psi_1force:out of stack memory',0,MA_ERR) 4782 4783 call psi_1toelectron() 4784 call electron_gen_hml(dcpl_mb(psi1(1)),dbl_mb(tmp(1))) 4785 4786 call psp_paw_overlap_fion(ispin, 4787 > dbl_mb(tmp(1)), 4788 > dcpl_mb(psi1(1)), 4789 > fion) 4790 4791 value = Dneall_m_pop_stack(tmp) 4792 if (.not.value) 4793 > call errquit('psi_1force:error popping stack',0,MA_ERR) 4794 end if 4795 4796 if (psp_U_psputerm()) then 4797 call f_psp_U_v_nonlocal(ispin, 4798 > neq, 4799 > dcpl_mb(psi1(1)), 4800 > fion, 4801 > occupation_on,dbl_mb(occ1(1)),.true.) 4802 end if 4803 return 4804 end 4805 4806* *********************************** 4807* * * 4808* * psi_1force_local * 4809* * * 4810* *********************************** 4811 subroutine psi_1force_local(fion) 4812 implicit none 4813 real*8 fion(3,*) 4814 4815#include "bafdecls.fh" 4816#include "errquit.fh" 4817#include "psi.fh" 4818 4819* **** local variables **** 4820 logical value 4821 integer r_grid(2),tmp(2) 4822 4823* **** external functions **** 4824 integer control_version 4825 external control_version 4826 4827c call electron_gen_psi_r(dcpl_mb(psi1(1))) 4828c call electron_gen_densities(dcpl_mb(psi1(1)), 4829c > dbl_mb(rho1(1)), 4830c > dcpl_mb(dng1(1))) 4831 4832 call f_vlocal(dcpl_mb(dng1(1)),fion) 4833 4834 if (control_version().eq.4) then 4835 value = BA_push_get(mt_dbl,(2*nfft3d),'tmp', 4836 > tmp(2),tmp(1)) 4837 value = value.and. 4838 > BA_push_get(mt_dbl,(6*nfft3d),'r_grid', 4839 > r_grid(2),r_grid(1)) 4840 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 4841 call dcopy(2*nfft3d,0.0d0,0,dbl_mb(tmp(1)),1) 4842 4843 call D3dB_rr_Sum(1,dbl_mb(rho1(1)), 4844 > dbl_mb(rho1(1)+(ispin-1)*2*nfft3d), 4845 > dbl_mb(tmp(1))) 4846 call lattice_r_grid(dbl_mb(r_grid(1))) 4847 call grad_v_lr_local(dbl_mb(r_grid(1)), 4848 > dbl_mb(tmp(1)), 4849 > fion) 4850 4851 value = BA_pop_stack(r_grid(2)) 4852 value = value.and.BA_pop_stack(tmp(2)) 4853 if (.not. value) call errquit('error popping stack memory',0, 4854 & MA_ERR) 4855 end if 4856 4857 return 4858 end 4859 4860* *********************************** 4861* * * 4862* * psi_1force_nonlocal * 4863* * * 4864* *********************************** 4865 subroutine psi_1force_nonlocal(fion) 4866 implicit none 4867 real*8 fion(3,*) 4868 4869#include "bafdecls.fh" 4870#include "errquit.fh" 4871#include "psi.fh" 4872 4873 call f_vnonlocal(ispin, 4874 > neq, 4875 > dcpl_mb(psi1(1)), 4876 > fion, 4877 > occupation_on,dbl_mb(occ1(1))) 4878 return 4879 end 4880 4881 4882* *********************************** 4883* * * 4884* * psi_1force_psp_U_v_nonlocal * 4885* * * 4886* *********************************** 4887 subroutine psi_1force_psp_U_v_nonlocal(fion) 4888 implicit none 4889 real*8 fion(3,*) 4890 4891#include "bafdecls.fh" 4892#include "errquit.fh" 4893#include "psi.fh" 4894 4895c **** external functions **** 4896 logical psp_U_psputerm 4897 external psp_U_psputerm 4898 4899 if (psp_U_psputerm()) then 4900 call f_psp_U_v_nonlocal(ispin, 4901 > neq, 4902 > dcpl_mb(psi1(1)), 4903 > fion, 4904 > occupation_on,dbl_mb(occ1(1)),.true.) 4905 end if 4906 4907 return 4908 end 4909 4910 4911 4912 4913 4914* *********************************** 4915* * * 4916* * psi_1ke_stress * 4917* * * 4918* *********************************** 4919 subroutine psi_1ke_stress(stress) 4920 implicit none 4921 real*8 stress(3,3) 4922 4923#include "bafdecls.fh" 4924#include "psi.fh" 4925 4926 call ke_euv(ispin,neq,dcpl_mb(psi1(1)),stress) 4927 return 4928 end 4929 4930* *********************************** 4931* * * 4932* * psi_1coulomb_stress * 4933* * * 4934* *********************************** 4935 subroutine psi_1coulomb_stress(stress) 4936 implicit none 4937 real*8 stress(3,3) 4938 4939#include "bafdecls.fh" 4940#include "psi.fh" 4941 4942 call coulomb_euv(dcpl_mb(dng1(1)),stress) 4943 return 4944 end 4945 4946* *********************************** 4947* * * 4948* * rho_1exc_stress * 4949* * * 4950* *********************************** 4951 subroutine rho_1exc_stress(stress) 4952 implicit none 4953 real*8 stress(3,3) 4954 4955#include "bafdecls.fh" 4956#include "psi.fh" 4957 4958* ***** rhoall common block **** 4959 integer rho1_all(2) 4960 integer rho2_all(2) 4961 common / rhoall_block / rho1_all,rho2_all 4962 4963* ***** local variables **** 4964 integer u,v,gga 4965 real*8 exc,pxc 4966 real*8 pi,scal,hm(3,3),tstress(3,3) 4967 4968* **** external functions **** 4969 integer control_gga 4970 real*8 rho_1exc,rho_1pxc,lattice_unitg,lattice_omega 4971 external control_gga 4972 external rho_1exc,rho_1pxc,lattice_unitg,lattice_omega 4973 4974* *** define hm **** 4975 pi = 4.0d0*datan(1.0d0) 4976 scal = 1.0d0/(2.0d0*pi) 4977 do v=1,3 4978 do u=1,3 4979 hm(u,v) = scal*lattice_unitg(u,v) 4980 end do 4981 end do 4982 4983* **** LDA part **** 4984 exc = rho_1exc() 4985 pxc = rho_1pxc() 4986 do v=1,3 4987 do u=1,3 4988 stress(u,v) = (exc-pxc)*hm(u,v) 4989 end do 4990 end do 4991 !write(*,*) "hm(1,1):",hm(1,1),1.0d0/hm(1,1) 4992 !write(*,*) "exc:",exc,pxc 4993 !write(*,*) "D:",stress(1,1) 4994 4995* **** PBE96 GGA part **** 4996* **** finished? 11/24/04 - still need to test *** 4997 gga = control_gga() 4998 if ((gga.ge.10).and.(gga.lt.100)) then 4999 call v_bwexc_euv(gga,2*nfft3d,ispin,dbl_mb(rho1_all(1)), 5000 > 1.0d0,1.0d0,tstress) 5001 do v=1,3 5002 do u=1,3 5003 stress(u,v) = stress(u,v) + tstress(u,v) 5004 end do 5005 end do 5006 end if 5007 5008 if (gga.eq.110) then 5009 call v_bwexc_euv(10,2*nfft3d,ispin,dbl_mb(rho1_all(1)), 5010 > 0.75d0,1.0d0,tstress) 5011 do v=1,3 5012 do u=1,3 5013 stress(u,v) = stress(u,v) + tstress(u,v) 5014 end do 5015 end do 5016 end if 5017 5018 if (gga.eq.112) then 5019 call v_bwexc_euv(12,2*nfft3d,ispin,dbl_mb(rho1_all(1)), 5020 > 0.75d0,1.0d0,tstress) 5021 do v=1,3 5022 do u=1,3 5023 stress(u,v) = stress(u,v) + tstress(u,v) 5024 end do 5025 end do 5026 end if 5027 5028 return 5029 end 5030 5031* *********************************** 5032* * * 5033* * rho_1semicore_stress * 5034* * * 5035* *********************************** 5036 subroutine rho_1semicore_stress(stress) 5037 implicit none 5038 real*8 stress(3,3) 5039 5040#include "bafdecls.fh" 5041#include "psi.fh" 5042 5043* **** not finished **** 5044 call semicore_euv(stress) 5045 5046 return 5047 end 5048 5049 5050 5051 5052* *********************************** 5053* * * 5054* * dng_1vlocal_stress * 5055* * * 5056* *********************************** 5057 5058 subroutine dng_1vlocal_stress(stress) 5059 implicit none 5060 real*8 stress(3,3) 5061 5062 5063#include "bafdecls.fh" 5064#include "psi.fh" 5065 5066 call v_local_euv(dcpl_mb(dng1(1)),stress) 5067 5068 return 5069 end 5070 5071* *********************************** 5072* * * 5073* * psi_1vnonlocal_stress * 5074* * * 5075* *********************************** 5076 subroutine psi_1vnonlocal_stress(stress) 5077 implicit none 5078 real*8 stress(3,3) 5079 5080#include "bafdecls.fh" 5081#include "psi.fh" 5082 5083 5084 5085* ***** local variables **** 5086 integer u,v 5087 real*8 evnl 5088 real*8 pi,scal,hm(3,3) 5089 5090* **** external functions **** 5091 real*8 psi_1vnl,lattice_unitg 5092 external psi_1vnl,lattice_unitg 5093 5094* *** define hm **** 5095 pi = 4.0d0*datan(1.0d0) 5096 scal = 1.0d0/(2.0d0*pi) 5097 do v=1,3 5098 do u=1,3 5099 hm(u,v) = scal*lattice_unitg(u,v) 5100 end do 5101 end do 5102 5103 call v_nonlocal_euv_2(ispin,neq,dcpl_mb(psi1(1)),stress) 5104 evnl = psi_1vnl() 5105 5106 do v=1,3 5107 do u=1,3 5108 stress(u,v) = stress(u,v) - evnl*hm(u,v) 5109 end do 5110 end do 5111 5112 return 5113 end 5114 5115 5116 5117 5118* *********************************** 5119* * * 5120* * psi_1Orb_Analysis * 5121* * * 5122* *********************************** 5123 subroutine psi_1Orb_Analysis(iunit) 5124 implicit none 5125 integer iunit 5126 5127#include "bafdecls.fh" 5128#include "psi.fh" 5129 5130c call Orb_Analysis(iunit,ispin,ne,dcpl_mb(psi1(1))) 5131 return 5132 end 5133 5134* *********************************** 5135* * * 5136* * psi_1Shml * 5137* * * 5138* *********************************** 5139 subroutine psi_1Shml(S0,S0hml) 5140 implicit none 5141 complex*16 S0(*) 5142 complex*16 S0hml(*) 5143 5144#include "bafdecls.fh" 5145#include "psi.fh" 5146 5147 integer ms,n,shift1,shift2 5148 5149 call electron_gen_hml(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 5150 do ms=1,ispin 5151 n = ne(ms) 5152 if (n.le.0) go to 30 5153 shift1 = 1 + (ms-1)*ne(1)*npack1 5154 shift2 = (ms-1)*ne(1)*ne(1) 5155 call DGEMM('N','N',2*npack1,n,n, 5156 > (1.0d0), 5157 > S0(shift1), 2*npack1, 5158 > dbl_mb(hml(1)+shift2), n, 5159 > (0.0d0), 5160 > S0hml(shift1), 2*npack1) 5161 30 continue 5162 end do 5163 return 5164 end 5165 5166 5167 5168* *********************************** 5169* * * 5170* * psi_1gen_hml * 5171* * * 5172* *********************************** 5173 subroutine psi_1gen_hml() 5174 implicit none 5175 5176#include "bafdecls.fh" 5177#include "psi.fh" 5178 5179 5180 call electron_gen_hml(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 5181 5182 return 5183 end 5184 5185 5186 5187 5188* *********************************** 5189* * * 5190* * psi_1gen_hml_g * 5191* * * 5192* *********************************** 5193 subroutine psi_1gen_hml_g() 5194 implicit none 5195 5196#include "bafdecls.fh" 5197#include "psi.fh" 5198 5199 5200 call electron_gen_hml_g(dcpl_mb(psi1(1)),dbl_mb(hml(1))) 5201 5202 return 5203 end 5204 5205 5206* *********************************** 5207* * * 5208* * psi_2gen_hml * 5209* * * 5210* *********************************** 5211 subroutine psi_2gen_hml() 5212 implicit none 5213 5214#include "bafdecls.fh" 5215#include "psi.fh" 5216 5217 5218 call electron_gen_hml(dcpl_mb(psi2(1)),dbl_mb(hml(1))) 5219 5220 return 5221 end 5222 5223* *********************************** 5224* * * 5225* * psi_hml_value * 5226* * * 5227* *********************************** 5228 real*8 function psi_hml_value(ms,i,j) 5229 implicit none 5230 integer ms 5231 integer i,j 5232 5233#include "bafdecls.fh" 5234#include "psi.fh" 5235 5236 psi_hml_value = dbl_mb(hml(1) + (i-1)+(j-1)*ne(ms) 5237 > + (ms-1)*ne(1)*ne(1)) 5238 return 5239 end 5240 5241* *********************************** 5242* * * 5243* * psi_eigenvalue * 5244* * * 5245* *********************************** 5246 real*8 function psi_eigenvalue(ms,i) 5247 implicit none 5248 integer ms 5249 integer i 5250 5251#include "bafdecls.fh" 5252#include "psi.fh" 5253 5254 real*8 sum 5255 5256 sum = dbl_mb(eig(1)+(i-1)+(ms-1)*ne(1)) 5257 psi_eigenvalue = sum 5258 5259 return 5260 end 5261 5262* *********************************** 5263* * * 5264* * psi_occupation * 5265* * * 5266* *********************************** 5267 real*8 function psi_occupation(ms,i) 5268 implicit none 5269 integer ms 5270 integer i 5271 5272#include "bafdecls.fh" 5273#include "psi.fh" 5274 5275 if (occupation_on) then 5276 psi_occupation = dbl_mb(occ1(1)+(i-1)+(ms-1)*ne(1)) 5277 else 5278 psi_occupation = 1.0d0 5279 end if 5280 return 5281 end 5282 5283* *********************************** 5284* * * 5285* * psi_1reverse_occupation * 5286* * * 5287* *********************************** 5288 subroutine psi_1reverse_occupation() 5289 implicit none 5290 5291#include "bafdecls.fh" 5292#include "psi.fh" 5293 5294 integer ms,i,indx1,indx2 5295 5296 do ms=1,ispin 5297 indx1 = occ1(1) + ne(ms) - 1 + (ms-1)*ne(1) 5298 indx2 = occ2(1) + (ms-1)*ne(1) 5299 do i=1,ne(ms) 5300 dbl_mb(indx2)=dbl_mb(indx1) 5301 indx1 = indx1 - 1 5302 indx2 = indx2 + 1 5303 end do 5304 end do 5305 call dcopy((ne(1)+ne(2)),dbl_mb(occ2(1)),1,dbl_mb(occ1(1)),1) 5306 return 5307 end 5308 5309* *********************************** 5310* * * 5311* * psi_1assending_occupation * 5312* * * 5313* *********************************** 5314* 5315* makes sure the occupation follows an assending order 5316* 5317 subroutine psi_1assending_occupation() 5318 implicit none 5319 5320#include "bafdecls.fh" 5321#include "psi.fh" 5322 5323 if (dbl_mb(occ1(1)).lt.dbl_mb(occ1(1)+ne(1)-1)) then 5324 call psi_1reverse_occupation() 5325 end if 5326 return 5327 end 5328 5329* *********************************** 5330* * * 5331* * psi_1desending_occupation * 5332* * * 5333* *********************************** 5334* 5335* makes sure the occupation follows a desending order 5336* 5337 subroutine psi_1desending_occupation() 5338 implicit none 5339 5340#include "bafdecls.fh" 5341#include "psi.fh" 5342 5343 if (dbl_mb(occ1(1)+ne(1)-1).lt.dbl_mb(occ1(1))) then 5344 call psi_1reverse_occupation() 5345 end if 5346 return 5347 end 5348 5349 5350 5351 5352 5353* *********************************** 5354* * * 5355* * psi_1define_occupation * 5356* * * 5357* *********************************** 5358 subroutine psi_1define_occupation(initial_alpha,use_hml) 5359 implicit none 5360 real*8 initial_alpha 5361 logical use_hml 5362 5363#include "bafdecls.fh" 5364#include "psi.fh" 5365 5366* **** local variables **** 5367 integer it,itmax 5368 parameter (itmax=50) 5369 5370 integer ms,nb,n,shift1,shift2,occ1_tag,ndiff 5371 real*8 e,x,kT,f,g,alpha,pi,f0 5372 real*8 ZZ,Z(2),Zlower,Zmid,Zupper,elower,emid,eupper 5373 real*8 flower,fmid,fupper,lmbda 5374 5375* **** external functions **** 5376 integer control_multiplicity 5377 real*8 control_TotalCharge,ion_TotalCharge_qm 5378 real*8 psi_occ_distribution,control_fractional_alpha 5379 external control_multiplicity 5380 external control_TotalCharge,ion_TotalCharge_qm 5381 external psi_occ_distribution,control_fractional_alpha 5382 5383 if (use_hml) then 5384 do ms=1,ispin 5385 shift1 = eig(1) + (ms-1)*ne(1) 5386 shift2 = hml(1) + (ms-1)*ne(1)*ne(1) 5387 do n=1,ne(ms) 5388 dbl_mb(shift1+n-1) = dbl_mb(shift2+(n-1)+(n-1)*ne(ms)) 5389 end do 5390 end do 5391 end if 5392 5393 smearfermi(1) = 0.0d0 5394 smearfermi(2) = 0.0d0 5395 smearcorrection = 0.0d0 5396 5397 if (occupation_on) then 5398 if (initial_alpha.lt.0.0d0) then 5399 alpha = control_fractional_alpha() 5400 else 5401 alpha = initial_alpha 5402 end if 5403 kT = smearkT 5404 ZZ = ion_TotalCharge_qm() - control_TotalCharge() 5405 if (dabs(ZZ).lt.1.0d-9) go to 98 5406 5407 if (ispin.eq.2) then 5408 ndiff = control_multiplicity() - 1 5409 Z(1) = 0.5d0*(ZZ+ndiff) 5410 Z(2) = 0.5d0*(ZZ-ndiff) 5411 else 5412 Z(1) = 0.5d0*ZZ 5413 Z(2) = 0.0d0 5414 end if 5415 5416 pi = 4.0d0*datan(1.0d0) 5417 !if (initial) alpha = 1.0d0 5418 if (smeartype.le.0) alpha = 0.0d0 5419 5420* **** outer loop over spins **** 5421 smearcorrection = 0.0d0 5422 do ms=1,ispin 5423 5424* **** find eupper and elower **** 5425 elower = 9.9d12 5426 eupper = -9.9d12 5427 shift1 = eig(1) + (ms-1)*ne(1) 5428 do n=1,ne(ms) 5429 e = dbl_mb(shift1) 5430 if (e.lt.elower) elower = e 5431 if (e.gt.eupper) eupper = e 5432 shift1 = shift1 + 1 5433 end do 5434 5435 5436* **** find fermi level **** 5437 Zlower = 0.0d0 5438 Zupper = 0.0d0 5439 shift1 = eig(1) + (ms-1)*ne(1) 5440 do n=1,ne(ms) 5441 e = dbl_mb(shift1) 5442 Zlower = Zlower 5443 > + psi_occ_distribution(smeartype,(e-elower)/kT) 5444 Zupper = Zupper 5445 > + psi_occ_distribution(smeartype,(e-eupper)/kT) 5446 shift1 = shift1 + 1 5447 end do 5448 5449 flower = Zlower - Z(ms) 5450 fupper = Zupper - Z(ms) 5451 5452 if (flower*fupper.ge.0.0d0) 5453 > call errquit( 5454 > 'psi_1define_occupation:Fermi energy not found',ms,0) 5455 5456 it = 0 5457 20 it = it + 1 5458 emid = 0.5d0*(elower + eupper) 5459 Zmid = 0.0d0 5460 shift1 = eig(1) + (ms-1)*ne(1) 5461 do n=1,ne(ms) 5462 e = dbl_mb(shift1) 5463 Zmid = Zmid + psi_occ_distribution(smeartype,(e-emid)/kT) 5464 shift1 = shift1 + 1 5465 end do 5466 fmid = Zmid - Z(ms) 5467 if (fmid.lt.0.0d0) then 5468 flower = fmid 5469 elower = emid 5470 else 5471 fupper = fmid 5472 eupper = emid 5473 end if 5474 if ( (dabs(fmid) .gt.1.0d-11) .and. 5475 > ((eupper-elower).gt.1.0d-11) .and. 5476 > (it.lt.itmax)) goto 20 5477 5478 5479 smearfermi(ms) = emid 5480 5481* **** determine filling and correction **** 5482 shift1 = eig(1) + (ms-1)*ne(1) 5483 shift2 = occ1(1) + (ms-1)*ne(1) 5484 do n=1,ne(ms) 5485 e = dbl_mb(shift1) 5486 x = (e - smearfermi(ms))/kT 5487 f = psi_occ_distribution(smeartype,x) 5488 f0 = dbl_mb(shift2) 5489 5490 dbl_mb(shift2) = (1.0d0-alpha)*f0 + alpha*f 5491 5492 if (smeartype.eq.1) then 5493 if ( (dbl_mb(shift2) .gt.1.0d-6) .and. 5494 > ( (1.0d0-dbl_mb(shift2)).gt.1.0d-6)) then 5495 smearcorrection = smearcorrection 5496 > + kT*( dbl_mb(shift2)*log(dbl_mb(shift2)) 5497 > + (1.0d0-dbl_mb(shift2))*log(1.0d0-dbl_mb(shift2)) ) 5498 end if 5499 else if (smeartype.eq.2) then 5500 smearcorrection 5501 > = smearcorrection 5502 > - kT*dexp(-x*x)/(4.0d0*dsqrt(pi)) 5503 end if 5504 5505 shift1 = shift1 + 1 5506 shift2 = shift2 + 1 5507 end do 5508 5509 end do !** ms*** 5510 if (ms.eq.1) smearcorrection=smearcorrection+smearcorrection 5511 5512 go to 99 5513 5514 98 continue 5515 smearcorrection = 0.0d0 5516 do ms=1,ispin 5517 shift2 = occ1(1) + (ms-1)*ne(1) 5518 call dcopy(ne(ms),0.0d0,0,dbl_mb(shift2),1) 5519 end do 5520 5521 99 continue 5522 5523 end if 5524 5525 return 5526 end 5527 5528 5529c set nwpw:fractional_smeartype 1 #0-none, 1-Fermi-Dirac, 2-Gaussian, 3-Hermite 5530 5531 real*8 function psi_occ_distribution(smeartype,e) 5532 implicit none 5533 integer smeartype 5534 real*8 e 5535 real*8 f 5536 5537* **** external functions **** 5538 real*8 util_erfc 5539 external util_erfc 5540 5541 if (smeartype.eq.1) then 5542 if (e.gt.30.0d0) then 5543 f = 0.0d0 5544 else if (e.lt.(-30.0d0)) then 5545 f = 1.0d0 5546 else 5547 f = 1.0d0/(1.0d0+dexp(e)) 5548 end if 5549 else if (smeartype.eq.2) then 5550 f = 0.5d0*util_erfc(e) 5551 else 5552 if (e.gt.0.0d0) then 5553 f = 0.0d0 5554 else 5555 f = 1.0d0 5556 end if 5557 end if 5558 psi_occ_distribution = f 5559 return 5560 end 5561 5562 real*8 function psi_smearfermi(ms) 5563 implicit none 5564 integer ms 5565#include "psi.fh" 5566 psi_smearfermi = smearfermi(ms) 5567 return 5568 end 5569 real*8 function psi_smearcorrection() 5570 implicit none 5571#include "psi.fh" 5572 psi_smearcorrection = smearcorrection 5573 return 5574 end 5575 5576 5577 5578 5579 5580* *********************************** 5581* * * 5582* * psi_virtual * 5583* * * 5584* *********************************** 5585 real*8 function psi_virtual(ms,i) 5586 implicit none 5587 integer ms 5588 integer i 5589 5590#include "bafdecls.fh" 5591#include "psi.fh" 5592 5593 psi_virtual=dbl_mb(eig_excited(1)+(i-1)+(ms-1)*ne_excited(1)) 5594 5595 return 5596 end 5597 5598* *********************************** 5599* * * 5600* * psi_hml * 5601* * * 5602* *********************************** 5603 real*8 function psi_hml(ms,i,j) 5604 implicit none 5605 integer ms 5606 integer i,j 5607 5608#include "bafdecls.fh" 5609#include "psi.fh" 5610 5611 psi_hml = dbl_mb(hml(1)-1 + i 5612 > + (j-1)*ne(ms) 5613 > + (ms-1)*ne(1)*ne(1)) 5614 5615 return 5616 end 5617 5618 5619* *********************************** 5620* * * 5621* * psi_iptr_hml * 5622* * * 5623* *********************************** 5624 integer function psi_iptr_hml(ms,i,j) 5625 implicit none 5626 integer ms 5627 integer i,j 5628 5629#include "bafdecls.fh" 5630#include "psi.fh" 5631 5632 psi_iptr_hml = (hml(1)-1 + i 5633 > + (j-1)*ne(ms) 5634 > + (ms-1)*ne(1)*ne(1)) 5635 5636 return 5637 end 5638 5639 5640* *********************************** 5641* * * 5642* * psi_spin_density * 5643* * * 5644* *********************************** 5645 subroutine psi_spin_density(en) 5646 implicit none 5647 real*8 en(2) 5648 5649#include "bafdecls.fh" 5650#include "psi.fh" 5651 5652* **** local variables **** 5653 integer ms,nx,ny,nz,n2ft3d 5654 real*8 scale,sumall 5655 5656* **** external functions **** 5657 real*8 lattice_omega 5658 external lattice_omega 5659 5660 call D3dB_n2ft3d(1,n2ft3d) 5661 !n2ft3d = 2*n2ft3d 5662 call D3dB_nx(1,nx) 5663 call D3dB_ny(1,ny) 5664 call D3dB_nz(1,nz) 5665 scale = lattice_omega()/dble(nx*ny*nz) 5666 5667* **** check total number of electrons **** 5668 en(2) = 0.0d0 5669 do ms =1,ispin 5670 call D3dB_r_dsum(1,dbl_mb(rho1(1)+(ms-1)*n2ft3d),sumall) 5671 en(ms) = sumall*scale 5672 end do 5673 5674 return 5675 end 5676 5677* *********************************** 5678* * * 5679* * psi_spin2 * 5680* * * 5681* *********************************** 5682 subroutine psi_spin2(Sab) 5683 implicit none 5684 real*8 Sab 5685 5686#include "bafdecls.fh" 5687#include "psi.fh" 5688 5689 call Calculate_psi_spin2(ispin,ne,npack1,dcpl_mb(psi1(1)), 5690 > occupation_on,dbl_mb(occ1(1)),Sab) 5691 return 5692 end 5693 5694* *********************************** 5695* * * 5696* * psi_1rotate2 * 5697* * * 5698* *********************************** 5699 subroutine psi_1rotate2() 5700 implicit none 5701 5702#include "bafdecls.fh" 5703#include "psi.fh" 5704 5705c* ***** local variables ***** 5706c integer ms,index,i,j,shift1,shift2 5707 5708 call Dneall_fmf_Multiply(0,dcpl_mb(psi1(1)),npack1, 5709 > dbl_mb(hml(1)),1.0d0, 5710 > dcpl_mb(psi2(1)),0.0d0) 5711 5712 5713 5714c !call dcopy(2*npack1*(ne(1)+ne(2)),0.0d0,0,dcpl_mb(psi2(1)),1) 5715c do ms=1,ispin 5716c if (ne(ms).le.0) go to 30 5717c shift1 = (ms-1)*ne(1) 5718c shift2 = (ms-1)*ne(1)*ne(1) 5719c 5720c call DGEMM('N','N',2*npack1,ne(ms),ne(ms), 5721c > (1.0d0), 5722c > dcpl_mb(psi1(1)+shift1*npack1),2*npack1, 5723c > dbl_mb(hml(1)+shift2),ne(ms), 5724c > (0.0d0), 5725c > dcpl_mb(psi2(1)+shift1*npack1),2*npack1) 5726cc do j=1,ne(ms) 5727cc do i=1,ne(ms) 5728cc index = (i-1) + (j-1)*ne(ms) + shift2 5729cc 5730cc call D3dB_cc_daxpy(1,dbl_mb(hml(1)+index), 5731cc > dcpl_mb(psi1(1)+(i-1+shift1)*nfft3d), 5732cc > dcpl_mb(psi2(1)+(j-1+shift1)*nfft3d)) 5733cc call Pack_cc_daxpy(1,dbl_mb(hml(1)+index), 5734cc > dcpl_mb(psi1(1)+(i-1+shift1)*npack1), 5735cc > dcpl_mb(psi2(1)+(j-1+shift1)*npack1)) 5736cc end do 5737cc end do 5738c 5739c 30 continue 5740c end do 5741 5742 return 5743 end 5744 5745* *********************************** 5746* * * 5747* * psi_2rotate1 * 5748* * * 5749* *********************************** 5750 subroutine psi_2rotate1() 5751 implicit none 5752 5753#include "bafdecls.fh" 5754#include "psi.fh" 5755 5756c* ***** local variables ***** 5757c integer ms,index,i,j,shift1,shift2 5758 5759 call Dneall_fmf_Multiply(0,dcpl_mb(psi2(1)),npack1, 5760 > dbl_mb(hml(1)),1.0d0, 5761 > dcpl_mb(psi1(1)),0.0d0) 5762 5763c do ms=1,ispin 5764c if (ne(ms).le.0) go to 30 5765c shift1 = (ms-1)*ne(1) 5766c shift2 = (ms-1)*ne(1)*ne(1) 5767c 5768c call DGEMM('N','N',2*npack1,ne(ms),ne(ms), 5769c > (1.0d0), 5770c > dcpl_mb(psi2(1)+shift1*npack1),2*npack1, 5771c > dbl_mb(hml(1)+shift2),ne(ms), 5772c > (0.0d0), 5773c > dcpl_mb(psi1(1)+shift1*npack1),2*npack1) 5774c 5775c 30 continue 5776c end do 5777 5778 return 5779 end 5780 5781 5782* *********************************** 5783* * * 5784* * psi_diagonalize_hml * 5785* * * 5786* *********************************** 5787 subroutine psi_diagonalize_hml() 5788 implicit none 5789 5790#include "bafdecls.fh" 5791#include "psi.fh" 5792 5793 5794c* ***** local variables **** 5795c logical value 5796c integer ms,shift1,shift2,ierr,i,j,indx 5797c integer tmp1(2) 5798 5799 5800 call Dneall_m_diagonalize(0,dbl_mb(hml(1)), 5801 > dbl_mb(eig(1)),.false.) 5802 5803c value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 5804c if (.not. value) 5805c > call errquit('out of stack memory in psi_diagonalize_hml',0, 5806c & MA_ERR) 5807 5808 5809c* ***** diagonalize the hamiltonian matrix ***** 5810c call dcopy((ne(1)+ne(2)),0.0d0,0,dbl_mb(eig(1)),1) 5811c do ms=1,ispin 5812c shift1 = (ms-1)*ne(1) 5813c shift2 = (ms-1)*ne(1)*ne(1) 5814c if (ne(ms).le.0) go to 30 5815 5816cc call eigen(ne(ms),ne(ms), 5817cc > dbl_mb(hml(1)+shift2), 5818cc > dbl_mb(eig(1)+shift1), 5819cc > dbl_mb(tmp1(1)),ierr) 5820c 5821c call DSYEV('V','U',ne(ms), 5822c > dbl_mb(hml(1)+shift2),ne(ms), 5823c > dbl_mb(eig(1)+shift1), 5824c > dbl_mb(tmp1(1)),2*ne(1)*ne(1), 5825c > ierr) 5826c 5827c call eigsrt(dbl_mb(eig(1)+shift1), 5828c > dbl_mb(hml(1)+shift2), 5829c > ne(ms),ne(ms)) 5830c 5831c 30 continue 5832c end do 5833c 5834c 5835c value = BA_pop_stack(tmp1(2)) 5836c if (.not. value) 5837c > call errquit('error popping stack in psi_diagonalize_hml',0, 5838c & MA_ERR) 5839 5840 return 5841 end 5842 5843* *********************************** 5844* * * 5845* * psi_diagonalize_hml_assending * 5846* * * 5847* *********************************** 5848 subroutine psi_diagonalize_hml_assending() 5849 implicit none 5850 5851#include "bafdecls.fh" 5852#include "psi.fh" 5853 5854 5855c* ***** local variables **** 5856c logical value 5857c integer ms,shift1,shift2,ierr 5858c integer tmp1(2) 5859 5860 call Dneall_m_diagonalize(0,dbl_mb(hml(1)),dbl_mb(eig(1)),.true.) 5861 5862c value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'tmp1',tmp1(2),tmp1(1)) 5863c if (.not. value) 5864c > call errquit( 5865c > 'out of stack memory in psi_diagonalize_hml_assending',0, 5866c > MA_ERR) 5867c 5868c 5869c* ***** diagonalize the hamiltonian matrix ***** 5870c call dcopy((ne(1)+ne(2)),0.0d0,0,dbl_mb(eig(1)),1) 5871c do ms=1,ispin 5872c shift1 = (ms-1)*ne(1) 5873c shift2 = (ms-1)*ne(1)*ne(1) 5874c if (ne(ms).le.0) go to 30 5875c 5876c call DSYEV('V','U',ne(ms), 5877c > dbl_mb(hml(1)+shift2),ne(ms), 5878c > dbl_mb(eig(1)+shift1), 5879c > dbl_mb(tmp1(1)),2*ne(1)*ne(1), 5880c > ierr) 5881c 5882c 30 continue 5883c end do 5884 5885 5886c value = BA_pop_stack(tmp1(2)) 5887c if (.not. value) 5888c > call errquit( 5889c > 'error popping stack in psi_diagonalize_hml_assending',0, 5890c > MA_ERR) 5891 5892 return 5893 end 5894 5895 5896 5897* *************************** 5898* * * 5899* * psi_error * 5900* * * 5901* *************************** 5902 real*8 function psi_error() 5903 implicit none 5904#include "errquit.fh" 5905 5906#include "bafdecls.fh" 5907#include "psi.fh" 5908 5909* ***** local variables **** 5910 logical value 5911 integer k,n 5912 real*8 error,sum,size 5913 integer tmp1(2) 5914 5915 value = BA_push_get(mt_dcpl,(npack1),'tmp1',tmp1(2),tmp1(1)) 5916 if (.not. value) 5917 > call errquit('out of stack memory in psi_error',0, MA_ERR) 5918 5919 5920 error = 0.0d0 5921 size = dble(ne(1)+ne(2)) 5922 do n=1, (neq(1)+neq(2)) 5923 do k=1,npack1 5924 dcpl_mb(tmp1(1)+k-1) = dcpl_mb(psi2(1)+k-1+(n-1)*npack1) 5925 > - dcpl_mb(psi1(1)+k-1+(n-1)*npack1) 5926 end do 5927c call D3dB_cc_dot(1,dcpl_mb(tmp1(1)),dcpl_mb(tmp1(1)),sum) 5928 call Pack_cc_dot(1,dcpl_mb(tmp1(1)),dcpl_mb(tmp1(1)),sum) 5929 5930 error = error + sum 5931 end do 5932 call D1dB_SumAll(error) 5933 error = dsqrt(error)/size 5934 5935 value = BA_pop_stack(tmp1(2)) 5936 if (.not. value) 5937 > call errquit('error popping stack memory in psi_error',0, MA_ERR) 5938 5939 5940 psi_error = error 5941 return 5942 end 5943 5944* *************************** 5945* * * 5946* * rho_error * 5947* * * 5948* *************************** 5949 real*8 function rho_error() 5950 implicit none 5951 5952#include "errquit.fh" 5953#include "bafdecls.fh" 5954#include "psi.fh" 5955 5956* ***** local variables **** 5957 logical value 5958 integer k,nx,ny,nz 5959 real*8 error,scale 5960 integer tmp1(2) 5961 5962 real*8 e1,e2 5963 common /eenergy_tmp_common/ e1,e2 5964 5965* ***** external functions ***** 5966 real*8 lattice_omega 5967 external lattice_omega 5968 5969 value = BA_push_get(mt_dbl,(2*nfft3d),'tmp1',tmp1(2),tmp1(1)) 5970 if (.not. value) 5971 > call errquit('out of stack memory in rho_error',0, MA_ERR) 5972 5973 5974 call D3dB_nx(1,nx) 5975 call D3dB_ny(1,ny) 5976 call D3dB_nz(1,nz) 5977 scale = lattice_omega() 5978 5979 scale = (scale)/dble(nx*ny*nz) 5980* scale = (scale)/dble(nx*ny*nz) 5981* scale = (scale*scale) 5982 5983!$OMP DO private(k) 5984 do k=1,(2*nfft3d) 5985 dbl_mb(tmp1(1)+k-1) = (dbl_mb(rho2(1)+k-1) 5986 > -dbl_mb(rho1(1)+k-1)) 5987 dbl_mb(tmp1(1)+k-1) = dbl_mb(tmp1(1)+k-1) 5988 > + (dbl_mb(rho2(1)+k-1+(ispin-1)*(2*nfft3d)) 5989 > -dbl_mb(rho1(1)+k-1+(ispin-1)*(2*nfft3d))) 5990 end do 5991!$OMP END DO 5992 call D3dB_rr_dot(1,dbl_mb(tmp1(1)),dbl_mb(tmp1(1)),e1) 5993 error = e1*scale 5994* error = dsqrt(error) 5995 5996 value = BA_pop_stack(tmp1(2)) 5997 if (.not. value) 5998 > call errquit('error popping stack memory in rho_error',0, MA_ERR) 5999 6000 6001 rho_error = error 6002 return 6003 end 6004 6005 6006* *************************** 6007* * * 6008* * psi_a_sum * 6009* * * 6010* *************************** 6011 real*8 function psi_a_sum(npack1,psi) 6012 implicit none 6013 integer npack1 6014 complex*16 psi(*) 6015 6016 integer k 6017 real*8 a,tmp 6018 6019 a = 0.0d0 6020 do k=1,npack1 6021 tmp = dble(psi(k)) 6022 a = a + tmp*tmp 6023 end do 6024 call D3dB_SumAll(a) 6025 6026 psi_a_sum = a 6027 return 6028 end 6029 6030 6031 6032 6033* *************************** 6034* * * 6035* * psi_b_sum * 6036* * * 6037* *************************** 6038 real*8 function psi_b_sum(npack1,psi) 6039 implicit none 6040 integer npack1 6041 complex*16 psi(*) 6042 6043 6044 integer k 6045 real*8 b,tmp 6046 6047 b = 0.0d0 6048 do k=1,npack1 6049 tmp = dimag(psi(k)) 6050 b = b + tmp*tmp 6051 end do 6052 call D3dB_SumAll(b) 6053 6054 psi_b_sum = b 6055 return 6056 end 6057 6058* ************************************** 6059* * * 6060* * psi_symm_project * 6061* * * 6062* ************************************** 6063 subroutine psi_a_project(npack1,psi) 6064 implicit none 6065 integer npack1 6066 complex*16 psi(*) 6067 integer k 6068 real*8 tmp 6069 do k=1,npack1 6070 tmp = dble(psi(k)) 6071 psi(k) = dcmplx(tmp,0.0d0) 6072 end do 6073 return 6074 end 6075 subroutine psi_b_project(npack1,psi) 6076 implicit none 6077 integer npack1 6078 complex*16 psi(*) 6079 integer k 6080 real*8 tmp 6081 do k=1,npack1 6082 tmp = dimag(psi(k)) 6083 psi(k) = dcmplx(0.0d0,tmp) 6084 end do 6085 return 6086 end 6087 6088 subroutine psi_symm_project(ispin,ne,npack1,psi1) 6089 implicit none 6090 integer ispin,ne(2),npack1 6091 complex*16 psi1(npack1,*) 6092 6093 integer i 6094 real*8 a,b 6095 real*8 psi_a_sum,psi_b_sum 6096 external psi_a_sum,psi_b_sum 6097 6098 do i=1,(ne(1)+ne(2)) 6099 a = psi_a_sum(npack1,psi1(1,i)) 6100 b = psi_b_sum(npack1,psi1(1,i)) 6101 if (a.ge.b) then 6102 call psi_a_project(npack1,psi1(1,i)) 6103 else 6104 call psi_b_project(npack1,psi1(1,i)) 6105 end if 6106 end do 6107 return 6108 end 6109 6110 subroutine psi_ab_gen_irrep_names(virtual) 6111 implicit none 6112 logical virtual 6113 6114#include "bafdecls.fh" 6115#include "errquit.fh" 6116#include "psi.fh" 6117 6118 integer irreps(2) 6119 common / ab_irrep / irreps 6120 6121 integer k,n,ptr,nn 6122 real*8 a,b,tmpa,tmpb 6123 6124 if (virtual) then 6125 ptr = psi1_excited(1) 6126 nn = ne_excited(1)+ne_excited(2) 6127 else 6128 ptr = psi1(1) 6129 nn = ne(1)+ne(2) 6130 end if 6131 6132 if (.not.BA_alloc_get(mt_int,nn, 6133 > 'irreps',irreps(2),irreps(1))) 6134 > call errquit('psi_ab_gen_irrep_names',0, MA_ERR) 6135 6136 do n=1,nn 6137 a = 0.0d0 6138 b = 0.0d0 6139 do k=1,npack1 6140 tmpa = dble( dcpl_mb(ptr+k-1+(n-1)*npack1)) 6141 tmpb = dimag(dcpl_mb(ptr+k-1+(n-1)*npack1)) 6142 a = a + tmpa*tmpa 6143 b = b + tmpb*tmpb 6144 end do 6145 call D3dB_SumAll(a) 6146 call D3dB_SumAll(b) 6147 6148 if ((b .lt. 1.0d-6).and.(a .gt. 1.0d-6)) then 6149 int_mb(irreps(1)+n-1) = 1 6150 else if ((a .lt. 1.0d-6).and.(b .gt. 1.0d-6)) then 6151 int_mb(irreps(1)+n-1) = -1 6152 else 6153 int_mb(irreps(1)+n-1) = 0 6154 end if 6155 end do 6156 6157 6158 return 6159 end 6160 6161 subroutine psi_ab_kill_irrep_names() 6162 implicit none 6163 6164#include "bafdecls.fh" 6165#include "errquit.fh" 6166 6167 integer irreps(2) 6168 common / ab_irrep / irreps 6169 6170 if (.not.BA_free_heap(irreps(2))) 6171 > call errquit('psi_ab_gen_irrep_names: error freeing heap', 6172 > 0, MA_ERR) 6173 6174 return 6175 end 6176 6177 6178 6179* ************************************** 6180* * * 6181* * psi_ab_irrep_name * 6182* * * 6183* ************************************** 6184 6185* This function resturns 6186* '[ag]' - if psi(n) is purely real 6187* '[au]' - if psi(n) is purely imaginary 6188* ' ' - if psi(n) is mixed 6189* 6190* Not psi_ab_gen_irrep_names needs to be called before this is used. 6191* 6192 character*4 function psi_ab_irrep_name(n) 6193 implicit none 6194 integer n 6195 6196#include "bafdecls.fh" 6197 6198 integer irreps(2) 6199 common / ab_irrep / irreps 6200 6201 character*4 abvalue 6202 6203 if (int_mb(irreps(1)+n-1).eq.1) then 6204 abvalue = '[ag]' 6205 else if (int_mb(irreps(1)+n-1).eq.-1) then 6206 abvalue = '[au]' 6207 else 6208 abvalue = ' ' 6209 end if 6210 6211 psi_ab_irrep_name = abvalue 6212 return 6213 end 6214 6215 6216* *************************** 6217* * * 6218* * psi1_crystal_dipole * 6219* * * 6220* *************************** 6221* 6222* Uses - electron_crystal_dipole 6223* 6224 subroutine psi1_crystal_dipole(dipole) 6225 implicit none 6226 real*8 dipole(3) 6227 6228#include "bafdecls.fh" 6229#include "psi.fh" 6230 6231 call Calculate_Resta_Dipole(.true.,ispin,ne,neq,npack1,nfft3d, 6232 > dcpl_mb(psi1(1)),dipole) 6233 6234 return 6235 end 6236 6237 6238 6239* *************************** 6240* * * 6241* * rho_dipole * 6242* * * 6243* *************************** 6244* 6245* Uses - Calculate_dipole (pspw/lib/psi/dipole.f) 6246* 6247 subroutine rho_dipole(dipole) 6248 implicit none 6249 real*8 dipole(3) 6250 6251#include "bafdecls.fh" 6252#include "psi.fh" 6253 6254 call Calculate_Dipole(ispin,ne,2*nfft3d,dbl_mb(rho1(1)),dipole) 6255 return 6256 end 6257 6258 6259* *************************** 6260* * * 6261* * psi_ispin * 6262* * * 6263* *************************** 6264 integer function psi_ispin() 6265 implicit none 6266 6267#include "bafdecls.fh" 6268#include "psi.fh" 6269 6270 psi_ispin = ispin 6271 return 6272 end 6273 6274 6275* *************************** 6276* * * 6277* * psi_ne * 6278* * * 6279* *************************** 6280 integer function psi_ne(ms) 6281 implicit none 6282 integer ms 6283 6284#include "bafdecls.fh" 6285#include "psi.fh" 6286 6287 psi_ne = ne(ms) 6288 return 6289 end 6290 6291* *************************** 6292* * * 6293* * psi_neq * 6294* * * 6295* *************************** 6296 integer function psi_neq(ms) 6297 implicit none 6298 integer ms 6299 6300#include "bafdecls.fh" 6301#include "psi.fh" 6302 6303 psi_neq = neq(ms) 6304 return 6305 end 6306 6307 6308 6309* *************************** 6310* * * 6311* * psi_cpmd_start * 6312* * * 6313* *************************** 6314 subroutine psi_cpmd_start() 6315 implicit none 6316 6317#include "bafdecls.fh" 6318#include "errquit.fh" 6319#include "psi.fh" 6320 6321 logical value 6322 6323 value = BA_alloc_get(mt_dbl,2*ispin*nfft3d,'rho0',rho0(2),rho0(1)) 6324 if (.not.value) 6325 > call errquit('psi_cpmd_start',0,MA_ERR) 6326 6327 call dcopy(2*ispin*nfft3d,dbl_mb(rho1(1)),1,dbl_mb(rho0(1)),1) 6328 return 6329 end 6330 6331* *************************** 6332* * * 6333* * psi_cpmd_end * 6334* * * 6335* *************************** 6336 subroutine psi_cpmd_end() 6337 implicit none 6338 6339#include "bafdecls.fh" 6340#include "errquit.fh" 6341#include "psi.fh" 6342 6343 if (.not.BA_free_heap(rho0(2))) 6344 > call errquit('psi_cpmd_end',0,MA_ERR) 6345 return 6346 end 6347 6348 6349* *************************** 6350* * * 6351* * psi_cpmd_step * 6352* * * 6353* *************************** 6354 subroutine psi_cpmd_step(dte) 6355 implicit none 6356 real*8 dte 6357 6358#include "bafdecls.fh" 6359#include "errquit.fh" 6360#include "psi.fh" 6361 6362 logical control_precondition 6363 external control_precondition 6364 integer control_ks_algorithm 6365 external control_ks_algorithm 6366 real*8 control_tole 6367 external control_tole 6368 6369 6370* **** psi2 = 2*psi1 - psi0 + dt*dt/fmass*Hpsi **** 6371c call electron_cpmd_update(dcpl_mb(psi0(1)), 6372c > dcpl_mb(psi1(1)), 6373c > dcpl_mb(psi2(1)), 6374c > dbl_mb(hml(1)), 6375c > dte) 6376c call Dneall_f_ortho(0,dcpl_mb(psi2(1)),npack1) 6377c write(*,*) "psi1 ortho:" 6378c call OrthoCheck_geo(ispin,ne,dcpl_mb(psi1(1))) 6379c* **** lagrange multiplier corrections **** 6380c if (pspw_SIC().or.occupation_on) then 6381c call psi_lmbda_sic(ispin,ne,(neq(1)+neq(2)),npack1, 6382c > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dte, 6383c > dbl_mb(lmd_cpmd(1)), 6384c > dbl_mb(tmp_L_cpmd(1)),ierr) 6385c else 6386c call psi_lmbda(ispin,ne,(neq(1)+neq(2)),npack1, 6387c > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dte, 6388c > dbl_mb(lmd_cpmd(1)), 6389c > dbl_mb(tmp_L_cpmd(1)),ierr) 6390c end if 6391c write(*,*) "psi2 ortho:" 6392c call OrthoCheck_geo(ispin,ne,dcpl_mb(psi2(1))) 6393 6394 call dcopy(2*ispin*nfft3d,dbl_mb(rho1(1)),1,dbl_mb(rho2(1)),1) 6395 call dscal(2*ispin*nfft3d, 2.0d0,dbl_mb(rho2(1)),1) 6396 call daxpy(2*ispin*nfft3d,-1.0d0, 6397 > dbl_mb(rho0(1)),1,dbl_mb(rho2(1)),1) 6398 call dcopy(2*ispin*nfft3d,dbl_mb(rho1(1)),1,dbl_mb(rho0(1)),1) 6399 6400 call psi_set_density(1,dbl_mb(rho2(1))) 6401 6402* **** diaganolize KS matrix **** 6403 call psi_KS_update(1, 6404 > control_ks_algorithm(), 6405 > control_precondition(), 6406 > control_tole()) 6407 6408 return 6409 end 6410 6411 6412 6413* *************************** 6414* * * 6415* * psi_initialize * 6416* * * 6417* *************************** 6418 6419 logical function psi_initialize() 6420 implicit none 6421 6422 6423#include "bafdecls.fh" 6424#include "btdb.fh" 6425#include "stdio.fh" 6426#include "util.fh" 6427#include "errquit.fh" 6428#include "psi.fh" 6429 6430* ***** rhoall common block **** 6431 integer rho1_all(2) 6432 integer rho2_all(2) 6433 common / rhoall_block / rho1_all,rho2_all 6434 6435* **** local variables **** 6436 integer MASTER,taskid 6437 parameter (MASTER=0) 6438 logical value,psi_nogrid 6439 integer nemax 6440 real*8 sum1,sum2,sum3 6441 integer hversion,hnfft(3),hispin,hne(2) 6442 real*8 hunita(3,3) 6443 integer rtdb,ind,vers 6444 integer control_rtdb,control_ngrid,control_symmetry 6445 external control_rtdb,control_ngrid,control_symmetry 6446 integer psi_get_version 6447 external psi_get_version 6448 character*50 filename 6449 character*50 control_input_psi 6450 external control_input_psi 6451 logical wvfnc_expander,Dneall_m_allocate,band_reformat_c_wvfnc 6452 external wvfnc_expander,Dneall_m_allocate,band_reformat_c_wvfnc 6453 logical psp_pawexist,control_print 6454 external psp_pawexist,control_print 6455 integer control_fractional_smeartype 6456 double precision control_fractional_kT 6457 external control_fractional_smeartype 6458 external control_fractional_kT 6459 logical control_ortho 6460 external control_ortho 6461 6462 6463 ne_excited(1) = 0 6464 ne_excited(2) = 0 6465 6466* **** reformat wavefunction if it is a band wavefunction **** 6467 vers = psi_get_version() 6468 if (vers.eq.5) then 6469 call Parallel_taskid(taskid) 6470 if (taskid.eq.MASTER) then 6471 value= band_reformat_c_wvfnc(1) 6472 end if 6473 end if 6474 pawexist = psp_pawexist() 6475 6476 6477* ***** get ispin,ne,neq,nfft3d,npack0,npack1 **** 6478 call psi_get_ne_occupation(ispin,ne,smearoccupation) 6479 call Dneall_neq(neq) 6480 call D3dB_nfft3d(1,nfft3d) 6481 call Pack_npack(1,npack1) 6482 call Pack_npack(0,npack0) 6483 nemax = ne(1)+ne(2) 6484 occupation_on = .false. 6485 if (smearoccupation.gt.0) occupation_on = .true. 6486 6487 6488* **** allocate memory **** 6489 value = BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 6490 > 'psi2',psi2(2),psi2(1)) 6491 value = value.and. 6492 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 6493 > 'psi1',psi1(2),psi1(1)) 6494 if (pawexist) then 6495 value = value.and. 6496 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 6497 > 'spsi1',spsi1(2),spsi1(1)) 6498 end if 6499 value = value.and. 6500 > BA_alloc_get(mt_dbl,4*nfft3d, 6501 > 'rho1',rho1(2),rho1(1)) 6502 value = value.and. 6503 > BA_alloc_get(mt_dbl,4*nfft3d, 6504 > 'rho2',rho2(2),rho2(1)) 6505 value = value.and. 6506 > BA_alloc_get(mt_dcpl,npack0, 6507 > 'dng1',dng1(2),dng1(1)) 6508 value = value.and. 6509 > BA_alloc_get(mt_dcpl,npack0, 6510 > 'dng2',dng2(2),dng2(1)) 6511c value = value.and. 6512c > BA_alloc_get(mt_dbl,(2*nemax*nemax),'hml',hml(2),hml(1)) 6513 value = value.and.Dneall_m_allocate(0,hml) 6514 6515 value = value.and. 6516 > BA_alloc_get(mt_dbl,(2*nemax),'eig',eig(2),eig(1)) 6517 6518 if (occupation_on) then 6519 value = value.and. 6520 > BA_alloc_get(mt_dbl,(2*nemax),'occ1',occ1(2),occ1(1)) 6521 value = value.and. 6522 > BA_alloc_get(mt_dbl,(2*nemax),'occ2',occ2(2),occ2(1)) 6523 smeartype = control_fractional_smeartype() 6524 smearkT = control_fractional_kT() 6525 end if 6526 6527 value = value.and. 6528 > BA_alloc_get(mt_dbl,4*nfft3d, 6529 > 'rho1_all',rho1_all(2),rho1_all(1)) 6530 value = value.and. 6531 > BA_alloc_get(mt_dbl,4*nfft3d, 6532 > 'rho2_all',rho2_all(2),rho2_all(1)) 6533 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 6534c call dcopy(4*nfft3d,0.0d0,0,dbl_mb(rho1_all(1)),1) 6535c call dcopy(4*nfft3d,0.0d0,0,dbl_mb(rho2_all(1)),1) 6536 call Parallel_shared_vector_zero(.true.,4*nfft3d, 6537 > dbl_mb(rho1_all(1))) 6538 call Parallel_shared_vector_zero(.true.,4*nfft3d, 6539 > dbl_mb(rho2_all(1))) 6540 6541* ***** read initial wavefunctions into psi1 **** 6542 rtdb = control_rtdb() 6543 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 6544 > mt_log,1,psi_nogrid)) 6545 > psi_nogrid = .true. 6546 6547 if (psi_nogrid) then 6548 6549 call psi_get_header(hversion,hnfft,hunita,hispin,hne) 6550 6551 if ( (hnfft(1).ne.control_ngrid(1)) .or. 6552 > (hnfft(2).ne.control_ngrid(2)) .or. 6553 > (hnfft(3).ne.control_ngrid(3)) ) then 6554 6555 hnfft(1) = control_ngrid(1) 6556 hnfft(2) = control_ngrid(2) 6557 hnfft(3) = control_ngrid(3) 6558 call Parallel_taskid(taskid) 6559 6560 call ga_sync() 6561 value = btdb_parallel(.false.) 6562 call ga_sync() 6563 if (taskid.eq.MASTER) then 6564 6565 filename = control_input_psi() 6566 6567 ind = index(filename,' ') - 1 6568 if (.not. btdb_cput(rtdb,'xpndr:old_wavefunction_filename', 6569 > 1,filename(1:ind))) 6570 > call errquit( 6571 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 6572 6573 if (.not. btdb_cput(rtdb,'xpndr:new_wavefunction_filename', 6574 > 1,filename(1:ind))) 6575 > call errquit( 6576 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 6577 6578 if (.not. btdb_put(rtdb,'xpndr:ngrid',mt_int,3,hnfft)) 6579 > call errquit( 6580 > 'wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 6581 6582 if (control_print(print_medium)) then 6583 write(luout,*) 6584 write(luout,*) "Grid is being converted:" 6585 write(luout,*) "------------------------" 6586 write(luout,*) 6587 write(luout,*) "To turn off automatic grid conversion:" 6588 write(luout,*) 6589 write(luout,*) "set nwpw:psi_nogrid .false." 6590 write(luout,*) 6591 endif 6592 value = wvfnc_expander(rtdb) 6593 6594 end if 6595 call ga_sync() 6596 value = btdb_parallel(.true.) 6597 value = .true. 6598 6599 end if 6600 6601 end if 6602 6603 call psi_read(ispin,ne,dcpl_mb(psi1(1)), 6604 > smearoccupation,dbl_mb(occ1(1))) 6605 6606 if (occupation_on) then 6607 call frac_occ_set(rtdb,ispin,ne,dbl_mb(occ1(1))) 6608 end if 6609 6610 call psi_history_read(ispin,ne, 6611 > dcpl_mb(psi1(1)), 6612 > dcpl_mb(psi2(1))) 6613 6614 6615* **** force inversion symmetry **** 6616 if (control_symmetry().eq.1) then 6617 call Parallel_taskid(taskid) 6618 if ((taskid.eq.MASTER).and. 6619 > (control_print(print_medium))) then 6620 write(luout,*) 6621 write(luout,*) 6622 > "Projecting wavefunctions to have inversion symmetry" 6623 write(luout,*) 6624 end if 6625 call psi_symm_project(ispin,neq,npack1,dcpl_mb(psi1(1))) 6626 end if 6627 6628* **** Ortho Check **** 6629 call psi_1ortho_check_fix() 6630c if (pawexist) then 6631c call psp_overlap_S(ispin,neq,dcpl_mb(psi1(1)),dcpl_mb(psi2(1))) 6632c call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6633c > dcpl_mb(psi1(1)), 6634c > dcpl_mb(psi2(1)), 6635c > sum2) 6636c else 6637c call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6638c > dcpl_mb(psi1(1)), 6639c > dcpl_mb(psi1(1)), 6640c > sum2) 6641c end if 6642c call D1dB_SumAll(sum2) 6643c 6644c 6645c sum1 = dble(ne(1) + ne(2)) 6646c if ((control_ortho()).and.(dabs(sum2-sum1).gt.1.0d-10)) then 6647c call Parallel_taskid(taskid) 6648c if (pawexist) then 6649c call Dneall_f_Sortho(0,dcpl_mb(psi1(1)), 6650c > dcpl_mb(psi2(1)),npack1) 6651c call psp_overlap_S(ispin,neq,dcpl_mb(psi1(1)), 6652c > dcpl_mb(psi2(1))) 6653c call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6654c > dcpl_mb(psi1(1)), 6655c > dcpl_mb(psi2(1)), 6656c > sum3) 6657c else 6658cc call Dneall_f_ortho(0,dcpl_mb(psi1(1)),npack1) 6659c call Dneall_f_GramSchmidt(0,dcpl_mb(psi1(1)),npack1) 6660c call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6661c > dcpl_mb(psi1(1)), 6662c > dcpl_mb(psi1(1)), 6663c > sum3) 6664c end if 6665c call D1dB_SumAll(sum3) 6666c if ((taskid.eq.MASTER).and.(control_print(print_medium))) 6667c > write(luout,*) 6668c > "Warning - Gram-Schmidt being performed on psi:", 6669c > sum1,sum2,sum3,dabs(sum2-sum1) 6670c 6671c 6672c end if 6673 6674 psi_initialize = value 6675 return 6676 end 6677 6678* *************************** 6679* * * 6680* * psi_1ortho_check_fix * 6681* * * 6682* *************************** 6683 subroutine psi_1ortho_check_fix() 6684 implicit none 6685 6686#include "bafdecls.fh" 6687#include "stdio.fh" 6688#include "util.fh" 6689#include "psi.fh" 6690 6691* **** local variables **** 6692 integer taskid,MASTER 6693 parameter (MASTER=0) 6694 6695 real*8 sum2,sum1,sum3 6696 6697* **** external functions **** 6698 logical control_ortho,control_print 6699 external control_ortho,control_print 6700 6701* **** Ortho Check **** 6702 if (pawexist) then 6703 call psp_overlap_S(ispin,neq,dcpl_mb(psi1(1)),dcpl_mb(psi2(1))) 6704 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6705 > dcpl_mb(psi1(1)), 6706 > dcpl_mb(psi2(1)), 6707 > sum2) 6708 else 6709 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6710 > dcpl_mb(psi1(1)), 6711 > dcpl_mb(psi1(1)), 6712 > sum2) 6713 end if 6714 call D1dB_SumAll(sum2) 6715 6716 6717 sum1 = dble(ne(1) + ne(2)) 6718 if ((control_ortho()).and.(dabs(sum2-sum1).gt.1.0d-10)) then 6719 call Parallel_taskid(taskid) 6720 if (pawexist) then 6721 call Dneall_f_Sortho(0,dcpl_mb(psi1(1)), 6722 > dcpl_mb(psi2(1)),npack1) 6723 call psp_overlap_S(ispin,neq,dcpl_mb(psi1(1)), 6724 > dcpl_mb(psi2(1))) 6725 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6726 > dcpl_mb(psi1(1)), 6727 > dcpl_mb(psi2(1)), 6728 > sum3) 6729 else 6730c call Dneall_f_ortho(0,dcpl_mb(psi1(1)),npack1) 6731 call Dneall_f_GramSchmidt(0,dcpl_mb(psi1(1)),npack1) 6732 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 6733 > dcpl_mb(psi1(1)), 6734 > dcpl_mb(psi1(1)), 6735 > sum3) 6736 end if 6737 call D1dB_SumAll(sum3) 6738 if ((taskid.eq.MASTER).and.(control_print(print_medium))) then 6739 write(luout,321) sum1,sum2,sum3,dabs(sum2-sum1) 6740 end if 6741 321 format(/" Warning - K.S. orbitals are not orthonormal. ", 6742 > "Applying Gram-Schmidt orthonormalization."/, 6743 > " - exact norm=",E12.6," norm=",E12.6, 6744 > " corrected norm=",E12.6, 6745 > " (error=",E12.6,")"/) 6746 6747 end if 6748 return 6749 end 6750 6751 6752 6753* *************************** 6754* * * 6755* * psi_tmp_write * 6756* * * 6757* *************************** 6758 subroutine psi_tmp_write() 6759 implicit none 6760 6761#include "bafdecls.fh" 6762#include "psi.fh" 6763 6764* ***** write psi1 wavefunctions **** 6765 call psi_write(ispin,ne,dcpl_mb(psi1(1)), 6766 > smearoccupation,dbl_mb(occ1(1))) 6767 6768 return 6769 end 6770 6771 6772 6773* *************************** 6774* * * 6775* * psi_finalize * 6776* * * 6777* *************************** 6778 6779 logical function psi_finalize(wpsi) 6780 implicit none 6781 logical wpsi 6782 6783#include "errquit.fh" 6784#include "bafdecls.fh" 6785#include "psi.fh" 6786 6787* ***** rhoall common block **** 6788 integer rho1_all(2) 6789 integer rho2_all(2) 6790 common / rhoall_block / rho1_all,rho2_all 6791 6792 6793* **** local variables **** 6794 logical value 6795 6796 logical Dneall_m_free 6797 external Dneall_m_free 6798 6799* ***** write psi1 wavefunctions **** 6800 if (wpsi) then 6801 call psi_write(ispin,ne,dcpl_mb(psi1(1)), 6802 > smearoccupation,dbl_mb(occ1(1))) 6803 call psi_history_write(ispin,ne,dcpl_mb(psi1(1))) 6804 end if 6805 6806 value = BA_free_heap(eig(2)) 6807 value = value.and.Dneall_m_free(hml) 6808 value = value.and.BA_free_heap(dng2(2)) 6809 value = value.and.BA_free_heap(dng1(2)) 6810 value = value.and.BA_free_heap(rho2(2)) 6811 value = value.and.BA_free_heap(rho1(2)) 6812 value = value.and.BA_free_heap(psi2(2)) 6813 value = value.and.BA_free_heap(psi1(2)) 6814 if (pawexist) then 6815 value = value.and.BA_free_heap(spsi1(2)) 6816 end if 6817 value = value.and.BA_free_heap(rho2_all(2)) 6818 value = value.and.BA_free_heap(rho1_all(2)) 6819 if (occupation_on) then 6820 value = value.and.BA_free_heap(occ2(2)) 6821 value = value.and.BA_free_heap(occ1(2)) 6822 end if 6823 6824 if (.not. value) 6825 > call errquit('psi_finalize: error freeing heap',0, MA_ERR) 6826 6827 psi_finalize = value 6828 return 6829 end 6830 6831 6832* *************************** 6833* * * 6834* * psi_ne_excited * 6835* * * 6836* *************************** 6837 integer function psi_ne_excited(ms) 6838 implicit none 6839 integer ms 6840 6841#include "bafdecls.fh" 6842#include "psi.fh" 6843 6844 psi_ne_excited = ne_excited(ms) 6845 return 6846 end 6847 6848* *************************** 6849* * * 6850* * psi_2epsi_gradients * 6851* * * 6852* *************************** 6853 logical function psi_2epsi_gradients(gn) 6854 implicit none 6855 integer gn 6856 6857#include "bafdecls.fh" 6858#include "btdb.fh" 6859#include "psi.fh" 6860#include "errquit.fh" 6861 6862* **** local variables **** 6863 logical value 6864 integer nfac,nex(2),nemax,nexmax,ii,n 6865 real*8 tsum,tsum0 6866 6867 nfac = 3 6868 if (gn.eq.2) then 6869 nfac = 9 6870 else if (gn.eq.3) then 6871 nfac = 19 6872 else 6873 nfac = gn 6874 end if 6875 nex(1) = ne(1)*nfac 6876 nex(2) = ne(2)*nfac 6877 nemax = ne(1) + ne(2) 6878 nexmax = nex(1) + nex(2) 6879 6880* **** allocate memory **** 6881 value = BA_alloc_get(mt_dcpl,npack1*(nexmax), 6882 > 'psi1_excited',psi1_excited(2),psi1_excited(1)) 6883 if (.not.value) 6884 > call errquit('psi_2epsi_gradients: out of heap memory',0,MA_ERR) 6885 6886 do ii=1,nfac 6887 do n=1,nemax 6888 call Pack_cc_multiplegradients(1,ii, 6889 > dcpl_mb(psi1(1)+(n-1)*npack1), 6890 > dcpl_mb(psi1_excited(1)+((n-1)+(ii-1)*nemax)*npack1)) 6891 call Pack_cc_dot(1, 6892 > dcpl_mb(psi1_excited(1)+((n-1)+(ii-1)*nemax)*npack1), 6893 > dcpl_mb(psi1_excited(1)+((n-1)+(ii-1)*nemax)*npack1), 6894 > tsum) 6895 call Pack_cc_dot(1, 6896 > dcpl_mb(psi1(1)+(n-1)*npack1), 6897 > dcpl_mb(psi1(1)+(n-1)*npack1), 6898 > tsum0) 6899 write(*,*) "ii,n,tsum=",ii,n,tsum,tsum0 6900 end do 6901 end do 6902 6903 call epsi_write(ispin,nex,dcpl_mb(psi1_excited(1))) 6904 value = BA_free_heap(psi1_excited(2)) 6905 if (.not.value) 6906 > call errquit('psi_2epsi_gradients:freeing heap',0,MA_ERR) 6907 6908 return 6909 end 6910 6911 6912* *************************** 6913* * * 6914* * epsi_initialize * 6915* * * 6916* *************************** 6917 6918 logical function epsi_initialize() 6919 implicit none 6920 6921#include "bafdecls.fh" 6922#include "btdb.fh" 6923#include "psi.fh" 6924#include "errquit.fh" 6925 6926* **** local variables **** 6927 integer MASTER,taskid 6928 parameter (MASTER=0) 6929 logical value,psi_nogrid 6930 integer nemax,ispin0 6931 6932 integer hversion,hnfft(3),hispin,hne(2) 6933 real*8 hunita(3,3),sum1,sum2 6934 integer rtdb,ind,nn 6935 integer control_rtdb,control_ngrid 6936 external control_rtdb,control_ngrid 6937 character*50 filename 6938 character*50 control_input_epsi 6939 external control_input_epsi 6940 logical wvfnc_expander 6941 external wvfnc_expander 6942 integer control_symmetry 6943 external control_symmetry 6944 6945 6946* ***** get ispin, and ne, and nfft3d **** 6947 call psi_get_ne_excited(ispin0,ne_excited) 6948 nemax = ne_excited(1) + ne_excited(2) 6949 nn = ne_excited(1)*ne_excited(1) + ne_excited(2)*ne_excited(2) 6950 6951* **** allocate memory **** 6952 value = BA_alloc_get(mt_dcpl,npack1*(nemax), 6953 > 'psi2_excited',psi2_excited(2),psi2_excited(1)) 6954 value = value.and. 6955 > BA_alloc_get(mt_dcpl,npack1*(nemax), 6956 > 'psi1_excited',psi1_excited(2),psi1_excited(1)) 6957 value = value.and. 6958 > BA_alloc_get(mt_dbl,(2*nemax),'eig_excited', 6959 > eig_excited(2),eig_excited(1)) 6960 value = value.and. 6961 > BA_alloc_get(mt_dbl,nn,'hml_excited', 6962 > hml_excited(2),hml_excited(1)) 6963 if (.not.value) 6964 > call errquit('epsi_initialize: out of heap memory',0,MA_ERR) 6965 6966 !call dcopy(2*npack1*nemax,0.0d0,0,dcpl_mb(psi2_excited(1)),1) 6967 !call dcopy(2*npack1*nemax,0.0d0,0,dcpl_mb(psi1_excited(1)),1) 6968 !call dcopy(2*nemax,0.0d0,0,dbl_mb(eig_excited(1)),1) 6969 !call dcopy(nn,0.0d0,0,dbl_mb(hml_excited(1)),1) 6970 call Parallel_shared_vector_zero(.true.,2*npack1*nemax, 6971 > dcpl_mb(psi2_excited(1))) 6972 call Parallel_shared_vector_zero(.true.,2*npack1*nemax, 6973 > dcpl_mb(psi1_excited(1))) 6974 call Parallel_shared_vector_zero(.true.,2*nemax, 6975 > dbl_mb(eig_excited(1))) 6976 call Parallel_shared_vector_zero(.true.,nn, 6977 > dbl_mb(hml_excited(1))) 6978 6979 6980* ***** read initial wavefunctions into psi1 **** 6981 rtdb = control_rtdb() 6982 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 6983 > mt_log,1,psi_nogrid)) 6984 > psi_nogrid = .true. 6985 6986 if (psi_nogrid) then 6987 6988 6989 filename = control_input_epsi() 6990 call psi_get_header_filename(filename, 6991 > hversion,hnfft,hunita,hispin,hne) 6992 6993 if ( (hnfft(1).ne.control_ngrid(1)) .or. 6994 > (hnfft(2).ne.control_ngrid(2)) .or. 6995 > (hnfft(3).ne.control_ngrid(3)) ) then 6996 6997 hnfft(1) = control_ngrid(1) 6998 hnfft(2) = control_ngrid(2) 6999 hnfft(3) = control_ngrid(3) 7000 call Parallel_taskid(taskid) 7001 value = btdb_parallel(.false.) 7002 if (taskid.eq.MASTER) then 7003 7004 7005 ind = index(filename,' ') - 1 7006 if (.not. btdb_cput(rtdb,'xpndr:old_wavefunction_filename', 7007 > 1,filename(1:ind))) 7008 > call errquit( 7009 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 7010 7011 if (.not. btdb_cput(rtdb,'xpndr:new_wavefunction_filename', 7012 > 1,filename(1:ind))) 7013 > call errquit( 7014 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 7015 7016 if (.not. btdb_put(rtdb,'xpndr:ngrid',mt_int,3,hnfft)) 7017 > call errquit( 7018 > 'wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 7019 7020 write(*,*) 7021 write(*,*) "Grid is being converted:" 7022 write(*,*) "------------------------" 7023 write(*,*) 7024 write(*,*) "To turn off automatic grid conversion:" 7025 write(*,*) 7026 write(*,*) "set nwpw:psi_nogrid .false." 7027 write(*,*) 7028 value = wvfnc_expander(rtdb) 7029 7030 end if 7031 value = btdb_parallel(.true.) 7032 7033 end if 7034 7035 end if 7036 7037* ***** read initial wavefunctions into psi1 **** 7038 call epsi_read(ispin0,ne_excited,dcpl_mb(psi1_excited(1))) 7039 7040 7041* **** force inversion symmetry **** 7042 if (control_symmetry().eq.1) then 7043 call Parallel_taskid(taskid) 7044 if (taskid.eq.MASTER) then 7045 write(*,*) 7046 write(*,*) 7047 > "Projecting virtual wavefunctions to have inversion symmetry" 7048 write(*,*) 7049 end if 7050 call psi_symm_project(ispin0,ne_excited,npack1, 7051 > dcpl_mb(psi1_excited(1))) 7052 end if 7053 7054 7055c* **** Ortho Check **** 7056c call Grsm_gg_trace(npack1,(ne_excited(1)+ne_excited(2)), 7057c > dcpl_mb(psi1_excited(1)), 7058c > dcpl_mb(psi1_excited(1)), 7059c > sum2) 7060c 7061c sum1 = dble(ne_excited(1) + ne_excited(2)) 7062c if (dabs(sum2-sum1).gt.1.0d-10) then 7063c 7064c call Parallel_taskid(taskid) 7065c call Grsm_g_MakeOrtho(npack1,ne_excited(1), 7066c > dcpl_mb(psi1_excited(1))) 7067c if (ispin.gt.1) then 7068c call Grsm_g_MakeOrtho(npack1,ne_excited(2), 7069c > dcpl_mb(psi1_excited(1) 7070c > +ne_excited(1)*npack1)) 7071c end if 7072c call Grsm_gg_trace(npack1,(ne_excited(1)+ne_excited(2)), 7073c > dcpl_mb(psi1_excited(1)), 7074c > dcpl_mb(psi1_excited(1)), 7075c > sum2) 7076c if (taskid.eq.MASTER) 7077c > write(*,*) "Warning - Gram-Schmidt being performed on epsi:", 7078c > dabs(sum2-sum1) 7079c 7080c end if 7081 7082 7083 epsi_initialize = value 7084 return 7085 end 7086 7087 7088 7089* *************************** 7090* * * 7091* * epsi_finalize * 7092* * * 7093* *************************** 7094 7095 logical function epsi_finalize(writepsi) 7096 implicit none 7097 logical writepsi 7098 7099#include "bafdecls.fh" 7100#include "errquit.fh" 7101#include "psi.fh" 7102 7103* **** local variables **** 7104 logical value 7105 7106* ***** write psi1 wavefunctions **** 7107 if (writepsi) 7108 > call epsi_write(ispin,ne_excited,dcpl_mb(psi1_excited(1))) 7109 7110 value = BA_free_heap(hml_excited(2)) 7111 value = value.and.BA_free_heap(eig_excited(2)) 7112 value = value.and.BA_free_heap(psi2_excited(2)) 7113 value = value.and.BA_free_heap(psi1_excited(2)) 7114 if (.not. value) 7115 > call errquit('epsi_finalize: error freeing heap',0, MA_ERR) 7116 7117 epsi_finalize = value 7118 return 7119 end 7120 7121 7122* *********************** 7123* * * 7124* * psi_Mulliken * 7125* * * 7126* *********************** 7127 7128 subroutine psi_Mulliken(rtdb) 7129 implicit none 7130 integer rtdb 7131 7132#include "bafdecls.fh" 7133#include "psi.fh" 7134#include "stdio.fh" 7135 7136 7137* **** Lubin Water Analysis **** 7138 call pspw_Lubin_water_analysis(rtdb,ispin,ne,2*nfft3d, 7139 > dbl_mb(rho1(1))) 7140 7141* **** Atom Analysis **** 7142 call pspw_atom_analysis(rtdb,ispin,2*nfft3d,dbl_mb(rho1(1))) 7143 7144* **** Mulliken Analysis **** 7145 call pspw_analysis(0,rtdb,ispin,ne,dcpl_mb(psi2(1)), 7146 > dbl_mb(eig(1))) 7147 7148 call pspw_gen_APC(ispin,ne,dcpl_mb(dng1(1))) 7149 call pspw_print_APC(luout) 7150 7151 call pspw_gen_Efield(rtdb,ispin,dbl_mb(rho1(1)),dcpl_mb(dng1(1))) 7152 7153 call pspw_gen_Efield_grad(rtdb,ispin,ne, 7154 > dcpl_mb(psi1(1)), 7155 > dcpl_mb(dng1(1))) 7156 7157 return 7158 end 7159 7160* *********************** 7161* * * 7162* * epsi_Mulliken * 7163* * * 7164* *********************** 7165 7166 subroutine epsi_Mulliken(rtdb) 7167 implicit none 7168 integer rtdb 7169 7170#include "bafdecls.fh" 7171#include "psi.fh" 7172 7173 call pspw_analysis(1,rtdb,ispin,ne_excited, 7174 > dcpl_mb(psi1_excited(1)), 7175 > dbl_mb(eig_excited(1))) 7176 return 7177 end 7178 7179 7180 7181 7182* *********************** 7183* * * 7184* * psi_DOS * 7185* * * 7186* *********************** 7187 7188 subroutine psi_DOS(rtdb) 7189 implicit none 7190 integer rtdb 7191 7192#include "btdb.fh" 7193#include "bafdecls.fh" 7194#include "psi.fh" 7195#include "errquit.fh" 7196 7197 7198* **** local variables **** 7199 logical value 7200 integer npoints,ii 7201 integer weight(2),nemax 7202 real*8 emin,emax,alpha 7203 character*255 filename 7204 7205 nemax = ne(1) 7206 value = BA_push_get(mt_dbl,(nemax),'weight',weight(2),weight(1)) 7207 if (.not. value) 7208 > call errquit('psi_dos:out of stack memory',0, MA_ERR) 7209 call dcopy(nemax,1.0d0,0,dbl_mb(weight(1)),1) 7210 7211 7212 if (.not.btdb_get(rtdb,'dos:alpha',mt_dbl,1,alpha)) then 7213 alpha = 0.05d0/27.2116d0 7214 end if 7215 7216 if (.not.btdb_get(rtdb,'dos:npoints',mt_int,1,npoints)) then 7217 npoints = 500 7218 end if 7219 7220 if (.not.btdb_get(rtdb,'dos:emin',mt_dbl,1,emin)) then 7221 emin = 99999.0d0 7222 do ii=1,ne(1)+ne(2) 7223 if (dbl_mb(eig(1)+ii-1).lt.emin) emin = dbl_mb(eig(1)+ii-1) 7224 end do 7225 emin = emin - 0.1d0 7226 end if 7227 7228 if (.not.btdb_get(rtdb,'dos:emax',mt_dbl,1,emax)) then 7229 emax = -99999.0d0 7230 do ii=1,ne(1)+ne(2) 7231 if (dbl_mb(eig(1)+ii-1).gt.emax) emax = dbl_mb(eig(1)+ii-1) 7232 end do 7233 emax = emax + 0.1d0 7234 end if 7235 7236* **** generate DENSITY OF STATES ***** 7237 if (ispin.eq.1) then 7238 filename = "smear_dos_both" 7239 call densityofstates(filename,.false., 7240 > dbl_mb(eig(1)),dbl_mb(weight(1)),ne(1), 7241 > 1.0d0,alpha,npoints,emin,emax) 7242 filename = "smear_fdos_both" 7243 call densityofstates(filename,.false., 7244 > dbl_mb(eig(1)),dbl_mb(weight(1)),ne(1), 7245 > 1.0d0,alpha,npoints,emin,emax) 7246 end if 7247 7248 if (ispin.eq.2) then 7249 filename = "smear_dos_alpha" 7250 call densityofstates(filename,.false., 7251 > dbl_mb(eig(1)),dbl_mb(weight(1)),ne(1), 7252 > 1.0d0,alpha,npoints,emin,emax) 7253 filename = "smear_dos_beta" 7254 call densityofstates(filename,.false., 7255 > dbl_mb(eig(1)+ne(1)),dbl_mb(weight(1)),ne(2), 7256 > -1.0d0,alpha,npoints,emin,emax) 7257 filename = "smear_fdos_alpha" 7258 call densityofstates(filename,.false., 7259 > dbl_mb(eig(1)),dbl_mb(weight(1)),ne(1), 7260 > 1.0d0,alpha,npoints,emin,emax) 7261 filename = "smear_fdos_beta" 7262 call densityofstates(filename,.false., 7263 > dbl_mb(eig(1)+ne(1)),dbl_mb(weight(1)),ne(2), 7264 > -1.0d0,alpha,npoints,emin,emax) 7265 end if 7266 7267 7268 value = BA_pop_stack(weight(2)) 7269 if (.not. value) 7270 > call errquit('psi_dos: error freeing stack',0, MA_ERR) 7271 7272 return 7273 end 7274 7275* *********************** 7276* * * 7277* * epsi_DOS * 7278* * * 7279* *********************** 7280 7281 subroutine epsi_DOS(rtdb) 7282 implicit none 7283 integer rtdb 7284 7285#include "btdb.fh" 7286#include "bafdecls.fh" 7287#include "psi.fh" 7288#include "errquit.fh" 7289 7290 7291* **** local variables **** 7292 logical value 7293 integer weight(2),npoints,ii 7294 real*8 emin,emax,alpha 7295 character*255 filename 7296 7297 value = BA_push_get(mt_dbl,(ne_excited(1)+ne_excited(2)), 7298 > 'weight',weight(2),weight(1)) 7299 if (.not. value) 7300 > call errquit('epsi_dos:out of stack memory',0, MA_ERR) 7301 call dcopy((ne_excited(1)+ne_excited(2)), 7302 > 1.0d0,0,dbl_mb(weight(1)),1) 7303 7304 if (.not.btdb_get(rtdb,'dos:alpha',mt_dbl,1,alpha)) then 7305 alpha = 0.05d0/27.2116d0 7306 end if 7307 7308 if (.not.btdb_get(rtdb,'dos:npoints',mt_int,1,npoints)) then 7309 npoints = 500 7310 end if 7311 7312 if (.not.btdb_get(rtdb,'dos:emin',mt_dbl,1,emin)) then 7313 emin = 99999.0d0 7314 do ii=1,ne_excited(1)+ne_excited(2) 7315 if (dbl_mb(eig_excited(1)+ii-1).lt.emin) 7316 > emin = dbl_mb(eig_excited(1)+ii-1) 7317 end do 7318 emin = emin - 0.1d0 7319 end if 7320 7321 if (.not.btdb_get(rtdb,'dos:emax',mt_dbl,1,emax)) then 7322 emax = -99999.0d0 7323 do ii=1,ne_excited(1)+ne_excited(2) 7324 if (dbl_mb(eig_excited(1)+ii-1).gt.emax) 7325 > emax = dbl_mb(eig_excited(1)+ii-1) 7326 end do 7327 emax = emax + 0.1d0 7328 end if 7329 7330* **** generate DENSITY OF STATES ***** 7331 if (ispin.eq.1) then 7332 filename = "smear_vdos_both" 7333 call densityofstates(filename,.false., 7334 > dbl_mb(eig_excited(1)),dbl_mb(weight(1)), 7335 > ne_excited(1), 7336 > 1.0d0,alpha,npoints,emin,emax) 7337 filename = "smear_dos_both" 7338 call densityofstates(filename,.true., 7339 > dbl_mb(eig_excited(1)),dbl_mb(weight(1)), 7340 > ne_excited(1), 7341 > 1.0d0,alpha,npoints,emin,emax) 7342 end if 7343 7344 if (ispin.eq.2) then 7345 filename = "smear_vdos_alpha" 7346 call densityofstates(filename,.false., 7347 > dbl_mb(eig_excited(1)),dbl_mb(weight(1)), 7348 > ne_excited(1), 7349 > 1.0d0,alpha,npoints,emin,emax) 7350 filename = "smear_vdos_beta" 7351 call densityofstates(filename,.false., 7352 > dbl_mb(eig_excited(1)+ne_excited(1)),dbl_mb(weight(1)), 7353 > ne_excited(2), 7354 > -1.0d0,alpha,npoints,emin,emax) 7355 filename = "smear_dos_alpha" 7356 call densityofstates(filename,.true., 7357 > dbl_mb(eig_excited(1)),dbl_mb(weight(1)), 7358 > ne_excited(1), 7359 > 1.0d0,alpha,npoints,emin,emax) 7360 filename = "smear_dos_beta" 7361 call densityofstates(filename,.true., 7362 > dbl_mb(eig_excited(1)+ne_excited(1)),dbl_mb(weight(1)), 7363 > ne_excited(2), 7364 > -1.0d0,alpha,npoints,emin,emax) 7365 end if 7366 7367 value = BA_pop_stack(weight(2)) 7368 if (.not. value) 7369 > call errquit('epsi_dos: error freeing stack',0, MA_ERR) 7370 7371 return 7372 end 7373 7374cccccccccccccccccccccccccccccccccccccccccccccccccccccc 7375 subroutine psi_polariz() 7376 implicit none 7377#include "psi.fh" 7378#include "bafdecls.fh" 7379#include "errquit.fh" 7380#include "util.fh" 7381 integer psirx(2),asize 7382 logical val 7383 asize=(ne(1)+ne(2))*nfft3d*2 7384 val=BA_push_get(mt_dbl,asize,"psir",psirx(2),psirx(1)) 7385 if (.not.val) then 7386 call errquit("psi_polariz stack empty",0,0) 7387 end if 7388 call berry_phase_pol(psi2,psirx) 7389 val=BA_pop_stack(psirx(2)) 7390 if (.not.val) then 7391 call errquit("psi_polariz failed pop stack",0,0) 7392 end if 7393 return 7394 end 7395ccccccccccccccccccccccccccccccccccccccccccccccccccc 7396 7397 7398 subroutine epsi_generate_kb_vnm(vnm) 7399 implicit none 7400 real*8 vnm(*) 7401 7402#include "bafdecls.fh" 7403#include "errquit.fh" 7404#include "psi.fh" 7405 7406* **** local variables **** 7407 logical value 7408 integer ms,i,j,vshift,ishift,neall(2) 7409 integer Gx(2),Gy(2),Gz(2),tmp(2),psii_ptr,psij_ptr 7410 real*8 vv(3) 7411 7412* **** external functions **** 7413 integer Dneall_mne_size,G_indx 7414 external Dneall_mne_size,G_indx 7415 7416 neall(1) = ne(1)+ne_excited(1) 7417 neall(2) = ne(2)+ne_excited(2) 7418 7419 value = BA_push_get(mt_dbl, nfft3d,'Gx',Gx(2),Gx(1)) 7420 value = value.and. 7421 > BA_push_get(mt_dbl, nfft3d,'Gy',Gy(2),Gy(1)) 7422 value = value.and. 7423 > BA_push_get(mt_dbl, nfft3d,'Gz',Gz(2),Gz(1)) 7424 value = value.and. 7425 > BA_push_get(mt_dbl, 2*npack1,'tmp',tmp(2),tmp(1)) 7426 if (.not. value) 7427 > call errquit('epsi_generate_kb_vnm:pushing stack',1,MA_ERR) 7428 7429 7430* **** define Gx,Gy and Gz in packed space **** 7431 call D3dB_t_Copy(1,dbl_mb(G_indx(1)),dbl_mb(Gx(1))) 7432 call D3dB_t_Copy(1,dbl_mb(G_indx(2)),dbl_mb(Gy(1))) 7433 call D3dB_t_Copy(1,dbl_mb(G_indx(3)),dbl_mb(Gz(1))) 7434 call Pack_t_pack(1,dbl_mb(Gx(1))) 7435 call Pack_t_pack(1,dbl_mb(Gy(1))) 7436 call Pack_t_pack(1,dbl_mb(Gz(1))) 7437 7438 7439 vshift = Dneall_mne_size(0,neall) 7440 !call dcopy(3*vshift,0.0d0,0,vnm,1) 7441 call Parallel_shared_vector_zero(.true.,3*vshift,vnm) 7442 7443 do ms=1,ispin 7444 ishift = (ms-1)*ne(1) 7445 7446 do j=1,neall(ms) 7447 if (j.le.ne(ms)) then 7448 psij_ptr =(j-1+ishift)*npack1 + psi1(1) 7449 else 7450 psij_ptr =(j-ne(ms)-1+ishift)*npack1 + psi1_excited(1) 7451 end if 7452 7453 do i=j,neall(ms) 7454 if (i.le.ne(ms)) then 7455 psii_ptr =(i-1+ishift)*npack1 + psi1(1) 7456 else 7457 psii_ptr =(i-ne(ms)-1+ishift)*npack1 + psi1_excited(1) 7458 end if 7459 7460 call Pack_tc_Mul(1, 7461 > dbl_mb(Gx(1)), 7462 > dcpl_mb(psij_ptr), 7463 > dbl_mb(tmp(1))) 7464 call Pack_cc_dot(1, 7465 > dcpl_mb(psii_ptr), 7466 > dbl_mb(tmp(1)), 7467 > vv(1)) 7468 call Pack_tc_Mul(1, 7469 > dbl_mb(Gy(1)), 7470 > dcpl_mb(psij_ptr), 7471 > dbl_mb(tmp(1))) 7472 call Pack_cc_dot(1, 7473 > dcpl_mb(psii_ptr), 7474 > dbl_mb(tmp(1)), 7475 > vv(2)) 7476 call Pack_tc_Mul(1, 7477 > dbl_mb(Gz(1)), 7478 > dcpl_mb(psij_ptr), 7479 > dbl_mb(tmp(1))) 7480 call Pack_cc_dot(1, 7481 > dcpl_mb(psii_ptr), 7482 > dbl_mb(tmp(1)), 7483 > vv(3)) 7484 7485 7486 call Dneall_mne_set_value(vv(1),0,neall,ms,i,j,vnm) 7487 call Dneall_mne_set_value(vv(2),0,neall,ms,i,j, 7488 > vnm(1+vshift)) 7489 call Dneall_mne_set_value(vv(3),0,neall,ms,i,j, 7490 > vnm(1+vshift+vshift)) 7491 if (i.ne.j) then 7492 call Dneall_mne_set_value(vv(1),0,neall,ms,j,i,vnm) 7493 call Dneall_mne_set_value(vv(2),0,neall,ms,j,i, 7494 > vnm(1+vshift)) 7495 call Dneall_mne_set_value(vv(3),0,neall,ms,j,i, 7496 > vnm(1+vshift+vshift)) 7497 end if 7498 end do 7499 end do 7500 end do 7501 7502 value = BA_pop_stack(tmp(2)) 7503 value = value.and.BA_pop_stack(Gz(2)) 7504 value = value.and.BA_pop_stack(Gy(2)) 7505 value = value.and.BA_pop_stack(Gx(2)) 7506 if (.not. value) 7507 > call errquit('epsi_generate_kb_vnm:popping stack',1,MA_ERR) 7508 7509 7510 return 7511 end 7512 7513 7514 7515* ********************************* 7516* * * 7517* * psi_1pressure_stress * 7518* * * 7519* ********************************* 7520 7521 subroutine psi_1pressure_stress(pressure,p1,p2,stress) 7522 implicit none 7523 real*8 pressure,p1,p2,stress(3,3) 7524 7525#include "bafdecls.fh" 7526#include "errquit.fh" 7527#include "psi.fh" 7528 7529* ***** rhoall common block **** 7530 integer rho1_all(2) 7531 integer rho2_all(2) 7532 common / rhoall_block / rho1_all,rho2_all 7533 7534* ***** external functions ***** 7535 integer electron_xcp_ptr 7536 external electron_xcp_ptr 7537 real*8 psi_1vnl,rho_1exc,rho_1pxc 7538 external psi_1vnl,rho_1exc,rho_1pxc 7539 7540 call cgsd_pressure_stress(ispin,neq, 7541 > dcpl_mb(psi1(1)), 7542 > dbl_mb(rho1_all(1)), 7543 > dcpl_mb(dng1(1)), 7544 > dbl_mb(electron_xcp_ptr()), 7545 > psi_1vnl(),rho_1exc(),rho_1pxc(), 7546 > pressure,p1,p2,stress) 7547 7548 7549 return 7550 end 7551 7552 7553 7554* *********************** 7555* * * 7556* * psi_MP2_energy * 7557* * * 7558* *********************** 7559 subroutine psi_MP2_energy(rtdb) 7560 implicit none 7561 integer rtdb 7562 7563#include "stdio.fh" 7564#include "bafdecls.fh" 7565#include "errquit.fh" 7566#include "psi.fh" 7567 7568* **** local variables **** 7569 integer taskid,MASTER 7570 parameter (MASTER=0) 7571 7572 logical oprint,value 7573 integer ms,a,b,r,s,n2ft3d,icount 7574 real*8 ea,eb,er,es 7575 real*8 e2,d2,tmp2 7576 integer vpsi_r(2),v1h(2),v2h(2) 7577 7578* **** allocate memory from heap **** 7579 n2ft3d = 2*nfft3d 7580 value = BA_alloc_get(mt_dcpl,nfft3d,'v1h',v1h(2),v1h(1)) 7581 value = value.and. 7582 > BA_alloc_get(mt_dcpl,nfft3d,'v2h',v2h(2),v2h(1)) 7583 value = value.and. 7584 > BA_alloc_get(mt_dbl,n2ft3d*(ne_excited(1)+ne_excited(2)), 7585 > 'vpsi_r',vpsi_r(2),vpsi_r(1)) 7586 if (.not. value) 7587 > call errquit('psi_MP2_energy: error allocating heap memory',0, 7588 & MA_ERR) 7589 7590 7591 call Parallel_taskid(taskid) 7592 oprint = (taskid.eq.MASTER) 7593 7594 icount = 0 7595 do ms=1,ispin 7596 do a=1,ne(ms)-1 7597 ea = dbl_mb(eig(1)+a-1+(ms-1)*ne(1)) 7598 do r=1,ne_excited(ms)-1 7599 er = dbl_mb(eig_excited(1)+r-1+(ms-1)*ne(1)) 7600 7601 do b=a+1,ne(ms) 7602 eb = dbl_mb(eig(1)+b-1+(ms-1)*ne(1)) 7603 do s=r+1,ne_excited(ms) 7604 es = dbl_mb(eig_excited(1)+s-1+(ms-1)*ne(1)) 7605 d2 = ea+eb-er-es 7606 tmp2 = 1.0d0 7607 e2 = e2 + tmp2*tmp2/d2 7608 icount = icount + 1 7609 if (oprint) 7610 > write(luout,*) "a,b,r,s, Esub=",a,b,r,s,(tmp2*tmp2/d2) 7611 end do 7612 end do 7613 end do 7614 end do 7615 end do 7616 if (ispin.eq.1) e2=e2+e2 7617 7618* **** deallocate memory from heap **** 7619 value = BA_free_heap(v1h(2)) 7620 > .and.BA_free_heap(v2h(2)) 7621 > .and.BA_free_heap(vpsi_r(2)) 7622 if (.not. value) 7623 > call errquit('psi_MP2_energy: error freeing heap memory',1, 7624 & MA_ERR) 7625 7626 7627 7628 if (oprint) then 7629 write(luout,*) "EMP2 = ", e2, " icount=",icount 7630 end if 7631 7632 return 7633 end 7634 7635 7636 7637* *********************** 7638* * * 7639* * psi_2q_integrals * 7640* * * 7641* *********************** 7642 subroutine psi_2q_integrals(rtdb) 7643 implicit none 7644 integer rtdb 7645 7646#include "stdio.fh" 7647#include "bafdecls.fh" 7648#include "errquit.fh" 7649#include "psi.fh" 7650 7651* **** local variables **** 7652 integer taskid,MASTER,taskid_j,pj 7653 parameter (MASTER=0) 7654 7655 logical oprint,value 7656 integer icount,iicount,version 7657 integer i,j,k,l,ij,kl 7658 integer nall,nnall,n4all,nx,ny,nz 7659 integer psiij 7660 integer ipackm(2),jpackm(2) 7661 integer h1_integrals(2),h2_integrals(2) 7662 integer vij(2),dnij(2),orbi(2),orbj(2) 7663 real*8 e1,e2,scal1,scal2,dv 7664 7665* **** external functions **** 7666 integer control_version 7667 external control_version 7668 real*8 lattice_omega 7669 external lattice_omega 7670 7671 call Parallel_taskid(taskid) 7672 call Parallel2d_taskid_j(taskid_j) 7673 oprint = (taskid.eq.MASTER) 7674 version = control_version() 7675 if (version.eq.4) then 7676 7677 if (oprint) then 7678 write(luout,*) 7679 > "== Generating One-Electron and Two-Electron Integrals ==" 7680 write(luout,*) 7681 end if 7682 call D3dB_nx(1,nx) 7683 call D3dB_ny(1,ny) 7684 call D3dB_nz(1,nz) 7685 scal1 = 1.0d0/dble(nx*ny*nz) 7686 scal2 = 1.0d0/lattice_omega() 7687 dv = scal1*lattice_omega() 7688 7689 nall = ne(1) + ne_excited(1) 7690 nnall = (nall*(nall+1))/2 7691 n4all = (nnall*(nnall+1))/2 7692 7693 value = BA_alloc_get(mt_dcpl,nfft3d,'vij',vij(2),vij(1)) 7694 value = value.and. 7695 > BA_alloc_get(mt_dcpl,nfft3d,'dnij',dnij(2),dnij(1)) 7696 value = value.and. 7697 > BA_alloc_get(mt_dcpl,nfft3d,'orbi',orbi(2),orbi(1)) 7698 value = value.and. 7699 > BA_alloc_get(mt_dcpl,nfft3d,'orbj',orbj(2),orbj(1)) 7700 value = value.and. 7701 > BA_alloc_get(mt_dbl,nnall,'h1_integrals', 7702 > h1_integrals(2),h1_integrals(1)) 7703 value = value.and. 7704 > BA_alloc_get(mt_dbl,n4all,'h2_integrals', 7705 > h2_integrals(2),h2_integrals(1)) 7706 value = value.and. 7707 > BA_alloc_get(mt_int,nnall,'ipackm', 7708 > ipackm(2),ipackm(1)) 7709 value = value.and. 7710 > BA_alloc_get(mt_int,nnall,'jpackm', 7711 > jpackm(2),jpackm(1)) 7712 if (.not.value) 7713 > call errquit('psi_2q_integrals:error allocating heap',0,MA_ERR) 7714 7715 call Parallel_shared_vector_zero(.true.,nnall, 7716 > dbl_mb(h1_integrals(1))) 7717 call Parallel_shared_vector_zero(.true.,n4all, 7718 > dbl_mb(h2_integrals(1))) 7719 7720 !**** generate 1e integrals - H1 only contains kinetic + ion-electron potentials **** 7721 call electron_gen_vall0() 7722 icount = 0 7723 do i=1,nall 7724 call psi_fetch_orbi_replicated(i,dcpl_mb(orbi(1))) 7725 call electron_get_H0psi_k_orb(dcpl_mb(orbi(1)), 7726 > dcpl_mb(vij(1))) 7727 7728 do j=i,nall 7729 call psi_fetch_orbj_indx(j,psiij,pj) 7730 if (pj.eq.taskid_j) then 7731 call Pack_cc_idot(1,dcpl_mb(psiij),dcpl_mb(vij(1)),e1) 7732 dbl_mb(h1_integrals(1)+icount) = e1 7733 end if 7734 int_mb(ipackm(1) + icount) = i 7735 int_mb(jpackm(1) + icount) = j 7736 icount = icount + 1 7737 end do 7738 call ga_sync() 7739 end do 7740 call Parallel_Vector_SumAll(nall,dbl_mb(h1_integrals(1))) 7741 7742 if (oprint) then 7743 write(luout,*) "begin_one_electron_integrals" 7744 icount = 0 7745 do i=1,nall 7746 do j=i,nall 7747 write(luout,'(I5,I5,F21.10)') 7748 > j,i,dbl_mb(h1_integrals(1)+icount) 7749 icount = icount + 1 7750 end do 7751 end do 7752 write(luout,*) "end_one_electron_integrals" 7753 end if 7754 7755 !*** generate 2e integrals *** 7756 iicount = 0 7757 do ij=1,nnall 7758 i = int_mb(ipackm(1)+ij-1) 7759 j = int_mb(jpackm(1)+ij-1) 7760 call psi_fetch_orbi_replicated(i,dcpl_mb(orbi(1))) 7761 call Pack_c_unpack(1, dcpl_mb(orbi(1))) 7762 call D3dB_cr_pfft3b(1,1,dcpl_mb(orbi(1))) 7763 call psi_fetch_orbi_replicated(j,dcpl_mb(orbj(1))) 7764 call Pack_c_unpack(1, dcpl_mb(orbj(1))) 7765 call D3dB_cr_pfft3b(1,1,dcpl_mb(orbj(1))) 7766 7767* **** generate dnij for Vij **** 7768 call D3dB_rr_Mul(1,dcpl_mb(orbi(1)), 7769 > dcpl_mb(orbj(1)), 7770 > dcpl_mb(dnij(1))) 7771 call D3dB_r_SMul1(1,scal2,dcpl_mb(dnij(1))) 7772 call coulomb2_v(dcpl_mb(dnij(1)),dcpl_mb(vij(1))) 7773 7774 7775 do kl=ij,nnall 7776 k = int_mb(ipackm(1)+kl-1) 7777 l = int_mb(jpackm(1)+kl-1) 7778 call psi_fetch_orbi_replicated(k,dcpl_mb(orbi(1))) 7779 call Pack_c_unpack(1, dcpl_mb(orbi(1))) 7780 call D3dB_cr_pfft3b(1,1,dcpl_mb(orbi(1))) 7781 call psi_fetch_orbi_replicated(l,dcpl_mb(orbj(1))) 7782 call Pack_c_unpack(1, dcpl_mb(orbj(1))) 7783 call D3dB_cr_pfft3b(1,1,dcpl_mb(orbj(1))) 7784 7785 call D3dB_rr_Mul(1,dcpl_mb(orbi(1)), 7786 > dcpl_mb(orbj(1)), 7787 > dcpl_mb(dnij(1))) 7788 call D3dB_r_SMul1(1,scal2,dcpl_mb(dnij(1))) 7789 7790 call D3dB_rr_dot(1,dcpl_mb(dnij(1)),dcpl_mb(vij(1)),e2) 7791 !e2 = 0.5d0*e2*dv 7792 e2 = e2*dv 7793 dbl_mb(h2_integrals(1)+iicount) = e2 7794 7795 !write(*,*) "i,j,k,l,e2=",i,j,k,l,e2 7796 iicount = iicount + 1 7797 end do 7798 end do 7799 7800 7801 if (oprint) then 7802 write(luout,*) 7803 write(luout,*) "begin_two_electron_integrals" 7804 iicount = 0 7805 do ij=1,nnall 7806 i = int_mb(ipackm(1)+ij-1) 7807 j = int_mb(jpackm(1)+ij-1) 7808 do kl=ij,nnall 7809 k = int_mb(ipackm(1)+kl-1) 7810 l = int_mb(jpackm(1)+kl-1) 7811 write(luout,'(I5,I5,I5,I5,F20.10)') 7812 > j,i,l,k,dbl_mb(h2_integrals(1)+iicount) 7813 iicount = iicount + 1 7814 end do 7815 end do 7816 write(luout,*) "end_two_electron_integrals" 7817 write(luout,*) 7818 end if 7819 7820 7821* **** deallocate memory from heap **** 7822 value = BA_free_heap(vij(2)) 7823 value = value.and.BA_free_heap(dnij(2)) 7824 value = value.and.BA_free_heap(orbi(2)) 7825 value = value.and.BA_free_heap(orbj(2)) 7826 value = value.and.BA_free_heap(h1_integrals(2)) 7827 value = value.and.BA_free_heap(h2_integrals(2)) 7828 value = value.and.BA_free_heap(ipackm(2)) 7829 value = value.and.BA_free_heap(jpackm(2)) 7830 if (.not. value) 7831 > call errquit('psi_2q_integrals:error freeing heap',1,MA_ERR) 7832 end if 7833 7834 return 7835 end 7836 7837 subroutine psi_fetch_orbi_replicated(i,orbi) 7838 implicit none 7839 integer i 7840 complex*16 orbi(*) 7841 7842#include "bafdecls.fh" 7843#include "errquit.fh" 7844#include "psi.fh" 7845 7846* **** local variables **** 7847 integer pi,q,taskid_j,psii_ptr 7848 7849 call Parallel2d_taskid_j(taskid_j) 7850 7851 if (i.le.ne(1)) then 7852 call Dneall_ntoqp(i,q,pi) 7853 call Pack_c_Zero(1,orbi) 7854 if (pi.eq.taskid_j) then 7855 psii_ptr = psi1(1) + (q-1)*npack1 7856 call Pack_c_Copy(1,dcpl_mb(psii_ptr),orbi) 7857 end if 7858 call D1dB_Vector_SumAll(2*npack1,orbi) 7859 else 7860 psii_ptr = psi1_excited(1) + (i-1-ne(1))*npack1 7861 call Pack_c_Copy(1,dcpl_mb(psii_ptr),orbi) 7862 end if 7863 7864 return 7865 end 7866 7867 subroutine psi_fetch_orbj_indx(j,orbj_indx,pj) 7868 implicit none 7869 integer j 7870 integer orbj_indx 7871 integer pj 7872 7873#include "bafdecls.fh" 7874#include "errquit.fh" 7875#include "psi.fh" 7876 7877* **** local variables **** 7878 integer q,np_j 7879 7880 call Parallel2d_np_j(np_j) 7881 if (j.le.ne(1)) then 7882 call Dneall_ntoqp(j,q,pj) 7883 orbj_indx = psi1(1) + (q-1)*npack1 7884 else 7885 orbj_indx = psi1_excited(1) + (j-1-ne(1))*npack1 7886 pj = mod(pj+1,np_j) 7887 end if 7888 7889 return 7890 end 7891 7892 7893