1C
2C     rt_tddft_params_init_common.F
3C
4C     Computes parameters which are universal to closed shell,
5C     openshell, etc and stores them in the params struct.  Also loads
6C     rtdb parameters into params.
7C
8      subroutine rt_tddft_init_common (params)
9      implicit none
10
11#include "errquit.fh"
12#include "mafdecls.fh"
13#include "stdio.fh"
14#include "global.fh"
15#include "rtdb.fh"
16#include "geom.fh"
17#include "util.fh"
18#include "cdft.fh"
19#include "matutils.fh"
20#include "zora.fh"
21#include "rt_tddft.fh"
22
23
24C     == In/out ==
25      type(rt_params_t), intent(inout) ::  params
26
27
28C     == Parameters ==
29      character(*), parameter :: pname = "rt_tddft_init_common: "
30
31
32C     == Variables ==
33      integer n
34      logical ok
35      integer me
36
37      double precision local_d_conv
38
39      double precision center(3)
40      logical dft_converged
41
42      integer filesize, memsize
43
44C     (tmp dip and quad integral matricies; all nbf_ao x nbf_ao)
45      integer g_dip_x, g_dip_y, g_dip_z
46      integer g_quad_xx, g_quad_xy, g_quad_xz
47      integer g_quad_yy, g_quad_yz, g_quad_zz
48
49      double precision elapsed
50      integer g_fock_ao_core
51
52
53      me = ga_nodeid ()
54
55
56C
57C     Hardcoded parameters.
58C
59      params%tol_zero = 1d-8
60
61
62C
63C     Get density tolerance--dictates precision of Fock matrix calculations.
64C
65      if (.not. rtdb_get (params%rtdb, "dft:d_conv",
66     $     mt_dbl, 1, local_d_conv)) then
67
68         call errquit (pname//"Failed to read d_conv", 0, 0)
69      else
70         params%d_conv = local_d_conv
71      endif
72
73
74
75C
76C     Check that DFT SCF converged from before.
77C
78      if (.not. rtdb_get (params%rtdb, "dft:converged",
79     $     mt_log, 1, dft_converged)) then
80
81         call rt_tddft_print_warning (
82     $        "Failed to find previous SCF convergence status."//
83     $        " Either you are loading vectors from a previous job"//
84     $        " or database is corrupt--be careful.")
85
86      else
87
88         if (.not. dft_converged) call rt_tddft_print_warning (
89     $        "SCF did not converge previously--it is quite likely"//
90     $        " that the starting density is not in the ground state")
91
92      endif
93
94
95C
96C     Store params passed in from SCF or in headers.
97C
98      params%nbf_ao = nbf_ao
99      params%ipol = ipol
100      params%ao_bas_han = ao_bas_han
101      params%cd_bas_han = cd_bas_han
102
103
104C
105C     Get values from rtdb and store in params struct.
106C
107      call rt_tddft_init_rtdb (params)
108
109
110
111C     XXX DISABLED FOR NOW
112#if 0
113C
114C     Disable *all* disk I/O for DFT routines (default behavior).
115C
116      if (params%nodisk) then
117
118C     (dont use disk for gridpts)
119         if(.not.rtdb_put(params%rtdb,'dft:largenode',mt_log,1,.true.))
120     $        call errquit(pname//"rtdb_put largenode failed",
121     $        0, RTDB_ERR)
122
123
124C     (direct)
125         if (.not. rtdb_put(params%rtdb, 'dft:direct',mt_log,1,.true.))
126     $        call errquit(pname//"rtdb_put direct failed",0,RTDB_ERR)
127
128         filesize = -1
129         memsize  = -1
130         if (.not.rtdb_put(params%rtdb,
131     $        'int2e:filesize',mt_int,1,filesize))
132     $        call errquit(pname//'failed to store filesize',
133     $        0, RTDB_ERR)
134
135         if (.not.rtdb_put(params%rtdb,
136     $        'int2e:memsize',mt_int,1,memsize))
137     $        call errquit(pname//'failed to store filesize',
138     $        0, RTDB_ERR)
139
140
141C     (no I/O)
142         if (.not. rtdb_put(params%rtdb, 'dft:noio', mt_int, 1, .true.))
143     $        call errquit(pname//"rtdb_put noio failed", 0, RTDB_ERR)
144
145         if (me.eq.0) then
146            write (luout, *) "*** All disk I/O has been disabled ***"
147         endif
148
149      else
150
151         call rt_tddft_print_warning (
152     $        "Disk I/O using DFT options; proceed with caution")
153
154      endif
155#endif
156
157
158
159C
160C     Dummy values XXX NEVER USED?
161C
162CXXX  [KAL]: add check to ensure none of these are used
163      params%g_vxc(1) = -999999
164      params%g_vxc(2) = -999999
165      params%g_vxc(3) = -999999
166      params%g_vxc(4) = -999999
167      params%g_xcinv  = -999999
168      params%rho_n = -99d0
169      params%lzorasf = .false.
170
171C
172C     Nuclear repulsion energy
173C
174      if (.not. geom_nuc_rep_energy (geom, params%Enuc))
175     $     call errquit (pname//'Failed to compute Enuc', 0, BASIS_ERR)
176
177
178C
179C     I/O logical flag (used in dft call)
180C     this might have been set in the override i/o block above.
181C
182      if (noio.eq.1) then
183         params%iolgc = .false.
184      else
185         params%iolgc = .true.
186      endif
187
188
189C
190C     Charge density fitting flag
191C
192      if (cd_bas_han.eq.-99999)then
193         params%iVcoul_opt = 0
194      else
195         params%iVcoul_opt = 1
196      endif
197
198
199C
200C     XC fitting flag
201C
202      if (xc_bas_han.eq.-99999)then
203         params%iVxc_opt = 0
204      else
205         params%iVxc_opt = 1
206      endif
207
208
209C
210C     Check parameters.
211C
212      call rt_tddft_init_check_params (params)
213
214
215C
216C     Calculate the total number of time steps.  Also, if number of
217C     prints, checks, snapshots, etc is greater than nt, set them to nt.
218C
219      if (params%nprints < 0)
220     $     call errquit (pname//"nprints must be >= 0", 0, 0)
221
222      if (params%nchecks < 0)
223     $     call errquit (pname//"nchecks must be >= 0", 0, 0)
224
225c$$$      if (params%nsnapshots < 0)
226c$$$     $     call errquit (pname//"nsnapshots must be >= 0", 0, 0)
227
228      if (params%nrestarts < 0)
229     $     call errquit (pname//"nrestarts must be >= 0", 0, 0)
230
231
232      params%nt = ceiling((params%tmax-params%tmin)/params%dt)
233
234      if (params%nprints > params%nt) params%nprints = params%nt
235      if (params%nchecks > params%nt) params%nchecks = params%nt
236c$$$      if (params%nsnapshots > params%nt) params%nsnapshots = params%nt
237      if (params%nrestarts > params%nt) params%nrestarts = params%nt
238
239
240C     hardcoded for now:
241      params%nestims = params%nt / 100  !i.e., every 100 steps
242
243
244
245C
246C     Store geom tag for "active" geometry.
247C
248      params%geom_active_handle = geom
249
250
251C
252C     Center of mass.
253C
254      if (.not.geom_center_of_mass(params%geom_active_handle,center))
255     &     call errquit (pname//'could not get center of mass',
256     $     0, GEOM_ERR)
257
258      params%origin(1) = center(1)
259      params%origin(2) = center(2)
260      params%origin(3) = center(3)
261
262
263C
264C     Number of exchange correlation energies.
265C
266      params%nExc = idecomp + 1
267
268
269C
270C     Number of atoms.
271C
272      if (.not. geom_ncent(params%geom_active_handle, params%natoms))
273     &     call errquit(pname//'geom_ncent failed',73, GEOM_ERR)
274
275
276
277C
278C     Initialize all integrals.
279C
280      call rt_tddft_init_fock2e (params)
281      call rt_tddft_init_dip_quad_ints (params)
282      if (cdfit) call rt_tddft_init_coulcdfit (params)
283
284
285C
286C     Initialize geometries.
287C
288      call rt_tddft_init_geoms (params)
289
290
291C
292C     Compute and store 1e (time-dependent) core part of the Fock
293C     matrix.  Note, the GA created here will be destroyed in
294C     rt_tddft_clean().
295C
296       if (params%prof) call prof_start (elapsed)
297
298C     xxx check ns_ao vs nbf_ao
299      if (.not. ga_create (mt_dbl, params%nbf_ao, params%nbf_ao,
300     $     "fock_ao_core", 0, 0, g_fock_ao_core))
301     $     call errquit (pname//"create fock_ao_core failed",0,GA_ERR)
302
303      call ga_zero (g_fock_ao_core)
304      call int_1e_ga(ao_bas_han, ao_bas_han,
305     $     g_fock_ao_core, 'kinetic', oskel)
306      call int_1e_ga(ao_bas_han, ao_bas_han,
307     $     g_fock_ao_core, 'potential', oskel)
308
309      params%g_fock_ao_core = g_fock_ao_core
310
311      if (params%prof) call prof_end(elapsed,"Fock core initialization")
312
313
314
315      if (params%static)
316     $     call rt_tddft_print_warning ("Static calculation--"//
317     $     "Fock matrix will be frozen")
318
319
320C
321C     Spatial complex absorbing potential (CAP).
322C
323      if (params%cap_active) then
324         if (.not. ga_create(mt_dcpl, nbf_ao, nbf_ao ,"zCAP",
325     $        0, 0, params%g_zcap))
326     $        call errquit (pname//"zcap create failed", 0, GA_ERR)
327
328         call rt_tddft_spatialcap (params, nbf_ao, params%g_zcap) !calc and store CAP
329      endif
330
331
332      end subroutine rt_tddft_init_common
333c $Id$
334