1      logical function dft_energy_gradient(rtdb)
2      implicit none
3#include "rtdb.fh"
4      integer rtdb
5c
6      logical nwdft, grad_dft, xc_chktau
7      external nwdft, grad_dft, xc_chktau
8      logical status
9c
10      status = rtdb_cput(rtdb,'dft:theory', 1, 'dft')
11      status = nwdft(rtdb)
12      if (status) status = grad_dft(rtdb)
13      dft_energy_gradient = status
14      call grid_cleanup(.false.)
15c
16      end
17      logical function sodft_energy_gradient(rtdb)
18      implicit none
19#include "rtdb.fh"
20      integer rtdb
21c
22      logical nwdft, grad_dft
23      external nwdft, grad_dft
24      logical status
25c
26      status = rtdb_cput(rtdb,'dft:theory', 1, 'sodft')
27      status = nwdft(rtdb)
28      if (status) status = grad_dft(rtdb)
29      sodft_energy_gradient = status
30      call grid_cleanup(.false.)
31c
32      end
33      logical function grad_dft(rtdb)
34*
35* $Id$
36*
37      implicit none
38#include "errquit.fh"
39#include "mafdecls.fh"
40#include "global.fh"
41#include "bas.fh"
42#include "geom.fh"
43#include "rtdb.fh"
44#include "util.fh"
45#include "schwarz.fh"
46#include "stdio.fh"
47#include "cdft.fh"
48c
49c     dftgrad module.
50c
51c     Context is '...:dftgrad' --> changed to '...:dft'   JAN
52c
53c     Assumes DFT has been completed, MO vectors stored
54c     and all information is still in the RTDB
55c
56      logical int_normalize
57      external int_normalize
58      logical int_norm_2c
59      external int_norm_2c
60c
61      integer rtdb              ! [input] database handle
62      integer nbases
63      logical converged, status
64      integer bases(3), plevel
65      character*80 theory
66c
67c-----------------------------------------------------------------------
68c
69c     Push context down to DFT
70c
71      status = rtdb_parallel(.true.) ! Broadcast reads to all processes
72      call util_print_push()
73      call util_print_rtdb_load(rtdb,'dft')
74      call ecce_print_module_entry ('dft')
75c
76      if (.not. rtdb_cget(rtdb, 'dft:theory', 1, theory))
77     $     call errquit('task:energy: theory not specified',0, RTDB_ERR)
78      if(theory .eq. 'dft')then
79         status = rtdb_get(rtdb, 'dft:converged', MT_LOG, 1, converged)
80      else if(theory .eq. 'sodft')then
81         status = rtdb_get(rtdb, 'sodft:converged', MT_LOG,1,converged)
82      endif
83      if (.not.( status .and. converged ))then
84         if ((ga_nodeid() .eq. 0)
85     &      .and. util_print('information', print_none) )then
86            write (luout,*)'status: ', status, '   converged: ',
87     &                     converged
88        endif
89        call errquit(
90     &       'dft grad_dft: no converged DFT wavefunction available',
91     &     0, UNKNOWN_ERR)
92      endif
93c
94c     Extract high level info from the data-base setting defaults
95c
96c     load DFT parameters into common; turn off DFT printing
97c
98      call dft_pstat_init(rtdb)
99c
100      call util_print_get_level(plevel)
101      call util_print_set_level(print_none)
102      call dft_rdinput(rtdb)
103      call util_print_set_level(plevel)
104c
105      if (XCFIT.and.CDFIT)then
106         nbases = 3
107         bases(1) = AO_bas_han
108         bases(2) = CD_bas_han
109         bases(3) = XC_bas_han
110      elseif((.not.XCFIT).and.CDFIT)then
111         nbases = 2
112         bases(1) = AO_bas_han
113         bases(2) = CD_bas_han
114      elseif((.not.CDFIT).and.XCFIT)then
115         nbases = 2
116         bases(1) = AO_bas_han
117         bases(2) = XC_bas_han
118      else
119         nbases = 1
120         bases(1) = AO_bas_han
121      endif
122      call dft_inpana(rtdb)
123c
124c     initialize for schwarz screening
125c
126      if (nbases .gt. 1) call int_app_set_no_texas(rtdb)
127      call int_init(rtdb, 1, ao_bas_han)
128      call schwarz_init (geom, ao_bas_han)
129      call int_terminate()
130c
131c     initialize for derivative integrals
132c
133      call intd_init(rtdb,nbases,bases)
134c
135      if (ga_nodeid() .eq. 0)then
136         if (util_print('information',print_default) )then
137            write(luout,*)
138            write(luout,*)
139            call util_print_centered(luout,
140     &          'NWChem DFT Gradient Module',40,.true.)
141            write(luout,*)
142            write(luout,*)
143            if (title .ne. ' ')then
144               call util_print_centered(luout, title, 40, .false.)
145               write(luout,*)
146               write(luout,*)
147            endif
148
149            if (ipol .eq. 1)then
150               write(luout,1) rcharge, 'closed shell'
151            else
152               write(luout,1) rcharge, 'open shell'
153            endif
154            call util_flush(luout)
155
156         endif
157
158 1       format(/
159     &        '  charge          = ', f6.2/
160     &        '  wavefunction    = ', a/)
161         if (util_print('debug',print_debug) )then
162            if (.not. geom_print(geom))
163     &         call errquit('grad_dft: geom_print ?',0, GEOM_ERR)
164            if (.not. bas_print(ao_bas_han))
165     &         call errquit('grad_dft: bas_print ?',0, BASIS_ERR)
166         endif
167      endif
168c
169c     Compute the gradients arising from CD fit and XC
170c
171      if(theory .eq. 'dft')call dft_gradients(rtdb)
172      if(theory .eq. 'sodft')call dft_gradients_so(rtdb)
173c
174      call schwarz_tidy ()
175c
176      if(CDFIT)then
177        if (.not. bas_destroy(cd_bas_han))
178     &     call errquit('grad_dft:not able to destroy CD_bas:',86,
179     &       BASIS_ERR)
180      endif
181      if(XCFIT)then
182        if (.not.bas_destroy(XC_bas_han))
183     &     call errquit('grad_dft:not able to destroy XC_bas:',86,
184     &       BASIS_ERR)
185      endif
186
187      call intd_terminate()
188      if (nbases .gt. 1) call int_app_unset_no_texas(rtdb)
189c
190c     Compute the rest of the gradients (1-e and 2-e)
191c
192      if(theory .eq. 'dft')call grad_force(rtdb, ao_bas_han, geom)
193      if(theory .eq. 'sodft')call grad_force_so(rtdb, ao_bas_han, geom)
194c
195      if (.not.((bas_destroy(ao_bas_han)) .and. (geom_destroy(geom)) ) )
196     &   call errquit
197     &   ('grad_dft:error destroying geom and ao_bas_han handles',911,
198     &       BASIS_ERR)
199      call dft_pstat_print
200c
201c     terminate integral scope
202c
203      call ecce_print_module_exit ('dft','ok')
204      call util_print_pop
205      grad_dft = .true.
206c
207      return
208      end
209