1c ... jochen: this comes from rohf_hessv2.F but has everything 2c "v2" replaced by "v3", for use with frequency-dependent response 3c and is called from the cphf_solve3 part of the cphf code 4c 5c ... jochen: further mods were made to accomodate the situation that 6c we might have damping in the response. That causes all quantities to 7c have an imaginary part, too 8 9 subroutine rohf_hessv3_ext( 10 & acc, 11 & g_x, 12 & g_ax, 13 & g_x_im, 14 & g_Ax_im, 15 & omega, 16 & limag, 17 & lifetime, 18 & gamwidth, 19 & ncomp) 20c 21c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 22c Date : 03-24-12 23c Note.- Equivalent to rohf_hessv3() (routine written by 24c J. Autschbach) but using optimizing routines. 25 26 implicit none 27#include "errquit.fh" 28#include "crohf.fh" 29#include "cscf.fh" 30#include "stdio.fh" 31#include "util.fh" 32#include "global.fh" 33c 34c $Id$ 35c 36c ... jochen: these two arrays now have two components: 37 integer g_x(2) ! [input] A-matrix elements for density matrix 38 integer g_ax(2) ! [output] Perturbed Fock operator 39c ... jochen: also, we might have imaginary components: 40 integer g_x_im(2) ! [input] A-matrix elements, Im 41 integer g_ax_im(2) ! [output] Perturbed Fock operator, Im 42 43 double precision acc, omega, gamwidth 44 logical limag, lifetime 45c 46 integer gtype,grow,gcol,growp,gcolp, ipm, ncomp 47 logical oprint, debug 48 external rohf_hessv_xx3_ext 49c 50c ================================================================ 51 52 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 53c debug = (.true. .and. ga_nodeid().eq.0) ! for code development 54 55 if (debug) write (6,*) 56 & 'rohf_hessv3: limag, omega, lifetime, gamwidth', 57 & limag, omega, lifetime, gamwidth 58c 59c Check for debug 60c 61 oprint= util_print('rohf_hessv2',print_debug) 62 if (crohf_init_flag.ne.1) 63 $ call errquit('rohf_hessv3: ROHF internal block invalid',0, 64 & UNKNOWN_ERR) 65c 66c ... jochen: use first component for the dimension checks. 67c the second component MUST have the same dimensions 68c otherwise there will be problems 69 call ga_inquire(g_x(1),gtype,grow,gcol) 70 if (grow.ne.crohf_vlen) 71 $ call errquit('rohf_hessv3: invalid vector length',0, 72 & UNKNOWN_ERR) 73 call ga_inquire(g_ax(1),gtype,growp,gcolp) 74 if (growp.ne.crohf_vlen) 75 $ call errquit('rohf_hessv3: invalid vector length',0, 76 & UNKNOWN_ERR) 77 if (gcol.ne.gcolp) 78 $ call errquit('rohf_hessv3: invalid no. of vectors',0, 79 & UNKNOWN_ERR) 80c 81c Call internal routine 82c 83 if (debug) write (6,*) 'calling rohf_hessv_xx3' 84 85 call rohf_hessv_xx3_ext( basis, geom, nbf, nmo, 86 $ nclosed, nopen, 87 $ pflg, 88 & g_movecs, 89 & oskel, noskew, 90 $ crohf_g_fcv, crohf_g_fpv, crohf_g_fcp, 91 $ acc, lshift, 92 & g_x, g_ax, 93 & g_x_im, g_Ax_im, 94 & omega, limag, 95 & lifetime, gamwidth, ncomp) 96 97 if (debug) write (6,*) 'back from rohf_hessv_xx3' 98c 99c Zap numbers much smaller than acc to ensure hard zeroes 100c remain unpolluted ... cannot use a threshold larger than the 101c integral accuracy since can break symmetry in non-abelian groups 102c Also must ensure that the threshold tends to zero to permit 103c tight convergence. 104c 105c ... jochen: screen components 106 do ipm = 1,ncomp 107 call ga_screen(g_ax(ipm), 108 & max(min(acc*acc,acc*0.01d0,1d-12),1d-16)) 109 if (lifetime) call ga_screen(g_ax_im(ipm), 110 & max(min(acc*acc,acc*0.01d0,1d-12),1d-16)) 111 enddo 112c 113 end 114 115 subroutine rohf_hessv_xx3_ext( 116 & basis, geom, 117 & nbf, nmo, nclosed, nopen, 118 $ pflg, 119 $ g_movecs, 120 & oskel, noskew, 121 & g_fcv, g_fpv, g_fcp, 122 $ acc, 123 & lshift, 124 & g_x, g_ax, 125 & g_x_im, g_Ax_im, 126 & omega, 127 & limag, 128 & lifetime, 129 & gamwidth, 130 & ncomp) 131c Note.- Improvements to this subroutine done by 132c Fredy W. Aquino, Northwestern University (Oct 2012) 133c Using rohf_hessv_2e2_opt() instead of old rohf_hessv_2e2() 134c rohf_hessv_2e3_opt() instead of old rohf_hessv_2e3() 135 136C $Id$ 137 implicit none 138#include "errquit.fh" 139#include "global.fh" 140#include "mafdecls.fh" 141#include "rtdb.fh" 142#include "bgj.fh" 143c 144 integer basis, geom 145 integer nbf, nmo, nclosed, nopen 146 integer pflg 147 integer g_movecs 148 logical oskel, noskew 149 integer g_fcv, g_fpv, g_fcp 150 double precision acc 151 double precision lshift 152c ... jochen: input arrays g_x and g_Ax have two components here 153 integer g_x(2), g_ax(2), vlen, nvec, g_tmp, gtype, ipm 154 integer g_x_im(2) ! [input] A-matrix elements, Im 155 integer g_ax_im(2) ! [output] Perturbed Fock operator, Im 156 157 double precision omega, gamwidth, 158 & wls, wlsim 159 logical limag, lifetime 160 integer ncomp 161 double precision omg(ncomp),gam(ncomp),coeffw 162 logical debug 163 external rohf_hessv_2e3_opt 164 165c 166c ================================================================= 167 168 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 169 170c debug = (.true. .and. ga_nodeid().eq.0) ! for code development 171c 172 if (debug) write (6,*) 'hessv3: omega =',omega 173 if (debug) write (6,*) 'hessv3: limag =',limag 174 if (debug) write (6,*) 175 & 'hessv3: lifetime, gamwidth =',lifetime, gamwidth 176c 177 do ipm = 1,ncomp 178 call ga_zero(g_Ax(ipm)) 179 if (lifetime) call ga_zero(g_Ax_im(ipm)) 180 end do 181 if (pflg.gt.2 .or. pflg.le.0) then 182 call errquit('rohf_hessv_xx: pflg invalid ', pflg, 183 & UNKNOWN_ERR) 184 endif 185c 186c ... jochen: to be consistent with the preconditioner, where 187c the level shift is added, we need to do the same thing here 188c and also add and subtract the frequency times 4 (it is times 189c 4 because of the factors of 4 in rohf_hessv_1e and in the 190c preconditioner) 191c During a response calculation, pflg is equal to 2 192c 193c what do we do here? Compare Gauss' paper Eqs. (32) and (135): 194c The lhs of the CPHF equations contain a term 195c (e_a - e_i -/+ omega) U_ai. First, we initialize g_Ax with 196c the term proportional to omega, then we add the delta-e term 197c (the e's are the orbital energies, calculated in hessv_1e as 198c the diagonal of the Fock matrix transformed to the MO basis) 199 200 if (pflg .gt. 0) then 201 coeffw=4.0d0 ! r-dft 202 omg(1)=-omega 203 omg(2)= omega 204 gam(1)=-gamwidth 205 gam(2)= gamwidth 206 if (.not.lifetime) then 207c no damping: initialize Ax with terms proportional omega 208 do ipm=1,ncomp 209 wls = lshift + coeffw * omg(ipm) 210 call ga_dadd(wls,g_x(ipm),0.d0,g_ax(ipm),g_ax(ipm)) 211 enddo ! end-loop-ipm 212 else ! lifetime 213c take care of damping here: Re and Im are coupled by gamwidth 214 do ipm=1,ncomp 215 wls = lshift + coeffw * omg(ipm) 216 wlsim = -coeffw * gam(ipm) 217 call ga_dadd(wls ,g_x(ipm), 218 & wlsim,g_x_im(ipm), 219 & g_ax(ipm)) 220 wls = coeffw * gam(ipm) 221 wlsim = lshift + coeffw * omg(ipm) 222 call ga_dadd(wls ,g_x(ipm), 223 & wlsim,g_x_im(ipm), 224 & g_ax_im(ipm)) 225 enddo ! end-lopp-ipm 226 endif ! .not.lifetime 227 call ga_sync() 228 if (debug) write (6,*) 'calling rohf_hessv_1e' 229 230c ============== debug g_ax ==================== START 231 if (debug) then 232 do ipm=1,ncomp 233 if (ga_nodeid().eq.0) 234 & write(*,*) '------- g_ax_re-0-c(',ipm,')------ START' 235 call ga_print(g_ax(ipm)) 236 if (ga_nodeid().eq.0) 237 & write(*,*) '------- g_ax_re-0-c(',ipm,')------ END' 238 enddo ! end-loop-ipm 239 if (lifetime) then 240 do ipm=1,2 241 if (ga_nodeid().eq.0) 242 & write(*,*) '------- g_ax_im-0-c(',ipm,')------ START' 243 call ga_print(g_ax_im(ipm)) 244 if (ga_nodeid().eq.0) 245 & write(*,*) '------- g_ax_im-0-c(',ipm,')------ END' 246 enddo ! end-loop-ipm 247 endif ! end-if-lifetime 248 endif ! end-if-debug 249c ============== debug g_ax ==================== END 250c 251c next: add (e_a - e_i) times A (also called U) matrix to Ax 252 do ipm=1,ncomp 253 call rohf_hessv_1e(basis,geom, ! in : handles 254 & nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell 255 $ g_fcv,g_fpv,g_fcp, ! in : densities 256 $ g_x(ipm), ! in : g_x 257 & g_ax(ipm)) ! out: 1e contrib to Ax 258 enddo ! end-loop-ipm 259 if (lifetime) then 260 do ipm=1,ncomp 261 call rohf_hessv_1e(basis,geom, ! in : handles 262 & nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell 263 $ g_fcv,g_fpv,g_fcp, ! in : densities 264 $ g_x_im(ipm), ! in : g_x_im 265 & g_ax_im(ipm)) ! out: 1e contrib to Ax_im 266 enddo ! end-loop-ipm 267 endif ! end-if-lifetime 268 endif ! pflg.gt.0 269c ============== debug g_ax ==================== START 270 if (debug) then 271 do ipm=1,ncomp 272 if (ga_nodeid().eq.0) 273 & write(*,*) '------- g_ax_re-0-d(',ipm,')------ START' 274 call ga_print(g_ax(ipm)) 275 if (ga_nodeid().eq.0) 276 & write(*,*) '------- g_ax_re-0-d(',ipm,')------ END' 277 enddo ! end-loop-ipm 278 if (lifetime) then 279 do ipm=1,2 280 if (ga_nodeid().eq.0) 281 & write(*,*) '------- g_ax_im-0-d(',ipm,')------ START' 282 call ga_print(g_ax_im(ipm)) 283 if (ga_nodeid().eq.0) 284 & write(*,*) '------- g_ax_im-0-d(',ipm,')------ END' 285 enddo ! end-loop-ipm 286 endif ! end-if-lifetime 287 endif ! end-if-debug 288c ============== debug g_ax ==================== END 289 if (pflg .gt. 1) then 290c 291c the next call basically uses the current guess for the solution 292c vector x (in g_x, which is the perturbed density matrix in the 293c MO basis) and calculates the perturbed Fock operator in the MO basis. 294c real and imaginary part of that Fock operator can be handled 295c separately here 296 297 if (ncomp.gt.1) then ! call 2e code for dynamic case 298 299 call rohf_hessv_2e3_opt( 300 & basis, ! in: basis handle 301 & geom, ! in: geom handle 302 & nbf, ! in: nr. basis functions 303 & nmo, ! in: nr. MOs 304 & nclosed, ! in: nr. double occupied MOs 305 & nopen, ! in: nr. single occupied MOs 306 $ g_movecs,! in: MO coefficients 307 & oskel, ! in: =.true. -> 308 & noskew, ! in: =.true. -> symmetric density matrix 309 & g_x, ! in: 310 & g_x_im, ! in: 311 & g_ax, ! in/out: Hessian product 312 & g_ax_im, ! in/out: Hessian product 313 & acc, ! in: accuracy Fock construction 314 & limag, ! in: =.true. -> imaginary component allowed 315 & lifetime)! in: =.true. -> RE-IM =.false -> RE 316 317 else ! call static 2e code 318 319 call rohf_hessv_2e2_opt( 320 & basis, ! in : basis handle 321 & geom, ! in : geom handle 322 & nbf, ! in : nr. basis functions 323 & nmo, ! in : nr. MOs vecs 324 & nclosed, ! in : nr. occupied MOs 325 & nopen, ! in : nr. open shells (unpaired e's) 326 $ g_movecs, ! in : MO vec coeffs 327 & oskel, ! in : 328 & noskew, ! in : symm density ? 329 & acc, ! in : accuracy Fock construction 330 & g_x(1), ! in : RE g_x 331 & g_x_im(1), ! in : IM g_x 332 & g_ax(1), ! in/ou : RE g_ax Hessian product 333 & g_ax_im(1),! in/ou : IM g_ax Hessian product 334 & lifetime) 335 336 endif ! ncomp 337 338 endif ! pflg.gt.1 339c 340 end 341 342 subroutine rohf_hessv3_cmplx( 343 & acc, ! in : accuracy 344 & g_z, ! in : z 345 & g_Az1, ! in : Az1 346 & nsub, ! in 347 & omega, ! in : 348 & limag, ! in : 349 & lifetime, ! in : 350 & gamwidth, ! in : 351 & ncomp, ! in : nr. components (+/-) 352 & iter_cphf) ! in: iteration nr. in cphf 353c 354c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 355c Date : 03-24-12 356c Note.- Modifying rohf_hessv3 to handle complex vars 357 358 implicit none 359#include "errquit.fh" 360#include "crohf.fh" 361#include "cscf.fh" 362#include "stdio.fh" 363#include "util.fh" 364#include "mafdecls.fh" 365#include "global.fh" 366c 367c $Id$ 368 integer ncomp 369 integer g_z(ncomp), ! [input] A-matrix elements for density matrix 370 & g_Az(ncomp) ! [output] Perturbed Fock operator 371 integer g_Az1, 372 & nsub 373 double precision acc, omega, gamwidth 374 logical limag, lifetime 375 integer gtype,grow,gcol, 376 & growp,gcolp,ipm,iter_cphf 377 logical oprint, debug 378 external rohf_hessv_xx3_cmplx 379 380 381 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 382 if (debug) write (6,*) 383 & 'rohf_hessv3: limag, omega, lifetime, gamwidth', 384 & limag, omega, lifetime, gamwidth 385 386c 387c Check for debug 388c 389 oprint= util_print('rohf_hessv2',print_debug) 390 if (crohf_init_flag.ne.1) 391 $ call errquit('rohf_hessv3: ROHF internal block invalid',0, 392 & UNKNOWN_ERR) 393 394c 395c ... jochen: use first component for the dimension checks. 396c the second component MUST have the same dimensions 397c otherwise there will be problems 398 call ga_inquire(g_z(1),gtype,grow,gcol) 399 400 if (grow.ne.crohf_vlen) 401 $ call errquit('rohf_hessv3_cmplx: invalid vector length-1', 402 & 0,UNKNOWN_ERR) 403 404 goto 177 ! skip-for-the moment 405 if (growp.ne.crohf_vlen) 406 $ call errquit('rohf_hessv3_cmplx: invalid vector length-2', 407 & 0,UNKNOWN_ERR) 408 if (gcol.ne.gcolp) 409 $ call errquit('rohf_hessv3_cmplx: invalid no. of vectors', 410 & 0,UNKNOWN_ERR) 411 177 continue 412c 413c Call internal routine 414c 415 if (debug) write (6,*) 'calling rohf_hessv_xx3_cmplx' 416 417 call rohf_hessv_xx3_cmplx( 418 & g_z, ! in : 419 & g_Az1, 420 & nsub, 421 & ncomp, ! in : nr. components 422 & basis, ! in : basis handle 423 & geom, ! in : geom hande 424 & nbf, ! in : nr. basis functions 425 & nmo, ! in : nr. MOs 426 $ nclosed, ! in : nr. occupied MOs 427 & nopen, ! in : nr. 1/2 occupied MOs (=0 for closed shells) 428 $ pflg, 429 & g_movecs, ! in : MO coefficients 430 & oskel, 431 & noskew, 432 $ crohf_g_fcv, ! in : closed-virtual density 433 & crohf_g_fpv, ! in : open-virtual density (not used) 434 & crohf_g_fcp, ! in : closed-open density (not used) 435 $ acc, 436 & lshift, 437 & omega, 438 & limag, 439 & lifetime, 440 & gamwidth, 441 & iter_cphf) ! in : iteration nr. in cphf cycle 442 443 if (debug) write (6,*) 'back from rohf_hessv_xx3_cmplx' 444 445c 446c Zap numbers much smaller than acc to ensure hard zeroes 447c remain unpolluted ... cannot use a threshold larger than the 448c integral accuracy since can break symmetry in non-abelian groups 449c Also must ensure that the threshold tends to zero to permit 450c tight convergence. 451c 452c ... jochen: screen components 453c ++++++++++++++++++++++++++++++++++++++++++ 454c =====> WARNING-UNFINISHED ROUTINE <======= 455c ++++++++++++++++++++++++++++++++++++++++++ 456c I need to adapt it for complex (later) --> zapping routine 457c do ipm = 1,ncomp 458c call ga_screen(g_ax(ipm), 459c & max(min(acc*acc,acc*0.01d0,1d-12),1d-16)) 460c if (lifetime) call ga_screen(g_ax_im(ipm), 461c & max(min(acc*acc,acc*0.01d0,1d-12),1d-16)) 462c enddo 463 464 end 465 466 subroutine rohf_hessv_xx3_cmplx( 467 & g_z, ! in : 468 & g_Az1, ! ou : product: A x z 469 & nsub, ! in : nr. subspace 470 & ncomp, ! in : nr. components 471 & basis, ! in : basis handle 472 & geom, ! in : geom hande 473 & nbf, ! in : nr. basis functions 474 & nmo, ! in : nr. MOs 475 & nclosed, ! in : nr. occupied MOs 476 & nopen, ! in : nr. 1/2 occupied MOs (=0 for closed shells) 477 $ pflg, 478 $ g_movecs,! in : MO coefficients 479 & oskel, 480 & noskew, 481 & g_fcv, ! in : closed-virtual density 482 & g_fpv, ! in : open-virtual density (not used) 483 & g_fcp, ! in : closed-open density (not used) 484 $ acc, 485 & lshift, 486 & omega, 487 & limag, 488 & lifetime, 489 & gamwidth, 490 & iter_cphf) 491c 492c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 493c Date : 03-24-12 494c -> modified from rohf_hessv_xx3() 495c $Id$ 496 497 implicit none 498#include "errquit.fh" 499#include "global.fh" 500#include "mafdecls.fh" 501#include "rtdb.fh" 502#include "bgj.fh" 503 integer ipm,ncomp,iter_cphf 504 integer ncompmx 505 parameter(ncompmx=2) 506 integer g_z(ncompmx), ! [input] A-matrix elements 507 & g_Az(ncompmx) ! [output] Perturbed Fock operator 508 integer g_Az1,nsub, 509 & alo(3),ahi(3), 510 & n1,maxsub, 511 & m1,m2,p1,p2, 512 & pretty 513 integer g_movecs 514 integer g_fcv, 515 & g_fpv, 516 & g_fcp 517 integer g_xreim,g_Axreim ! scratch GA arrays 518 integer basis,geom 519 integer nbf,nmo,nvir, 520 & nclosed,nopen,nreim,n 521 integer pflg 522 double precision acc,lshift 523 integer vlen, nvec,g_tmp,gtype 524 double precision omega,gamwidth, 525 & omg(ncompmx), 526 & gam(ncompmx), 527 & wls,wlsim 528 double complex wls_cmplx 529 logical oskel, noskew 530 logical limag,lifetime 531 logical debug 532c DIM/QM JEM 533 integer g_x(ncomp), g_x_im(ncomp) 534 integer g_Ax(ncomp), g_Ax_im(ncomp) 535 logical ldimqm 536c 537 external rohf_hessv_2e3_opt_cmplx, 538 & rohf_hessv_2e2_opt_cmplx, 539 & getreorim,getreorim1, 540 & conv2complex3,conv2complex4 541c 542c ================================================================= 543 544 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 545 546 if (debug) write (6,*) 'hessv3: omega =',omega 547 if (debug) write (6,*) 'hessv3: limag =',limag 548 if (debug) write (6,*) 549 & 'hessv3: lifetime, gamwidth =',lifetime, gamwidth 550 551c ----- Clear sub-block ---- START 552 call ga_inquire(g_Az1,gtype,n1,maxsub) 553 call ga_inquire(g_z(1),gtype,n,nvec) ! get (n1,nvec) 554c 555 alo(1)=1 556 ahi(1)=n1 557 alo(2)=nsub+1 558 ahi(2)=nsub+nvec 559 call nga_zero_patch(g_Az1,alo,ahi) 560c ----- Clear sub-block ---- END 561c 562 if (pflg.gt.2 .or. pflg.le.0) then 563 call errquit('rohf_hessv_xx: pflg invalid ', pflg, 564 & UNKNOWN_ERR) 565 endif 566c 567c ... jochen: to be consistent with the preconditioner, where 568c the level shift is added, we need to do the same thing here 569c and also add and subtract the frequency times 4 (it is times 570c 4 because of the factors of 4 in rohf_hessv_1e and in the 571c preconditioner) 572c During a response calculation, pflg is equal to 2 573c 574c what do we do here? Compare Gauss' paper Eqs. (32) and (135): 575c The lhs of the CPHF equations contain a term 576c (e_a - e_i -/+ omega) U_ai. First, we initialize g_Ax with 577c the term proportional to omega, then we add the delta-e term 578c (the e's are the orbital energies, calculated in hessv_1e as 579c the diagonal of the Fock matrix transformed to the MO basis) 580 if (pflg .gt. 0) then 581 omg(1)=-omega 582 omg(2)= omega 583 gam(1)=-gamwidth 584 gam(2)= gamwidth 585c take care of damping here: Re and Im are coupled by gamwidth 586 m1=1 587 m2=n 588 p1=nsub+1 589 p2=nsub+nvec 590 do ipm=1,ncomp 591 wls = lshift + 4.0d0 * omg(ipm) 592 wlsim = -4d0 * gam(ipm) 593 wls_cmplx=dcmplx(wls,-wlsim) 594 call ga_copy_patch('n',g_z(ipm),1 ,n ,1 ,nvec, 595 & g_Az1 ,m1,m2,p1,p2) 596 call ga_scale_patch(g_Az1,m1,m2,p1,p2,wls_cmplx) 597 m1=m1+n 598 m2=m2+n 599 enddo ! end-lopp-ipm 600 601 if (debug) then 602 m1=1 603 m2=n 604 p1=nsub+1 605 p2=nsub+nvec 606 do ipm=1,ncomp 607 if (ga_nodeid().eq.0) 608 & write(*,*) '------- g_Az1-c(',ipm,')------ START' 609 call ga_print(g_Az1) 610 if (ga_nodeid().eq.0) 611 & write(*,*) '------- g_Az1-c(',ipm,')------ END' 612 m1=m1+n 613 m2=m2+n 614 enddo ! end-loop-ipm 615 endif ! end-if-debug 616c 617c next: add (e_a - e_i) times A (also called U) matrix to Ax 618 call ga_inquire(g_z(1),gtype,n,nvec) ! get (n,nvec) 619 if (.not. ga_create(MT_DBL,n,nvec, 620 & 'hessv_xx3_cmplx: g_xreim',0,0,g_xreim)) 621 $ call errquit('hessv_xx3_cmplx: failed allocating g_xreim', 622 & n,GA_ERR) 623 if (.not. ga_create(MT_DBL,n,nvec, 624 & 'hessv_xx3_cmplx: g_xreim',0,0,g_Axreim)) 625 $ call errquit('hessv_xx3_cmplx: failed allocating g_xreim', 626 & n,GA_ERR) 627 628 nvir = nmo - nclosed - nopen 629 do nreim=1,2 ! loop in RE,IM 630 do ipm=1,ncomp 631 call getreorim(g_xreim, ! out : real or im arr 632 & g_z(ipm), ! in : = complx(g_xre,g_xim) 633 & nvir, ! in : nr. virtual MOs 634 & nclosed, ! in : nr. occupied MOs 635 & nreim) ! in : =1 -> re =2 -> im 636 call getreorim1(g_Axreim,! out : real or im arr 637 & g_Az1, ! in : = complx(g_xre,g_xim) 638 & nsub, ! in : subblock index 639 & ipm, ! in : = 1,2 to access slctd component 640 & nvir, ! in : nr. virtual MOs 641 & nclosed, ! in : nr. occupied MOs 642 & nreim) ! in : =1 -> re =2 -> im 643 call rohf_hessv_1e( 644 & basis,geom, ! in : handles 645 & nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell 646 $ g_fcv,g_fpv,g_fcp, ! in : densities 647 & g_xreim, 648 & g_Axreim) 649c ++++++++ update g_Az +++++++++ START 650 call conv2complex4(g_Az1, ! out: = history matrix complex 651 & g_Axreim,! in : real arr 652 & nsub, ! in : subblock index 653 & ipm, ! in : = 1,2 to access slctd component 654 & nvir, ! in : nr. virtual MOs 655 & nclosed, ! in : nr. occupied MOs 656 & nreim) ! in : =1 -> re =2 -> im 657c ++++++++ update g_Az +++++++++ END 658 enddo ! end-loop-ipm 659 enddo ! end-loop-nreim 660 if (.not. ga_destroy(g_xreim)) call errquit 661 & ('hessv_xx3_cmplx: g_xreim',0, GA_ERR) 662 if (.not. ga_destroy(g_Axreim)) call errquit 663 & ('hessv_xx3_cmplx: g_xreim',0, GA_ERR) 664 endif ! pflg.gt.0 665c ============== debug g_ax ==================== START 666 if (debug) then 667 if (ga_nodeid().eq.0) 668 & write(*,*) '------- g_Az1-d------ START' 669 call ga_print(g_Az1) 670 if (ga_nodeid().eq.0) 671 & write(*,*) '------- g_Az1-d------ END' 672 endif ! end-if-debug 673c ============== debug g_ax ==================== END 674 if (pflg .gt. 1) then 675c 676c the next call basically uses the current guess for the solution 677c vector x (in g_x, which is the perturbed density matrix in the 678c MO basis) and calculates the perturbed Fock operator in the MO basis. 679c real and imaginary part of that Fock operator can be handled 680c separately here 681 682 if (ncomp.gt.1) then ! call 2e code for dynamic case 683 call rohf_hessv_2e3_opt_cmplx( 684 & g_z, ! in : 685 & g_Az1, ! out: history of g_Az 686 & nsub, ! in : 687 & basis, ! in : basis handle 688 & geom, ! in : geom handle 689 & nbf, ! in : nr. basis functions 690 & nmo, ! in : nr. MOs 691 & nclosed, ! in : nr. double occupied MOs 692 & nopen, ! in : nr. single occupied MOs 693 $ g_movecs, ! in : MO coefficients 694 & oskel, ! in : =.true. -> 695 & noskew, ! in : =.true. -> symmetric density matrix 696 & acc, ! in : accuracy Fock construction 697 & limag, ! in : =.true. -> imaginary component allowed 698 & lifetime) ! in : =.true. -> RE-IM =.false -> RE 699 else ! call static 2e code 700 call rohf_hessv_2e2_opt_cmplx( 701 & g_z(1), ! in : 702 & g_Az1, ! out: history of g_Az 703 & nsub, ! in : 704 & basis, ! in : basis handle 705 & geom, ! in : geom handle 706 & nbf, ! in : nr. basis functions 707 & nmo, ! in : nr. MOs vecs 708 & nclosed, ! in : nr. occupied MOs 709 & nopen, ! in : nr. open shells (unpaired e's) 710 $ g_movecs,! in : MO vec coeffs 711 & oskel, ! in : 712 & noskew, ! in : symm density ? 713 & acc, ! in : accuracy Fock construction 714 & lifetime) 715 endif ! ncomp 716 717 endif ! pflg.gt.1 718c 719c DIM/QM JEM 720 if (.not.rtdb_get(bgj_get_rtdb_handle(), 'dimqm:lrsp', MT_LOG, 721 $ 1, ldimqm)) ldimqm = .false. 722 if(ldimqm) then 723c Transforming complex arrays to 2 real arrays to fit DIM structure 724c This could be cleaned up with a substantial rework of dimqm_addop 725 do ipm = 1, ncomp 726 if (.not. ga_create(MT_DBL,n,nvec,'dimqm: g_x', 727 $ 0,0,g_x(ipm))) 728 $ call errquit('dimqm: failed allocating g_x',n,GA_ERR) 729 if (.not. ga_create(MT_DBL,n,nvec,'dimqm: g_x', 730 $ 0,0,g_x_im(ipm))) 731 $ call errquit('dimqm: failed allocating g_x_im',n,GA_ERR) 732 if (.not. ga_create(MT_DBL,n,nvec,'dimqm: g_Ax', 733 $ 0,0,g_Ax(ipm))) 734 $ call errquit('dimqm: failed allocating g_Ax',n,GA_ERR) 735 if (.not. ga_create(MT_DBL,n,nvec,'dimqm: g_Ax_im', 736 $ 0,0,g_Ax_im(ipm))) 737 $ call errquit('dimqm: failed allocating g_Ax_im',n,GA_ERR) 738 call getreorim(g_x(ipm), ! out : real or im arr 739 & g_z(ipm), ! in : = complx(g_xre,g_xim) 740 & nvir, ! in : nr. virtual MOs 741 & nclosed, ! in : nr. occupied MOs 742 & 1) ! in : =1 -> re =2 -> im 743 call getreorim(g_x_im(ipm), ! out : real or im arr 744 & g_z(ipm), ! in : = complx(g_xre,g_xim) 745 & nvir, ! in : nr. virtual MOs 746 & nclosed, ! in : nr. occupied MOs 747 & 2) ! in : =1 -> re =2 -> im 748 call getreorim1(g_Ax(ipm), ! out : real or im arr 749 & g_Az1, ! in : = complx(g_xre,g_xim) 750 & nsub, ! in : subblock index 751 & ipm, ! in : = 1,2 to access slctd component 752 & nvir, ! in : nr. virtual MOs 753 & nclosed, ! in : nr. occupied MOs 754 & 1) ! in : =1 -> re =2 -> im 755 call getreorim1(g_Ax_im(ipm), ! out : real or im arr 756 & g_Az1, ! in : = complx(g_xre,g_xim) 757 & nsub, ! in : subblock index 758 & ipm, ! in : = 1,2 to access slctd component 759 & nvir, ! in : nr. virtual MOs 760 & nclosed, ! in : nr. occupied MOs 761 & 2) ! in : =1 -> re =2 -> im 762 end do 763 call dimqm_addop(g_x, g_x_im, g_Ax, g_Ax_im, 764 & ncomp, limag, lifetime) 765 do ipm = 1,ncomp 766 call conv2complex4(g_Az1, ! out: = history matrix complex 767 & g_Ax(ipm), ! in : real arr 768 & nsub, ! in : subblock index 769 & ipm, ! in : = 1,2 to access slctd component 770 & nvir, ! in : nr. virtual MOs 771 & nclosed, ! in : nr. occupied MOs 772 & 1) ! in : =1 -> re =2 -> im 773 call conv2complex4(g_Az1, ! out: = history matrix complex 774 & g_Ax_im(ipm), ! in : real arr 775 & nsub, ! in : subblock index 776 & ipm, ! in : = 1,2 to access slctd component 777 & nvir, ! in : nr. virtual MOs 778 & nclosed, ! in : nr. occupied MOs 779 & 2) ! in : =1 -> re =2 -> im 780 if (.not. ga_destroy(g_x(ipm))) call errquit 781 & ('dimqm: failed destroying g_x',ipm, GA_ERR) 782 if (.not. ga_destroy(g_x_im(ipm))) call errquit 783 & ('dimqm: failed destroying g_x_im',ipm, GA_ERR) 784 if (.not. ga_destroy(g_Ax(ipm))) call errquit 785 & ('dimqm: failed destroying g_Ax',ipm, GA_ERR) 786 if (.not. ga_destroy(g_Ax_im(ipm))) call errquit 787 & ('dimqm: failed destroying g_Ax_im',ipm, GA_ERR) 788 789 end do 790 791 end if 792c 793 end 794 795 subroutine rohf_hessv_2e3_opt_cmplx( 796 & g_z, 797 & g_Az1, ! in: (n1,maxsub) history of Az matrix (large matrix) 798 & nsub, ! in: point to (n1,nvec) block to be updated in g_Az1 799 & basis, ! in: basis handle 800 & geom, ! in: geom handle 801 & nbf, ! in: nr. basis functions 802 & nmo, ! in: nr. MOs 803 & nclosed, ! in: nr. double occupied MOs 804 & nopen, ! in: nr. single occupied MOs 805 $ g_movec, ! in: MO coefficients 806 & oskel, ! in: =.true. -> 807 & noskew, ! in: =.true. -> symmetric density matrix 808 & acc, ! in: accuracy Fock construction 809 & limag, ! in: =.true. -> imaginary component allowed 810 & lifetime)! in: =.true. -> RE-IM =.false -> RE 811c 812c Purpose: Optimization of rohf_hessv_2e3() 813c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 814c Date : 03-28-12 815c Note1.- Modifying rohf_hessv_2e3, reducing computation of 816c two-electron integrals by putting together 817c symmetric+antisymmetric perturbed densities: g_dens(ipm) ipm=1,2 818c and doing a single call to shell_fock_build() when using 819c the routine shell_fock_build2(). 820c Note2.- size(g_Az1)=(n1,maxsub) n1=ncomp*n maxsumb=maxiter*nvec 821c ncomp=2 (+/-) n=nvir*nocc maxiter=10 (usually) nvec=3 (x,y,z) 822c nsub, will point to next (n1,nvec) block to be updated 823c --> The purpose of using g_Az1 is to reduce usage of memory. 824c by doing it I skip using g_Az(ipm) ipm=2 (n1,nvec) complex matrix. 825 826 implicit none 827#include "errquit.fh" 828#include "mafdecls.fh" 829#include "global.fh" 830#include "util.fh" 831#include "cscfps.fh" 832#include "rtdb.fh" 833#include "bgj.fh" 834#include "stdio.fh" 835#include "case.fh" 836c 837c Return the ROHF orbital 2e-Hessian vector product, g_ax = A * g_x 838c 839c ... jochen: modified version of rohf_hessv_2e2 which keeps track 840c of two sets of input vectors that couple via the density matrix. 841c one could likely save some memory here by re-using temp arrays 842c ... jochen: Also made modifications to calculate imaginary terms due 843c to finite lifetime damping 844ccccccccccccccc This code does NOT work for open shell!!!!!ccccccccccccccccc 845 integer g_z(2),g_Az(2) 846 integer g_Az1,nsub 847 integer m1,m2,p1,p2,pretty 848 integer basis, geom ! basis & geom handle 849 integer nclosed,nvir,nopen! Basis size and occupation 850 integer nbf,nmo ! No. of linearly dependent MOs 851 integer g_movec ! MO coefficients 852 logical oskel 853 854 double precision acc ! Accuracy of "Fock" construction 855 logical limag ! imaginary perturbation? 856 integer voff,xoff 857 integer nmul,gtype, 858 & ivec,vlen,nvec 859 integer g_tmp,g_tmp1, 860 & g_dens(4),g_fock(4), 861 & g_xreim(2) 862 double precision tol2e 863 logical odebug,oprint,noskew 864 integer dims(3),chunk(3), 865 & alo(3),ahi(3), 866 & blo(2),bhi(2) 867 integer ga_create_atom_blocked 868 external ga_create_atom_blocked, 869 & get_undosymm_fock,update_ax_fock, 870 & shell_fock_build2,getreorim, 871 & update_gz_reorim, 872 & update_gz_reorim1 873 double precision one,zero,mone,four,half,mhalf,two,mtwo 874 double precision itol_floor, itol_ceil 875 parameter(itol_floor=1.d-15, itol_ceil=1.d-3) 876 parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0) 877 parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0) 878 integer ipm ! counter for density matrix components 879 character*(255) cstemp 880 integer g_pmats(2),g_pmata(2),g_h1mat(2), 881 & cnt,ind,indx(2,2), 882 & npol,nset,nblock,ncomp 883 double precision tenm6,coef(2,2) 884 parameter (tenm6 = 1d-6) 885 logical lifetime,debug 886 data npol /1/ ! for restricted calculations 887 data indx /1,2, ! indx(1,1),indx(1,2) 888 & 3,4/ ! indx(2,1),indx(2,2) 889 ncomp=2 ! using two components 890 call ga_inquire(g_z(1),gtype,vlen,nvec) ! out: nvec,vlen 891 do ipm = 1,ncomp 892 if (.not. ga_create(MT_DBL,vlen,nvec, 893 & 'hessv_2e3_opt_cmplx: g_xreim',0,0,g_xreim(ipm))) 894 $ call errquit('rhessv_2e3_opt_cmplx: failed alloc g_xreim', 895 & nvec,GA_ERR) 896 enddo ! end-loop-ipm 897 nmul=1 898 if (npol.eq. 2) nmul=2 899 nset =1 ! for RE 900 nblock=2 ! for RE 901 if (lifetime) then 902 nset =2 ! for RE-IM 903 nblock=4 ! for RE-IM 904 endif 905 906 if (oskel) 907 $ call errquit('rohf_h2e3: no way',0, UNKNOWN_ERR) 908 909 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 910 911 if (debug) then 912 if (ga_nodeid().eq.0) then 913 write(*,2001) limag,lifetime,nset,nblock 914 2001 format('(limag,lifetime,nset,nblock)=(', 915 & L1,',',L1,',',i3,',',i3,')') 916 endif 917 endif ! end-if-debug 918 919 oprint= util_print('rohf_hessv2',print_debug) 920 921 if (nopen.ne.0) call errquit 922 $ ('rohf_h2e3: does not work for open shells',nopen, 923 & UNKNOWN_ERR) 924 odebug = util_print('rohf_hessv', print_debug) 925 tol2e = min(max(acc,itol_floor),itol_ceil) 926 nvir = nmo - nclosed - nopen 927 voff = nclosed + nopen + 1 928 929 dims(1) = nbf 930 dims(2) = nbf 931 chunk(1) = dims(1) 932 chunk(2) = -1 933 934c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---START 935 do ipm = 1,ncomp 936 write(cstemp,'(a,i1)') 'pmats_',ipm 937 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 938 & g_pmats(ipm))) call 939 & errquit('rohf_h2e3: nga_create failed '//cstemp(1:7), 940 & 0,GA_ERR) 941 call ga_zero(g_pmats(ipm)) 942 write(cstemp,'(a,i1)') 'pmata_',ipm 943 write(cstemp,'(a,i1)') 'h1mat_',ipm 944 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 945 & g_h1mat(ipm))) call 946 & errquit('rohf_h2e3: nga_create failed '//cstemp(1:7), 947 & 0,GA_ERR) 948 call ga_zero(g_h1mat(ipm)) 949 enddo ! end-loop-ipm 950c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---END 951 dims(1) = nvec 952 dims(2) = nbf 953 dims(3) = nbf 954 chunk(1) = dims(1) 955 chunk(2) = -1 956 chunk(3) = -1 957 call ga_sync() 958 do ipm = 1,nblock ! =2 or 4 959c ... allocate g_dens=[g_dens_re g_dens_im] 960 if (.not. nga_create (MT_DBL,3,dims,'CPKS dens',chunk, 961 & g_dens(ipm))) 962 & call errquit('rohf_h2e3: could not allocate g_dens',555, 963 & GA_ERR) 964 call ga_zero(g_dens(ipm)) 965c ... allocate g_fock=[g_fock_re g_fock_im] 966 if (.not. nga_create (MT_DBL,3,dims,'Fockv',chunk, 967 & g_fock(ipm))) 968 & call errquit('rohf_h2e3: could not allocate g_fock',555, 969 & GA_ERR) 970 call ga_zero(g_fock(ipm)) 971 enddo ! end-loop-ipm 972 if (debug) then 973 if (ga_nodeid().eq.0) 974 & write(*,*) 'BEF get_dens_reorim-RE' 975 endif ! end-if-debug 976c ---- Copy g_z --> g_x_reim ------ START 977 do ipm=1,ncomp 978 call ga_zero(g_xreim(ipm)) 979 call getreorim(g_xreim(ipm),! out : real or im arr 980 & g_z(ipm), ! in : = complx(g_xre,g_xim) 981 & nvir, ! in : nr. virtual MOs 982 & nclosed, ! in : nr. occupied MOs 983 & 1) ! in : =1 -> re =2 -> im 984 enddo ! end-loop-ipm 985c ---- Copy g_z --> g_x_reim ------ END 986 call get_dens_reorim( 987 & g_dens, ! in/ou: perturbed density matrix 988 & 1, ! in : =1 1st block RE 989 & g_xreim,! in : RE 990 & g_movec,! in : MO coefficients 991 & nbf, ! in : nr. basis functions 992 & nmo, ! in : nr. MOs 993 & 1, ! in : for ipol=1 -> istart=1 994 & nclosed,! in : nr. occupied MOs 995 & nvir, ! in : nr. virtual MOs 996 & nvec, ! in : nr. directions (x,y,z) 997 & npol, ! in : nr. polarizations 998 & limag, ! in : = .true. imaginary allowed 999 & g_pmats,! in : scratch GA array 1000 & g_pmata,! in : scratch GA array - NOT USED 1001 & g_h1mat)! in : scratch GA array 1002 1003 if (lifetime) then 1004c ---- Copy g_z --> g_x_reim ------ START 1005 do ipm=1,ncomp 1006 call ga_zero(g_xreim(ipm)) 1007 call getreorim(g_xreim(ipm),! out : real or im arr 1008 & g_z(ipm), ! in : = complx(g_xre,g_xim) 1009 & nvir, ! in : nr. virtual MOs 1010 & nclosed, ! in : nr. occupied MOs 1011 & 2) ! in : =1 -> re =2 -> im 1012 enddo ! end-loop-ipm 1013c ---- Copy g_z --> g_x_reim ------ END 1014 1015 call get_dens_reorim( 1016 & g_dens, ! in/ou: perturbed density matrix 1017 & 2, ! in : =2 2nd block IM 1018 & g_xreim,! in : IM 1019 & g_movec,! in : MO coefficients 1020 & nbf, ! in : nr. basis functions 1021 & nmo, ! in : nr. MOs 1022 & 1, ! in : for ipol=1 -> istart=1 1023 & nclosed,! in : nr. occupied MOs 1024 & nvir, ! in : nr. virtual MOs 1025 & nvec, ! in : nr. directions (x,y,z) 1026 & npol, ! in : nr. polarizations 1027 & limag, ! in : = .true. imaginary allowed 1028 & g_pmats,! in : scratch GA array 1029 & g_pmata,! in : scratch GA array - NOT USED 1030 & g_h1mat)! in : scratch GA array 1031 1032 endif ! end-if-lifetime 1033 do ipm = 1,ncomp 1034 if (.not.ga_destroy(g_xreim(ipm))) call errquit( 1035 & 'rohf_hessv3: ga_destroy failed g_xreim',0,GA_ERR) 1036 if (.not.ga_destroy(g_h1mat(ipm))) call errquit( 1037 & 'rohf_hessv3: ga_destroy failed g_h1mat',0,GA_ERR) 1038 enddo ! end-loop-ipm 1039 1040 if (debug) then 1041 do ipm=1,nblock 1042 if (ga_nodeid().eq.0) 1043 & write(*,*) '--------g_dens(',ipm,')-------START' 1044 call ga_print(g_dens(ipm)) 1045 if (ga_nodeid().eq.0) 1046 & write(*,*) '--------g_dens(',ipm,')-------END' 1047 enddo ! end-loop-ipm 1048 endif ! end-if-debug 1049 1050 call shell_fock_build2(g_fock, ! out: Fock matrices 1051 & g_dens, ! in : density matrices 1052 & geom, ! in : geom handle 1053 & basis, ! in : basis handle 1054 & nbf, ! in : nr. basis functions 1055 & nvec, ! in : nr. vecs (x,y,z) 1056 & npol, ! in : npol=1 for R-DFT =2 for U-DFT 1057 & ncomp, ! in : nr. components 1058 & nblock, ! in : nr. of g_dens,g_fock blocks 1059 & .true., ! in : =.true. for symm dens 1060 & tol2e, ! in : 1061 & debug) ! in : = .true. -> debugging printouts 1062 1063 if (debug) then 1064 do ipm=1,nblock 1065 if (ga_nodeid().eq.0) 1066 & write(*,*) '------- g_fock-0(',ipm,')------ START' 1067 call ga_print(g_fock(ipm)) 1068 if (ga_nodeid().eq.0) 1069 & write(*,*) '------- g_fock-0(',ipm,')------ END' 1070 enddo ! end-loop-ipm 1071 endif ! end-if-debug 1072 1073 call get_undosymm_fock( 1074 & g_fock, ! in/ou: fock matrix 1075 & nset, ! in : =1 g_x is real, =2 g_x is complex (g_x_re,g_x_im) 1076 & nvec, ! in : nr. directions (x,y,z) 1077 & nbf, ! in : nr. basis functions 1078 & npol, ! in : nr. polarizations 1079 & nmul, ! in : =1 npol=1 =2 npol=2 (acc. JK terms) 1080 & g_pmats, ! in : scratch GA array 1081 & limag) ! in : =.true. imaginary comp. exists 1082 1083c ------- Remove GA arrays: 1084 do ipm = 1,ncomp 1085 if (.not.ga_destroy(g_pmats(ipm))) call errquit( 1086 & 'rohf_hessv3: ga_destroy failed g_pmats',0,GA_ERR) 1087 enddo ! end-loop-ipm 1088 1089 if (debug) then 1090 do ipm=1,nblock 1091 if (ga_nodeid().eq.0) 1092 & write(*,*) '------- g_fock-unsym(',ipm,')------ START' 1093 call ga_print(g_fock(ipm)) 1094 if (ga_nodeid().eq.0) 1095 & write(*,*) '------- g_fock-unsym(',ipm,')------ END' 1096 enddo ! end-loop-ipm 1097 endif ! end-if-debug 1098 1099 if (debug) then 1100 if (ga_nodeid().eq.0) 1101 & write(*,*) '------- g_Az1-0------ START' 1102 call ga_print(g_Az1) 1103 if (ga_nodeid().eq.0) 1104 & write(*,*) '------- g_Az1-0------ END' 1105 endif ! end-if-debug 1106c 1107c start loop over components of perturbing field 1108c 1109 g_tmp = ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp') 1110 g_tmp1= ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp1') 1111 alo(2) = 1 1112 ahi(2) = nbf 1113 alo(3) = 1 1114 ahi(3) = nbf 1115 blo(1) = 1 1116 bhi(1) = nbf 1117 blo(2) = 1 1118 bhi(2) = nclosed 1119 do cnt=1,nset 1120 do ivec = 1, nvec 1121 alo(1) = ivec 1122 ahi(1) = ivec 1123 do ipm = 1,ncomp ! loop over Fock matrix components +/- here 1124 ind=indx(ipm,cnt) 1125 1126 if (debug) then 1127 if (ga_nodeid().eq.0) then 1128 write(*,117) cnt,ivec,ipm,ind 1129 117 format('XX:(cnt,ivec,ipm,ind)=(', 1130 & i3,',',i3,',',i3,',',i3,')') 1131 endif 1132 endif ! end-if-debug 1133c 1134c P = 4(ij|kl) - (ik|jl) - (il|kj) 1135c ij,kl 1136c 1137c K = (ik|jl) + (il|kj) 1138c ij,kl 1139c 1140c cv cv pv cp 1141c Z = 2P.[D ] + P.[D + D ] 1142c 1143c pv cv cp pv 1144c Z = 0.5d0*Z + 0.5*K.[D - D ] 1145c 1146c cp cv cp pv 1147c Z = 0.5d0*Z - 0.5*K.[D - D ] 1148c 1149c Add the Fock matrices together overwriting the density 1150c matrices to form the results above 1151c 1152c Closed-Virtual bit 1153 if (debug) then 1154 if (ga_nodeid().eq.0) 1155 & write(*,*) '--------- g_fck -------- START' 1156 call ga_print(g_fock(ind)) 1157 if (ga_nodeid().eq.0) 1158 & write(*,*) '--------- g_fck -------- END' 1159 if (ga_nodeid().eq.0) 1160 & write(*,*) '--------- g_vecs -------- START' 1161 call ga_print(g_movec) 1162 if (ga_nodeid().eq.0) 1163 & write(*,*) '--------- g_vecs -------- END' 1164 endif ! end-if-debug 1165 1166 call ga_zero(g_tmp) 1167 call nga_matmul_patch('n','n',four,zero, 1168 & g_fock(ind),alo,ahi, 1169 & g_movec ,blo,bhi, 1170 & g_tmp ,blo,bhi) 1171 1172 if (debug) then 1173 if (ga_nodeid().eq.0) 1174 & write(*,*) '--------- FnnCno -------- START' 1175 call ga_print(g_tmp) 1176 if (ga_nodeid().eq.0) 1177 & write(*,*) '--------- FnnCno -------- END' 1178 endif ! end-if-debug 1179 1180 call ga_zero(g_tmp1) 1181 call ga_matmul_patch('t','n',one,zero, 1182 $ g_movec,voff,nmo ,1,nbf, ! MO coefficients 1183 $ g_tmp ,1 ,nbf ,1,nclosed, ! result from step 1 1184 $ g_tmp1 ,1 ,nvir,1,nclosed) ! vir-occ Fock matrix 1185 1186 if (debug) then 1187 if (ga_nodeid().eq.0) then 1188 write(*,3701) cnt,ivec,ipm 11893701 format('----- CvnFnnCno(',i3,',',i3,',',i3,')------START') 1190 endif 1191 call ga_print(g_tmp1) 1192 if (ga_nodeid().eq.0) then 1193 write(*,3702) cnt,ivec,ipm 11943702 format('----- CvnFnnCno(',i3,',',i3,',',i3,')------END') 1195 endif 1196 endif ! end-if-debug 1197 1198 if (cnt.eq.1) then 1199 1200 if (debug) then 1201 if (ga_nodeid().eq.0) 1202 & write(*,*) '------- g_Az1-re-BEF(',ipm,')------ START' 1203 call ga_print(g_Az1) 1204 if (ga_nodeid().eq.0) 1205 & write(*,*) '------- g_Az1-re-BEF(',ipm,')------ END' 1206 endif ! end-if-debug 1207 1208c Note.- The operation below does: 1209c g_ax_re= g_ax_re + 4 [4 C^T F C] --> I am not sure if this is right. 1210 1211 call update_gz_reorim1(g_Az1, ! out: = complx(g_xre,g_xim) 1212 & g_tmp1, ! in : real arr 1213 & 1, ! in : =1 -> re =2 -> im 1214 & nsub, ! in : index to sub-block in g_z 1215 & ipm, ! in : = 1 or 2 index for component 1216 & vlen, ! in : = nocc*nvir 1217 & four, ! in : scaling factor 1218 & nvir, 1219 & nclosed, 1220 & ivec) 1221 1222 if (debug) then 1223 if (ga_nodeid().eq.0) 1224 & write(*,*) '------- g_Az1-re-AFT(',ipm,')------ START' 1225 call ga_print(g_Az1) 1226 if (ga_nodeid().eq.0) 1227 & write(*,*) '------- g_Az1-re-AFT(',ipm,')------ END' 1228 endif ! end-if-debug 1229 1230 else if (cnt.eq.2) then 1231 1232 if (debug) then 1233 if (ga_nodeid().eq.0) 1234 & write(*,*) '------- g_Az1-im-BEF(',ipm,')------ START' 1235 call ga_print(g_Az1) 1236 if (ga_nodeid().eq.0) 1237 & write(*,*) '------- g_Az1-im-BEF(',ipm,')------ END' 1238 endif ! end-if-debug 1239 1240 call update_gz_reorim1(g_Az1, ! out: = complx(g_xre,g_xim) 1241 & g_tmp1, ! in : real arr 1242 & 2, ! in : =1 -> re =2 -> im 1243 & nsub, ! in : index to sub-block in g_z 1244 & ipm, ! in : = 1 or 2 index for component 1245 & vlen, ! in : = nocc*nvir 1246 & four, ! in : scaling factor 1247 & nvir, 1248 & nclosed, 1249 & ivec) 1250 1251 if (debug) then 1252 if (ga_nodeid().eq.0) 1253 & write(*,*) '------- g_Az1-im-AFT(',ipm,')------ START' 1254 call ga_print(g_Az1) 1255 if (ga_nodeid().eq.0) 1256 & write(*,*) '------- g_Az1-im-AFT(',ipm,')------ END' 1257 endif ! end-if-debug 1258 1259 endif ! end-if-cnt 1260 enddo ! end-loop-ipm lopp in +/- components 1261 enddo ! end-loop-ivec-loop in field directions 1262 enddo ! end-loop-cnt 1263 if (debug) then 1264 do ipm=1,ncomp 1265 if (ga_nodeid().eq.0) 1266 & write(*,*) '------- g_Az1-1(',ipm,')------ START' 1267 call ga_print(g_Az1) 1268 if (ga_nodeid().eq.0) 1269 & write(*,*) '------- g_Az1-1(',ipm,')------ END' 1270 enddo ! end-loop-ipm 1271 endif ! end-if-debug 1272 1273 do ipm = 1,nblock 1274 if (.not. ga_destroy(g_dens(ipm))) call errquit( 1275 & 'rohf_hessv3: ga_destroy failed g_dens',0,GA_ERR) 1276 if (.not. ga_destroy(g_fock(ipm))) call errquit 1277 & ('rohf_hessv3: ga_destroy failed g_fock',0,GA_ERR) 1278 enddo ! end-loop-ipm 1279 if (.not.ga_destroy(g_tmp)) call errquit( 1280 & 'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR) 1281 if (.not.ga_destroy(g_tmp1)) call errquit( 1282 & 'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR) 1283 return 1284 end 1285 1286 subroutine rohf_hessv_2e3_opt( 1287 & basis, ! in: basis handle 1288 & geom, ! in: geom handle 1289 & nbf, ! in: nr. basis functions 1290 & nmo, ! in: nr. MOs 1291 & nclosed, ! in: nr. double occupied MOs 1292 & nopen, ! in: nr. single occupied MOs 1293 $ g_movec, ! in: MO coefficients 1294 & oskel, ! in: =.true. -> 1295 & noskew, ! in: =.true. -> symmetric density matrix 1296 & g_x_re, ! in: 1297 & g_x_im, ! in: 1298 & g_ax_re, ! in/out: Hessian product 1299 & g_ax_im, ! in/out: Hessian product 1300 & acc, ! in: accuracy Fock construction 1301 & limag, ! in: =.true. -> imaginary component allowed 1302 & lifetime)! in: =.true. -> RE-IM =.false -> RE 1303c 1304c Purpose: Optimization of rohf_hessv_2e3() 1305c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1306c Date : 03-28-12 1307c Note.- Modifying rohf_hessv_2e3, reducing computation of 1308c two-electron integrals by putting together 1309c symmetric+antisymmetric perturbed densities: g_dens(ipm) ipm=1,2 1310c and doing a single call to shell_fock_build() when using 1311c the routine shell_fock_build2(). 1312 1313 implicit none 1314#include "errquit.fh" 1315#include "mafdecls.fh" 1316#include "global.fh" 1317#include "util.fh" 1318#include "cscfps.fh" 1319#include "rtdb.fh" 1320#include "bgj.fh" 1321#include "stdio.fh" 1322#include "case.fh" 1323c 1324c Return the ROHF orbital 2e-Hessian vector product, g_ax = A * g_x 1325c 1326c ... jochen: modified version of rohf_hessv_2e2 which keeps track 1327c of two sets of input vectors that couple via the density matrix. 1328c one could likely save some memory here by re-using temp arrays 1329c ... jochen: Also made modifications to calculate imaginary terms due 1330c to finite lifetime damping 1331ccccccccccccccc This code does NOT work for open shell!!!!!ccccccccccccccccc 1332 integer basis, geom ! basis & geom handle 1333 integer nclosed,nvir,nopen! Basis size and occupation 1334 integer nbf,nmo ! No. of linearly dependent MOs 1335 integer g_movec ! MO coefficients 1336 logical oskel 1337 integer g_x_re(2), ! 1338 & g_x_im(2) ! Argument 1339 integer g_ax_re(2), ! Hessian product 1340 & g_ax_im(2) ! Hessian product 1341 double precision acc ! Accuracy of "Fock" construction 1342 logical limag ! imaginary perturbation? 1343 integer voff,xoff 1344 integer ivec,nvec,nmul,gtype,vlen 1345 integer g_tmp,g_tmp1,g_dens(4),g_fock(4) 1346 double precision tol2e 1347 logical odebug,oprint,noskew 1348 integer dims(3),chunk(3), 1349 & alo(3),ahi(3), 1350 & blo(2),bhi(2) 1351 integer ga_create_atom_blocked 1352 external ga_create_atom_blocked, 1353 & get_undosymm_fock,update_ax_fock, 1354 & shell_fock_build2 1355 double precision one,zero,mone,four,half,mhalf,two,mtwo 1356 double precision itol_floor, itol_ceil 1357 parameter(itol_floor=1.d-15, itol_ceil=1.d-3) 1358 parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0) 1359 parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0) 1360 integer ipm ! counter for density matrix components 1361 character*(255) cstemp 1362 integer g_pmats(2),g_pmata(2),g_h1mat(2), 1363 & cnt,ind,indx(2,2), 1364 & npol,nset,nblock,ncomp 1365 double precision tenm6,coef(2,2) 1366 parameter (tenm6 = 1d-6) 1367 logical lifetime,debug 1368 data npol /1/ ! for restricted calculations 1369 data indx /1,2, ! indx(1,1),indx(1,2) 1370 & 3,4/ ! indx(2,1),indx(2,2) 1371 call ga_inquire(g_x_re(1),gtype,vlen,nvec) ! out: nvec,vlen 1372 1373 ncomp=2 ! using two components 1374 nmul=1 1375 if (npol.eq. 2) nmul=2 1376 nset =1 ! for RE 1377 nblock=2 ! for RE 1378 if (lifetime) then 1379 nset =2 ! for RE-IM 1380 nblock=4 ! for RE-IM 1381 endif 1382 1383 if (oskel) 1384 $ call errquit('rohf_h2e3: no way',0, UNKNOWN_ERR) 1385 1386 debug = (.false. .and. ga_nodeid().eq.0) ! for code development 1387 1388 if (debug) then 1389 if (ga_nodeid().eq.0) then 1390 write(*,2001) limag,lifetime,nset,nblock 1391 2001 format('(limag,lifetime,nset,nblock)=(', 1392 & L1,',',L1,',',i3,',',i3,')') 1393 endif 1394 endif ! end-if-debug 1395 1396 oprint= util_print('rohf_hessv2',print_debug) 1397 1398 if (nopen.ne.0) call errquit 1399 $ ('rohf_h2e3: does not work for open shells',nopen, 1400 & UNKNOWN_ERR) 1401 odebug = util_print('rohf_hessv', print_debug) 1402 tol2e = min(max(acc,itol_floor),itol_ceil) 1403 nvir = nmo - nclosed - nopen 1404 voff = nclosed + nopen + 1 1405 1406 dims(1) = nbf 1407 dims(2) = nbf 1408 chunk(1) = dims(1) 1409 chunk(2) = -1 1410c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---START 1411 do ipm = 1,ncomp 1412 write(cstemp,'(a,i1)') 'pmats_',ipm 1413 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 1414 & g_pmats(ipm))) call 1415 & errquit('rohf_h2e3: nga_create failed '//cstemp(1:7), 1416 & 0,GA_ERR) 1417 call ga_zero(g_pmats(ipm)) 1418 write(cstemp,'(a,i1)') 'pmata_',ipm 1419 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 1420 & g_pmata(ipm))) call 1421 & errquit('rohf_h2e3: nga_create failed '//cstemp(1:7), 1422 & 0,GA_ERR) 1423 call ga_zero(g_pmata(ipm)) 1424 write(cstemp,'(a,i1)') 'h1mat_',ipm 1425 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 1426 & g_h1mat(ipm))) call 1427 & errquit('rohf_h2e3: nga_create failed '//cstemp(1:7), 1428 & 0,GA_ERR) 1429 call ga_zero(g_h1mat(ipm)) 1430 enddo ! end-loop-ipm 1431c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---END 1432 dims(1) = nvec 1433 dims(2) = nbf 1434 dims(3) = nbf 1435 chunk(1) = dims(1) 1436 chunk(2) = -1 1437 chunk(3) = -1 1438 do ipm = 1,nblock ! =2 or 4 1439c ... allocate g_dens=[g_dens_re g_dens_im] 1440 if (.not. nga_create (MT_DBL,3,dims,'CPKS dens',chunk, 1441 & g_dens(ipm))) 1442 & call errquit('rohf_h2e3: could not allocate g_dens',555, 1443 & GA_ERR) 1444 call ga_zero(g_dens(ipm)) 1445c ... allocate g_fock=[g_fock_re g_fock_im] 1446 if (.not. nga_create (MT_DBL,3,dims,'Fockv',chunk, 1447 & g_fock(ipm))) 1448 & call errquit('rohf_h2e3: could not allocate g_fock',555, 1449 & GA_ERR) 1450 call ga_zero(g_fock(ipm)) 1451 enddo ! end-loop-ipm 1452 1453 if (debug) then 1454 if (ga_nodeid().eq.0) 1455 & write(*,*) 'BEF get_dens_reorim-RE' 1456 endif ! end-if-debug 1457 1458 call get_dens_reorim( 1459 & g_dens, ! in/ou: perturbed density matrix 1460 & 1, ! in : =1 1st block RE 1461 & g_x_re, ! in : RE 1462 & g_movec,! in : MO coefficients 1463 & nbf, ! in : nr. basis functions 1464 & nmo, ! in : nr. MOs 1465 & 1, ! in : for ipol=1 -> istart=1 1466 & nclosed,! in : nr. occupied MOs 1467 & nvir, ! in : nr. virtual MOs 1468 & nvec, ! in : nr. directions (x,y,z) 1469 & npol, ! in : nr. polarizations 1470 & limag, ! in : = .true. imaginary allowed 1471 & g_pmats,! in : scratch GA array 1472 & g_pmata,! in : scratch GA array 1473 & g_h1mat)! in : scratch GA array 1474 1475 if (debug) then 1476 if (ga_nodeid().eq.0) 1477 & write(*,*) 'BEF get_dens_reorim-IM' 1478 endif ! end-if-debug 1479 1480 if (lifetime) then 1481 1482 call get_dens_reorim( 1483 & g_dens, ! in/ou: perturbed density matrix 1484 & 2, ! in : =2 2nd block IM 1485 & g_x_im, ! in : IM 1486 & g_movec,! in : MO coefficients 1487 & nbf, ! in : nr. basis functions 1488 & nmo, ! in : nr. MOs 1489 & 1, ! in : for ipol=1 -> istart=1 1490 & nclosed,! in : nr. occupied MOs 1491 & nvir, ! in : nr. virtual MOs 1492 & nvec, ! in : nr. directions (x,y,z) 1493 & npol, ! in : nr. polarizations 1494 & limag, ! in : = .true. imaginary allowed 1495 & g_pmats,! in : scratch GA array 1496 & g_pmata,! in : scratch GA array 1497 & g_h1mat)! in : scratch GA array 1498 1499 endif ! end-if-lifetime 1500 1501 if (debug) then 1502 do ipm=1,nblock 1503 if (ga_nodeid().eq.0) 1504 & write(*,*) '--------g_dens(',ipm,')-------START' 1505 call ga_print(g_dens(ipm)) 1506 if (ga_nodeid().eq.0) 1507 & write(*,*) '--------g_dens(',ipm,')-------END' 1508 enddo ! end-loop-ipm 1509 endif ! end-if-debug 1510 1511 call shell_fock_build2(g_fock, ! out: Fock matrices 1512 & g_dens, ! in : density matrices 1513 & geom, ! in : geom handle 1514 & basis, ! in : basis handle 1515 & nbf, ! in : nr. basis functions 1516 & nvec, ! in : nr. vecs (x,y,z) 1517 & npol, ! in : npol=1 for R-DFT =2 for U-DFT 1518 & ncomp, ! in : nr. components 1519 & nblock, ! in : nr. of g_dens,g_fock blocks 1520 & .true., ! in : =.true. for symm dens 1521 & tol2e, ! in : 1522 & debug) ! in : = .true. -> debugging printouts 1523 1524 if (debug) then 1525 do ipm=1,nblock 1526 if (ga_nodeid().eq.0) 1527 & write(*,*) '------- g_fock-0(',ipm,')------ START' 1528 call ga_print(g_fock(ipm)) 1529 if (ga_nodeid().eq.0) 1530 & write(*,*) '------- g_fock-0(',ipm,')------ END' 1531 enddo ! end-loop-ipm 1532 endif ! end-if-debug 1533 1534 call get_undosymm_fock( 1535 & g_fock, ! in/ou: fock matrix 1536 & nset, ! in : =1 g_x is real, =2 g_x is complex (g_x_re,g_x_im) 1537 & nvec, ! in : nr. directions (x,y,z) 1538 & nbf, ! in : nr. basis functions 1539 & npol, ! in : nr. polarizations 1540 & nmul, ! in : =1 npol=1 =2 npol=2 (acc. JK terms) 1541 & g_pmats, ! in : scratch GA array 1542 & limag) ! in : =.true. imaginary comp. exists 1543 1544c ------- Remove GA arrays: 1545 do ipm = 1,ncomp 1546 if (.not.ga_destroy(g_pmats(ipm))) call errquit( 1547 & 'rohf_hessv3: ga_destroy failed g_pmats',0,GA_ERR) 1548 if (.not.ga_destroy(g_pmata(ipm))) call errquit( 1549 & 'rohf_hessv3: ga_destroy failed g_pmata',0,GA_ERR) 1550 if (.not.ga_destroy(g_h1mat(ipm))) call errquit( 1551 & 'rohf_hessv3: ga_destroy failed g_h1mat',0,GA_ERR) 1552 enddo ! end-loop-ipm 1553 1554 if (debug) then 1555 do ipm=1,nblock 1556 if (ga_nodeid().eq.0) 1557 & write(*,*) '------- g_fock-unsym(',ipm,')------ START' 1558 call ga_print(g_fock(ipm)) 1559 if (ga_nodeid().eq.0) 1560 & write(*,*) '------- g_fock-unsym(',ipm,')------ END' 1561 enddo ! end-loop-ipm 1562 endif ! end-if-debug 1563 1564 if (debug) then 1565 do ipm=1,ncomp 1566 if (ga_nodeid().eq.0) 1567 & write(*,*) '------- g_ax_re-0(',ipm,')------ START' 1568 call ga_print(g_ax_re(ipm)) 1569 if (ga_nodeid().eq.0) 1570 & write(*,*) '------- g_ax_re-0(',ipm,')------ END' 1571 enddo ! end-loop-ipm 1572 if (lifetime) then 1573 do ipm=1,ncomp 1574 if (ga_nodeid().eq.0) 1575 & write(*,*) '------- g_ax_im-0(',ipm,')------ START' 1576 call ga_print(g_ax_im(ipm)) 1577 if (ga_nodeid().eq.0) 1578 & write(*,*) '------- g_ax_im-0(',ipm,')------ END' 1579 enddo ! end-loop-ipm 1580 endif ! end-if-lifetime 1581 endif ! end-if-debug 1582c 1583c start loop over components of perturbing field 1584c 1585 g_tmp = ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp') 1586 g_tmp1= ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp1') 1587 alo(2) = 1 1588 ahi(2) = nbf 1589 alo(3) = 1 1590 ahi(3) = nbf 1591 blo(1) = 1 1592 bhi(1) = nbf 1593 blo(2) = 1 1594 bhi(2) = nclosed 1595 do cnt=1,nset 1596 do ivec = 1, nvec 1597 alo(1) = ivec 1598 ahi(1) = ivec 1599 do ipm = 1,ncomp ! loop over Fock matrix components +/- here 1600 ind=indx(ipm,cnt) 1601 1602 if (debug) then 1603 if (ga_nodeid().eq.0) then 1604 write(*,117) cnt,ivec,ipm,ind 1605 117 format('XX:(cnt,ivec,ipm,ind)=(', 1606 & i3,',',i3,',',i3,',',i3,')') 1607 endif 1608 endif ! end-if-debug 1609c 1610c P = 4(ij|kl) - (ik|jl) - (il|kj) 1611c ij,kl 1612c 1613c K = (ik|jl) + (il|kj) 1614c ij,kl 1615c 1616c cv cv pv cp 1617c Z = 2P.[D ] + P.[D + D ] 1618c 1619c pv cv cp pv 1620c Z = 0.5d0*Z + 0.5*K.[D - D ] 1621c 1622c cp cv cp pv 1623c Z = 0.5d0*Z - 0.5*K.[D - D ] 1624c 1625c Add the Fock matrices together overwriting the density 1626c matrices to form the results above 1627c 1628c Closed-Virtual bit 1629 if (debug) then 1630 if (ga_nodeid().eq.0) 1631 & write(*,*) '--------- g_fck -------- START' 1632 call ga_print(g_fock(ind)) 1633 if (ga_nodeid().eq.0) 1634 & write(*,*) '--------- g_fck -------- END' 1635 if (ga_nodeid().eq.0) 1636 & write(*,*) '--------- g_vecs -------- START' 1637 call ga_print(g_movec) 1638 if (ga_nodeid().eq.0) 1639 & write(*,*) '--------- g_vecs -------- END' 1640 endif ! end-if-debug 1641 call ga_zero(g_tmp) 1642 call nga_matmul_patch('n','n',four,zero, 1643 & g_fock(ind),alo,ahi, 1644 & g_movec ,blo,bhi, 1645 & g_tmp ,blo,bhi) 1646 if (debug) then 1647 if (ga_nodeid().eq.0) 1648 & write(*,*) '--------- FnnCno -------- START' 1649 call ga_print(g_tmp) 1650 if (ga_nodeid().eq.0) 1651 & write(*,*) '--------- FnnCno -------- END' 1652 endif ! end-if-debug 1653 1654 call ga_zero(g_tmp1) 1655 call ga_matmul_patch('t','n',one,zero, 1656 $ g_movec,voff,nmo ,1,nbf, ! MO coefficients 1657 $ g_tmp ,1 ,nbf ,1,nclosed, ! result from step 1 1658 $ g_tmp1 ,1 ,nvir,1,nclosed) ! vir-occ Fock matrix 1659 1660 if (debug) then 1661 if (ga_nodeid().eq.0) then 1662 write(*,3701) cnt,ivec,ipm 16633701 format('----- CvnFnnCno(',i3,',',i3,',',i3,')------START') 1664 endif 1665 call ga_print(g_tmp1) 1666 if (ga_nodeid().eq.0) then 1667 write(*,3702) cnt,ivec,ipm 16683702 format('----- CvnFnnCno(',i3,',',i3,',',i3,')------END') 1669 endif 1670 endif ! end-if-debug 1671 1672 if (cnt.eq.1) then 1673 1674 if (debug) then 1675 if (ga_nodeid().eq.0) 1676 & write(*,*) '--------- g_ax-re-BEF-------- START' 1677 call ga_print(g_ax_re(ipm)) 1678 if (ga_nodeid().eq.0) 1679 & write(*,*) '--------- g_ax-re-BEF-------- END' 1680 endif ! end-if-debug 1681 1682c Note.- The operation below does: 1683c g_ax_re= g_ax_re + 4 [4 C^T F C] --> I am not sure if this is right. 1684 call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed, 1685 $ g_ax_re(ipm),1,ivec,four,'+') 1686 1687 if (debug) then 1688 if (ga_nodeid().eq.0) 1689 & write(*,*) '--------- g_ax-re-AFT-------- START' 1690 call ga_print(g_ax_re(ipm)) 1691 if (ga_nodeid().eq.0) 1692 & write(*,*) '--------- g_ax-re-AFT-------- END' 1693 endif ! end-if-debug 1694 1695 else if (cnt.eq.2) then 1696 1697 if (debug) then 1698 if (ga_nodeid().eq.0) 1699 & write(*,*) '--------- g_ax-im-BEF-------- START' 1700 call ga_print(g_ax_im(ipm)) 1701 if (ga_nodeid().eq.0) 1702 & write(*,*) '--------- g_ax-im-BEF-------- END' 1703 endif ! end-if-debug 1704 1705 call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed, 1706 $ g_ax_im(ipm),1,ivec,four,'+') 1707 1708 if (debug) then 1709 if (ga_nodeid().eq.0) 1710 & write(*,*) '--------- g_ax-im-AFT-------- START' 1711 call ga_print(g_ax_im(ipm)) 1712 if (ga_nodeid().eq.0) 1713 & write(*,*) '--------- g_ax-im-AFT-------- END' 1714 endif ! end-if-debug 1715 1716 endif ! end-if-cnt 1717 enddo ! end-loop-ipm lopp in +/- components 1718 enddo ! end-loop-ivec-loop in field directions 1719 enddo ! end-loop-cnt 1720 1721 if (debug) then 1722 do ipm=1,ncomp 1723 if (ga_nodeid().eq.0) 1724 & write(*,*) '------- g_ax_re-1(',ipm,')------ START' 1725 call ga_print(g_ax_re(ipm)) 1726 if (ga_nodeid().eq.0) 1727 & write(*,*) '------- g_ax_re-1(',ipm,')------ END' 1728 enddo ! end-loop-ipm 1729 if (lifetime) then 1730 do ipm=1,ncomp 1731 if (ga_nodeid().eq.0) 1732 & write(*,*) '------- g_ax_im-1(',ipm,')------ START' 1733 call ga_print(g_ax_im(ipm)) 1734 if (ga_nodeid().eq.0) 1735 & write(*,*) '------- g_ax_im-1(',ipm,')------ END' 1736 enddo ! end-loop-ipm 1737 endif ! end-if-lifetime 1738 endif ! end-if-debug 1739 1740 do ipm = 1,nblock 1741 if (.not. ga_destroy(g_dens(ipm))) call errquit( 1742 & 'rohf_hessv3: ga_destroy failed g_dens',0,GA_ERR) 1743 if (.not. ga_destroy(g_fock(ipm))) call errquit 1744 & ('rohf_hessv3: ga_destroy failed g_fock',0,GA_ERR) 1745 enddo ! end-loop-ipm 1746 if (.not.ga_destroy(g_tmp)) call errquit( 1747 & 'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR) 1748 if (.not.ga_destroy(g_tmp1)) call errquit( 1749 & 'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR) 1750 return 1751 end 1752 1753 subroutine get_dens_reorim( 1754 & g_dens, ! out : perturbed density matrix 1755 & cnt, ! in/ou: counter of g_dens, =1 or 2 1756 & g_x, ! in : 1757 & g_movec,! in : MO coefficients 1758 & nbf, ! in : nr. basis functions 1759 & nmo, ! in : nr. MOs 1760 & istart, ! in : shift nocc-nvirt block 1761 & nocc, ! in : nr. occupied MOs 1762 & nvir, ! in : nr. virtual MOs 1763 & nvec, ! in : nr. directions (x,y,z) 1764 & ipol, ! in : nr. polarizations 1765 & limag, ! in : = .true. imaginary allowed 1766 & g_pmats,! in : scratch GA array 1767 & g_pmata,! in : scratch GA array - NOT USED 1768 & g_h1mat)! in : scratch GA array 1769c 1770c Purpose: Calculate perturbed density matrix 1771c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1772c Date : 03-15-12 1773 1774 implicit none 1775#include "errquit.fh" 1776#include "mafdecls.fh" 1777#include "global.fh" 1778#include "util.fh" 1779#include "cscfps.fh" 1780#include "rtdb.fh" 1781#include "bgj.fh" 1782#include "stdio.fh" 1783#include "case.fh" 1784 integer g_x(2) ! Argument: g_x_re or g_x_im 1785 integer g_dens(*) ! size= 2 RE only or 4 RE+IM 1786 integer nbf,nmo,nocc,nvir 1787 integer g_movec ! MO coefficients 1788 integer dims(3),chunk(3), 1789 & alo(3),ahi(3), 1790 & blo(2),bhi(2) 1791 integer ivec,nvec,ipm,ncomp, 1792 & ipol, ! =1 for Alpha =2 for Beta 1793 & shift,cnt,ind,indx(2,2) 1794 character*(255) cstemp 1795 integer g_pmats(2),g_pmata(2),g_h1mat(2) 1796 integer istart,iend 1797 double precision tenm6,coef(2,2) 1798 parameter (tenm6 = 1d-6) 1799 logical limag,debug 1800 double precision one, zero, mone, 1801 & four, half, mhalf, two, mtwo 1802 data indx /1,2, ! indx(1,1),indx(2,1) 1803 & 3,4/ ! indx(1,2),indx(2,2) 1804 data ncomp/2/ 1805 parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0) 1806 parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0) 1807 external ga_vec_to_mat, 1808 & CalcPerturbedTDPmat1_opt 1809 1810 debug=.false. ! allow debugging printouts 1811 1812c ----- Construct coeffs for P(S),P(A) ------- START 1813 coef(1,1)= 0.5d0 1814 coef(1,2)= 0.5d0 1815 coef(2,1)=-0.5d0 1816 coef(2,2)= 0.5d0 1817 if (limag) then 1818 coef(1,1)= 0.5d0 1819 coef(1,2)=-0.5d0 1820 coef(2,1)=-0.5d0 1821 coef(2,2)=-0.5d0 1822 endif ! end-if-limag 1823 if (debug) then 1824 if (ga_nodeid().eq.0) then 1825 write(*,10) limag, 1826 & coef(1,1),coef(1,2), 1827 & coef(2,1),coef(2,2) 1828 10 format('(limag,coeff)=(',L1,',',f15.8,',', 1829 & f15.8,',',f15.8,',',f15.8,')') 1830 endif 1831 endif 1832c ----- Construct coeffs for P(S),P(A) ------- END 1833 alo(2) = 1 1834 ahi(2) = nbf 1835 alo(3) = 1 1836 ahi(3) = nbf 1837 blo(1) = 1 1838 bhi(1) = nbf 1839 blo(2) = 1 1840 bhi(2) = nbf 1841 shift=(ipol-1)*nvec 1842 iend = istart + nocc*nvir - 1 1843 do ivec = 1, nvec 1844c 1845c Compute CV, PV & CP "densities" from argument vector 1846c 1847c ... jochen: skip this part and place a subroutine call instead. 1848c it calculates the perturbed density matrix in the AO basis. 1849c I keep this source code here for reference; it is left 1850c unmodified from the version of rohf_hessv2 that this 1851c subroutine was created from. 1852 do ipm = 1,ncomp 1853 call ga_zero(g_h1mat(ipm)) 1854 call ga_copy_patch('n', ! Reshape vector into matrix 1855 $ g_x(ipm) ,istart,iend,ivec,ivec, 1856 $ g_h1mat(ipm),1 ,nvir,1 ,nocc) 1857 1858 enddo ! end-loop-ipm 1859 1860 if (debug) then 1861 do ipm=1,2 1862 if (ga_nodeid().eq.0) 1863 & write(*,*) '----------g_h1mat(',ipm,')-----START' 1864 call ga_print(g_h1mat(ipm)) 1865 if (ga_nodeid().eq.0) 1866 & write(*,*) '----------g_h1mat(',ipm,')-----END' 1867 enddo ! end-loop-ip 1868 endif ! end-if-debug 1869 1870 if (debug) then 1871 if (ga_nodeid().eq.0) 1872 & write(*,*) 'In rohf_hessv_2e3: BEF CalcPerturbedTDPmat1' 1873 if (ga_nodeid().eq.0) then 1874 write(*,667) 2,nbf,nocc,nvir,nmo,limag 1875 667 format('(ncomp,nbf,nclosed,nvir,nmo,limag)=(', 1876 & i3,',',i3,',',i3,',',i3,',',i3,',',L1,')') 1877 endif 1878 endif ! end-if -debug 1879 1880 call CalcPerturbedTDPmat1_opt( 1881 & 2, ! in : nr. components to calculate 1882 & g_pmats, ! out: density matrix symmetrized 1883 & g_pmata, ! out: density matrix antisymmetrized 1884 & g_h1mat, ! in : perturbed MO coefficients 1885 & g_movec, ! in : unperturbed MO coefficients 1886 & nbf, ! in : nr. AOs 1887 & nocc, ! in : nr. occupied MOs 1888 & nvir, ! in : nr. virtual MOs 1889 & nmo, ! in : nr. MOs 1890 & .false., ! in : = .true. calc. (symm,antisymm)=(pmats,pmata) 1891 & .false., ! in : = .true. static response, dynamic otherwise 1892 & limag, ! in : = .true. if amat is imaginary instead of real 1893 & .false.) ! in : = .true. if amat contains occ-occ 1894 1895 if (debug) then 1896 if (ga_nodeid().eq.0) 1897 & write(*,*) '---- g_pmats-1-------- START' 1898 call ga_print(g_pmats(1)) 1899 if (ga_nodeid().eq.0) 1900 & write(*,*) '---- g_pmats-1-------- END' 1901 if (ga_nodeid().eq.0) 1902 & write(*,*) '---- g_pmats-2-------- START' 1903 call ga_print(g_pmats(2)) 1904 if (ga_nodeid().eq.0) 1905 & write(*,*) '---- g_pmats-2-------- END' 1906 endif ! end-if-debug 1907 1908c next 2 lines for debugging only, to force uncoupled CPKS 1909c call ga_zero(g_pmata(1)) 1910c call ga_zero(g_pmata(2)) 1911c 1912c Instead of P(+) and P(-) which are both non-symmetric for 1913c non-zero frequency 1914c we will work with a symmetrized (S) and an antisymmetrized (A) 1915c component, calculate F(S) and F(A), respectively, and construct 1916c the Fock operators F(+/-) afterwards from F(S) +/- F(A). 1917c If it works for the skew-symmetric density matrix of NMR then 1918c it should work for this problem here, too 1919c note: here is one of those ominous scalings by 1/4 1920c needed to get the correct final results 1921 call ga_scale(g_pmats(1),0.25d0) 1922 call ga_scale(g_pmats(2),0.25d0) 1923 alo(1) = shift+ivec 1924 ahi(1) = alo(1) 1925c we need to take care here of the symmetry of the density 1926c matrices depending on whether the perturbation is real 1927c or purely imaginary. 1928c 1929c this works for real, symmetric, perturbations 1930c calculate P(S) = [ P(+) + P(-)]/2 1931c calculate P(A) = [-P(+) + P(-)]/2 (wrong results 1932c with opposite sign ...) 1933 do ipm=1,ncomp 1934 ind=indx(ipm,cnt) 1935 1936 if (debug) then 1937 if (ga_nodeid().eq.0) then 1938 write(*,11) cnt,ipm,ind,ivec,coef(ipm,1),coef(ipm,2) 1939 11 format('check-ind: (cnt,ipm,ind,ivec)=(', 1940 & i3,',',i3,',',i3,',',i3,')', 1941 & 'coeff12=(',f15.8,',',f15.8,')') 1942 endif 1943 endif ! end-if-debug 1944 1945 call nga_add_patch(coef(ipm,1),g_pmats(1) ,blo,bhi, 1946 & coef(ipm,2),g_pmats(2) ,blo,bhi, 1947 & g_dens(ind),alo,ahi) 1948 if (debug) then 1949 if (ga_nodeid().eq.0) then 1950 write(*,*) '---g_dens-acc(',ind,',',ivec,')-------START' 1951 endif 1952 call ga_print(g_dens(ind)) 1953 if (ga_nodeid().eq.0) 1954 & write(*,*) '----g_dens-acc(',ind,',',ivec,')-------END' 1955 endif ! end-if-debug 1956 1957 enddo ! end-loop-ipm 1958 enddo ! end-loop-ivec 1959 return 1960 end 1961