1c
2c     Routines for calculating XC functional second derivatives
3c     by finite difference
4c
5c     BGJ - 1/02
6c
7c     $Id$
8c
9c     Finite difference step size
10c
11#define STEP_SIZE 0.0001
12c
13      subroutine xc_setup_fd(tol_rho, rho, delrho, qwght, nq, ipol, GC,
14     &     l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc,
15     &     i_qwght_copy)
16c
17      implicit none
18#include "errquit.fh"
19c
20      integer nq, ipol, l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat,
21     &     i_pfunc, i_qwght_copy
22      double precision tol_rho, rho(nq,*), delrho(nq,3,*), qwght(nq)
23      logical GC
24c
25c     Local variables
26c
27      integer npert, len_prho, len_pAmat, len_prho_tot, len_pAmat_tot
28      integer i
29c
30#include "mafdecls.fh"
31c
32      npert = 2
33      len_prho = 2*3*nq
34      len_pAmat = npert*npert*nq
35      if (GC) then
36         npert = 5
37         len_prho_tot = len_prho + 2*6*nq
38         len_pAmat_tot = npert*npert*nq
39      else
40         len_prho_tot = len_prho
41         len_pAmat_tot = len_pAmat
42      endif
43      if (.not.
44     &     MA_Alloc_Get(MT_DBL,len_prho_tot+len_pAmat_tot+4*nq*npert,
45     &     'xc_fd',l_storage,i_prho))
46     &     call errquit('xc_setup_fd: cannot allocate storage',0,
47     &       MA_ERR)
48      call dfill(len_prho_tot+len_pAmat_tot+4*nq*npert, 0d0,
49     &     dbl_mb(i_prho), 1)
50      i_pAmat = i_prho + len_prho_tot
51      i_pdelrho = i_prho + len_prho
52      i_pCmat = i_pAmat + len_pAmat
53      i_pfunc = i_pAmat + len_pAmat_tot
54      i_qwght_copy = i_pfunc + 2*nq*npert
55c
56c     Copy weights
57c
58      do i = 1, 2*npert
59         call dcopy(nq, qwght, 1, dbl_mb(i_qwght_copy+(i-1)*nq), 1)
60      enddo
61c
62c     Set up perturbed density values
63c
64      call xc_pert_rho(tol_rho, rho, delrho, nq, ipol, GC,
65     &     dbl_mb(i_prho), dbl_mb(i_pdelrho))
66c
67      return
68      end
69c
70      subroutine xc_pert_rho(tol_rho, rho, delrho, nq, ipol, npert,
71     &     prho, pdelrho)
72c
73      implicit none
74c
75      integer nq, ipol, npert
76      double precision tol_rho, rho(nq,*), delrho(nq,3,*),
77     &     prho(nq,3,2,npert), pdelrho(nq,3,2,2,npert)
78c
79c     Local variables
80c
81      double precision h
82      parameter (h = STEP_SIZE)
83      logical GC, DoA, DoB
84      integer i, ipert
85      double precision ra, drax, dray, draz, Gaa, Gab, Gbb,
86     &     mod_drhoa, mod_drhob, cos_theta, x, mod_drhob_p,
87     &     cos_theta_p, sin_theta_p
88c
89      GC = npert.eq.5
90c
91c     Initial copy of density values - copy as unrestricted
92c
93      call dfill(nq*3, 0d0, prho, 1)
94      call dfill(nq*6, 0d0, pdelrho, 1)
95      if (ipol.eq.1) then
96         do i = 1, nq
97            if (rho(i,1).gt.tol_rho) then
98               ra = rho(i,1) * 0.5d0
99               prho(i,1,1,1) = rho(i,1)
100               prho(i,2,1,1) = ra
101               prho(i,3,1,1) = ra
102            endif
103         enddo
104         if (GC) then
105            do i = 1, nq
106               drax = delrho(i,1,1) * 0.5d0
107               dray = delrho(i,2,1) * 0.5d0
108               draz = delrho(i,3,1) * 0.5d0
109               pdelrho(i,1,1,1,1) = drax
110               pdelrho(i,2,1,1,1) = dray
111               pdelrho(i,3,1,1,1) = draz
112               pdelrho(i,1,2,1,1) = drax
113               pdelrho(i,2,2,1,1) = dray
114               pdelrho(i,3,2,1,1) = draz
115            enddo
116         endif
117      else
118         do i = 1, nq
119            if (rho(i,1).gt.tol_rho) then
120               prho(i,1,1,1) = rho(i,1)
121               prho(i,2,1,1) = rho(i,2)
122               prho(i,3,1,1) = rho(i,3)
123            endif
124         enddo
125         if (GC) then
126            call dcopy(nq*6, delrho, 1, pdelrho, 1)
127         endif
128      endif
129c
130c     Copy to all other perturbation locations
131c
132      call dcopy(nq*3, prho, 1, prho(1,1,2,1), 1)
133      if (GC) then
134         call dcopy(nq*6, pdelrho, 1, pdelrho(1,1,1,2,1), 1)
135      endif
136      do ipert = 2, npert
137         do i = 1, 2
138            call dcopy(nq*3, prho, 1, prho(1,1,i,ipert), 1)
139            if (GC) then
140               call dcopy(nq*6, pdelrho, 1, pdelrho(1,1,1,i,ipert), 1)
141            endif
142         enddo
143      enddo
144c
145c     Perturb the density parameter values - we aren't concerned if
146c     some density values go negative as a result of this since the
147c     functional implementation routines map negative densities to
148c     a zero value for the functional and its derivatives (e.g. like
149c     analytic continuation)
150c
151      do ipert = 1, 2
152         do i = 1, nq
153            if (prho(i,1,1,ipert).gt.tol_rho) then
154               prho(i,1,1,ipert) = prho(i,1,1,ipert) + h
155               prho(i,ipert+1,1,ipert) = prho(i,ipert+1,1,ipert) + h
156               prho(i,1,2,ipert) = prho(i,1,2,ipert) - h
157               prho(i,ipert+1,2,ipert) = prho(i,ipert+1,2,ipert) - h
158            endif
159         enddo
160      enddo
161      if (GC) then
162c
163c     Perturb gamma values
164c
165         do i = 1, nq
166            Gaa = pdelrho(i,1,1,1,3)*pdelrho(i,1,1,1,3)
167     &          + pdelrho(i,2,1,1,3)*pdelrho(i,2,1,1,3)
168     &          + pdelrho(i,3,1,1,3)*pdelrho(i,3,1,1,3)
169            Gab = pdelrho(i,1,1,1,3)*pdelrho(i,1,2,1,3)
170     &          + pdelrho(i,2,1,1,3)*pdelrho(i,2,2,1,3)
171     &          + pdelrho(i,3,1,1,3)*pdelrho(i,3,2,1,3)
172            Gbb = pdelrho(i,1,2,1,3)*pdelrho(i,1,2,1,3)
173     &          + pdelrho(i,2,2,1,3)*pdelrho(i,2,2,1,3)
174     &          + pdelrho(i,3,2,1,3)*pdelrho(i,3,2,1,3)
175            mod_drhoa = sqrt(Gaa)
176            mod_drhob = sqrt(Gbb)
177c
178c     Ensure perturbed gammas are always > 0, otherwise take derivative as 0
179c
180            DoA = Gaa.gt.h
181            DoB = DoA.and.Gbb.gt.h
182            if (DoB) then
183               cos_theta = Gab/(mod_drhoa*mod_drhob)
184            else
185               cos_theta = 1d0
186            endif
187c
188c     Now that we have the vital parameters clear out locations
189c     in the perturbed density gradient array - we only need three
190c     locations to construct the proper perturbed gamma values
191c
192            pdelrho(i,1,1,1,3) = 0d0
193            pdelrho(i,2,1,1,3) = 0d0
194            pdelrho(i,3,1,1,3) = 0d0
195            pdelrho(i,1,2,1,3) = 0d0
196            pdelrho(i,2,2,1,3) = 0d0
197            pdelrho(i,3,2,1,3) = 0d0
198            pdelrho(i,1,1,1,4) = 0d0
199            pdelrho(i,2,1,1,4) = 0d0
200            pdelrho(i,3,1,1,4) = 0d0
201            pdelrho(i,1,2,1,4) = 0d0
202            pdelrho(i,2,2,1,4) = 0d0
203            pdelrho(i,3,2,1,4) = 0d0
204            pdelrho(i,1,1,1,5) = 0d0
205            pdelrho(i,2,1,1,5) = 0d0
206            pdelrho(i,3,1,1,5) = 0d0
207            pdelrho(i,1,2,1,5) = 0d0
208            pdelrho(i,2,2,1,5) = 0d0
209            pdelrho(i,3,2,1,5) = 0d0
210            if (DoA) then
211c
212c     Perturb Gaa up
213c
214               x = mod_drhoa*cos_theta
215               pdelrho(i,1,1,1,3) = x
216               pdelrho(i,2,1,1,3) = sqrt(Gaa-x*x+h)
217               pdelrho(i,1,2,1,3) = mod_drhob
218            endif
219            if (DoB) then
220c
221c     Perturb Gab up
222c
223               x = (Gab+h)/mod_drhob
224               pdelrho(i,1,1,1,4) = x
225               pdelrho(i,2,1,1,4) = sqrt(Gaa-x*x)
226               pdelrho(i,1,2,1,4) = mod_drhob
227c
228c     Perturb Gbb up
229c
230               mod_drhob_p = sqrt(Gbb+h)
231               cos_theta_p = (mod_drhob/mod_drhob_p)*cos_theta
232               sin_theta_p = sqrt(1d0-cos_theta_p*cos_theta_p)
233               pdelrho(i,1,1,1,5) = mod_drhoa*cos_theta_p
234               pdelrho(i,2,1,1,5) = mod_drhoa*sin_theta_p
235               pdelrho(i,1,2,1,5) = mod_drhob_p
236            endif
237            if (DoA) then
238c
239c     Perturb Gaa back
240c
241               x = mod_drhoa*cos_theta
242               pdelrho(i,1,1,2,3) = x
243               pdelrho(i,2,1,2,3) = sqrt(Gaa-x*x-h)
244               pdelrho(i,1,2,2,3) = mod_drhob
245            endif
246            if (DoB) then
247c
248c     Perturb Gab back
249c
250               x = (Gab-h)/mod_drhob
251               pdelrho(i,1,1,2,4) = x
252               pdelrho(i,2,1,2,4) = sqrt(Gaa-x*x)
253               pdelrho(i,1,2,2,4) = mod_drhob
254c
255c     Perturb Gbb back
256c
257               mod_drhob_p = sqrt(Gbb-h)
258               cos_theta_p = (mod_drhob/mod_drhob_p)*cos_theta
259               sin_theta_p = sqrt(1d0-cos_theta_p*cos_theta_p)
260               pdelrho(i,1,1,2,5) = mod_drhoa*cos_theta_p
261               pdelrho(i,2,1,2,5) = mod_drhoa*sin_theta_p
262               pdelrho(i,1,2,2,5) = mod_drhob_p
263            endif
264         enddo
265      endif
266c
267      return
268      end
269c
270      subroutine xc_make_fd(Amat2, Cmat2, nq, GC, pAmat, pCmat)
271c
272      implicit none
273c
274      integer nq
275      double precision Amat2(nq,*), Cmat2(nq,*),
276     &     pAmat(nq,2,2,*), pCmat(nq,3,2,*)
277      logical GC
278c
279c     Local variables
280c
281      double precision h, r2h
282      parameter (h = STEP_SIZE)
283      integer i, npert, ipert
284c
285#include "dft2drv.fh"
286c
287      if (GC) then
288         npert = 5
289      else
290         npert = 2
291      endif
292      r2h = 0.5d0*h
293c
294c     Construct finite differences in the temporary arrays
295c
296      do ipert = 1, npert
297         do i = 1, nq
298            pAmat(i,1,1,ipert) = (pAmat(i,1,1,ipert)
299     &                           -pAmat(i,1,2,ipert))*r2h
300            pAmat(i,2,1,ipert) = (pAmat(i,2,1,ipert)
301     &                           -pAmat(i,2,2,ipert))*r2h
302         enddo
303         if (GC) then
304            do i = 1, nq
305               pCmat(i,1,1,ipert) = (pCmat(i,1,1,ipert)
306     &                              -pCmat(i,1,2,ipert))*r2h
307               pCmat(i,2,1,ipert) = (pCmat(i,2,1,ipert)
308     &                              -pCmat(i,2,2,ipert))*r2h
309               pCmat(i,3,1,ipert) = (pCmat(i,3,1,ipert)
310     &                              -pCmat(i,3,2,ipert))*r2h
311            enddo
312         endif
313      enddo
314c
315c     Now scatter to Amat2 and Cmat2
316c
317      call dcopy(nq, pAmat(1,1,1,1), 1, Amat2(1,D2_RA_RA), 1)
318      call dcopy(nq, pAmat(1,2,1,1), 1, Amat2(1,D2_RA_RB), 1)
319c     Or could be as below, etc., etc.
320c      call dcopy(nq, pAmat(1,1,1,2), 1, Amat2(1,D2_RA_RB), 1)
321      call dcopy(nq, pAmat(1,2,1,2), 1, Amat2(1,D2_RB_RB), 1)
322      if (GC) then
323         call dcopy(nq, pCmat(1,D1_GAA,1,1), 1, Cmat2(1,D2_RA_GAA), 1)
324c     Or could be as below, etc., etc.
325c         call dcopy(nq, pAmat(1,1,1,3), 1, Cmat2(1,D2_RA_GAA), 1)
326         call dcopy(nq, pCmat(1,D1_GAA,1,2), 1, Cmat2(1,D2_RB_GAA), 1)
327         call dcopy(nq, pCmat(1,D1_GAA,1,3), 1, Cmat2(1,D2_GAA_GAA), 1)
328c
329         call dcopy(nq, pCmat(1,D1_GAB,1,1), 1, Cmat2(1,D2_RA_GAB), 1)
330         call dcopy(nq, pCmat(1,D1_GAB,1,2), 1, Cmat2(1,D2_RB_GAB), 1)
331         call dcopy(nq, pCmat(1,D1_GAB,1,3), 1, Cmat2(1,D2_GAA_GAB), 1)
332         call dcopy(nq, pCmat(1,D1_GAB,1,4), 1, Cmat2(1,D2_GAB_GAB), 1)
333c
334         call dcopy(nq, pCmat(1,D1_GBB,1,1), 1, Cmat2(1,D2_RA_GBB), 1)
335         call dcopy(nq, pCmat(1,D1_GBB,1,2), 1, Cmat2(1,D2_RB_GBB), 1)
336         call dcopy(nq, pCmat(1,D1_GBB,1,3), 1, Cmat2(1,D2_GAA_GBB), 1)
337         call dcopy(nq, pCmat(1,D1_GBB,1,4), 1, Cmat2(1,D2_GAB_GBB), 1)
338         call dcopy(nq, pCmat(1,D1_GBB,1,5), 1, Cmat2(1,D2_GBB_GBB), 1)
339      endif
340c
341      return
342      end
343