1c ---------------------------------------------------------------------- 2c Subroutine used for calculating and adding in the DIM operator in 3c frequency dependent UHF calculations. This routine is highly experimental, 4c and is still being worked on. Adapted from dimqm_addop. This routine 5c works for the FD case with no damping. (Imaginary terms not needed) 6c 7c The dimqm_addop_uhf_damp subroutine is a copy to have the lifetime 8c effects taken into account. 9c 10c Author: Jeff Becca, jbb5516@psu.edu, 2017 11c ---------------------------------------------------------------------- 12 13 Subroutine dimqm_addop_uhf(g_x_r, g_x_i, g_Ax_r, g_Ax_i, 14 $ ncomp, limag, lifetime, g_dens_r, g_dens_i) 15 16c Called from: file - uhf_hessv2_ext.F 17c subroutine - uhf_hessv_2e3 18c 19c Subroutines called from: dimqm_EqmE.F, dimqm_f2d.F dim_fock_xc.F 20c 21c Calculates and adds the frequency-dependent DIM potential to the 22c response Fock matricies (real and imaginary). Requires knowledge 23c of both the real and imaginary vectors simultaneously which is 24c why this has to be located here, unlike the static routine. 25c 26 implicit none 27#include "errquit.fh" 28#include "stdio.fh" 29#include "rtdb.fh" 30#include "mafdecls.fh" 31#include "global.fh" 32#include "dimqm_constants.fh" 33#include "dimqm.fh" 34#include "geom.fh" 35#include "crohf.fh" 36#include "cscf.fh" 37c 38c Input Variables 39 integer g_x_r(2) ! A matrix handle (real) [IN] 40 integer g_x_i(2) ! A matrix handle (imaginary) [IN] 41 integer g_Ax_r(2) ! F matrix handle (real) [IN/OUT] 42 integer g_Ax_i(2) ! F matrix handle (imaginary) [IN/OUT] 43 integer ncomp ! num of components (+/-) [IN] 44 logical limag ! Imaginary perturbation? [IN] 45 logical lifetime ! Damping or no damping 46 integer g_dens_r(2) ! Perturbed pmat [IN] 47 integer g_dens_i(2) ! perturbed pmat IMAG [IN] 48c 49c Local variables 50 integer g_tmp1, g_tmp2, g_dcv 51c integer l_dimxyz, k_dimxyz 52 double precision dimxyz(3, nDIM) 53c integer l_muind, k_muind 54 double precision muind(3, nDIM, 2) 55 integer dims(3), chunk(3) 56 character*(255) cstemp 57 integer g_pmats(2), g_pmata(2), g_h1mat(2) 58 integer g_tmpwork(2) 59 integer ipm 60c integer g_dens_r(2) 61c integer g_dens_i(2) 62 integer alo(3), ahi(3) 63 integer blo(2), bhi(2) 64 integer g_dens_comp_r 65 integer g_dens_comp_i 66 integer xend 67 double precision pre_factor 68 double precision muold(3, nDIM, 2) 69 70 double precision dx_r, dy_r, dz_r 71 double precision dx_i, dy_i, dz_i 72 double precision dsum, rmax 73 external dsum 74 integer i3, ivec, n 75c integer l_fld, k_fld 76 double precision fld(3, nDIM, 2) 77 integer g_dim_r(2) 78 integer g_dim_i(2) 79 integer nvir, voff, xoff 80 integer ga_create_atom_blocked 81 external ga_create_atom_blocked 82 character*(1) direction(3) 83 character*(1) pm(2) 84 data direction /'x', 'y', 'z'/ 85 data pm /'+', '-'/ 86 logical firsttime 87c double precision screen(nDIM) 88 double precision dimErr(3,2,2) 89 double precision calcErr 90 external calcErr 91 integer id 92c 93 id = ga_nodeid() 94 if (ldebug .and. id .eq. 0) then 95 write(luout,*) "Start dimqm_addop_uhf" 96 end if 97 nvir = nmo - nclosed - nopen 98 i3 = nDIM * 3 99 g_tmp1 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp1') 100 g_tmp2 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp2') 101 102 dims(1) = nbf 103 dims(2) = nbf 104 chunk(1) = dims(1) 105 chunk(2) = -1 106 107c 108c Allocate new arrays 109c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:fld',l_fld,k_fld)) 110c $ call errquit('malloc dimrsp:fld failed',1,MA_ERR) 111c 112c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:muind', 113c $ l_muind,k_muind)) 114c $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 115c 116c if(.not.ma_push_get(mt_dbl,i3,'dimrsp:xyz',l_dimxyz,k_dimxyz)) 117c $ call errquit('malloc dimrsp:xyz failed',1,MA_ERR) 118c 119 if(.not.rtdb_get(dimqm_rtdb,'dimpar:coords', mt_dbl, i3, dimxyz)) 120 $ call errquit('get dimpar:coords failed', 1, RTDB_ERR) 121c 122 g_dens_comp_r = ga_create_atom_blocked(geom,basis, 123 $ 'real density matrix comp') 124 if (lifetime) then 125 g_dens_comp_i = ga_create_atom_blocked(geom,basis, 126 $ 'imag density matrix comp') 127 end if 128c 129c 130c call dimqm_screening(dimqm_rtdb, geom, basis, dbl_mb(k_dimxyz), 131c $ screen) 132c screen = ONE 133c 134c ============================= 135c Solve for induced dipoles +/- 136c ============================= 137c 138c Investigate the shape of g_dens_r 139c write(luout,*)'g_dens_r+' 140c call ga_print(g_dens_r(1)) 141c write(luout,*)'g_dens_r-' 142c call ga_print(g_dens_r(2)) 143c set dimension variables for patching 144 alo(2) = 1 145 ahi(2) = nbf 146 alo(3) = 1 147 ahi(3) = nbf 148 blo(1) = 1 149 bhi(1) = nbf 150 blo(2) = 1 151 bhi(2) = nbf 152c Loop over perturbations 153 do n = 1, 3 154 do ipm = 1, 2 155 call ga_zero(g_dens_comp_r) 156 if (lifetime) call ga_zero(g_dens_comp_i) 157 alo(1) = n 158 ahi(1) = n 159c 160c Copy current perturbation into g_dens_comp 161 call nga_copy_patch('N',g_dens_r(ipm), alo, ahi, 162 $ g_dens_comp_r, blo, bhi) 163 if (lifetime) then 164 call nga_copy_patch('N',g_dens_i(ipm), alo, ahi, 165 $ g_dens_comp_i, blo, bhi) 166 end if 167 muind = ZERO 168 fld = ZERO 169 firsttime = .false. 170 if(.not.rtdb_get(dimqm_rtdb, 171 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 172 $ mt_dbl, i3, muold(:,:,1))) then 173 if(id.eq.0 .and. ldebug) 174 $ write(luout,*) "First cycle, no old dipoles!" 175 muold = ZERO 176 firsttime = .true. 177 dimqm_seeded = .false. 178c xyz_seeded(3*(n-1)+ipm) = .false. 179 if(dimtol0 < 1.0d-4 .and. .not. dimqm_noseed) then 180 dimtolxyz(ipm*3 - 1 + n) = 1.0d-4 181 if(id.eq.0 .and. ldebug) then 182 write(luout,*) "Requested tolerance below 1.0d-4" 183 write(luout,*) "Setting "//direction(n)//pm(ipm)// 184 $ " dir tolerance to 1.0d-4 to start" 185 end if 186 end if 187 else 188 if(.not.rtdb_get(dimqm_rtdb, 189 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 190 $ mt_dbl, i3, muold(:,:,2))) 191 $ call errquit('get dimqm:muold failed',1,RTDB_ERR) 192 end if 193c Set convergence tolerance 194c dimtol = dimtolxyz(ipm*3 - 1 + n) 195c dimqm_seeded = xyz_seeded(ipm*3 - 1 + n) 196c dimtol = 1.0d-7 197c dimqm_noseed = .true. 198c call dfill(i3*2, ZERO, dbl_mb(k_muind), 1) 199c call dfill(i3*2, ZERO, dbl_mb(k_fld), 1) 200c 201c Real portion of E-Field 202c write(luout,*) "REAL" 203c call ga_print(g_dens_comp_r) 204 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_r, geom, basis, 205 $ fld(:,:,1), dimxyz) 206c 207c Imaginary portion of E-Field 208c write(luout,*) "IMAG" 209c call ga_print(g_dens_comp_i) 210 if (lifetime) then 211 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_i, geom, basis, 212 $ fld(:,:,2), dimxyz) 213 end if 214c 215c Solve for induced dipoles 216 call dimqm_f2d(dimqm_rtdb, fld, muind, muold, dimxyz, 2, 217 $ direction(n), pm(ipm),.false.) 218c 219c Write induced dipoles to RTDB 220 dx_r = SUM(muind(1,:,1)) 221 dy_r = SUM(muind(2,:,1)) 222 dz_r = SUM(muind(3,:,1)) 223 dx_i = SUM(muind(1,:,2)) 224 dy_i = SUM(muind(2,:,2)) 225 dz_i = SUM(muind(3,:,2)) 226 if(id.eq.0.and.ldebug) then 227 write(luout,*) "Total induced dipole moment for "// 228 $ direction(n)//pm(ipm)//" perturbation" 229 write(luout,*) "X:", dx_r, dx_i 230 write(luout,*) "Y:", dy_r, dy_i 231 write(luout,*) "Z:", dz_r, dz_i 232 write(luout,*) '' 233 end if 234 dimErr(n, ipm, 1) = calcErr(i3, muold(:,:,1), muind(:,:,1)) 235 dimErr(n, ipm, 2) = calcErr(i3, muold(:,:,2), muind(:,:,2)) 236 if(id.eq.0.and.ldebug) then 237 write(luout,*) "Max error in real dipoles:", 238 $ dimErr(n, ipm, 1) 239 write(luout,*) "Max error in imag dipoles:", 240 $ dimErr(n, ipm, 2) 241 end if 242c if(dimErr(n, ipm, 1)/dimtol < HUNDRED 243c $ .and. dimErr(n, ipm, 2)/dimtol < HUNDRED 244c $ .and. .not. xyz_seeded(ipm*3 - 1 + n) 245c $ .and. .not. firsttime) then 246c xyz_seeded(ipm*3 - 1 + n) = .true. 247c write(luout,*) "Error within 10^2 of", dimtol, "for "// 248c $ direction(n)//pm(ipm)//" dir" 249c write(luout,*) "Setting current "//direction(n)//pm(ipm)// 250c $ " dir as seed" 251c write(luout,*)"Reverting tolerance back to", dimtol0 252c dimtolxyz(ipm*3 - 1 + n) = dimtol0 253c end if 254 if(.not.rtdb_put(dimqm_rtdb, 255 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 256 $ mt_dbl, i3, muind(:,:,1))) 257 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 258 if(.not.rtdb_put(dimqm_rtdb, 259 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 260 $ mt_dbl, i3, muind(:,:,2))) 261 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 262 end do ! ipm = 1, 2 263 end do ! ivec = 1, 3 264c if(MAXVAL(dimErr) <= 1.0d-4) then 265c write(luout,*) "Dipole error below 1d-4" 266c write(luout,*) "Shutting down DIM" 267c dimqm_on = .false. 268c end if 269c 270c Destroy GAs we don't need anymore 271 if (.not. ga_destroy(g_dens_comp_r)) call errquit 272 $ ('addop: dens_comp_r GA?',0, GA_ERR) 273 if (lifetime) then 274 if (.not. ga_destroy(g_dens_comp_i)) call errquit 275 $ ('addop: dens_comp_i GA?',0, GA_ERR) 276 end if 277c do ipm = 1,2 278c if (.not. ga_destroy(g_dens_r(ipm))) call errquit 279c $ ('addop: dens_r GA?',0, GA_ERR) 280c if (lifetime) then 281c if (.not. ga_destroy(g_dens_i(ipm))) call errquit 282c $ ('addop: dens_i GA?',0, GA_ERR) 283c endif 284c end do 285c 286c Deallocate l_fld, l_muind, l_dimxyz 287c if (.not. ma_chop_stack(l_fld)) call errquit 288c $ ('addop: fld MA?', 0, MA_ERR) 289c 290c ==================================================== 291c Solve for DIM potential, both real and imaginary S/A 292c ==================================================== 293c 294 dims(1) = 3 295 dims(2) = nbf 296 dims(3) = nbf 297 chunk(1) = dims(1) 298 chunk(2) = -1 299 chunk(3) = -1 300c 301c Real + 302 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r+',chunk, 303 & g_dim_r(1))) 304 & call errquit('addop: could not allocate g_dim_r+',1,GA_ERR) 305 call ga_zero(g_dim_r(1)) 306 call fock_dim(geom, nbf, basis, 3, g_dim_r(1), 1, 1) 307 call ga_symmetrize(g_dim_r(1)) 308c 309c Real - 310 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r-',chunk, 311 & g_dim_r(2))) 312 & call errquit('addop: could not allocate g_dim_r-',1,GA_ERR) 313 call ga_zero(g_dim_r(2)) 314 call fock_dim(geom, nbf, basis, 3, g_dim_r(2), 2, 1) 315 call ga_antisymmetrize(g_dim_r(2)) 316 if (lifetime) then 317c 318c Imaginary + 319 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i+',chunk, 320 & g_dim_i(1))) 321 & call errquit('addop: could not allocate g_dim_i+',1,GA_ERR) 322 call ga_zero(g_dim_i(1)) 323 call fock_dim(geom, nbf, basis, 3, g_dim_i(1), 1, 2) 324 call ga_symmetrize(g_dim_i(1)) 325c 326c Imaginary - 327 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i-',chunk, 328 & g_dim_i(2))) 329 & call errquit('addop: could not allocate g_dim_i-',1,GA_ERR) 330 call ga_zero(g_dim_i(2)) 331 call fock_dim(geom, nbf, basis, 3, g_dim_i(2), 2, 2) 332 call ga_antisymmetrize(g_dim_i(2)) 333 end if 334c 335c ====================================== 336c Undo the symmetrization to recover +/- 337c ====================================== 338 blo(1) = nbf 339 blo(2) = nbf 340 chunk(1) = blo(1) 341 chunk(2) = -1 342 343 do ipm = 1, ncomp 344 write(cstemp, '(a,i1)') 'g_tmpwork_',ipm 345 if (.not.nga_create(MT_DBL,2,blo,cstemp(1:11),chunk, 346 $ g_tmpwork(ipm))) call errquit('dim_addop_ufh: 347 $ nga_create failed '//cstemp(1:11),0,GA_ERR) 348 call ga_zero(g_tmpwork(ipm)) 349 enddo 350c reset blo and bhi for future use 351 blo(1) = 1 352 bhi(1) = nbf 353 blo(2) = 1 354 bhi(2) = nbf 355c 356 do ivec = 1, 3 357 alo(1) = ivec 358 ahi(1) = ivec 359c ************ 360c Real portion 361c ************ 362c TODO: I think the g_pmats here are being used just as temp arrays, 363c but I really do not know for sure. 364c all uses of g_pmats are being switched to g_tmpwork and the arrays 365c are allocated here. 366 call nga_copy_patch('N',g_dim_r(1),alo,ahi,g_tmpwork(1),blo,bhi) 367 call nga_copy_patch('N',g_dim_r(2),alo,ahi,g_tmpwork(2),blo,bhi) 368 369c 370c it might be necessary to use 0.5 here instead of 1.0 371c (note: that turned out NOT to be the case after some testing) 372 pre_factor = 1.0d0 373 call ga_sync() 374 if (.not.limag) then 375c real perturbation: 376 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 377 & pre_factor, g_tmpwork(2), blo, bhi, 378 & g_dim_r(1), alo, ahi) 379 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 380 & -pre_factor, g_tmpwork(2), blo, bhi, 381 & g_dim_r(2), alo, ahi) 382 else 383c imaginary perturbation: 384 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 385 & pre_factor, g_tmpwork(2), blo, bhi, 386 & g_dim_r(1), alo, ahi) 387 call nga_add_patch (-pre_factor, g_tmpwork(1), blo, bhi, 388 & pre_factor, g_tmpwork(2), blo, bhi, 389 & g_dim_r(2), alo, ahi) 390 end if ! if .not.limag 391 if (lifetime) then 392c ***************** 393c Imaginary portion 394c ***************** 395 call nga_copy_patch('N',g_dim_i(1),alo,ahi,g_tmpwork(1),blo,bhi) 396 call nga_copy_patch('N',g_dim_i(2),alo,ahi,g_tmpwork(2),blo,bhi) 397c 398c it might be necessary to use 0.5 here instead of 1.0 399c (note: that turned out NOT to be the case after some testing) 400 pre_factor = 1.0d0 401 call ga_sync() 402 if (.not.limag) then 403c real perturbation: 404 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 405 & pre_factor, g_tmpwork(2), blo, bhi, 406 & g_dim_i(1), alo, ahi) 407 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 408 & -pre_factor, g_tmpwork(2), blo, bhi, 409 & g_dim_i(2), alo, ahi) 410 else 411c imaginary perturbation: 412 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 413 & pre_factor, g_tmpwork(2), blo, bhi, 414 & g_dim_i(1), alo, ahi) 415 call nga_add_patch (-pre_factor, g_tmpwork(1), blo, bhi, 416 & pre_factor, g_tmpwork(2), blo, bhi, 417 & g_dim_i(2), alo, ahi) 418 end if ! if .not.limag 419 end if ! lifetime 420 enddo ! ivec = 1,nvec 421 422c Deallocate arrays no longer needed 423 do ipm = 1, ncomp 424 if (.not.ga_destroy(g_tmpwork(ipm))) call errquit('dim_addop 425 $ _uhf: g_tmpwork GA?', 0, GA_ERR) 426 enddo 427 428100 continue 429c 430c ==================================== 431c Add DIM potential to the Fock matrix 432c ==================================== 433c 434c call ga_print(g_movecs) 435 436 g_dcv = ga_create_atom_blocked(geom, basis, 'rohf_h2e3: dcv') 437 xoff = 1 438 voff = nclosed + nopen + 1 439 xend = nvir*nclosed 440 do ivec = 1, 3 ! Loop over perturbations 441 alo(1) = ivec 442 ahi(1) = ivec 443 do ipm = 1, ncomp! Loop over +/- 444c We only add the + direction of the DIM potential to both +/- of the Fock matrix 445c Real Portion 446 call nga_copy_patch('N',g_dim_r(ipm),alo,ahi,g_dcv,blo,bhi) 447 call ga_scale(g_dcv, four) 448 call ga_matmul_patch('n', 'n', two, zero, 449 $ g_dcv, 1, nbf, 1, nbf, 450 $ g_movecs, 1, nbf, 1, nclosed, 451 $ g_tmp1, 1, nbf, 1, nclosed) 452 call ga_sync() 453 call ga_matmul_patch('t', 'n', one, zero, 454 $ g_movecs, voff, nmo, 1, nbf, 455 $ g_tmp1, 1, nbf, 1, nclosed, 456 $ g_tmp2, 1, nvir, 1, nclosed) 457 call ga_sync() 458 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_r(ipm), 459 $ xoff, ivec, four, '+') 460c 461c Imaginary Portion 462 if (lifetime) then 463 call nga_copy_patch('N',g_dim_i(ipm),alo,ahi,g_dcv,blo,bhi) 464 call ga_scale(g_dcv, two) 465 call ga_matmul_patch('n', 'n', two, zero, 466 $ g_dcv, 1, nbf, 1, nbf, 467 $ g_movecs, 1, nbf, 1, nclosed, 468 $ g_tmp1, 1, nbf, 1, nclosed) 469 call ga_sync() 470 call ga_matmul_patch('t', 'n', one, zero, 471 $ g_movecs, voff, nmo, 1, nbf, 472 $ g_tmp1, 1, nbf, 1, nclosed, 473 $ g_tmp2, 1, nvir, 1, nclosed) 474 call ga_sync() 475c 476c ******NOTE JEM****** 477c We can remove the - sign if we apply the sign change discussed in aoresponse and rohf_hessv_xx3 478c ******************** 479c 480 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_i(ipm), 481 $ xoff, ivec, four, '+') 482 end if !lifetime 483 end do !ipm = 1, 2 484 end do !ivec = 1, 3 485c ======== 486c Clean up 487c ======= 488 lfirst = .false. 489 do ipm = 1,2 490c if (.not. ga_destroy(g_pmats(ipm))) call errquit 491c $ ('addop: pmats GA?', 0, GA_ERR) 492c if (.not. ga_destroy(g_pmata(ipm))) call errquit 493c $ ('addop: pmata GA?', 0, GA_ERR) 494c if (.not. ga_destroy(g_h1mat(ipm))) call errquit 495c $ ('addop: h1mat GA?', 0, GA_ERR) 496 if (.not. ga_destroy(g_dim_r(ipm))) call errquit 497 $ ('addop: dim_r GA?', 0, GA_ERR) 498 if (lifetime) then 499 if (.not. ga_destroy(g_dim_i(ipm))) call errquit 500 $ ('addop: dim_i GA?', 0, GA_ERR) 501 end if 502 enddo ! ipm = 1,2 503c 504 if (.not. ga_destroy(g_tmp1)) call errquit 505 $ ('addop: tmp1 GA?', 0, GA_ERR) 506 if (.not. ga_destroy(g_tmp2)) call errquit 507 $ ('addop: tmp2 GA?', 0, GA_ERR) 508 if (.not. ga_destroy(g_dcv)) call errquit 509 $ ('addop: dcv GA?',0, GA_ERR) 510 write(luout,*)'end of dimqm_addop_uhf' 511c 512 end subroutine dimqm_addop_uhf 513 514 515c ---------------------------------------------------------------------- 516c Subroutine used for calculating and adding in the DIM operator in 517c frequency dependent UHF calculations. This routine is highly experimental, 518c and is still being worked on. Adapted from dimqm_addop. This routine 519c works for the FD case with damping. 520c 521c Modeled after dimqm_addop_uhf 522c 523c Author: Jeff Becca, jbb5516@psu.edu, 2017 524c currently not used 525c ---------------------------------------------------------------------- 526 527 528 Subroutine dimqm_addop_uhf_damp(g_Ax_r, g_Ax_i, 529 $ ncomp, limag, lifetime, g_dens_r, g_dens_i) 530 531c Called from: Nothing currently, but will be uhf_hessv2_ext. 532c 533c Subroutines called from: dimqm_EqmE.F 534c dimqm_f2d.F dim_fock_xc.F 535c 536c Calculates and adds the frequency-dependent DIM potential to the 537c response Fock matricies (real and imaginary). Requires knowledge 538c of both the real and imaginary vectors simultaneously which is 539c why this has to be located here, unlike the static routine. 540c 541 implicit none 542#include "errquit.fh" 543#include "stdio.fh" 544#include "rtdb.fh" 545#include "mafdecls.fh" 546#include "global.fh" 547#include "dimqm_constants.fh" 548#include "dimqm.fh" 549#include "geom.fh" 550#include "crohf.fh" 551#include "cscf.fh" 552c 553c Input Variables 554 integer g_Ax_r(2) ! F matrix handle (real) [IN/OUT] 555 integer g_Ax_i(2) ! F matrix handle (imaginary) [IN/OUT] 556 integer ncomp ! num of components (+/-) [IN] 557 logical limag ! Imaginary perturbation? [IN] 558 logical lifetime ! Damping or no damping 559 integer g_dens_r(2) ! Perturbed pmat [IN] 560 integer g_dens_i(2) ! perturbed pmat IMAG [IN] 561c 562c Local variables 563 integer g_tmp1, g_tmp2, g_dcv 564c integer l_dimxyz, k_dimxyz 565 double precision dimxyz(3, nDIM) 566c integer l_muind, k_muind 567 double precision muind(3, nDIM, 2) 568 integer dims(3), chunk(3) 569 character*(255) cstemp 570c integer g_pmats(2), g_pmata(2), g_h1mat(2) 571 integer g_tmpwork(2) 572 integer ipm 573c integer g_dens_r(2) 574c integer g_dens_i(2) 575 integer alo(3), ahi(3) 576 integer blo(2), bhi(2) 577 integer g_dens_comp_r 578 integer g_dens_comp_i 579 integer xend 580 double precision pre_factor 581 double precision muold(3, nDIM, 2) 582 583 double precision dx_r, dy_r, dz_r 584 double precision dx_i, dy_i, dz_i 585 double precision dsum, rmax 586 external dsum 587 integer i3, ivec, n 588c integer l_fld, k_fld 589 double precision fld(3, nDIM, 2) 590 integer g_dim_r(2) 591 integer g_dim_i(2) 592 integer nvir, voff, xoff 593 integer ga_create_atom_blocked 594 external ga_create_atom_blocked 595 character*(1) direction(3) 596 character*(1) pm(2) 597 data direction /'x', 'y', 'z'/ 598 data pm /'+', '-'/ 599 logical firsttime 600c double precision screen(nDIM) 601 double precision dimErr(3,2,2) 602 double precision calcErr 603 external calcErr 604 integer id 605c 606 id = ga_nodeid() 607 if (ldebug .and. id .eq. 0) then 608 write(luout,*) "Start dimqm_addop_uhf" 609 end if 610 nvir = nmo - nclosed - nopen 611 i3 = nDIM * 3 612 g_tmp1 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp1') 613 g_tmp2 = ga_create_atom_blocked(geom, basis, 'dim_addop: tmp2') 614 615 dims(1) = nbf 616 dims(2) = nbf 617 chunk(1) = dims(1) 618 chunk(2) = -1 619 620c 621c Allocate new arrays 622c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:fld',l_fld,k_fld)) 623c $ call errquit('malloc dimrsp:fld failed',1,MA_ERR) 624c 625c if(.not.ma_push_get(mt_dbl,i3*2,'dimrsp:muind', 626c $ l_muind,k_muind)) 627c $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 628c 629c if(.not.ma_push_get(mt_dbl,i3,'dimrsp:xyz',l_dimxyz,k_dimxyz)) 630c $ call errquit('malloc dimrsp:xyz failed',1,MA_ERR) 631c 632 if(.not.rtdb_get(dimqm_rtdb,'dimpar:coords', mt_dbl, i3, dimxyz)) 633 $ call errquit('get dimpar:coords failed', 1, RTDB_ERR) 634c 635 g_dens_comp_r = ga_create_atom_blocked(geom,basis, 636 $ 'real density matrix comp') 637 if (lifetime) then 638 g_dens_comp_i = ga_create_atom_blocked(geom,basis, 639 $ 'imag density matrix comp') 640 end if 641c 642c 643c call dimqm_screening(dimqm_rtdb, geom, basis, dbl_mb(k_dimxyz), 644c $ screen) 645c screen = ONE 646c 647c ============================= 648c Solve for induced dipoles +/- 649c ============================= 650c 651c Investigate the shape of g_dens_r 652 write(luout,*)'g_dens_r+' 653c call ga_print(g_dens_r(1)) 654c write(luout,*)'g_dens_r-' 655c call ga_print(g_dens_r(2)) 656c set dimension variables for patching 657 alo(2) = 1 658 ahi(2) = nbf 659 alo(3) = 1 660 ahi(3) = nbf 661 blo(1) = 1 662 bhi(1) = nbf 663 blo(2) = 1 664 bhi(2) = nbf 665c Loop over perturbations 666 do n = 1, 3 667 do ipm = 1, 2 668 call ga_zero(g_dens_comp_r) 669 if (lifetime) call ga_zero(g_dens_comp_i) 670 alo(1) = n 671 ahi(1) = n 672c 673 write(luout,*)'debug 1' 674c Copy current perturbation into g_dens_comp 675 call nga_copy_patch('N',g_dens_r(ipm), alo, ahi, 676 $ g_dens_comp_r, blo, bhi) 677 write(luout,*)'debug 2' 678 if (lifetime) then 679 call nga_copy_patch('N',g_dens_i(ipm), alo, ahi, 680 $ g_dens_comp_i, blo, bhi) 681 write(luout,*)'debug 3' 682 end if 683 muind = ZERO 684 fld = ZERO 685 firsttime = .false. 686 if(.not.rtdb_get(dimqm_rtdb, 687 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 688 $ mt_dbl, i3, muold(:,:,1))) then 689 if(id.eq.0 .and. ldebug) 690 $ write(luout,*) "First cycle, no old dipoles!" 691 muold = ZERO 692 firsttime = .true. 693 dimqm_seeded = .false. 694c xyz_seeded(3*(n-1)+ipm) = .false. 695 if(dimtol0 < 1.0d-4 .and. .not. dimqm_noseed) then 696 dimtolxyz(ipm*3 - 1 + n) = 1.0d-4 697 if(id.eq.0) then 698 write(luout,*) "Requested tolerance below 1.0d-4" 699 write(luout,*) "Setting "//direction(n)//pm(ipm)// 700 $ " dir tolerance to 1.0d-4 to start" 701 end if 702 end if 703 else 704 if(.not.rtdb_get(dimqm_rtdb, 705 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 706 $ mt_dbl, i3, muold(:,:,2))) 707 $ call errquit('get dimqm:muold failed',1,RTDB_ERR) 708 end if 709 write(luout,*)'debug 4' 710c Set convergence tolerance 711c dimtol = dimtolxyz(ipm*3 - 1 + n) 712c dimqm_seeded = xyz_seeded(ipm*3 - 1 + n) 713c dimtol = 1.0d-7 714c dimqm_noseed = .true. 715c call dfill(i3*2, ZERO, dbl_mb(k_muind), 1) 716c call dfill(i3*2, ZERO, dbl_mb(k_fld), 1) 717c 718c Real portion of E-Field 719c write(luout,*) "REAL" 720c call ga_print(g_dens_comp_r) 721 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_r, geom, basis, 722 $ fld(:,:,1), dimxyz) 723c 724c Imaginary portion of E-Field 725c write(luout,*) "IMAG" 726c call ga_print(g_dens_comp_i) 727 if (lifetime) then 728 call dimqm_EqmE(dimqm_rtdb, g_dens_comp_i, geom, basis, 729 $ fld(:,:,2), dimxyz) 730 end if 731c 732c Solve for induced dipoles 733c TODO: this is always called in with 2, which makes it go to complex 734c solver. Should this happen for FD not damped case? 735 call dimqm_f2d(dimqm_rtdb, fld, muind, muold, dimxyz, 2, 736 $ direction(n), pm(ipm),.false.) 737c 738c Write induced dipoles to RTDB 739 dx_r = SUM(muind(1,:,1)) 740 dy_r = SUM(muind(2,:,1)) 741 dz_r = SUM(muind(3,:,1)) 742 dx_i = SUM(muind(1,:,2)) 743 dy_i = SUM(muind(2,:,2)) 744 dz_i = SUM(muind(3,:,2)) 745 if(id.eq.0) then 746 write(luout,*) "Total induced dipole moment for "// 747 $ direction(n)//pm(ipm)//" perturbation" 748 write(luout,*) "X:", dx_r, dx_i 749 write(luout,*) "Y:", dy_r, dy_i 750 write(luout,*) "Z:", dz_r, dz_i 751 write(luout,*) '' 752 end if 753 dimErr(n, ipm, 1) = calcErr(i3, muold(:,:,1), muind(:,:,1)) 754 dimErr(n, ipm, 2) = calcErr(i3, muold(:,:,2), muind(:,:,2)) 755 if(id.eq.0) then 756 write(luout,*) "Max error in real dipoles:", 757 $ dimErr(n, ipm, 1) 758 write(luout,*) "Max error in imag dipoles:", 759 $ dimErr(n, ipm, 2) 760 end if 761c if(dimErr(n, ipm, 1)/dimtol < HUNDRED 762c $ .and. dimErr(n, ipm, 2)/dimtol < HUNDRED 763c $ .and. .not. xyz_seeded(ipm*3 - 1 + n) 764c $ .and. .not. firsttime) then 765c xyz_seeded(ipm*3 - 1 + n) = .true. 766c write(luout,*) "Error within 10^2 of", dimtol, "for "// 767c $ direction(n)//pm(ipm)//" dir" 768c write(luout,*) "Setting current "//direction(n)//pm(ipm)// 769c $ " dir as seed" 770c write(luout,*)"Reverting tolerance back to", dimtol0 771c dimtolxyz(ipm*3 - 1 + n) = dimtol0 772c end if 773 if(.not.rtdb_put(dimqm_rtdb, 774 $ 'dimqm:muind_'//direction(n)//'_r'//pm(ipm), 775 $ mt_dbl, i3, muind(:,:,1))) 776 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 777 if(.not.rtdb_put(dimqm_rtdb, 778 $ 'dimqm:muind_'//direction(n)//'_i'//pm(ipm), 779 $ mt_dbl, i3, muind(:,:,2))) 780 $ call errquit('put dimqm:muind_p failed',1,RTDB_ERR) 781 end do ! ipm = 1, 2 782 end do ! ivec = 1, 3 783c if(MAXVAL(dimErr) <= 1.0d-4) then 784c write(luout,*) "Dipole error below 1d-4" 785c write(luout,*) "Shutting down DIM" 786c dimqm_on = .false. 787c end if 788c 789c Destroy GAs we don't need anymore 790 if (.not. ga_destroy(g_dens_comp_r)) call errquit 791 $ ('addop: dens_comp_r GA?',0, GA_ERR) 792 if (lifetime) then 793 if (.not. ga_destroy(g_dens_comp_i)) call errquit 794 $ ('addop: dens_comp_i GA?',0, GA_ERR) 795 end if 796c do ipm = 1,2 797c if (.not. ga_destroy(g_dens_r(ipm))) call errquit 798c $ ('addop: dens_r GA?',0, GA_ERR) 799c if (lifetime) then 800c if (.not. ga_destroy(g_dens_i(ipm))) call errquit 801c $ ('addop: dens_i GA?',0, GA_ERR) 802c endif 803c end do 804c 805c Deallocate l_fld, l_muind, l_dimxyz 806c if (.not. ma_chop_stack(l_fld)) call errquit 807c $ ('addop: fld MA?', 0, MA_ERR) 808c 809c ==================================================== 810c Solve for DIM potential, both real and imaginary S/A 811c ==================================================== 812c 813 dims(1) = 3 814 dims(2) = nbf 815 dims(3) = nbf 816 chunk(1) = dims(1) 817 chunk(2) = -1 818 chunk(3) = -1 819c 820c Real + 821 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r+',chunk, 822 & g_dim_r(1))) 823 & call errquit('addop: could not allocate g_dim_r+',1,GA_ERR) 824 call ga_zero(g_dim_r(1)) 825 call fock_dim(geom, nbf, basis, 3, g_dim_r(1), 1, 1) 826 call ga_symmetrize(g_dim_r(1)) 827c 828c Real - 829 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r-',chunk, 830 & g_dim_r(2))) 831 & call errquit('addop: could not allocate g_dim_r-',1,GA_ERR) 832 call ga_zero(g_dim_r(2)) 833 call fock_dim(geom, nbf, basis, 3, g_dim_r(2), 2, 1) 834 call ga_antisymmetrize(g_dim_r(2)) 835 if (lifetime) then 836c 837c Imaginary + 838 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i+',chunk, 839 & g_dim_i(1))) 840 & call errquit('addop: could not allocate g_dim_i+',1,GA_ERR) 841 call ga_zero(g_dim_i(1)) 842 call fock_dim(geom, nbf, basis, 3, g_dim_i(1), 1, 2) 843 call ga_symmetrize(g_dim_i(1)) 844c 845c Imaginary - 846 if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i-',chunk, 847 & g_dim_i(2))) 848 & call errquit('addop: could not allocate g_dim_i-',1,GA_ERR) 849 call ga_zero(g_dim_i(2)) 850 call fock_dim(geom, nbf, basis, 3, g_dim_i(2), 2, 2) 851 call ga_antisymmetrize(g_dim_i(2)) 852 end if 853c 854c ====================================== 855c Undo the symmetrization to recover +/- 856c ====================================== 857 blo(1) = nbf 858 blo(2) = nbf 859 chunk(1) = blo(1) 860 chunk(2) = -1 861 862 do ipm = 1, ncomp 863 write(cstemp, '(a,i1)') 'g_tmpwork_',ipm 864 if (.not.nga_create(MT_DBL,2,blo,cstemp(1:11),chunk, 865 $ g_tmpwork(ipm))) call errquit('dim_addop_ufh: 866 $ nga_create failed '//cstemp(1:11),0,GA_ERR) 867 call ga_zero(g_tmpwork(ipm)) 868 enddo 869 write(luout,*)'debug 5' 870c reset blo and bhi for future use 871 blo(1) = 1 872 bhi(1) = nbf 873 blo(2) = 1 874 bhi(2) = nbf 875c 876 do ivec = 1, 3 877 alo(1) = ivec 878 ahi(1) = ivec 879c ************ 880c Real portion 881c ************ 882c TODO: I think the g_pmats here are being used just as temp arrays, 883c but I really do not know for sure. 884c all uses of g_pmats are being switched to g_tmpwork and the arrays 885c are allocated here. 886 call nga_copy_patch('N',g_dim_r(1),alo,ahi,g_tmpwork(1),blo,bhi) 887 call nga_copy_patch('N',g_dim_r(2),alo,ahi,g_tmpwork(2),blo,bhi) 888 889 write(luout,*)'debug 6' 890c 891c it might be necessary to use 0.5 here instead of 1.0 892c (note: that turned out NOT to be the case after some testing) 893 pre_factor = 1.0d0 894 call ga_sync() 895 if (.not.limag) then 896c real perturbation: 897 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 898 & pre_factor, g_tmpwork(2), blo, bhi, 899 & g_dim_r(1), alo, ahi) 900 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 901 & -pre_factor, g_tmpwork(2), blo, bhi, 902 & g_dim_r(2), alo, ahi) 903 write(luout,*)'debug 7' 904 else 905c imaginary perturbation: 906 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 907 & pre_factor, g_tmpwork(2), blo, bhi, 908 & g_dim_r(1), alo, ahi) 909 call nga_add_patch (-pre_factor, g_tmpwork(1), blo, bhi, 910 & pre_factor, g_tmpwork(2), blo, bhi, 911 & g_dim_r(2), alo, ahi) 912 end if ! if .not.limag 913 if (lifetime) then 914c ***************** 915c Imaginary portion 916c ***************** 917 call nga_copy_patch('N',g_dim_i(1),alo,ahi,g_tmpwork(1),blo,bhi) 918 call nga_copy_patch('N',g_dim_i(2),alo,ahi,g_tmpwork(2),blo,bhi) 919c 920c it might be necessary to use 0.5 here instead of 1.0 921c (note: that turned out NOT to be the case after some testing) 922 pre_factor = 1.0d0 923 call ga_sync() 924 if (.not.limag) then 925c real perturbation: 926 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 927 & pre_factor, g_tmpwork(2), blo, bhi, 928 & g_dim_i(1), alo, ahi) 929 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 930 & -pre_factor, g_tmpwork(2), blo, bhi, 931 & g_dim_i(2), alo, ahi) 932 write(luout,*)'debug 8' 933 else 934c imaginary perturbation: 935 call nga_add_patch (pre_factor, g_tmpwork(1), blo, bhi, 936 & pre_factor, g_tmpwork(2), blo, bhi, 937 & g_dim_i(1), alo, ahi) 938 call nga_add_patch (-pre_factor, g_tmpwork(1), blo, bhi, 939 & pre_factor, g_tmpwork(2), blo, bhi, 940 & g_dim_i(2), alo, ahi) 941 end if ! if .not.limag 942 end if ! lifetime 943 enddo ! ivec = 1,nvec 944 945 write(luout,*)'debug 9' 946c Deallocate arrays no longer needed 947 do ipm = 1, ncomp 948 if (.not.ga_destroy(g_tmpwork(ipm))) call errquit('dim_addop 949 $ _uhf: g_tmpwork GA?', 0, GA_ERR) 950 enddo 951 952100 continue 953c 954c ==================================== 955c Add DIM potential to the Fock matrix 956c ==================================== 957c 958c call ga_print(g_movecs) 959c write(luout,*)'debug 10' 960 961 g_dcv = ga_create_atom_blocked(geom, basis, 'rohf_h2e3: dcv') 962 xoff = 1 963 voff = nclosed + nopen + 1 964 xend = nvir*nclosed 965 do ivec = 1, 3 ! Loop over perturbations 966 alo(1) = ivec 967 ahi(1) = ivec 968 do ipm = 1, ncomp! Loop over +/- 969c We only add the + direction of the DIM potential to both +/- of the Fock matrix 970c Real Portion 971 call nga_copy_patch('N',g_dim_r(ipm),alo,ahi,g_dcv,blo,bhi) 972 call ga_scale(g_dcv, four) 973 call ga_matmul_patch('n', 'n', two, zero, 974 $ g_dcv, 1, nbf, 1, nbf, 975 $ g_movecs, 1, nbf, 1, nclosed, 976 $ g_tmp1, 1, nbf, 1, nclosed) 977 call ga_sync() 978 call ga_matmul_patch('t', 'n', one, zero, 979 $ g_movecs, voff, nmo, 1, nbf, 980 $ g_tmp1, 1, nbf, 1, nclosed, 981 $ g_tmp2, 1, nvir, 1, nclosed) 982 call ga_sync() 983 write(luout,*)'debug 10.5' 984 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_r(ipm), 985 $ xoff, ivec, four, '+') 986c 987 write(luout,*)'debug 11' 988c Imaginary Portion 989 if (lifetime) then 990 call nga_copy_patch('N',g_dim_i(ipm),alo,ahi,g_dcv,blo,bhi) 991 call ga_scale(g_dcv, four) 992 call ga_matmul_patch('n', 'n', two, zero, 993 $ g_dcv, 1, nbf, 1, nbf, 994 $ g_movecs, 1, nbf, 1, nclosed, 995 $ g_tmp1, 1, nbf, 1, nclosed) 996 call ga_sync() 997 call ga_matmul_patch('t', 'n', one, zero, 998 $ g_movecs, voff, nmo, 1, nbf, 999 $ g_tmp1, 1, nbf, 1, nclosed, 1000 $ g_tmp2, 1, nvir, 1, nclosed) 1001 call ga_sync() 1002c 1003c ******NOTE JEM****** 1004c We can remove the - sign if we apply the sign change discussed in aoresponse and rohf_hessv_xx3 1005c ******************** 1006c 1007 call ga_mat_to_vec(g_tmp2, 1, nvir, 1, nclosed, g_Ax_i(ipm), 1008 $ xoff, ivec, four, '+') 1009 end if !lifetime 1010 end do !ipm = 1, 2 1011 end do !ivec = 1, 3 1012c ======== 1013c Clean up 1014c ======= 1015 lfirst = .false. 1016 do ipm = 1,2 1017c if (.not. ga_destroy(g_pmats(ipm))) call errquit 1018c $ ('addop: pmats GA?', 0, GA_ERR) 1019c if (.not. ga_destroy(g_pmata(ipm))) call errquit 1020c $ ('addop: pmata GA?', 0, GA_ERR) 1021c if (.not. ga_destroy(g_h1mat(ipm))) call errquit 1022c $ ('addop: h1mat GA?', 0, GA_ERR) 1023 if (.not. ga_destroy(g_dim_r(ipm))) call errquit 1024 $ ('addop: dim_r GA?', 0, GA_ERR) 1025 if (lifetime) then 1026 if (.not. ga_destroy(g_dim_i(ipm))) call errquit 1027 $ ('addop: dim_i GA?', 0, GA_ERR) 1028 end if 1029 enddo ! ipm = 1,2 1030c 1031 if (.not. ga_destroy(g_tmp1)) call errquit 1032 $ ('addop: tmp1 GA?', 0, GA_ERR) 1033 if (.not. ga_destroy(g_tmp2)) call errquit 1034 $ ('addop: tmp2 GA?', 0, GA_ERR) 1035 if (.not. ga_destroy(g_dcv)) call errquit 1036 $ ('addop: dcv GA?',0, GA_ERR) 1037 write(luout,*)'end of dimqm_addop_uhf' 1038c 1039 end subroutine dimqm_addop_uhf_damp 1040