1C> \ingroup task
2C> @{
3C>
4C> \brief The main driver for analytical Hessian calculations
5C>
6C> This routine checks the level of theory requested and if no
7C> analytical Hessians are available it will issue the evaluation of
8C> numerical Hessians. The resulting Hessian is stored on a file,
9C> rather than in the RTDB. The name of the Hessian file is stored
10C> on the RTDB instead of the Hessian itself.
11C>
12C> \return Return .true. if the Hessian was calculated successfully,
13C> and .false. otherwise.
14C>
15      logical function task_hessian(rtdb)
16*
17* $Id$
18*
19      implicit none
20#include "errquit.fh"
21#include "rtdb.fh"
22#include "mafdecls.fh"
23#include "global.fh"
24#include "inp.fh"
25#include "util.fh"
26#include "stdio.fh"
27      integer rtdb !< [Input] The RTDB handle
28c
29c     Generic NWChem interface to compute the analytic hessian.
30c
31c     RTDB input parameters
32c     ---------------------
33c     task:theory (string) - name of (QM) level of theory to use
34c
35c     RTDB output parameters no for analytic hessian at the moment.
36c     ----------------------
37c     task:hessian file name - that has a lower triangular
38C                              (double precision) array
39c                              derivative w.r.t. geometry cart. coords.
40c     task:status (logical)  - T/F for success/failure
41c     task:cputime (real)    - cpu time to execute the task
42c     task:walltime (real)   - wall time to execute the task
43c
44c     Also returns status through the function value
45c
46c     If the method does not have analytic derivatives automatically call
47c     the numerical derivative routine (not true in my case at the time).
48c
49      logical task_hessian_doit
50      external task_hessian_doit
51c
52      logical task_bsse
53      logical bsse_hessian
54      external bsse_hessian
55c
56      double precision cpu, wall, delta_pass
57      logical status
58      logical lcgmin
59      character*80 prefix
60      character*(nw_max_path_len) filehess
61c
62      call ecce_print_module_entry('task hessian')
63c
64      task_hessian = .false.
65c
66      cpu  = util_cpusec()
67      wall = util_wallsec()
68c
69      if (.not.rtdb_get(rtdb, 'dft:cgmin', mt_log, 1, lcgmin))
70     &    lcgmin = .false.
71c
72c     Right now only have a QM component
73c     but the calculation might be called with *doit as
74c     the others task_*
75c
76      if (.not. rtdb_get(rtdb,'bsse',mt_log,1,task_bsse))
77     $     task_bsse = .false.
78c
79      if(.not.task_bsse) then
80        status= task_hessian_doit(rtdb)
81        if(.not.status)
82     $    call
83     $  errquit('task_gradient: error task_hessian_doit',0,UNKNOWN_ERR)
84      else
85        status = bsse_hessian(rtdb)
86        if(.not.status)
87     $    call
88     $  errquit('task_hessian: error call bsse hessian',911,UNKNOWN_ERR)
89      endif
90c
91      if(.not.status)
92     $    call errquit('task_hessian: error task_hessian_doit',0,0)
93c
94      cpu  = util_cpusec() - cpu
95      wall = util_wallsec() - wall
96c
97      task_hessian = status
98      if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, task_hessian))
99     $     call errquit('task_hessian: failed to store status',0,
100     &       RTDB_ERR)
101c
102      call util_file_name('hess',  .false., .false.,filehess)
103      if (.not.rtdb_cput(rtdb,'task:hessian file name',1,filehess))
104     &    call errquit('task_hessian: failed to store filename',911,
105     &       RTDB_ERR)
106      if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cpu))
107     $    call errquit('task_gradient: failed storing cputime',0,
108     &       RTDB_ERR)
109      if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, wall))
110     $    call errquit('task_gradient: failed storing walltime',0,
111     &       RTDB_ERR)
112c
113      call ecce_print1('cpu time', mt_dbl, cpu, 1)
114      call ecce_print1('wall time', mt_dbl, wall, 1)
115      if (task_hessian) then
116         call ecce_print_module_exit('task hessian', 'ok')
117      else
118         call ecce_print_module_exit('task hessian', 'failed')
119      endif
120      end
121C> @}
122c
123      logical function task_hessian_doit(rtdb)
124c
125      implicit none
126#include "errquit.fh"
127#include "rtdb.fh"
128#include "mafdecls.fh"
129#include "global.fh"
130#include "inp.fh"
131#include "util.fh"
132      integer rtdb
133c
134c     Generic NWChem interface to compute the analytic hessian.
135c
136c     RTDB input parameters
137c     ---------------------
138c     task:theory (string) - name of (QM) level of theory to use
139c     task:numerical (logical) - optional - if true use numerical
140c         differentiation. if
141c     task:analytic  (logical) - force analytic hessian
142c
143c     RTDB output parameters no for analytic hessian at the moment.
144c     ----------------------
145c     task:hessian file name - that has a lower triangular
146C                              (double precision) array
147c                              derivative w.r.t. geometry cart. coords.
148c     task:status (logical)  - T/F for success/failure
149c     task:cputime (real)    - cpu time to execute the task
150c     task:walltime (real)   - wall time to execute the task
151c
152c     Also returns status through the function value
153c
154c     If the method does not have analytic derivatives automatically call
155c     the numerical derivative routine (not true in my case at the time).
156c
157      logical stpr_gen_hess_at
158      external stpr_gen_hess_at
159      logical  scf,hess_anal,hess_check
160      external scf,hess_anal,hess_check
161c
162      logical task_bsse
163c
164      integer ecce_old_print, ecce_junk_print
165      logical status, ignore, numerical, analytic
166      double precision delta_pass
167      double precision default_delta
168      character*80 prefix
169      character*32 theory
170c
171c     vdw contrib
172      logical disp
173      logical xc_chkdispauto
174      external xc_chkdispauto
175c
176      call ecce_print_module_entry('task hessian')
177c
178c     Prevent BSSE calculations
179      if (.not. rtdb_get(rtdb,'bsse',mt_log,1,task_bsse))
180     $    task_bsse = .false.
181c
182      if (task_bsse) then
183        if (.not. rtdb_put(rtdb,'bsse',mt_log,1,.false.))
184     $    call errquit('task_hessian_doit:rtdb_put failed',911,RTDB_ERR)
185      endif
186c
187c     Check for vdw (disp) calculations
188      disp=.false.
189      if (.not.rtdb_get(rtdb, 'dft:disp', mt_log, 1, disp))
190     &      disp=.false.
191c
192      if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory))
193     $     call errquit('task:hessian: theory not specified',0,
194     &       RTDB_ERR)
195      if (.not. rtdb_get(rtdb,'task:numerical',mt_log,1,numerical))
196     &    numerical = .false.
197      if (.not. rtdb_get(rtdb,'task:analytic',mt_log,1,analytic))
198     &    analytic = .true.
199c
200      ignore = rtdb_delete(rtdb, 'task:numerical')
201      ignore = rtdb_delete(rtdb, 'task:hessian file name')
202      ignore = rtdb_delete(rtdb,'task:hessian')
203      if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .false.))
204     $     call errquit('task_hessian: failed to invalidate status',0,
205     &       RTDB_ERR)
206
207      prefix = theory ! Most common
208c
209c  Make sure that we don't do analytic for anything other than
210c  scf and dft.
211c
212      if ((.not.inp_compare(.false.,theory,'scf')).and.
213     $    (.not.inp_compare(.false.,theory,'dft'))) then
214        numerical = .true.
215        analytic  = .false.
216      end if
217c
218c     == make sure the theory supports analytic hessians ==
219      if (.not.hess_check(rtdb)) then
220          numerical = .true.
221          analytic  = .false.
222      endif
223c
224c     Actually do the deed
225c
226      if (numerical) then
227        if (ga_nodeid().eq.0 .and.
228     $      util_print('task_hessian', print_low)) then
229          write(6,*)
230          write(6,*)
231          call util_print_centered(6,
232     $        'NWChem Finite-difference Hessian ',
233     $        40,.true.)
234          write(6,*)
235          write(6,*)
236        endif
237c
238c define default delta value
239c
240        default_delta = 0.001d00  ! should it be 0.0025
241        if (inp_compare(.false.,theory,'dft'))
242     &       default_delta = 0.01d00
243        if (inp_compare(.false.,theory,'pspw')   .or.
244     &      inp_compare(.false.,theory,'band')   .or.
245     &      inp_compare(.false.,theory,'paw'))
246     &       default_delta = 0.01d00
247c
248c     check for change of finite difference delta
249c
250        if (.not.rtdb_get(rtdb,'stpr_gen:delta',MT_DBL,1,delta_pass))
251     &      delta_pass = default_delta
252c
253        call ecce_print_control(0, ecce_old_print) ! disable ECCE printing
254        status = stpr_gen_hess_at(rtdb,delta_pass)
255        call ecce_print_control(ecce_old_print,ecce_junk_print) ! re-enable ECCE printing
256
257      else  ! analytic hessian
258c
259c       == calculate analytic hessian ==
260        if (ga_nodeid().eq.0 .and.
261     $      util_print('task_hessian', print_low)) then
262          write(6,*)
263          write(6,*)
264          call util_print_centered(6,
265     $        'NWChem Analytic Hessian ',
266     $        40,.true.)
267          write(6,*)
268        endif
269c
270        if (theory .eq. 'scf' .or. theory .eq. 'dft') then
271          status = hess_anal(rtdb)
272        else
273          call errquit('task_hessian: unknown analytic theory',0,
274     &       INPUT_ERR)
275        endif
276c
277c       vdw contribution to analytic hessian
278c       for numerical hessians, this contribution is added via the gradients
279        if(theory.eq.'dft') then
280          if(disp.or.xc_chkdispauto()) call xc_vdw_to_hessian(rtdb)
281        endif ! vdw
282      endif  ! numerical or analytic
283c
284      if (task_bsse) then
285        if (.not. rtdb_put(rtdb,'bsse',mt_log,1,task_bsse))
286     $   call errquit('task_hessian_doit:rtdb_put failed',911,RTDB_ERR)
287      endif
288c
289      task_hessian_doit = status
290      end
291