1 subroutine dimqm_lclfld(g_dipel, omega, lifetime, g_dipel_i) 2c 3c 4c 5c Called from aoresponse_driver.F 6c 7 implicit none 8#include "errquit.fh" 9#include "stdio.fh" 10#include "rtdb.fh" 11#include "mafdecls.fh" 12#include "global.fh" 13#include "dimqm_constants.fh" 14#include "dimqm.fh" 15#include "geom.fh" 16#include "crohf.fh" 17#include "cscf.fh" 18c 19c Input Variables 20 integer g_dipel ! Global array handle to the dipole matrix 21 logical lifetime !damping or not 22 integer g_dipel_i !currently defined in dimqm.fh, might change 23 double precision omega ! freq value 24c 25c Local variables 26 integer g_tmp1, g_tmp2, g_dcv, g_dim_temp 27 integer dims(3), chunk(3) 28 integer alo(3), ahi(3) 29 integer blo(2), bhi(2) 30 integer clo(3), chi(3) 31 integer xend 32 integer icmplx 33 34 double precision dx_r, dy_r, dz_r 35 double precision dx_i, dy_i, dz_i 36 integer l_dimxyz, k_dimxyz 37 integer l_muind, k_muind 38c integer l_muold, k_muold 39 double precision dsum 40 external dsum 41 integer i3, ivec, n 42 integer l_fld, k_fld 43 integer g_dim_r(2) 44 integer g_dim_i(2) 45 integer g_temp(2) 46 integer nvir, voff, xoff 47 integer ga_create_atom_blocked 48 external ga_create_atom_blocked 49 integer ipm 50 double precision pre_factor 51 character*(1) direction(3) 52 character*(1) pm(2) 53 logical seed_save 54 data direction /'x', 'y', 'z'/ 55 data pm /'+', '-'/ 56c 57 nvir = nmo - nclosed - nopen 58 voff = nclosed + nopen + 1 59 i3 = nDIM * 3 60 icmplx = 1 61 if(omega > ZERO) icmplx = 2 62 63c write(luout,*) "omega:", omega 64c write(luout,*) "icmplx:", icmplx 65 66 67c g_tmp1 = ga_create_atom_blocked(geom, basis, 'dim_lclfld: tmp1') 68c g_tmp2 = ga_create_atom_blocked(geom, basis, 'dim_lclfld: tmp2') 69c TODO: work out a way for FD without damping to run correctly, 70c although, the case itself makes little sense 71c if(icmplx > 1 ) then 72 if(lifetime) then 73 alo(1) = nbf 74 alo(2) = -1 75 alo(3) = -1 76 ahi(1) = nbf 77 ahi(2) = nbf 78 ahi(3) = 3 79 if (.not.nga_create(MT_DBL,3,ahi,'e-dipole-i',alo,g_dipel_i)) 80 $ call errquit('lclfld: nga_create failed g_dipel_i',0,GA_ERR) 81 call ga_zero(g_dipel_i) 82 end if 83 84 alo(2) = 1 85 ahi(2) = nbf 86 alo(3) = 1 87 ahi(3) = nbf 88 blo(1) = 1 89 bhi(1) = nbf 90 blo(2) = 1 91 bhi(2) = nbf 92 93 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:fld',l_fld,k_fld)) 94 $ call errquit('malloc dimrsp:fld failed',1,MA_ERR) 95c 96 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:muind', 97 $ l_muind,k_muind)) 98 $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 99 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:muold', 100 $ l_muold,k_muold)) 101 $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 102c 103 if(.not.ma_push_get(mt_dbl,i3,'dimrsp:xyz',l_dimxyz,k_dimxyz)) 104 $ call errquit('malloc dimrsp:xyz failed',1,MA_ERR) 105c 106 if(.not.rtdb_get(dimqm_rtdb,'dimpar:coords',mt_dbl,i3, 107 $ dbl_mb(k_dimxyz))) 108 $ call errquit('get dimpar:coords failed', 1, RTDB_ERR) 109 110c Loop over xyz dsrections 111 seed_save = dimqm_noseed 112 dimqm_noseed = .true. 113 do n = 1, 3 114c do ipm = 1, icmplx ! +/- direction if it's complex 115c 116c ============================= 117c Solve for dipoles 118c ============================= 119c 120c Zero arrays 121 call dfill(i3*icmplx, ZERO, dbl_mb(k_muind), 1) 122 call dfill(i3*icmplx, ZERO, dbl_mb(k_muold), 1) 123 call dfill(i3*icmplx, ZERO, dbl_mb(k_fld), 1) 124c 125c Ones for current cartesian direction 126 call dfill(nDIM, ONE, dbl_mb(k_fld+n-1), 3) 127c 128c Calculate induced dipoles 129 call dimqm_f2d(dimqm_rtdb, dbl_mb(k_fld), 130 $ dbl_mb(k_muind), dbl_mb(k_muold), 131 $ dbl_mb(k_dimxyz), icmplx, 's', ' ',.false.) 132c 133 if(icmplx.eq.2 .and. lifetime) then 134 if(.not.rtdb_put(dimqm_rtdb, 135 $ 'dimqm:muind_'//direction(n)//'_r'//pm(1), 136 $ mt_dbl, i3, dbl_mb(k_muind))) 137 $ call errquit('put dimqm:muind_r failed',1,RTDB_ERR) 138 if(.not.rtdb_put(dimqm_rtdb, 139 $ 'dimqm:muind_'//direction(n)//'_i'//pm(1), 140 $ mt_dbl, i3, dbl_mb(k_muind+i3))) 141 $ call errquit('put dimqm:muind_i failed',1,RTDB_ERR) 142 else 143 if(.not.rtdb_put(dimqm_rtdb, 144 $ 'dimqm:muind_'//direction(n), 145 $ mt_dbl, i3, dbl_mb(k_muind))) 146 $ call errquit('put dimqm:muind_r failed',1,RTDB_ERR) 147 end if 148c end do 149 end do ! ivec = 1, 3 150c 151 if (.not. ma_chop_stack(l_fld)) call errquit('lclfld: MA?', 0, 152 $ MA_ERR) 153c 154c ==================================================== 155c Solve for DIM potential 156c ==================================================== 157c 158 dims(1) = 3 159 dims(2) = nbf 160 dims(3) = nbf 161 chunk(1) = dims(1) 162 chunk(2) = -1 163 chunk(3) = -1 164c 165 if((icmplx.eq.2) .and. lifetime) then 166c Real + 167 if (.not. nga_create (MT_DBL, 3, dims, 'lclfld:dim_r+',chunk, 168 & g_dim_r(1))) 169 & call errquit('lclfld: could not allocate g_dim_r+',1,GA_ERR) 170 call ga_zero(g_dim_r(1)) 171 call fock_dim(geom, nbf, basis, 3, g_dim_r(1), 1, 1) 172 call ga_symmetrize(g_dim_r(1)) 173c! Real - 174 ! if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_r-',chunk, 175 ! & g_dim_r(2))) 176 ! & call errquit('addop: could not allocate g_dim_r-',1,GA_ERR) 177 ! call ga_zero(g_dim_r(2)) 178 ! call fock_dim(geom, nbf, basis, 3, g_dim_r(2), 2, 1) 179 ! call ga_antisymmetrize(g_dim_r(2)) 180c Imaginary + 181 if (.not. nga_create (MT_DBL, 3, dims, 'lclfld:dim_i+',chunk, 182 & g_dim_i(1))) 183 & call errquit('lclfld: could not allocate g_dim_i+',1,GA_ERR) 184 call ga_zero(g_dim_i(1)) 185 call fock_dim(geom, nbf, basis, 3, g_dim_i(1), 1, 2) 186 call ga_symmetrize(g_dim_i(1)) 187c! Imaginary - 188 ! if (.not. nga_create (MT_DBL, 3, dims, 'addop:dim_i-',chunk, 189 ! & g_dim_i(2))) 190 ! & call errquit('addop: could not allocate g_dim_i-',1,GA_ERR) 191 ! call ga_zero(g_dim_i(2)) 192 ! call fock_dim(geom, nbf, basis, 3, g_dim_i(2), 2, 2) 193 ! call ga_antisymmetrize(g_dim_i(2)) 194 ! 195 else 196 if (.not. nga_create (MT_DBL, 3, dims, 'lclfld:dim_r',chunk, 197 & g_dim_r(1))) 198 & call errquit('lclfld: could not allocate g_dim_r',1,GA_ERR) 199 call ga_zero(g_dim_r) 200 call fock_dim(geom, nbf, basis, 3, g_dim_r(1), 0, 1) 201 call ga_symmetrize(g_dim_r(1)) 202c! ====================================== 203c! Undo the symmetrization to recover +/- 204c! ====================================== 205c! 206 ! dims(1) = nbf 207 ! dims(2) = nbf 208 ! chunk(1) = dims(1) 209 ! chunk(2) = -1 210 ! 211 ! if (.not.nga_create(MT_DBL,2,dims,'gtemp',chunk, 212 ! & g_temp(1))) call 213 ! & errquit('dim_addop: nga_create failed gtemp', 214 ! & 0,GA_ERR) 215 ! call ga_zero(g_temp(1)) 216 ! if (.not.nga_create(MT_DBL,2,dims,'gtemp',chunk, 217 ! & g_temp(2))) call 218 ! & errquit('dim_addop: nga_create failed gtemp', 219 ! & 0,GA_ERR) 220 ! call ga_zero(g_temp(2)) 221 ! 222 ! do ivec = 1, 3 223 ! alo(1) = ivec 224 ! ahi(1) = ivec 225c! ************ 226c! Real portion 227c! ************ 228 ! call nga_copy_patch ('N',g_dim_r(1),alo,ahi,g_temp(1),blo,bhi) 229 ! call nga_copy_patch ('N',g_dim_r(2),alo,ahi,g_temp(2),blo,bhi) 230!c 231!c it might be necessary to use 0.5 here instead of 1.0 232!c (note: that turned out NOT to be the case after some testing) 233! pre_factor = 1.0d0 234 ! call ga_sync() 235 ! call nga_add_patch (pre_factor, g_temp(1), blo, bhi, 236 ! & pre_factor, g_temp(2), blo, bhi, 237 ! & g_dim_r(1), alo, ahi) 238 ! call nga_add_patch (pre_factor, g_temp(1), blo, bhi, 239 ! & -pre_factor, g_temp(2), blo, bhi, 240 ! & g_dim_r(2), alo, ahi) 241c! ***************** 242c! Imaginary portion 243c! ***************** 244 ! call nga_copy_patch ('N',g_dim_i(1),alo,ahi,g_temp(1),blo,bhi) 245 ! call nga_copy_patch ('N',g_dim_i(2),alo,ahi,g_temp(2),blo,bhi) 246!c 247c! it might be necessary to use 0.5 here instead of 1.0 248c! (note: that turned out NOT to be the case after some testing) 249 ! pre_factor = 1.0d0 250 ! call ga_sync() 251c! real perturbation: 252 ! call nga_add_patch (pre_factor, g_temp(1), blo, bhi, 253 ! & pre_factor, g_temp(2), blo, bhi, 254 ! & g_dim_i(1), alo, ahi) 255 ! call nga_add_patch (pre_factor, g_temp(1), blo, bhi, 256 ! & -pre_factor, g_temp(2), blo, bhi, 257 ! & g_dim_i(2), alo, ahi) 258 ! enddo ! ivec = 1,nvec 259 end if 260c 261c ======================================== 262c Add DIM local field to the dipole matrix 263c ======================================== 264c start: some debug stuff ------------------- 265c call ga_print(g_dim_r(1)) 266cc call ga_print(g_dim_r(2)) 267c call ga_print(g_dipel) 268c 269c if (icmplx > 1) then 270c call ga_print(g_dim_i(1)) 271cc call ga_print(g_dim_i(2)) 272c endif 273c write(luout,*)'end of g_dim print' 274c end: some debug stuff -------------------- 275 276c 277 xoff = 1 278 xend = nvir*nclosed 279c call ga_print(g_movecs) 280 call ga_sync() 281 clo(1) = 1 282 chi(1) = nbf 283 clo(2) = 1 284 chi(2) = nbf 285!jbecca START: commenting this stuff out just to test 286! do ivec = 1, 3 ! Loop over perturbations 287! alo(1) = ivec 288! ahi(1) = ivec 289! clo(3) = ivec 290! chi(3) = ivec 291! call nga_copy_patch('N',g_dim_r(1),alo,ahi,g_dcv,blo,bhi) 292! call nga_add_patch(ONE, g_dipel, clo, chi, 293! $ ONE, g_dcv, blo, bhi, 294! $ g_dipel, clo, chi) 295! if(icmplx > 1) then 296! call nga_copy_patch('N',g_dim_i(1),alo,ahi,g_dcv,blo,bhi) 297! call nga_add_patch(ONE, g_dipel_i, clo, chi, 298! $ ONE, g_dcv, blo, bhi, 299! $ g_dipel_i, clo, chi) 300! 301!c call ga_zero(g_dipel_i) 302! end if 303! end do !ivec = 1, 3 304!jbecca END 305!jbecca START: It seems that ga_create_atom_blocked does not work 306! for this case. Making this by hand now. 307 if (.not. ga_create(MT_DBL,nbf,nbf,'dim_temp_pert',chunk(1), 308 $ chunk(2),g_dim_temp)) call errquit('lclfld: 309 $ nga_create failed g_dim_temp',0,GA_ERR) 310 311 do ivec = 1, 3 ! loop over pert 312 alo(1) = ivec 313 ahi(1) = ivec 314 clo(3) = ivec 315 chi(3) = ivec 316 call ga_zero(g_dim_temp) 317 call nga_copy_patch('N',g_dim_r(1),alo,ahi,g_dim_temp,blo,bhi) 318 call nga_add_patch(ONE, g_dipel, clo, chi, 319 $ ONE, g_dim_temp, blo, bhi, 320 $ g_dipel, clo, chi) 321 322c if (icmplx > 1 ) then 323 if (lifetime) then 324 call ga_zero(g_dim_temp) 325 call nga_copy_patch('N',g_dim_i(1),alo,ahi,g_dim_temp, 326 $ blo,bhi) 327 call nga_add_patch(ONE, g_dipel_i, clo, chi, 328 $ ONE, g_dim_temp, blo, bhi, 329 $ g_dipel_i, clo, chi) 330 endif 331 enddo !ivec 332 333 if (.not. ga_destroy(g_dim_temp)) 334 & call errquit('lclfld: g_dim_temp destroy?',0,GA_ERR) 335!jbecca END 336c call ga_print(g_dim_r(1)) 337c ======== 338c Clean up 339c ======= 340 if (.not. ga_destroy(g_dim_r(1))) 341 & call errquit('lclfld: g_dim_r destroy?',0, GA_ERR) 342 if(icmplx > 1 .and. lifetime) then 343 ! if (.not. ga_destroy(g_dim_r(2))) call errquit('addop: GA?',0, 344 ! & GA_ERR) 345 if (.not. ga_destroy(g_dim_i(1))) 346 & call errquit('lclfld:g_dim_i destroy?',0,GA_ERR) 347 ! if (.not. ga_destroy(g_dim_i(2))) call errquit('addop: GA?',0, 348 ! & GA_ERR) 349 ! if (.not. ga_destroy(g_temp(1))) call errquit('addop: GA?',0, 350 ! & GA_ERR) 351 ! if (.not. ga_destroy(g_temp(2))) call errquit('addop: GA?',0, 352 ! & GA_ERR) 353 end if 354!jbecca START 355 dimqm_noseed = seed_save 356c call ga_print(g_dipel_i) 357 end subroutine dimqm_lclfld 358 359 360 361 subroutine dimqm_addlclfld(g_dip, omega) 362c----------------------------------------------------------------------- 363c This subroutine is used to create and add in the local field 364c operator to the U (or A) matrix before solving the linear equations 365c in response routines. This ensures that the QM system's response is 366c also dependent on any local fields that exist. 367c 368c Called from: nothing right now 369c 370c Calls: subroutine dimqm_f2d 371c 372c Author: Jeff Becca, jbb5516@psu.edu, 2017 373c----------------------------------------------------------------------- 374 375 implicit none 376#include "errquit.fh" 377#include "stdio.fh" 378#include "rtdb.fh" 379#include "mafdecls.fh" 380#include "global.fh" 381#include "dimqm_constants.fh" 382#include "dimqm.fh" 383#include "geom.fh" 384#include "crohf.fh" 385#include "cscf.fh" 386 387c--------------------------- Input variables --------------------------- 388 integer g_dip !H10, used to initialize rhs 389 integer omega !frequency 390 391c--------------------------- Local variables --------------------------- 392 integer icmplx, i3, ivec, n 393 integer alo(3), ahi(3), dims(3), chunk(3) 394 integer l_fld, k_fld, l_muind, k_muind, l_dimxyz, k_dimxyz 395 logical seed_save 396 character*(1) direction(3) 397 character*(1) pm(2) 398 data direction /'x', 'y', 'z'/ 399 data pm /'+', '-'/ 400c---------------------------- Routine START ---------------------------- 401c TODO: Update this for the case of complex local fields. Currently 402c only caring about the real part of the local field. 403 i3 = nDIM * 3 404 icmplx = 1 405 if (omega > ZERO) icmplx = 2 406 407c get some DIM stuff needed for setting up field 408 409 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:fld',l_fld,k_fld)) 410 $ call errquit('malloc dimrsp:fld failed',1,MA_ERR) 411c 412 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:muind', 413 $ l_muind,k_muind)) 414 $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 415 if(.not.ma_push_get(mt_dbl,i3*icmplx,'dimrsp:muold', 416 $ l_muold,k_muold)) 417 $ call errquit('malloc dimrsp:muind failed',1,MA_ERR) 418c 419 if(.not.ma_push_get(mt_dbl,i3,'dimrsp:xyz',l_dimxyz,k_dimxyz)) 420 $ call errquit('malloc dimrsp:xyz failed',1,MA_ERR) 421c 422 if(.not.rtdb_get(dimqm_rtdb,'dimpar:coords',mt_dbl,i3, 423 $ dbl_mb(k_dimxyz))) 424 $ call errquit('get dimpar:coords failed', 1, RTDB_ERR) 425 426 seed_save = dimqm_noseed 427 dimqm_noseed = .true. 428 429 do n = 1, 3 430 431c zero arrays 432 call dfill(i3*icmplx, ZERO, dbl_mb(k_muind), 1) 433 call dfill(i3*icmplx, ZERO, dbl_mb(k_muold), 1) 434 call dfill(i3*icmplx, ZERO, dbl_mb(k_fld), 1) 435 436c Create unit field in current cartesian direction for field 437 call dfill(nDIM, ONE, dbl_mb(k_fld+n-1), 3) 438 439c Calculate the induced dipoles from the unit field 440 call dimqm_f2d(dimqm_rtdb, dbl_mb(k_fld), 441 $ dbl_mb(k_muind), dbl_mb(k_muold), 442 $ dbl_mb(k_dimxyz), icmplx, 's', ' ',.false.) 443 444c Store dipoles 445 if (icmplx.eq.1) then 446 if(.not.rtdb_put(dimqm_rtdb, 447 $ 'dimqm:muind_'//direction(n), 448 $ mt_dbl, i3, dbl_mb(k_muind))) 449 $ call errquit('put dimqm:muindr_r failed',1,RTDB_ERR) 450 else 451 if(.not.rtdb_put(dimqm_rtdb, 452 $ 'dimqm:muind_'//direction(n)//'_r'//pm(1), 453 $ mt_dbl, i3, dbl_mb(k_muind))) 454 $ call errquit('put dimqm:muind_r failed',1,RTDB_ERR) 455 if(.not.rtdb_put(dimqm_rtdb, 456 $ 'dimqm:muind_'//direction(n)//'_i'//pm(1), 457 $ mt_dbl, i3, dbl_mb(k_muind+i3))) 458 $ call errquit('put dimqm:muind_i failed',1,RTDB_ERR) 459 end if 460 end do ! ivec = 1, 3 461 462 if (.not. ma_chop_stack(l_fld)) call errquit('lclfld: MA?', 0, 463 $ MA_ERR) 464c 465c ==================================================== 466c Solve for DIM potential 467c ==================================================== 468c 469 dims(1) = 3 470 dims(2) = nbf 471 dims(3) = nbf 472 chunk(1) = dims(1) 473 chunk(2) = -1 474 chunk(3) = -1 475 end subroutine dimqm_addlclfld 476