1C
2C     zfock_cs_build.F
3C
4C     Builds closed shell complex Fock matrix from a complex density
5C     matrix.
6C
7C     A comment on notation.  The "z" in g_zfock means this GA is
8C     complex data type.  The corresponding real and imaginary parts are
9C     called g_fockre and g_fockim, which are both real data types, so
10C     no "z".
11C
12C     All quantities are in the atomic orbital (AO) basis.
13C
14
15      subroutine zfock_cs_build (params, g_zdens, energies, g_zfock)
16      implicit none
17
18#include "errquit.fh"
19#include "mafdecls.fh"
20#include "stdio.fh"
21#include "global.fh"
22#include "msgids.fh"
23#include "util.fh"
24#include "cdft.fh"
25#include "matutils.fh"
26#include "zora.fh"
27#include "case.fh"
28#include "rt_tddft.fh"
29
30
31C     == Inputs ==
32      type(rt_params_t), intent(in) :: params
33      integer, intent(in)           :: g_zdens
34
35
36C     == Outputs ==
37      type(rt_energies_t), intent(out) :: energies
38      integer, intent(in)              :: g_zfock
39
40
41C     == Parameters ==
42      character(*), parameter :: pname = "zfock_cs_build: "
43
44
45C     == Variables ==
46      logical need_exactexch
47      logical need_dftxc
48      integer g_densre          !real part of dens mat
49      integer g_densim          !imag part of dens mat
50      integer g_fockre          !real part of fock mat
51      integer g_fockim          !imag part of fock mat
52      integer g_v1, g_v2        !potentials--results of each piece of F
53      double precision ener1, ener2 !energies-results of each piece of F
54      double precision dftExc(2)
55      logical status
56      double precision elapsed, elapsed_tot
57      integer g_zfock_mo, g_zcapmo ! for MO damping
58
59C     == External routines ==
60      logical xc_gotxc
61      external xc_gotxc
62      external sandbox_complex
63
64
65      call rt_tddft_cs_confirm (params,'zfock_cs_build.F')
66
67
68C     == Initializations ==
69      if (params%prof) call prof_start (elapsed_tot)
70
71
72      if (params%nExc .ne. 1)
73     $     call errquit (pname//"might not work for nExc /= 1",0,0)
74
75      if (oskel)
76     $     call errquit (pname//"not working with symmetry", 0, 0)
77
78      if (nbf_ao .ne. params%nbf_ao) call errquit (pname//
79     $     "nbf_ao in params /= nbf_ao in cdft header", 0, 0)
80
81
82C
83C     If asked for (usually a "static" calculation), dont build but
84C     instead use stored Fock matrix.
85C
86      if (params%skip_fock) then
87         call ga_copy (params%g_zfock_ao_init(1), g_zfock)
88         call rt_tddft_print_notice ("Static CS Fock matrix")
89         return
90      endif
91
92
93      call ga_zero (g_zfock)
94      need_exactexch = abs(xfac(1)).gt.1d-8
95      need_dftxc = xc_gotxc()
96
97
98      if (.not. ga_create(mt_dbl, nbf_ao, nbf_ao,
99     $     "Re[Dens]", 0, 0, g_densre))
100     $     call errquit ("couldnt create g_densre", 0, GA_ERR)
101      if (.not. ga_duplicate(g_densre, g_densim, "Im[Dens]"))
102     $     call errquit ("couldnt duplicate g_densre", 0, GA_ERR)
103      if (.not. ga_duplicate(g_densre, g_fockre, "Re[F]"))
104     $     call errquit ("couldnt duplicate g_densre", 0, GA_ERR)
105      if (.not. ga_duplicate(g_densre, g_fockim, "Im[F]"))
106     $     call errquit ("couldnt duplicate g_densre", 0, GA_ERR)
107      if (.not. ga_duplicate(g_densre, g_v1, "V1"))
108     $     call errquit ("couldnt duplicate g_densre", 0, GA_ERR)
109      if (.not. ga_duplicate(g_densre, g_v2, "V2"))
110     $     call errquit ("couldnt duplicate g_densre", 0, GA_ERR)
111
112
113C     == Extract real and imag parts of density matrix ==
114      if (params%prof) call prof_start (elapsed)
115      call ga_zero (g_densre)
116      call ga_zero (g_densim)
117      call convert_z2d (g_zdens, g_densre, g_densim)
118      if (params%prof) call prof_end (elapsed, "Fock CS z2d")
119
120
121C
122C     == Compute complex Fock matrix ==
123C
124C     For each piece we compute the energy and potential, then
125C     accumulate the result in the real or imag part of F.  Note that
126C     the only piece that involves the imag part of the density matrix
127C     is the exact exchange.  We also only call the DFT XC routine if we
128C     need it, i.e., anything but pure Hartree-Fock.
129C
130      energies%core = 0d0
131      energies%coul = 0d0
132      energies%xc(1) = 0d0
133      energies%xc(2) = 0d0
134
135
136      call ga_zero (g_fockre)
137      call ga_zero (g_fockim)
138
139
140C     == Standard core (kinetic+potential) ==
141      call zfock_cs_core (params, g_densre, g_v1)
142      call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
143C      energies%core = ga_ddot (g_densre, g_v1)
144
145
146C     == Scalar ZORA core ==
147      if (params%lzorasf) then
148C         call zfock_cs_core_scalarzora (params, g_densre, g_v1)
149         call ga_add (1d0, params%g_zora_sf(1), 1d0, g_fockre, g_fockre)
150C         energies%core = energies%core + ga_ddot (g_densre, g_v1)
151      endif
152
153      energies%core = ga_ddot (g_densre, g_fockre)
154
155
156C     == Standard coul + standard exch ==
157      if ((need_exactexch).and.(.not.cam_exch).and.(.not.cdfit)) then
158         call zfock_cs_coul_exchre(params, g_densre, ener1, ener2, g_v1)
159         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
160         energies%coul = ener1
161         energies%xc(1) = ener2
162         energies%xc(2) = 0d0
163
164         call zfock_cs_exchim (params, g_densim, ener1, g_v1)
165         call ga_add (1d0, g_v1, 1d0, g_fockim, g_fockim)
166         energies%xc(1) = energies%xc(1) + ener1
167
168
169C     == Standard coul + CAM exchange ==
170      elseif ((need_exactexch).and.(cam_exch).and.(.not.cdfit)) then
171         call zfock_cs_coul (params, g_densre, ener1, g_v1)
172         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
173         energies%coul = ener1
174
175         call zfock_cs_exch (params, g_densre, g_densim, .true.,
176     $        ener1, g_v1, g_v2)
177         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
178         call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim)
179         energies%xc(1) = ener1
180         energies%xc(2) = 0d0
181
182
183C     == CD fit coul + standard exch ==
184      elseif ((need_exactexch).and.(.not. cam_exch).and.(cdfit)) then
185         call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1)
186         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
187         energies%coul = ener1
188
189         call zfock_cs_exch (params, g_densre, g_densim, .false.,
190     $        ener1, g_v1, g_v2)
191         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
192         call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim)
193         energies%xc(1) = ener1
194
195
196C     == CD fit coul + CAM exch ==
197      elseif ((need_exactexch).and.(cam_exch).and.(cdfit)) then
198         call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1)
199         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
200         energies%coul = ener1
201
202         call zfock_cs_exch (params, g_densre, g_densim, .true.,
203     $        ener1, g_v1, g_v2)
204         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
205         call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim)
206         energies%xc(1) = ener1
207         energies%xc(2) = 0d0
208
209
210C     == Standard coul + no exchange ==
211      elseif ((.not.need_exactexch).and.(.not.cdfit)) then
212         energies%xc(1) = 0d0
213         energies%xc(2) = 0d0
214
215         call zfock_cs_coul (params, g_densre, ener1, g_v1)
216         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
217         energies%coul = ener1
218
219
220C     == CD fit coul + no exchange ==
221      elseif ((.not.need_exactexch).and.(cdfit)) then
222         energies%xc(1) = 0d0
223         energies%xc(2) = 0d0
224
225         call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1)
226         call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre)
227         energies%coul = ener1
228
229      else
230         call errquit (pname//"Unsupported combination of flags", 0, 0)
231      endif
232
233
234C     == DFT XC ==
235      if (need_dftxc) then
236         dftExc(1) = 0d0
237         dftExc(2) = 0d0
238         call zfock_cs_dftxc (params, g_densre, dftExc, g_v1)
239         call ga_dadd (1d0, g_v1, 1d0, g_fockre, g_fockre)
240         energies%xc(1) = energies%xc(1) + dftExc(1)
241         energies%xc(2) = energies%xc(2) + dftExc(2)
242      endif
243
244
245C
246C     Fudge symmetry
247C
248      if (oskel) call errquit (pname//"not working with oskel", 0, 0)
249
250
251
252C     == Pack real and imag parts of Fock matrix back into g_zfock ==
253C     No need to repack density matrix, as its value should
254C     not have changed. Perhaps double check this??
255
256      if (params%prof) call prof_start (elapsed)
257      call convert_d2z (1d0, g_fockre, 1d0, g_fockim, g_zfock)
258      if (params%prof) call prof_end (elapsed, "Fock CS d2z")
259
260
261
262C
263C     Apply molecular orbital complex absorbing potential (MO CAP)
264C
265      if (params%mocap_active) then
266         if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
267     $        "Fock MO tmp" , 0, 0, g_zfock_mo))
268     $        call errquit(pname//"GA alloc failed",0, GA_ERR)
269
270         if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
271     $        "Fock MO tmp" , 0, 0, g_zcapmo))
272     $        call errquit(pname//"GA alloc failed",0, GA_ERR)
273
274         call canorg_trans (params,"F","AO->MO", g_zfock, g_zfock_mo)
275         call rt_tddft_mocap (params, g_zfock_mo, g_zcapmo)
276         call ga_zadd (z1, g_zfock_mo, z1, g_zcapmo, g_zfock_mo)
277         call canorg_trans (params,"F","MO->AO", g_zfock_mo, g_zfock)
278         call rt_tddft_print_notice ("Applied closed shell MO CAP")
279
280         if (.not. ga_destroy(g_zfock_mo))
281     $        call errquit(pname//"GA destroy failed",0, GA_ERR)
282
283         if (.not. ga_destroy(g_zcapmo))
284     $        call errquit(pname//"GA destroy failed",0, GA_ERR)
285      endif
286
287
288C     == Spatial complex absorbing boundary conditions (CAP) ==
289      if (params%cap_active) then
290         call ga_zadd (z1, g_zfock, z1, params%g_zcap, g_zfock)
291         call rt_tddft_print_notice ("Applied spatial CAP to "//
292     $        "restricted Fock matrix")
293      endif
294
295
296C
297C     Extra checks, disabled for speed.
298C
299      if (params%checklvl .ge. 2) then
300         if (.not. mat_is_hermitian (g_zfock, params%tol_zero))
301     $        call errquit ("F not hermitian in fock builder", 0, 0)
302         if (.not. mat_is_symmetric (g_fockim, "A", params%tol_zero))
303     $        call errquit ("Im[F] not antisym in fock builder", 0, 0)
304         if (energies%xc(2) > 1d-8)
305     $        call errquit (pname//"Exc(2) /= 0?", 0, 0)
306      endif
307
308      status = .true.
309      status=status.and.ga_destroy(g_densre)
310      status=status.and.ga_destroy(g_densim)
311      status=status.and.ga_destroy(g_fockre)
312      status=status.and.ga_destroy(g_fockim)
313      status=status.and.ga_destroy(g_v1)
314      status=status.and.ga_destroy(g_v2)
315
316      if (.not.status)
317     $     call errquit (pname//"couldnt free arrays", 0, GA_ERR)
318
319      if (params%prof) call prof_end (elapsed_tot,
320     $     "Fock CS total build")
321
322      end subroutine
323
324c $Id$
325