1      integer function rt_tddft_find_itreference (params)
2#include "errquit.fh"
3#include "global.fh"
4#include "stdio.fh"
5#include "rt_tddft.fh"
6
7C     == Inputs ==
8      type(rt_params_t), intent(in) :: params
9
10
11C     == Parameters ==
12      character(*) ,parameter  :: pname = "rt_tddft_find_itreference: "
13      double precision, parameter :: thresh = 1d-9
14
15C     == Variables ==
16      integer it, itreference
17      double precision tt
18      logical found
19
20      found = .false.
21      itreference = -99
22
23      do it = 1, params%nt
24         tt = params%tmin + (it-1)*params%dt
25
26         if ( abs(tt - params%viz_treference) .lt. thresh) then
27            if (found) then
28               call errquit (pname//"multiple values found",0,0)
29            endif
30
31            itreference = it
32            found = .true.
33         endif
34      enddo
35
36      if (.not. found) then
37         call errquit (pname//"failed to find value",0,0)
38      endif
39
40      call ga_sync ()
41
42      if (itreference < 1)
43     $     call errquit (pname//"bad itreference",0,0)
44
45      rt_tddft_find_itreference = itreference
46
47      end function
48
49
50
51      logical function rt_tddft_at_snap (params, tt)
52#include "errquit.fh"
53#include "global.fh"
54#include "stdio.fh"
55#include "rt_tddft.fh"
56
57C     == Inputs ==
58      type(rt_params_t), intent(in) :: params
59      double precision, intent(in)  :: tt
60
61C     == Parameters ==
62      character(*) ,parameter  :: pname = "rt_tddft_at_snap: "
63      double precision, parameter :: thresh = 1d-9
64
65C     == Variables ==
66      logical at_snap
67
68
69C
70C     if visualization is active, only make snapshots between the
71C     specified range.
72C
73      if ((tt.le.params%viz_tend).and.(tt.ge.params%viz_tstart)) then
74         at_snap = .true.
75      else
76         at_snap = .false.
77      endif
78
79
80C
81C     also if this is the "reference" density matrix, make a snapshot
82C
83      if ( abs(tt - params%viz_treference) .lt. thresh) then
84         at_snap = .true.
85      endif
86
87      at_snap = at_snap .and. params%viz_active
88
89      rt_tddft_at_snap = at_snap
90
91      end function
92
93
94
95C====================================================================
96C
97C     Generate a filename for density matrix dump (snapshot) with
98C     with time step index.
99C
100      subroutine rt_tddft_snapshot_fname (params, it, tt, fout)
101      implicit none
102
103#include "errquit.fh"
104#include "global.fh"
105#include "stdio.fh"
106#include "rt_tddft.fh"
107
108
109C     == Inputs ==
110      type(rt_params_t), intent(in) :: params
111      integer, intent(in)           :: it
112      double precision, intent(in)  :: tt
113
114
115C     == Outputs ==
116      character(50), intent(out)    :: fout  !hardcoded str len; fix if changing fname format
117
118
119C     == Parameters ==
120      character(*) ,parameter  :: pname = "rt_tddft_snapshot_fname: "
121
122      if (params%nt .gt. 99999999)
123     $     call errquit (pname//"nt too big; fix formatting", 0, 0)
124
125C      write (fout, "(a,i0.10,a,e12.6,a)")
126C     $        "ptot_ao_re.", it, ".", tt, "-au"
127
128      write (fout, "(a,i0.10)")
129     $        "ptot_ao_re.", it
130
131      end subroutine
132
133
134C====================================================================
135      subroutine rt_tddft_snapshot_fname_cube (params, it, tt, fout)
136      implicit none
137
138#include "errquit.fh"
139#include "global.fh"
140#include "stdio.fh"
141#include "rt_tddft.fh"
142
143C     == Inputs ==
144      type(rt_params_t), intent(in) :: params
145      integer, intent(in)           :: it
146      double precision, intent(in)  :: tt
147
148C     == Outputs ==
149      character(50), intent(out)    :: fout  !note hardcoded from fname call
150
151
152C     == Parameters ==
153      character(*) ,parameter :: pname="rt_tddft_snapshot_fname_cube: "
154
155      if (params%nt .gt. 99999999)
156     $     call errquit (pname//"nt too big; fix formatting", 0, 0)
157
158      if (params%viz_subgs) then
159         write (fout, "(a,i0.10,a)")
160     $        "density_subgs.", it, ".cube"
161      else
162         write (fout, "(a,i0.10,a)")
163     $        "density.", it, ".cube"
164      endif
165
166      end subroutine
167
168
169
170C====================================================================
171C
172C     Dumps the supplied real part of the dens mat in AO basis to file;
173C     we later read in this snapshot to make density plots, etc.
174C
175      subroutine rt_tddft_snapshot_save (params, it, tt, g_densre_ao)
176      implicit none
177
178#include "errquit.fh"
179#include "mafdecls.fh"
180#include "global.fh"
181#include "stdio.fh"
182#include "rtdb.fh"
183#include "matutils.fh"
184#include "rt_tddft.fh"
185
186
187C     == Inputs ==
188      type(rt_params_t), intent(in) :: params
189      integer, intent(in)           :: it
190      double precision, intent(in)  :: tt
191      integer, intent(in)           :: g_densre_ao  !re part of *total* dens mat in AO basis
192
193
194C     == Parameters ==
195      character(*), parameter :: pname = "rt_tddft_snapshot_save: "
196
197
198C     == Variables ==
199      character(255) fname
200      character(50) fname_tail  !note hardcoded str length; update if changing fname format
201      double precision elapsed
202
203      if (params%prof) call prof_start (elapsed)
204
205      call rt_tddft_snapshot_fname (params, it, tt, fname_tail)
206      call util_file_name (fname_tail, .false., .false., fname)
207
208      if (.not. dmat_io_dump (g_densre_ao, trim(fname)))
209     $     call errquit (pname//"failed to dump densao_re",0,0)
210
211      if (params%prof) call prof_end (elapsed, "Saving snapshot")
212      end subroutine
213
214
215C====================================================================
216C
217C     Convert all stored density matrix snapshots (real part in AO
218C     basis) to charge densities .cube files using dplot.  Note: dplot
219C     will clean up after itself, which includes shutting down the
220C     integrals, etc.  Call this at the very end.
221C
222      subroutine rt_tddft_snapshot_dplot (params, subgs)
223      implicit none
224
225#include "errquit.fh"
226#include "mafdecls.fh"
227#include "global.fh"
228#include "stdio.fh"
229#include "rtdb.fh"
230#include "matutils.fh"
231#include "rt_tddft.fh"
232
233C     == Inputs ==
234      type(rt_params_t), intent(in) :: params
235      logical, intent(in)           :: subgs  !if true, subtract ground state density
236
237
238C     == Parameters ==
239      character(*), parameter :: pname = "rt_tddft_snapshot_dplot: "
240
241
242C     == Variables ==
243      character(255) fname_dens, fname_dens_gs, fname_cube
244      character(50) fname_tail  !note hardcoded str length; update if changing fname format
245      character(50) fname_tail_cube
246      integer it, itreference
247      double precision tt
248      integer isnap
249      integer me
250      integer num_basis,bas_handles(2)
251      double precision elapsed
252      character(64) title
253
254
255C     == External ==
256      logical, external :: dplot
257C      logical, external :: rt_tddft_atpoint
258      logical, external :: rt_tddft_at_snap
259      integer, external :: rt_tddft_find_itreference
260
261
262      if (params%prof) call prof_start (elapsed)
263
264      me = ga_nodeid ()
265
266      if (me.eq.0) then
267         call util_print_centered (luout,
268     $        "Post-processing of density matrix snapshots", 20, .true.)
269
270         write (luout, *) ""
271         write (luout, *) ""
272      endif
273
274
275C     (determine index of density matrix we want to subtract)
276      itreference = rt_tddft_find_itreference (params)
277
278C
279C     If desired, set up dplot to subtract ground state: rho = rho(t) -
280C     rho(t=0):
281C
282      if (subgs) then
283
284C         call rt_tddft_snapshot_fname(params,1, params%tmin, fname_tail)  !this is hardcoded for t=0
285
286         call rt_tddft_snapshot_fname (params, itreference,
287     $        params%viz_treference, fname_tail)
288
289         call util_file_name(fname_tail, .false.,.false., fname_dens_gs)
290
291         if (.not.rtdb_cput(params%rtdb,"dplot:File_Mat2",
292     $        1,fname_dens_gs))
293     $        call errquit(pname//"Write failed to rtdb",0,RTDB_ERR)
294         call rt_tddft_print_notice (
295     $        "Subtracting reference density matrix")
296      else
297         call rt_tddft_print_notice (
298     $        "Not subtracting reference density matrix")
299      endif
300
301c     int cleanup to avoid clashes in dplot
302            call int_terminate
303
304C
305C     Loop over all times and convert each stored dens mat snapshot to
306C     cube.
307C
308      isnap = 0
309      it = 1
310      do while (it .lt. params%nt)
311         tt = params%tmin + (it-1)*params%dt
312
313C         if (rt_tddft_atpoint(it, params%nt, params%nsnapshots)) then
314         if (rt_tddft_at_snap (params, tt)) then
315
316            isnap = isnap + 1
317
318            if (me.eq.0) then
319               write (luout, "(a,i0,a,f12.3,a)")
320     $              "Postprocessing snapshot ",
321     $              isnap, ", t = ", tt, " au ..."
322            endif
323
324C     (title of the slice is the time in atomic units)
325            write (title, "(a,1es14.7,a,1es14.7,a)")
326     $           "t = ", tt, " au = ", tt*au2fs, " fs"
327
328            call rt_tddft_snapshot_fname (params, it, tt, fname_tail)
329            call util_file_name (fname_tail,.false.,.false., fname_dens)
330
331            call rt_tddft_snapshot_fname_cube (params, it, tt,
332     $           fname_tail_cube)
333
334C     changed to match new dplot api (no leading ./perm)
335c$$$            call util_file_name (trim(fname_tail_cube),
336c$$$     $        .false., .false., fname_cube)
337            fname_cube = trim(fname_tail_cube)
338
339C     (load dplot params in rtdb)
340            if (.not.rtdb_cput(params%rtdb,'dplot:Title',1,title))
341     &           call errquit(pname//"Write failed to rtdb",0,RTDB_ERR)
342
343            if (.not.rtdb_cput(params%rtdb,"dplot:File_Out" ,
344     $           1,fname_cube))
345     $           call errquit(pname//"Write failed to rtdb",0,RTDB_ERR)
346
347            if (.not.rtdb_cput(params%rtdb,"dplot:File_Mat1",
348     $           1,fname_dens))
349     $           call errquit(pname//"Write failed to rtdb",0,RTDB_ERR)
350
351C     (call dplot to convert)
352            if (.not. dplot (params%rtdb))
353     $           call errquit (pname//"dplot call failed", 0, 0)
354         endif
355
356         it = it + 1
357      enddo
358
359c     int cleanup to avoid clashes in dplot
360      bas_handles(1)=params%ao_bas_han
361      num_basis=1
362      if(params%iVcoul_opt.eq.1) then
363         num_basis=2
364         bas_handles(2)=params%cd_bas_han
365      endif
366      call int_init(params%rtdb,num_basis,bas_handles)
367
368      if (params%prof) call prof_end (elapsed,
369     $     "Converting snapshots to cube files")
370      end subroutine
371c $Id$
372