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