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