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