1C 2C 3C rt_tddft_spatial_potential.F 4C 5C Compute spatial potential projected onto nbf_ao x nbf_ao matrix. 6C Currently this reads in complex absorbing potential parameters 7C from the rtdb. 8C 9C Note, that for spin-orbit calcns this will be called from 10C the openshell branch (thus all dimensions here are nbf_ao *not* 11C ns_ao). 12C 13C Outputs a real-valued nao x nao GA. 14C 15C 16 subroutine rt_tddft_spatial_potential (params, nao, g_v) 17 implicit none 18 19#include "errquit.fh" 20#include "mafdecls.fh" 21#include "stdio.fh" 22#include "global.fh" 23#include "msgids.fh" 24#include "util.fh" 25#include "cdft.fh" 26#include "geomP.fh" 27#include "geom.fh" 28#include "bas.fh" 29#include "rtdb.fh" 30#include "matutils.fh" 31#include "rt_tddft.fh" 32 33 34C == Inputs == 35 type(rt_params_t), intent(in) :: params 36 integer, intent(in) :: nao 37 38 39C == Outputs == 40 integer, intent(in) :: g_v !real-valued AO basis potential 41 42 43C == Parameters == 44 character(len=*), parameter::pname="rt_tddft_spatial_potential: " 45 46 47C == Variables == 48 integer :: rtdb 49 integer :: me, np, n, m, m0, istart, iend, id 50 integer :: nq(3), nqtot, nqlocal 51 integer :: pid 52 double precision :: rval, rvalfull 53 double precision :: cap_gammamax, cap_width 54 double precision :: r 55 double precision :: r1, r2 ,r3 56 double precision :: qmax_au(3), qmax_ang(3) 57 double precision :: qmin_au(3), qmin_ang(3) 58 double precision :: dq(3) 59 double precision :: elapsed 60 double precision :: uniform_wght , val 61 double precision :: x, y, z, fx, fy, fz 62 integer :: lxyz, ixyz, lpot, ipot, iwgt, lwgt 63 integer :: iq, ix, iy, iz, i, j 64 65 integer, parameter :: pot_msg = 29348, ovl_msg = 89723 66 67C xc_eval_basis() stuff 68 integer :: ncontrset 69 integer :: lbas_cent_info, ibas_cent_info 70 integer :: lbas_cset_info, ibas_cset_info 71 integer :: ldocset, idocset 72 integer :: lniz, iniz 73 integer :: nxyz_atom, lxyz_atom, ixyz_atom 74 integer :: lcharge, icharge, ltags, itags 75 integer :: l_rchi_atom, i_rchi_atom 76 integer :: l_rq, i_rq 77 integer :: lchi_ao, ichi_ao 78 integer :: ld1chi_ao, id1chi_ao !gradients of basis functions 79 80 character :: ctag 81 character*16 :: tag 82 83 integer :: iovl_ao, lovl_ao, ipot_ao, lpot_ao 84 integer :: ia, ic 85 86 double precision :: atom_x(ncenters), atom_y(ncenters), 87 $ atom_z(ncenters) 88 double precision :: atom_rstart(ncenters), atom_rend(ncenters) 89 90 logical :: lreduced_memory 91 logical :: lcap_print 92 93 character*64 :: tagstr 94 95 rtdb = params%rtdb 96 97 if (params%prof) call prof_start (elapsed) 98 99 me = ga_nodeid() 100 np = ga_nnodes() 101 102C 103C Rectangular grid 104C 105C (read params from rtdb) 106 if (.not. rtdb_get (rtdb, "rt_tddft:cap_potmax", 107 $ mt_dbl, 1, cap_gammamax)) call errquit (pname// 108 $ "failed to read potmax from rtdb", 0, RTDB_ERR) ! called cap_gammamax 109 110 if (.not. rtdb_get (rtdb, "rt_tddft:cap_nq", 111 $ mt_int, 3, nq)) call errquit (pname// 112 $ "failed to read nq from rtdb", 0, RTDB_ERR) 113 114 if (.not. rtdb_get (rtdb, "rt_tddft:cap_qmin", 115 $ mt_dbl, 3, qmin_au)) call errquit (pname// 116 $ "failed to read qmin from rtdb", 0, RTDB_ERR) !note: stored internally in au 117 118 if (.not. rtdb_get (rtdb, "rt_tddft:cap_qmax", 119 $ mt_dbl, 3, qmax_au)) call errquit (pname// 120 $ "failed to read qmax from rtdb", 0, RTDB_ERR) 121 122 123 if (.not.rtdb_get(params%rtdb, "rt_tddft:cap_reduce_memory", 124 $ mt_log, 1, lreduced_memory)) 125 $ lreduced_memory = .false. 126 127 if (lreduced_memory) 128 $ call errquit (pname// 129 $ "reduced memory method not available yet", 0, 0) 130 131 132C Default: do not print CAP 133 if (.not. rtdb_get (rtdb, "rt_tddft:cap_print", 134 $ mt_log, 1, lcap_print)) lcap_print = .false. 135 136 137 138C Convert lengths to angstroms (for outputs) and compute grid spacing. 139 do id = 1, 3 140 qmin_ang(id) = qmin_au(id) * au2ang 141 qmax_ang(id) = qmax_au(id) * au2ang 142 dq(id) = (qmax_au(id) - qmin_au(id)) / dble(nq(id)) 143 enddo 144 nqtot = nq(1)*nq(2)*nq(3) 145 146 147C 148C Now divide grid up over processors (broken in z-direction) 149C 150C m is the number of z-grid points on this proc 151C m0 is the number of z-grid points on processor 0 (has extra) 152C 153 pid = 3 !parallel over the z grid points !XXX SIMPLY CHANGING THIS MIGHT NOT WORK 154 n = nq(pid) 155 156 m0 = n/np + mod (n, np) 157 158C call rt_tddft_calc_1d_partitioning (n, m, istart, iend) 159 160 if (me.eq.0) then 161 m = m0 162 else 163 m = n/np 164 endif 165 166 if (me.eq.0) then 167 istart = 1 168 iend = m0 169 else 170 istart = m0 + 1 + (me-1)*m 171 iend = istart + m - 1 172 endif 173 174 nqlocal = nq(1) * nq(2) * m 175 176 if (me.eq.0) then 177 write (luout, *) "" 178 call util_print_centered (luout, 179 $ "Spatial grid-based complex absorbing potential", 180 $ 40,.true.) 181 write (luout, *) "" 182 183 write (luout, "(1x,a,i0,a,i0,a,i0,a,i0)") 184 $ "Spatial grid points : ", 185 $ nq(1), " x ", nq(2), " x ", nq(3), " = ", nqtot 186 write (luout, *) "" 187 188 write (luout, "(1x,a,3f12.4)") " : ", 189 $ qmin_ang(1), qmax_ang(1), dq(1) 190 write (luout, "(1x,a,3f12.4)") "Grid geometry [A] : ", 191 $ qmin_ang(2), qmax_ang(2), dq(2) 192 write (luout, "(1x,a,3f12.4)") " : ", 193 $ qmin_ang(3), qmax_ang(3), dq(3) 194 write (luout, *) "" 195 196 write (luout, "(1x,a,3f12.4)") " : ", 197 $ qmin_au(1), qmax_au(1), dq(1) 198 write (luout, "(1x,a,3f12.4)") "Grid geometry [au] : ", 199 $ qmin_au(2), qmax_au(2), dq(2) 200 write (luout, "(1x,a,3f12.4)") " : ", 201 $ qmin_au(3), qmax_au(3), dq(3) 202 write (luout, *) "" 203 204 call util_flush (luout) 205 endif 206 call ga_sync () 207 208 write (luout, "(a,i6,a,i0,a)") 209 $ "Proc ", me, ": ", nqlocal, " grid points" 210 211 212C 213C Allocate local parts of the grid. 214C 215 if (.not. ma_push_get (mt_dbl, 3*nqlocal, "grid", lxyz, ixyz)) 216 $ call errquit (pname//"cannot alloc grid", 0, MA_ERR) 217 218 if (.not. ma_push_get (mt_dbl, nqlocal, "potential", lpot, ipot)) 219 $ call errquit (pname//"cannot alloc potential", 0, MA_ERR) 220 221 if (.not. ma_push_get (mt_dbl, nqlocal, "weight", lwgt, iwgt)) 222 $ call errquit (pname//"cannot alloc weight", 0, MA_ERR) 223 224 225C (first misc basis set info req'd) 226 if (.not.bas_numcont(ao_bas_han, ncontrset)) 227 $ call errquit(pname//"bas_numcont failed",0, BASIS_ERR) 228 229 if (.not.ma_push_get(mt_int, 3*ncenters, "bas_cent_info", 230 & lbas_cent_info, ibas_cent_info)) 231 & call errquit(pname//"cannot allocate bas_cent_info",0, 232 & MA_ERR) 233 234 if (.not.ma_push_get(mt_int, 6*ncontrset, 'bas_cset_info', 235 & lbas_cset_info, ibas_cset_info)) 236 & call errquit(pname//"cannot allocate bas_cset_info",0, 237 & MA_ERR) 238 239 call xc_make_basis_info(ao_bas_han, int_mb(ibas_cent_info), 240 & int_mb(ibas_cset_info), ncenters) 241 242 if (.not.ma_push_get(mt_log, ncontrset, 'docset', 243 & ldocset, idocset)) 244 & call errquit(pname//'cannot allocate ccdocset', 245 . ncontrset, MA_ERR) 246 do i=1,ncontrset 247 log_mb(idocset+i-1)=.true. 248 enddo 249 250 if(.not.ma_push_get(MT_int, ncenters, 'iniz', 251 & lniz, iniz)) 252 & call errquit(pname//"iniz",0, MA_ERR) 253 do i= 1, ncenters 254 int_mb(iniz+i-1)=1 255 enddo 256 257 nxyz_atom = 3*ncenters 258 if (.not.ma_push_get(mt_dbl,nxyz_atom,'xyz_atom', 259 $ lxyz_atom,ixyz_atom)) 260 & call errquit(pname//'cannot allocate xyz_atom',0, MA_ERR) 261 262 if (.not.ma_push_get(mt_dbl,ncenters,'charge',lcharge,icharge)) 263 & call errquit(pname//'cannot allocate charge',0, MA_ERR) 264 265 if (.not.ma_push_get(mt_Byte,ncenters*16,'tags',ltags,itags)) 266 & call errquit(pname//'cannot allocate tags',0, MA_ERR) 267 268 if (.not. geom_cart_get(geom, ncenters, byte_mb(itags), 269 & dbl_mb(ixyz_atom), dbl_mb(icharge))) 270 & call errquit(pname//'geom_cart_get failed', 0, GEOM_ERR) 271 272 273C (get geometry information) 274 call ga_sync() 275 call util_flush(luout) 276 if (me.eq.0) then 277 write(luout,*)"" 278 call util_print_centered (luout, 279 $ "Atom-centered CAP geometry in atomic units", 280 $ 40,.true.) 281 282 write(luout,*)"" 283 write(luout,*) 284 $ "Atom tag ", 285 $ "x y z rstart rend" 286 write(luout,*) "-------------------------------------------"// 287 $ "--------------------------------------------" 288 endif 289 290 do ic = 1, ncenters 291 atom_x(ic) = dbl_mb(ixyz_atom + 3*(ic-1) + 0) 292 atom_y(ic) = dbl_mb(ixyz_atom + 3*(ic-1) + 1) 293 atom_z(ic) = dbl_mb(ixyz_atom + 3*(ic-1) + 2) 294 295 tag = "" 296 do ia = 1, 16 297 ctag = byte_mb(itags + 16*(ic-1) + ia - 1) 298 tag = trim(tag) // ctag 299 enddo 300 301C 302C Rstart and Rend for each atom type. 303C 304C if (trim(tag) .ne. "bqH") then 305 tagstr = "rt_tddft:cap_rstart_"//trim(tag) ! e.g. "rt_tddft:cap_rstart_O" 306 if (.not. rtdb_get (rtdb, trim(tagstr), mt_dbl, 307 $ 1, atom_rstart(ic))) 308 $ call errquit (pname// 309 $ "failed to read rstart for "//trim(tag), 0, RTDB_ERR) 310 311 tagstr = "rt_tddft:cap_rend_"//trim(tag) ! e.g. "rt_tddft:cap_rend_O" 312 if (.not. rtdb_get (rtdb, trim(tagstr), mt_dbl, 313 $ 1, atom_rend(ic))) 314 $ call errquit (pname// 315 $ "failed to read rend for "//trim(tag), 0, RTDB_ERR) 316 317C Convert Rstart and Rend from angstroms to atomic units 318C Not needed, since we now store internally in au. 319c$$$ atom_rstart(ic) = atom_rstart(ic) * 1.889725989d0 320c$$$ atom_rend(ic) = atom_rend(ic) * 1.889725989d0 321 322 if (me.eq.0) then 323 write(luout,"(1x,i6,4x,a,5f12.4)") 324 $ ic, tag, atom_x(ic), atom_y(ic), atom_z(ic), 325 $ atom_rstart(ic), atom_rend(ic) 326 endif 327 enddo 328 329 if (me.eq.0) then 330 write(luout,*)"" 331 endif 332 333 334C 335C Radial CAP from each atom 336C 337 iq = 0 338 uniform_wght = dq(1) * dq(2) * dq(3) 339 340 do iz = istart, iend !note this only loops over slice local to this process 341 z = qmin_au(3) + (iz-1)*dq(3) 342 343 do iy = 1, nq(2) 344 y = qmin_au(2) + (iy-1)*dq(2) 345 346 do ix = 1, nq(1) 347 x = qmin_au(1) + (ix-1)*dq(1) 348 349C (product of radial hann functions centered on each atom) 350 rvalfull = 1d0 351 do ic = 1, ncenters 352 353 tag = "" 354 do ia = 1, 16 355 ctag = byte_mb(itags + 16*(ic-1) + ia - 1) 356 tag = trim(tag) // ctag 357 enddo 358 359 cap_width = atom_rend(ic) - atom_rstart(ic) 360 361 r = sqrt((x-atom_x(ic))**2 + 362 $ (y-atom_y(ic))**2 + (z-atom_z(ic))**2) 363 364 if (r .lt. atom_rstart(ic)) then 365 rval = 0d0 366 elseif ((r.gt.atom_rstart(ic)) .and. 367 $ (r.le.atom_rend(ic))) then 368 rval = 1d0 * 369 $ sin((dpi/2d0*(r-atom_rstart(ic))) 370 $ / (cap_width))**2 371 else 372 rval = 1d0 373 endif 374 375 rvalfull = rvalfull*rval 376 enddo 377 378 379C Normalize cap such that max value is gamma_max 380 rvalfull = cap_gammamax * rvalfull 381 382 383C (store in local MA) 384 iq = iq + 1 385 dbl_mb(ixyz + 3*(iq-1)+0) = x 386 dbl_mb(ixyz + 3*(iq-1)+1) = y 387 dbl_mb(ixyz + 3*(iq-1)+2) = z 388 dbl_mb(ipot + iq - 1) = rvalfull 389 dbl_mb(iwgt + iq - 1) = uniform_wght 390 391C 392C Print CAP stdout. Ugly secret option for now. 393C 394 if (lcap_print) then 395 if ((abs(x) .lt. dq(1))) then 396 write(luout,"(3es22.12e3,a)") 397 $ y, z, rvalfull, " # CAP slice x=0" 398 endif 399 call util_flush(luout) 400 401 if ((abs(y) .lt. dq(2))) then 402 write(luout,"(3es22.12e3,a)") 403 $ x, z, rvalfull, " # CAP slice y=0" 404 endif 405 call util_flush(luout) 406 407 if ((abs(z) .lt. dq(3))) then 408 write(luout,"(3es22.12e3,a)") 409 $ x, y, rvalfull, " # CAP slice z=0" 410 endif 411 call util_flush(luout) 412 endif 413 414 enddo 415 enddo 416 enddo 417 418 419 if (params%prof) call prof_end (elapsed, 420 $ "Computing spatial potential on the grid") 421 422 423C 424C Evalulate this potential over the AO basis for the slice local to 425C this process. 426C 427C (modified from dft_frozemb.F) 428 429 if (params%prof) call prof_start (elapsed) 430 431C (now compute basis functions over the grid) 432 if(.not.ma_push_get(mt_dbl, ncenters, 'rchi_atom', 433 & l_rchi_atom,i_rchi_atom)) 434 & call errquit(pname//"failed to allocate rchi_atom, "// 435 $ "not enough memory?",0, MA_ERR) 436 437 if(.not.ma_push_get(mt_dbl, nqlocal*ncenters, 'rq', 438 & l_rq,i_rq)) 439 & call errquit(pname//"failed to allocate rq, "// 440 $ "not enough memory?",0, MA_ERR) 441 442 if (.not.ma_push_get(mt_dbl, nqlocal*nao, 443 & 'chi_ao', lchi_ao, ichi_ao)) 444 & call errquit(pname//"failed to allocate chi_ao, "// 445 $ "not enough memory?",0, MA_ERR) 446 447 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 448 & dbl_mb(ixyz), dbl_mb(ixyz_atom), nqlocal, ncenters) 449 450C 451C Compute the basis functions on the grid. 452C 453 call xc_eval_basis(ao_bas_han, 0, dbl_mb(ichi_ao), 454 & 0d0, 0d0, 0d0, dbl_mb(i_rq), 455 & dbl_mb(ixyz), dbl_mb(ixyz_atom), nqlocal, ncenters, 456 & int_mb(iniz), log_mb(idocset), 457 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 458 459 460 if (params%prof) call prof_end (elapsed, 461 $ "Computing basis functions over the grid") 462 463 464C 465C Compute the basis functions on the grid. This version also 466C computes gradients of the basis functions on the grid. Note that 467C maxder = 1. I do not currently use this, but I've left it here 468C for future use. If using this, comment out above block. 469C 470c$$$ if (.not.ma_push_get(mt_dbl, nqlocal*nao*3 471c$$$ & 'd1chi_ao', ld1chi_ao, id1chi_ao)) 472c$$$ & call errquit(pname//'d1chi',0, MA_ERR) 473c$$$ 474c$$$ call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 475c$$$ & dbl_mb(id1chi_ao), 0d0, 0d0, dbl_mb(i_rq), 476c$$$ & dbl_mb(ixyz), dbl_mb(ixyz_atom), nqlocal, ncenters, 477c$$$ & int_mb(iniz), log_mb(idocset), 478c$$$ & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 479 480 481 482C 483C Now integrate over the basis to compute the AO potential and 484C overlap (to diagnose the quality of the spatial grid) 485C 486 487C Each process has its own nao x nao pot_ao and ovl_ao, which we 488C populate independently (each proc gets a subsection of the grid 489C points), then accumulate at the end. 490 491 if (params%prof) call prof_start (elapsed) 492 493 if (.not.ma_push_get(mt_dbl, nao * nao, 494 & "pot_ao", lpot_ao, ipot_ao)) 495 & call errquit(pname//"failed to allocate pot_ao, "// 496 $ "not enough memory?",0, MA_ERR) 497 498 if (.not.ma_push_get(mt_dbl, nao * nao, 499 & "ovl_ao", lovl_ao, iovl_ao)) 500 & call errquit(pname//"failed to allocate ovl_ao, "// 501 $ "not enough memory?",0, MA_ERR) 502 503 504C Turned on by adonay 505C (print projection of AO basis on grid to file) 506c$$$ call rt_tddft_spatial_potential_dump (rtdb, nao, 507c$$$ $ nqlocal, dbl_mb(ixyz), dbl_mb(ichi_ao), dbl_mb(id1chi_ao)) 508C 509C above turned on by adonay 510C added by adonay 06/29/17 511c$$$ call rt_tddft_spatial_potential_read (rtdb, nao, 512c$$$ $ nqlocal, dbl_mb(ixyz), dbl_mb(ichi_ao), dbl_mb(id1chi_ao)) 513C 514 515C (integrate over AO basis for this slice of the grid) 516 call rt_tddft_spatial_potential_integrate (rtdb, nao, 517 $ nqlocal, dbl_mb(ichi_ao), dbl_mb(ipot), dbl_mb(iwgt), 518 $ dbl_mb(ipot_ao), dbl_mb(iovl_ao)) 519 520 521 if (params%prof) call prof_end (elapsed, 522 $ "Integrating over the basis functions") 523 524 525C (accumulate) 526 call ga_sync () 527 call ga_dgop (pot_msg, dbl_mb(ipot_ao), 528 $ nao*nao, "+") 529 call ga_dgop (ovl_msg, dbl_mb(iovl_ao), 530 $ nao*nao, "+") 531 call ga_sync () 532 533C (print total potential and overlap to screen) 534 call rt_tddft_spatial_potential_print (rtdb, nao, 535 $ nqtot, dbl_mb(ipot_ao), dbl_mb(iovl_ao)) 536 537C 538C Real-valued output 539C 540 if (me.eq.0) 541 $ call ga_put (g_v, 1, nao, 1, nao, dbl_mb(ipot_ao), nao) 542 543C 544C Clean up 545C 546 if (.not. ma_chop_stack (lxyz)) 547 $ call errquit (pname//"chop failed", 0, MA_ERR) 548 549 550 end subroutine 551 552 553 554C============================================================ 555C 556C Compute integral: < mu(g) | V(g) | nu(g) > 557C 558C Ripped from dft_frozemb.F : acc_fock() 559C 560 subroutine rt_tddft_spatial_potential_integrate (rtdb, nao, nq, 561 $ chi_ao, pot, wgt, pot_ao, ovl_ao) 562 implicit none 563 564#include "errquit.fh" 565#include "mafdecls.fh" 566#include "stdio.fh" 567#include "global.fh" 568#include "msgids.fh" 569#include "util.fh" 570#include "cdft.fh" 571#include "geomP.fh" 572#include "geom.fh" 573#include "bas.fh" 574#include "matutils.fh" 575#include "rt_tddft.fh" 576 577 integer, intent(in) :: rtdb 578C type(rt_params_t), intent(in) :: params 579 integer, intent(in) :: nao, nq 580 double precision, intent(in) :: chi_ao(nq, nao), wgt(nq), pot(nq) 581 double precision, intent(inout) :: pot_ao(nao,nao) 582 double precision, intent(inout) :: ovl_ao(nao,nao) 583 584 character(len=*), parameter :: pname = 585 $ "rt_tddft_spatial_potential_integrate: " 586 587 integer :: i, j, k 588 integer :: me 589 590 me = ga_nodeid() 591 592 593 do i = 1, nao 594 do j = 1, nao 595 pot_ao(i,j) = 0d0 596 ovl_ao(i,j) = 0d0 597 do k = 1, nq 598 pot_ao(i,j) = pot_ao(i,j) + 599 $ chi_ao(k,i)*wgt(k)*chi_ao(k,j)*pot(k) 600 ovl_ao(i,j) = ovl_ao(i,j) + 601 $ chi_ao(k,i)*wgt(k)*chi_ao(k,j) 602 enddo 603 enddo 604 enddo 605 606 end subroutine 607 608 609 610C============================================================ 611C 612C Dump grid and projection of AO basis functions to file 613C 614 subroutine rt_tddft_spatial_potential_dump (rtdb, nao, nq, 615 $ xyz, chi_ao, d1chi_ao) 616 implicit none 617 618#include "errquit.fh" 619#include "mafdecls.fh" 620#include "stdio.fh" 621#include "global.fh" 622#include "msgids.fh" 623#include "util.fh" 624#include "cdft.fh" 625#include "geomP.fh" 626#include "geom.fh" 627#include "bas.fh" 628#include "matutils.fh" 629#include "rt_tddft.fh" 630 631 integer, intent(in) :: rtdb 632C type(rt_params_t), intent(in) :: params 633 integer, intent(in) :: nao, nq 634 double precision, intent(in) :: chi_ao(nq, nao), xyz(3, nq) 635 double precision, intent(in) :: d1chi_ao(nq, 3, nao) 636 637 character(len=*), parameter :: pname = 638 $ "rt_tddft_spatial_potential_dump: " 639 640 integer :: me, iq, ibf 641 integer :: ios, unitno 642 integer :: ix, iy, iz 643 double precision :: x, y, z 644 character(len=100) :: fname 645 646 character(len=*), parameter :: outfmt="(3es15.6e3)" 647C character(len=*), parameter :: outfmt="(3e20.10)" 648 649 650 character(len=100) :: fnamex, fnamey, fnamez 651 integer :: unitx, unity, unitz 652 653 654 me = ga_nodeid() 655 656C 657C Dump basis functions (on grid) to file 658C XXX DOESNT WORK IN PARALLEL YET 659C 660 if (me.eq.0) then 661 662 unitno = 34957 663 664 call util_file_name ("basis_grid.dat", .false., .false., fname) 665 666c$$$ Turned off by Adonay 667 668 unitx = 349834 669 call util_file_name ("basis_grid_der_x.dat", 670 $ .false., .false., fnamex) 671 unity = 12928 672 call util_file_name ("basis_grid_der_y.dat", 673 $ .false., .false., fnamey) 674 unitz = 74555 675 call util_file_name ("basis_grid_der_z.dat", 676 $ .false., .false., fnamez) 677C turned on by adonay 678 679 open (unitno, status="replace", file=fname, 680 $ iostat=ios, action="write") 681 if (ios .ne. 0) then 682 call errquit (pname//"failed to open: "//trim(fname),0,0) 683 endif 684 685 open (unitx, status="replace", file=fnamex, 686 $ iostat=ios, action="write") 687 if (ios .ne. 0) then 688 call errquit (pname//"failed to open: "//trim(fnamex),0,0) 689 endif 690 691 open (unity, status="replace", file=fnamey, 692 $ iostat=ios, action="write") 693 if (ios .ne. 0) then 694 call errquit (pname//"failed to open: "//trim(fnamey),0,0) 695 endif 696 697 open (unitz, status="replace", file=fnamez, 698 $ iostat=ios, action="write") 699 if (ios .ne. 0) then 700 call errquit (pname//"failed to open: "//trim(fnamez),0,0) 701 endif 702 703 704 write (unitno, fmt="(a)", iostat=ios) 705 $ "# Projection of AO basis on spatial grid" 706 write (unitx, fmt="(a)", iostat=ios) 707 $ "# Projection of x-derivative of AO basis on spatial grid" 708 write (unity, fmt="(a)", iostat=ios) 709 $ "# Projection of y-derivative of AO basis on spatial grid" 710 write (unitz, fmt="(a)", iostat=ios) 711 $ "# Projection of z-derivative of AO basis on spatial grid" 712 713 write (unitno, fmt="(a,i0)", iostat=ios) 714 $ "# nq = ", nq 715 write (unitx, fmt="(a,i0)", iostat=ios) 716 $ "# nq = ", nq 717 write (unity, fmt="(a,i0)", iostat=ios) 718 $ "# nq = ", nq 719 write (unitz, fmt="(a,i0)", iostat=ios) 720 $ "# nq = ", nq 721 722 723 write (unitno, fmt="(a,i0)", iostat=ios) 724 $ "# nbf = ", nao 725 write (unitx, fmt="(a,i0)", iostat=ios) 726 $ "# nbf = ", nao 727 write (unity, fmt="(a,i0)", iostat=ios) 728 $ "# nbf = ", nao 729 write (unitz, fmt="(a,i0)", iostat=ios) 730 $ "# nbf = ", nao 731 732 write (unitno, fmt="(a)", iostat=ios, advance="no") 733 $ "# x [au] y [au] z [au]" 734 write (unitx, fmt="(a)", iostat=ios, advance="no") 735 $ "# x [au] y [au] z [au]" 736 write (unity, fmt="(a)", iostat=ios, advance="no") 737 $ "# x [au] y [au] z [au]" 738 write (unitz, fmt="(a)", iostat=ios, advance="no") 739 $ "# x [au] y [au] z [au]" 740 741 do ibf = 1, nao 742 write (unitno, fmt="(i15)", iostat=ios, advance="no") ibf 743 write (unitx, fmt="(i15)", iostat=ios, advance="no") ibf 744 write (unity, fmt="(i15)", iostat=ios, advance="no") ibf 745 write (unitz, fmt="(i15)", iostat=ios, advance="no") ibf 746 enddo 747 write(unitno,*) "" 748 write(unitx,*) "" 749 write(unity,*) "" 750 write(unitz,*) "" 751 752 753c$$$ write (unitno, fmt="(a,i0,1es22.12e3,1es22.12e3)", iostat=ios) 754c$$$ $ "x ", nq(1), dq(1), qmin_ang(1) 755c$$$ 756c$$$ write (unitno, fmt="(a,i0,1es22.12e3,1es22.12e3)", iostat=ios) 757c$$$ $ "y ", nq(2), dq(2), qmin_ang(2) 758c$$$ 759c$$$ write (unitno, fmt="(a,i0,1es22.12e3,1es22.12e3)", iostat=ios) 760c$$$ $ "z ", nq(3), dq(3), qmin_ang(3) 761 762 do iq = 1, nq 763 764 x = xyz(1, iq) 765 y = xyz(2, iq) 766 z = xyz(3, iq) 767 768C if (abs(z) < 0.0001d0) then 769 770 write (unitno, fmt=outfmt, advance="no") 771 $ x, y, z 772 write (unitx, fmt=outfmt, advance="no") 773 $ x, y, z 774 write (unity, fmt=outfmt, advance="no") 775 $ x, y, z 776 write (unitz, fmt=outfmt, advance="no") 777 $ x, y, z 778 779c$$$ $ dbl_mb(ixyz + 3*(iq - 1) + 0), !x 780c$$$ $ dbl_mb(ixyz + 3*(iq - 1) + 1), !y 781c$$$ $ dbl_mb(ixyz + 3*(iq - 1) + 2) !z 782 783 do ibf = 1, nao 784 write (unitno, fmt=outfmt, advance="no") 785 $ chi_ao(iq, ibf) 786 write (unitx, fmt=outfmt, advance="no") 787 $ d1chi_ao(iq, 1, ibf) 788 write (unity, fmt=outfmt, advance="no") 789 $ d1chi_ao(iq, 2, ibf) 790 write (unitz, fmt=outfmt, advance="no") 791 $ d1chi_ao(iq, 3, ibf) 792 enddo 793 write (unitno, *) "" 794 write (unitx, *) "" 795 write (unity, *) "" 796 write (unitz, *) "" 797 798C endif 799 800 enddo 801 802c$$$ if (ios .ne. 0) then 803c$$$ call errquit(pname//"failed to write to: "//trim(fname),0,0) 804c$$$ endif 805c$$$ 806 close (unitno, iostat=ios) 807 if (ios .ne. 0) then 808 call errquit (pname//"failed to close: "//trim(fname),0,0) 809 endif 810 811 close (unitx, iostat=ios) 812 if (ios .ne. 0) then 813 call errquit (pname//"failed to close: "//trim(fnamex),0,0) 814 endif 815 close (unity, iostat=ios) 816 if (ios .ne. 0) then 817 call errquit (pname//"failed to close: "//trim(fnamey),0,0) 818 endif 819 close (unitz, iostat=ios) 820 if (ios .ne. 0) then 821 call errquit (pname//"failed to close: "//trim(fnamez),0,0) 822 endif 823 824 endif 825 826C call halt () 827 828 end subroutine 829 830 831 832 833 834 835C============================================================ 836C 837C 838 subroutine rt_tddft_spatial_potential_print (rtdb, nao, nq, 839 $ pot_ao, ovl_ao) 840 implicit none 841 842#include "errquit.fh" 843#include "mafdecls.fh" 844#include "stdio.fh" 845#include "global.fh" 846#include "msgids.fh" 847#include "util.fh" 848#include "cdft.fh" 849#include "geomP.fh" 850#include "geom.fh" 851#include "bas.fh" 852#include "matutils.fh" 853#include "rt_tddft.fh" 854 855 integer, intent(in) :: rtdb 856C type(rt_params_t), intent(in) :: params 857 integer, intent(in) :: nao, nq 858 double precision, intent(inout) :: pot_ao(nao,nao) 859 double precision, intent(inout) :: ovl_ao(nao,nao) 860 861 character(len=*), parameter :: pname = 862 $ "rt_tddft_spatial_potential_print: " 863 864 integer :: i, j, k 865 integer :: me 866 integer :: icen 867 character*16 icen_tag 868 double precision icen_loc(3), loc_ang(3) 869 double precision icen_charge 870 double precision :: intgrd, intgrd_max 871 872 873 me = ga_nodeid() 874 875C (check grid projection quality) 876 intgrd = 0d0 877 do i = 1, nao 878 intgrd = intgrd + abs(ovl_ao(i,i)) !abs needed? 879 enddo 880 881 intgrd = intgrd / dble(nao) !on-diagonal should be 1.0, so divide by nao to ideally get 1.0 882 intgrd_max = 1d0 !ideal case 883 884 885C 886C Print projected potential and overlap for diagnostic purposes. 887C 888 if (me.eq.0) then 889 write (luout, *) "" 890 call util_print_centered (luout, 891 $ "Projection of grid-based potential onto AO basis", 892 $ 40,.true.) 893 894 write (luout, *) "" 895 write (luout, "(1x,a,i0)") "Spatial grid points : ", nq 896 write (luout, "(1x,a,i0)") "AO basis functions : ", nao 897 write (luout, "(1x,a,1f10.4,a)") "Overall integral : ", 898 $ intgrd, " (ideal 1.0)" 899 900 write (luout, *) "" 901 write (luout, *) 902 $ " On-diagonal elements (overlap should be 1.0)" 903 write (luout, *) "Function Atom "// 904 $ "Element Overlap Potential" 905 write (luout, *) 906 $ "-------------------------------"// 907 $ "------------------------------------------" 908 endif 909 910 911C (doesnt use params struct values, gets ao_bas_hand and geom from common blocks) 912 do i = 1, nao 913 if (.not. bas_bf2ce (ao_bas_han, i, icen)) 914 $ call errquit (pname//"bas_bf2ce failed", 0, 0) 915 916c$$$ icen = -99 917 918 919C (note this acts on full active geom, specified by the handle 920C stored in params) 921 922c$$$ if (.not. geom_cent_get (params%geom_active_handle, icen, 923c$$$ $ icen_tag, icen_loc, icen_charge)) 924c$$$ $ call errquit (pname//"geom_cent_get active failed",0,0) 925 926 if (.not. geom_cent_get (geom, icen, 927 $ icen_tag, icen_loc, icen_charge)) 928 $ call errquit (pname//"geom_cent_get active failed",0,0) 929 930 if (me.eq.0) then 931 write (luout, "(i11, i9, 4x, a, 1f10.2, 1es22.12e3)") 932 $ i, icen, icen_tag, ovl_ao(i,i), pot_ao(i,i) 933 endif 934 935 enddo 936 937 end subroutine 938 939 940