1 subroutine dimqm_addop_complx(g_z, g_Az, ncomp, limag, lifetime) 2c ================================================================== 3c Calculates and adds the frequency-dependent DIM potential to the 4c response Fock matricies (real and imaginary). Requires knowledge 5c of both the real and imaginary vectors simultaneously which is 6c why this has to be located here, unlike the static routine. 7c 8c THIS ROUTINE IS CURRENTLY NOT USED 9c 10c ================================================================== 11 implicit none 12#include "errquit.fh" 13#include "stdio.fh" 14#include "rtdb.fh" 15#include "mafdecls.fh" 16#include "global.fh" 17#include "dimqm_constants.fh" 18#include "dimqm.fh" 19#include "geom.fh" 20#include "crohf.fh" 21#include "cscf.fh" 22c 23c Input Variables 24 integer ncomp 25 integer g_z(ncomp) ! A matrix handle (complex) [IN] 26 integer g_Az(ncomp) ! F matrix handle (complex) [IN/OUT] 27 logical limag ! Imaginary perturbation? [IN] 28 logical lifetime ! Damping or no damping [IN] 29c 30c Local variables 31 integer g_tmp1, g_tmp2, g_dcv 32c integer l_dimxyz, k_dimxyz 33 double precision dimxyz(3, nDIM) 34c integer l_muind, k_muind 35 double precision muind(3, nDIM, 2) 36 integer dims(3), chunk(3) 37 character*(255) cstemp 38 integer g_pmats(2), g_pmata(2), g_h1mat(2) 39 integer ipm 40 integer g_dens_r(2) 41 integer g_dens_i(2) 42 integer alo(3), ahi(3) 43 integer blo(2), bhi(2) 44 integer g_dens_comp_r 45 integer g_dens_comp_i 46 integer xend 47 double precision pre_factor 48 double precision muold(3, nDIM, 2) 49 50 double precision dx_r, dy_r, dz_r 51 double precision dx_i, dy_i, dz_i 52 double precision dsum 53 external dsum 54 integer i3, ivec, n 55c integer l_fld, k_fld 56 double precision fld(3, nDIM, 2) 57 integer g_dim_r(2) 58 integer g_dim_i(2) 59 integer nvir, voff, xoff 60 integer nmul, nset, nblock 61 integer ga_create_atom_blocked 62 external ga_create_atom_blocked 63 character*(1) direction(3) 64 character*(1) pm(2) 65 data direction /'x', 'y', 'z'/ 66 data pm /'+', '-'/ 67 logical firsttime 68c double precision screen(nDIM) 69 double precision dimErr(3,2,2) 70 double precision calcErr 71 external calcErr 72 integer id 73c 74 id = ga_nodeid() 75 if (ldebug .and. id .eq. 0) then 76 write(luout,*) "Start dimqm_addop" 77 end if 78 79 nvir = nmo - nclosed - nopen 80 voff = nclosed + nopen + 1 81 nmul = 1 82 if (npol .eq. 2) numl = 2 83 if (lifetime) then 84 nset = 2 85 nblock = 4 86 else 87 nset = 1 88 nblock = 2 89 end if 90 91 i3 = nDIM * 3 92 g_tmp1 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp1') 93 g_tmp2 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp2') 94 95 dims(1) = nbf 96 dims(2) = nbf 97 chunk(1) = dims(1) 98 chunk(2) = -1 99c 100c Allocate GAs 101 do ipm = 1,ncomp 102 write(cstemp,'(a,i1)') 'pmats_',ipm 103 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 104 & g_pmats(ipm))) call 105 & errquit('dim_addop: nga_create failed '//cstemp(1:7), 106 & 0,GA_ERR) 107 call ga_zero(g_pmats(ipm)) 108c write(cstemp,'(a,i1)') 'pmata_',ipm 109c if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 110c & g_pmata(ipm))) call 111c & errquit('dim_addop: nga_create failed '//cstemp(1:7), 112c & 0,GA_ERR) 113c call ga_zero(g_pmata(ipm)) 114 write(cstemp,'(a,i1)') 'h1mat_',ipm 115 if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk, 116 & g_h1mat(ipm))) call 117 & errquit('dim_addop: nga_create failed '//cstemp(1:7), 118 & 0,GA_ERR) 119 call ga_zero(g_h1mat(ipm)) 120 enddo 121 122 dims(1) = nvec 123 dims(2) = nbf 124 dims(3) = nbf 125 chunk(1) = dims(1) 126 chunk(2) = -1 127 chunk(3) = -1 128 do ipm = 1, nblock 129 if (.not. nga_create (MT_DBL, 3, dims, 'CPKS dens_r',chunk, 130 & g_dens_r(ipm))) 131 & call errquit('dim_addop: could not allocate g_dens_r',555, 132 & GA_ERR) 133 call ga_zero(g_dens_r(ipm)) 134 if (.not. nga_create (MT_DBL, 3, dims, 'CPKS dens_i',chunk, 135 & g_dens_i(ipm))) 136 & call errquit('dim_addop: could not allocate g_dens_i',555, 137 & GA_ERR) 138 call ga_zero(g_dens_i(ipm)) 139 end do 140 141 alo(2) = 1 142 ahi(2) = nbf 143 alo(3) = 1 144 ahi(3) = nbf 145 blo(1) = 1 146 bhi(1) = nbf 147 blo(2) = 1 148 bhi(2) = nbf 149c 150 do ivec = 1, 3 151c 152c ========================== 153c Real Portion of QM Density 154c ========================== 155c 156 do ipm = 1,2 157 call ga_zero(g_h1mat(ipm)) 158 call ga_vec_to_mat(g_h1mat(ipm), 1, nvir, 1, nclosed, 159 $ g_x_r(ipm), 1, ivec) ! g_h1mat now holds A-matrix for 160 call ga_zero(g_pmats(ipm)) 161 call ga_zero(g_pmata(ipm)) 162 end do ! ipm = 1,2 163 call ga_sync() 164c note: the last argument tells it not to use an occ-occ 165c block to build the density marix. 166 call CalcPerturbedTDPmat1 167 & (2, g_pmats, g_pmata, g_h1mat, g_movecs, nbf, nclosed, 168 & nvir, nmo, .false., .false., 169 & limag, .false.) ! density matrix -> pmats 170 171c call ga_zero(g_pmata(1)) 172c call ga_zero(g_pmata(2)) 173 call ga_sync() 174 call ga_scale(g_pmats(1),0.25d0) 175 call ga_scale(g_pmats(2),0.25d0) 176c 177 alo(1) = ivec 178 ahi(1) = ivec 179 if (.not.limag) then 180c this works for real, symmetric, perturbations 181c calculate P(S) = [P(+) + P(-)]/2 182 call nga_add_patch (0.5d0, g_pmats(1), blo, bhi, 183 & 0.5d0, g_pmats(2), blo, bhi, 184 & g_dens_r(1), alo, ahi) 185c caluclate P(A) = [-P(+) + P(-)]/2 (wrong results with opposite sign ...) 186 call nga_add_patch (-0.5d0, g_pmats(1), blo, bhi, 187 & 0.5d0, g_pmats(2), blo, bhi, 188 & g_dens_r(2), alo, ahi) 189 else 190c 191c this here is for imaginary, antisymmetric, perturbations 192c calculate P(S) = [P(+) - P(-)]/2 193 call nga_add_patch (0.5d0, g_pmats(1), blo, bhi, 194 & -0.5d0, g_pmats(2), blo, bhi, 195 & g_dens_r(1), alo, ahi) 196c caluclate P(A) = -[P(+) + P(-)]/2 ! sign needs to be determined 197 call nga_add_patch (-0.5d0, g_pmats(1), blo, bhi, 198 & -0.5d0, g_pmats(2), blo, bhi, 199 & g_dens_r(2), alo, ahi) 200 end if 201 call ga_sync() 202c 203c =============================== 204c Imaginary Portion of QM Density 205c =============================== 206c 207 if (lifetime) then 208 do ipm = 1,2 209 call ga_zero(g_h1mat(ipm)) 210 call ga_vec_to_mat(g_h1mat(ipm), 1, nvir, 1, nclosed, 211 $ g_x_i(ipm), 1, ivec) ! g_h1mat now holds A-matrix for 212 call ga_zero(g_pmats(ipm)) 213 call ga_zero(g_pmata(ipm)) 214 end do ! ipm = 1,2 215 call ga_sync() 216c note: the last argument tells it not to use an occ-occ 217c block to build the density marix. 218 call CalcPerturbedTDPmat1 219 & (2, g_pmats, g_pmata, g_h1mat, g_movecs, nbf, nclosed, 220 & nvir, nmo, .false., .false., 221 & limag, .false.) ! density matrix -> pmats 222c write(luout,*) "pmata" 223c call ga_print(g_pmata(1)) 224c call ga_print(g_pmata(2)) 225 226 call ga_zero(g_pmata(1)) 227 call ga_zero(g_pmata(2)) 228 call ga_sync() 229 call ga_scale(g_pmats(1),0.25d0) 230 call ga_scale(g_pmats(2),0.25d0) 231c 232 alo(1) = ivec 233 ahi(1) = ivec 234 if (.not.limag) then 235c calculate P(S) = [P(+) + P(-)]/2 236 call nga_add_patch (0.5d0, g_pmats(1), blo, bhi, 237 & 0.5d0, g_pmats(2), blo, bhi, 238 & g_dens_i(1), alo, ahi) 239c caluclate P(A) = [-P(+) + P(-)]/2 (wrong results with opposite sign ...) 240 call nga_add_patch (-0.5d0, g_pmats(1), blo, bhi, 241 & 0.5d0, g_pmats(2), blo, bhi, 242 & g_dens_i(2), alo, ahi) 243 else 244 call nga_add_patch (0.5d0, g_pmats(1), blo, bhi, 245 & -0.5d0, g_pmats(2), blo, bhi, 246 & g_dens_i(1), alo, ahi) 247c caluclate P(A) = -[P(+) + P(-)]/2 ! sign needs to be determined 248c 249 call nga_add_patch (-0.5d0, g_pmats(1), blo, bhi, 250 & -0.5d0, g_pmats(2), blo, bhi, 251 & g_dens_i(2), alo, ahi) 252 end if ! limag 253 end if ! lifetime 254 call ga_sync() 255 end do ! ivec = 1, 3 256c 257c Allocate new arrays 258c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:fld',l_fld,k_fld)) 259c $ call errquit('malloc dimrsp:fld failed',1,MA_ERR) 260c 261c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:muind', 262c $ l_muind,k_muind)) 263c $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 264c 265c if(.not.ma_push_get(mt_dbl,i3,'dimrsp:xyz',l_dimxyz,k_dimxyz)) 266c $ call errquit('malloc dimrsp:xyz failed',1,MA_ERR) 267c 268 if(.not.rtdb_get(dimqm_rtdb,'dimpar:coords', mt_dbl, i3, dimxyz)) 269 $ call errquit('get dimpar:coords failed', 1, RTDB_ERR) 270c 271 g_dens_comp_r = ga_create_atom_blocked(geom,basis, 272 $ 'real density matrix comp') 273 if (lifetime) then 274 g_dens_comp_i = ga_create_atom_blocked(geom,basis, 275 $ 'imag density matrix comp') 276 end if 277c 278c 279c call dimqm_screening(dimqm_rtdb, geom, basis, dbl_mb(k_dimxyz), 280c $ screen) 281c screen = ONE 282c 283c ============================= 284c Solve for induced dipoles +/- 285c ============================= 286c 287c Loop over perturbations 288 do n = 1, 3 289 do ipm = 1, 2 290 call ga_zero(g_dens_comp_r) 291 if (lifetime) call ga_zero(g_dens_comp_i) 292 alo(1) = n 293 ahi(1) = n 294c 295c Copy current perturbation into g_dens_comp 296 call nga_copy_patch('N',g_dens_r(ipm), alo, ahi, 297 $ g_dens_comp_r, blo, bhi) 298 if (lifetime) then 299 call nga_copy_patch('N',g_dens_i(ipm), alo, ahi, 300 $ g_dens_comp_i, blo, bhi) 301 end if 302 muind = ZERO 303 fld = ZERO 304 firsttime = .false. 305 if(.not.rtdb_get(dimqm_rtdb, 306 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 307 $ mt_dbl, i3, muold(:,:,1))) then 308 if(id.eq.0) write(luout,*) "First cycle, no old dipoles!" 309 muold = ZERO 310 firsttime = .true. 311 dimqm_seeded = .false. 312c xyz_seeded(3*(n-1)+ipm) = .false. 313 if(dimtol0 < 1.0d-4 .and. .not. dimqm_noseed) then 314 dimtolxyz(ipm*3 - 1 + n) = 1.0d-4 315 if(id.eq.0) then 316 write(luout,*) "Requested tolerance below 1.0d-4" 317 write(luout,*) "Setting "//direction(n)//pm(ipm)// 318 $ " dir tolerance to 1.0d-4 to start" 319 end if 320 end if 321 else 322 if(.not.rtdb_get(dimqm_rtdb, 323 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 324 $ mt_dbl, i3, muold(:,:,2))) 325 $ call errquit('get dimqm:muold failed',1,RTDB_ERR) 326 end if 327c Set convergence tolerance 328c dimtol = dimtolxyz(ipm*3 - 1 + n) 329c dimqm_seeded = xyz_seeded(ipm*3 - 1 + n) 330c dimtol = 1.0d-7 331c dimqm_noseed = .true. 332c call dfill(i3*2, ZERO, dbl_mb(k_muind), 1) 333c call dfill(i3*2, ZERO, dbl_mb(k_fld), 1) 334c 335c Real portion of E-Field 336c write(luout,*) "REAL" 337c call ga_print(g_dens_comp_r) 338 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_r, geom, basis, 339 $ fld(:,:,1), dimxyz) 340c 341c Imaginary portion of E-Field 342c write(luout,*) "IMAG" 343c call ga_print(g_dens_comp_i) 344 if (lifetime) then 345 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_i, geom, basis, 346 $ fld(:,:,2), dimxyz) 347 end if 348c 349c Solve for induced dipoles 350 call dimqm_f2d(dimqm_rtdb, fld, muind, muold, dimxyz, 2, 351 $ direction(n), pm(ipm)) 352c 353c Write induced dipoles to RTDB 354 dx_r = SUM(muind(1,:,1)) 355 dy_r = SUM(muind(2,:,1)) 356 dz_r = SUM(muind(3,:,1)) 357 dx_i = SUM(muind(1,:,2)) 358 dy_i = SUM(muind(2,:,2)) 359 dz_i = SUM(muind(3,:,2)) 360 if(id.eq.0) then 361 write(luout,*) "Total induced dipole moment for "// 362 $ direction(n)//pm(ipm)//" perturbation" 363 write(luout,*) "X:", dx_r, dx_i 364 write(luout,*) "Y:", dy_r, dy_i 365 write(luout,*) "Z:", dz_r, dz_i 366 write(luout,*) '' 367 end if 368 dimErr(n, ipm, 1) = calcErr(i3, muold(:,:,1), muind(:,:,1)) 369 dimErr(n, ipm, 2) = calcErr(i3, muold(:,:,2), muind(:,:,2)) 370 if(id.eq.0) then 371 write(luout,*) "Max error in real dipoles:", 372 $ dimErr(n, ipm, 1) 373 write(luout,*) "Max error in imag dipoles:", 374 $ dimErr(n, ipm, 2) 375 end if 376c if(dimErr(n, ipm, 1)/dimtol < HUNDRED 377c $ .and. dimErr(n, ipm, 2)/dimtol < HUNDRED 378c $ .and. .not. xyz_seeded(ipm*3 - 1 + n) 379c $ .and. .not. firsttime) then 380c xyz_seeded(ipm*3 - 1 + n) = .true. 381c write(luout,*) "Error within 10^2 of", dimtol, "for "// 382c $ direction(n)//pm(ipm)//" dir" 383c write(luout,*) "Setting current "//direction(n)//pm(ipm)// 384c $ " dir as seed" 385c write(luout,*)"Reverting tolerance back to", dimtol0 386c dimtolxyz(ipm*3 - 1 + n) = dimtol0 387c end if 388 if(.not.rtdb_put(dimqm_rtdb, 389 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 390 $ mt_dbl, i3, muind(:,:,1))) 391 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 392 if(.not.rtdb_put(dimqm_rtdb, 393 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 394 $ mt_dbl, i3, muind(:,:,2))) 395 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 396 end do ! ipm = 1, 2 397 end do ! ivec = 1, 3 398c if(MAXVAL(dimErr) <= 1.0d-4) then 399c write(luout,*) "Dipole error below 1d-4" 400c write(luout,*) "Shutting down DIM" 401c dimqm_on = .false. 402c end if 403c 404c Destroy GAs we don't need anymore 405 if (.not. ga_destroy(g_dens_comp_r)) call errquit 406 $ ('addop: dens_comp_r GA?',0, GA_ERR) 407 if (lifetime) then 408 if (.not. ga_destroy(g_dens_comp_i)) call errquit 409 $ ('addop: dens_comp_i GA?',0, GA_ERR) 410 end if 411 do ipm = 1,2 412 if (.not. ga_destroy(g_dens_r(ipm))) call errquit 413 $ ('addop: dens_r GA?',0, GA_ERR) 414 if (lifetime) then 415 if (.not. ga_destroy(g_dens_i(ipm))) call errquit 416 $ ('addop: dens_i GA?',0, GA_ERR) 417 end if 418 end do 419c 420c Deallocate l_fld, l_muind, l_dimxyz 421c if (.not. ma_chop_stack(l_fld)) call errquit 422c $ ('addop: fld MA?', 0, MA_ERR) 423c 424c ==================================================== 425c Solve for DIM potential, both real and imaginary S/A 426c ==================================================== 427c 428 dims(1) = 3 429 dims(2) = nbf 430 dims(3) = nbf 431 chunk(1) = dims(1) 432 chunk(2) = -1 433 chunk(3) = -1 434c 435c Real + 436 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r+',chunk, 437 & g_dim_r(1))) 438 & call errquit('addop: could not allocate g_dim_r+',1,GA_ERR) 439 call ga_zero(g_dim_r(1)) 440 call fock_dim(geom, nbf, basis, 3, g_dim_r(1), 1, 1) 441 call ga_symmetrize(g_dim_r(1)) 442c 443c Real - 444 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r-',chunk, 445 & g_dim_r(2))) 446 & call errquit('addop: could not allocate g_dim_r-',1,GA_ERR) 447 call ga_zero(g_dim_r(2)) 448 call fock_dim(geom, nbf, basis, 3, g_dim_r(2), 2, 1) 449 call ga_antisymmetrize(g_dim_r(2)) 450 if (lifetime) then 451c 452c Imaginary + 453 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i+',chunk, 454 & g_dim_i(1))) 455 & call errquit('addop: could not allocate g_dim_i+',1,GA_ERR) 456 call ga_zero(g_dim_i(1)) 457 call fock_dim(geom, nbf, basis, 3, g_dim_i(1), 1, 2) 458 call ga_symmetrize(g_dim_i(1)) 459c 460c Imaginary - 461 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i-',chunk, 462 & g_dim_i(2))) 463 & call errquit('addop: could not allocate g_dim_i-',1,GA_ERR) 464 call ga_zero(g_dim_i(2)) 465 call fock_dim(geom, nbf, basis, 3, g_dim_i(2), 2, 2) 466 call ga_antisymmetrize(g_dim_i(2)) 467 end if 468c 469c ====================================== 470c Undo the symmetrization to recover +/- 471c ====================================== 472c 473 do ivec = 1, 3 474 alo(1) = ivec 475 ahi(1) = ivec 476c ************ 477c Real portion 478c ************ 479 call nga_copy_patch ('N',g_dim_r(1),alo,ahi,g_pmats(1),blo,bhi) 480 call nga_copy_patch ('N',g_dim_r(2),alo,ahi,g_pmats(2),blo,bhi) 481c 482c it might be necessary to use 0.5 here instead of 1.0 483c (note: that turned out NOT to be the case after some testing) 484 pre_factor = 1.0d0 485 call ga_sync() 486 if (.not.limag) then 487c real perturbation: 488 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 489 & pre_factor, g_pmats(2), blo, bhi, 490 & g_dim_r(1), alo, ahi) 491 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 492 & -pre_factor, g_pmats(2), blo, bhi, 493 & g_dim_r(2), alo, ahi) 494 else 495c imaginary perturbation: 496 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 497 & pre_factor, g_pmats(2), blo, bhi, 498 & g_dim_r(1), alo, ahi) 499 call nga_add_patch (-pre_factor, g_pmats(1), blo, bhi, 500 & pre_factor, g_pmats(2), blo, bhi, 501 & g_dim_r(2), alo, ahi) 502 end if ! if .not.limag 503 if (lifetime) then 504c ***************** 505c Imaginary portion 506c ***************** 507 call nga_copy_patch ('N',g_dim_i(1),alo,ahi,g_pmats(1),blo,bhi) 508 call nga_copy_patch ('N',g_dim_i(2),alo,ahi,g_pmats(2),blo,bhi) 509c 510c it might be necessary to use 0.5 here instead of 1.0 511c (note: that turned out NOT to be the case after some testing) 512 pre_factor = 1.0d0 513 call ga_sync() 514 if (.not.limag) then 515c real perturbation: 516 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 517 & pre_factor, g_pmats(2), blo, bhi, 518 & g_dim_i(1), alo, ahi) 519 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 520 & -pre_factor, g_pmats(2), blo, bhi, 521 & g_dim_i(2), alo, ahi) 522 else 523c imaginary perturbation: 524 call nga_add_patch (pre_factor, g_pmats(1), blo, bhi, 525 & pre_factor, g_pmats(2), blo, bhi, 526 & g_dim_i(1), alo, ahi) 527 call nga_add_patch (-pre_factor, g_pmats(1), blo, bhi, 528 & pre_factor, g_pmats(2), blo, bhi, 529 & g_dim_i(2), alo, ahi) 530 end if ! if .not.limag 531 end if ! lifetime 532 enddo ! ivec = 1,nvec 533100 continue 534c 535c ==================================== 536c Add DIM potential to the Fock matrix 537c ==================================== 538c 539 g_dcv = ga_create_atom_blocked(geom, basis, 'rohf_h2e3: dcv') 540 xoff = 1 541 voff = nclosed + nopen + 1 542 xend = nvir*nclosed 543 do ivec = 1, 3 ! Loop over perturbations 544 alo(1) = ivec 545 ahi(1) = ivec 546 do ipm = 1, 2! Loop over +/- 547c We only add the + direction of the DIM potential to both +/- of the Fock matrix 548c Real Portion 549 call nga_copy_patch('N',g_dim_r(ipm),alo,ahi,g_dcv,blo,bhi) 550 call ga_scale(g_dcv, two) 551 call ga_matmul_patch('n', 'n', two, zero, 552 $ g_dcv, 1, nbf, 1, nbf, 553 $ g_movecs, 1, nbf, 1, nclosed, 554 $ g_tmp1, 1, nbf, 1, nclosed) 555 call ga_sync() 556 call ga_matmul_patch('t', 'n', one, zero, 557 $ g_movecs, voff, nmo, 1, nbf, 558 $ g_tmp1, 1, nbf, 1, nclosed, 559 $ g_tmp2, 1, nvir, 1, nclosed) 560 call ga_sync() 561 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_r(ipm), 562 $ xoff, ivec, four, '+') 563c 564c Imaginary Portion 565 if (lifetime) then 566 call nga_copy_patch('N',g_dim_i(ipm),alo,ahi,g_dcv,blo,bhi) 567 call ga_scale(g_dcv, two) 568 call ga_matmul_patch('n', 'n', two, zero, 569 $ g_dcv, 1, nbf, 1, nbf, 570 $ g_movecs, 1, nbf, 1, nclosed, 571 $ g_tmp1, 1, nbf, 1, nclosed) 572 call ga_sync() 573 call ga_matmul_patch('t', 'n', one, zero, 574 $ g_movecs, voff, nmo, 1, nbf, 575 $ g_tmp1, 1, nbf, 1, nclosed, 576 $ g_tmp2, 1, nvir, 1, nclosed) 577 call ga_sync() 578c 579c ******NOTE JEM****** 580c We can remove the - sign if we apply the sign change discussed in aoresponse and rohf_hessv_xx3 581c ******************** 582c 583 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_i(ipm), 584 $ xoff, ivec, four, '+') 585 end if !lifetime 586 end do !ipm = 1, 2 587 end do !ivec = 1, 3 588c ======== 589c Clean up 590c ======= 591 lfirst = .false. 592 do ipm = 1,2 593 if (.not. ga_destroy(g_pmats(ipm))) call errquit 594 $ ('addop: pmats GA?', 0, GA_ERR) 595 if (.not. ga_destroy(g_pmata(ipm))) call errquit 596 $ ('addop: pmata GA?', 0, GA_ERR) 597 if (.not. ga_destroy(g_h1mat(ipm))) call errquit 598 $ ('addop: h1mat GA?', 0, GA_ERR) 599 if (.not. ga_destroy(g_dim_r(ipm))) call errquit 600 $ ('addop: dim_r GA?', 0, GA_ERR) 601 if (lifetime) then 602 if (.not. ga_destroy(g_dim_i(ipm))) call errquit 603 $ ('addop: dim_i GA?', 0, GA_ERR) 604 end if 605 enddo ! ipm = 1,2 606c 607 if (.not. ga_destroy(g_tmp1)) call errquit 608 $ ('addop: tmp1 GA?', 0, GA_ERR) 609 if (.not. ga_destroy(g_tmp2)) call errquit 610 $ ('addop: tmp2 GA?', 0, GA_ERR) 611 if (.not. ga_destroy(g_dcv)) call errquit 612 $ ('addop: dcv GA?',0, GA_ERR) 613c 614 end subroutine dimqm_addop_cmplx 615