1!  Copyright (C) 2018 Andre Severo Pereira Gomes, Christoph Jacob, Lucas Visscher and collaborators
2!
3!  This file is part of Embed, a program implementing the Frozen Density Embedding (FDE) framework
4!
5!  This Source Code Form is subject to the terms of the Mozilla Public
6!  License, v. 2.0. If a copy of the MPL was not distributed with this
7!  file, You can obtain one at http://mozilla.org/MPL/2.0/.
8!
9
10module fde_nadd_derv
11
12      use fde_types
13
14      use fde_max_block_length
15      use fde_xcfun_interface
16
17      public fde_initialize_nadd_functionals
18      public fde_print_nadd_functionals
19
20      public get_fde_nadd_fun_derv
21      public fde_calc_xc_nadd_fun_derv
22      public fde_calc_ke_nadd_fun_derv
23
24      public fde_set_skip_nadd_xc
25      public fde_set_skip_nadd_ke
26      public fde_set_nadd_kef
27      public fde_set_nadd_xcf
28      public fde_set_nadd_all_alda
29
30      public fde_xc_fun_is_lda
31      public fde_xc_fun_is_gga
32      public fde_ke_fun_is_lda
33      public fde_ke_fun_is_gga
34      public is_fdersp_alda
35
36      logical,save :: skip_kin_contribution = .false.
37      logical,save :: skip_xc_contribution  = .false.
38      logical,save :: is_fdersp_alda        = .false.
39
40      type(functional), save :: fde_kef,       fde_xcf
41      type(functional), save :: fde_kef_alda,  fde_xcf_alda
42      type(functional), save :: fde_kef_xalda, fde_xcf_xalda
43
44! the default (orbital-free) kinetic energy functional for now is the thomas-fermi one
45      character(len=80), save :: kefun_fde='kin_tf'
46      character(len=80), save :: xcfun_fde='lda'
47!      character(len=80), save :: kefun_fde='kin_pw91'
48!      character(len=80), save :: xcfun_fde='pbe'
49
50      character(len=80), save :: kefun_fde_alda='kin_tf'
51      character(len=80), save :: xcfun_fde_alda='lda'
52
53      logical, save :: fde_set_functionals(20)
54      logical, save :: functionals_initialized = .false.
55
56      real(8) :: hfx_out, mu_out, beta_out
57
58! derv_lenght here is set to the appropriate parameter from xc_derv
59! looking at xcint, this would be done differently if using openrsp...
60      integer :: derv_length = nr_nonzero_derv
61
62      public fde_new_derv_result
63      public fde_delete_derv_result
64
65      interface fde_new_derv_result
66         module procedure fde_new_derv_result_partial
67      end interface
68
69      interface fde_delete_derv_result
70         module procedure fde_delete_derv_result_partial
71      end interface
72
73      contains
74
75
76!   ----------------------------------------------------------------------------
77      subroutine fde_new_derv_result_partial(npoints,derv)
78!   ----------------------------------------------------------------------------
79         real(kind=8), allocatable :: derv(:,:)
80         integer :: npoints
81
82         allocate(derv(npoints, 0:derv_length))
83         derv = 0.0d0
84      end subroutine fde_new_derv_result_partial
85
86
87!   ----------------------------------------------------------------------------
88      subroutine fde_del_derv_result_partial(derv)
89!   ----------------------------------------------------------------------------
90         real(kind=8), allocatable :: derv(:,:)
91
92         deallocate(derv)
93      end subroutine fde_del_derv_result_partial
94
95
96!   ----------------------------------------------------------------------------
97      subroutine calc_xcke_energy_density(gf,  derv)
98!   ----------------------------------------------------------------------------
99         type(grid_function), intent(in)  :: gf
100         real(kind=8),        intent(out) :: derv(gf%npoints,0:derv_length)
101         real(kind=8) :: dervtmp(max_block_length, 0:derv_length)
102
103         integer :: i
104         integer :: order = 0
105         integer :: bllen = 1
106         real(kind=8) :: r_0(1)
107         real(kind=8) :: z_0(1)
108         real(kind=8) :: w(1)
109
110         dervtmp = 0.0d0
111         w(1) = 1.0d0
112
113         do i = 1, gf%npoints
114            r_0(1) = gf%n(i)
115            z_0(1) = gf%gn(i,1)*gf%gn(i,1) &
116                   + gf%gn(i,2)*gf%gn(i,2) &
117                   + gf%gn(i,3)*gf%gn(i,3)
118
119            call get_xcke_fun_derv(order,        &
120                                   bllen,        &
121                                   w,            &
122                                   r_0,          &
123                                   z_0,          &
124                                   dervtmp)
125
126            derv(i,:) = dervtmp(1,:)
127         enddo
128      end subroutine calc_xcke_energy_density
129
130
131!   ----------------------------------------------------------------------------
132      subroutine fde_delete_derv_result_partial(derv)
133!   ----------------------------------------------------------------------------
134         real(kind=8), allocatable :: derv(:,:)
135
136         deallocate(derv)
137      end subroutine fde_delete_derv_result_partial
138
139
140!   ----------------------------------------------------------------------------
141      function fde_xc_fun_is_lda()
142!   ----------------------------------------------------------------------------
143      logical :: fde_xc_fun_is_lda
144
145         fde_xc_fun_is_lda = (xc_get_type(fde_xcf%id) == 0)
146      end function
147
148
149!   ----------------------------------------------------------------------------
150      function fde_xc_fun_is_gga()
151!   ----------------------------------------------------------------------------
152      logical :: fde_xc_fun_is_gga
153
154         fde_xc_fun_is_gga = (xc_get_type(fde_xcf%id) == 1)
155      end function
156
157
158!   ----------------------------------------------------------------------------
159      function fde_ke_fun_is_lda()
160!   ----------------------------------------------------------------------------
161      logical :: fde_ke_fun_is_lda
162
163         fde_ke_fun_is_lda = (xc_get_type(fde_kef%id) == 0)
164      end function
165
166
167!   ----------------------------------------------------------------------------
168      function fde_ke_fun_is_gga()
169!   ----------------------------------------------------------------------------
170      logical :: fde_ke_fun_is_gga
171
172         fde_ke_fun_is_gga = (xc_get_type(fde_kef%id) == 1)
173      end function
174
175
176!   ----------------------------------------------------------------------------
177         subroutine fde_set_skip_nadd_xc(value)
178!   ----------------------------------------------------------------------------
179            logical :: value
180            skip_xc_contribution  = value
181         end subroutine
182
183
184!   ----------------------------------------------------------------------------
185         subroutine fde_set_skip_nadd_ke(value)
186!   ----------------------------------------------------------------------------
187            logical :: value
188            skip_kin_contribution  = value
189         end subroutine
190
191
192!   ----------------------------------------------------------------------------
193         subroutine fde_initialize_nadd_functionals
194!   ----------------------------------------------------------------------------
195            integer :: i, order = 5
196            logical :: parallel_fde
197
198            parallel_fde = .false.
199
200        if (.not.functionals_initialized) then
201            functionals_initialized = .true.
202
203            call parse_functional(kefun_fde,      fde_kef,      hfx_out, mu_out, beta_out, .false.)
204            call parse_functional(xcfun_fde,      fde_xcf,      hfx_out, mu_out, beta_out, .false.)
205            call parse_functional(kefun_fde_alda, fde_kef_alda, hfx_out, mu_out, beta_out, .false.)
206            call parse_functional(xcfun_fde_alda, fde_xcf_alda, hfx_out, mu_out, beta_out, .false.)
207
208
209!            ! we set them now so that we can verify whether or not they
210            call xcfun_set_functional(fde_xcf     , order, parallel_fde)
211            call xcfun_set_functional(fde_kef,      order, parallel_fde)
212            call xcfun_set_functional(fde_xcf_alda, order, parallel_fde)
213            call xcfun_set_functional(fde_kef_alda, order, parallel_fde)
214        endif
215
216         end subroutine fde_initialize_nadd_functionals
217
218
219!   ----------------------------------------------------------------------------
220         subroutine fde_set_nadd_kef(id)
221!   ----------------------------------------------------------------------------
222            character(len=80) :: id
223
224            write(kefun_fde,'(A)') id
225         end subroutine fde_set_nadd_kef
226
227
228!   ----------------------------------------------------------------------------
229         subroutine fde_set_nadd_xcf(id)
230!   ----------------------------------------------------------------------------
231            character(len=80) :: id
232
233            write(xcfun_fde,'(A)') id
234         end subroutine fde_set_nadd_xcf
235
236
237!   ----------------------------------------------------------------------------
238         subroutine fde_set_nadd_all_alda
239!   ----------------------------------------------------------------------------
240             is_fdersp_alda = .true.
241         end subroutine fde_set_nadd_all_alda
242
243
244!   ----------------------------------------------------------------------------
245         subroutine fde_print_nadd_functionals(unit)
246!   ----------------------------------------------------------------------------
247            integer :: unit
248
249            write(unit,'(A,A20)') ' * Non-additive functional : (KE) ',kefun_fde
250            write(unit,'(A,A20)') '                             (XC) ',xcfun_fde
251            write(unit,'(A)') ' '
252            write(unit,'(A,A20)') '   if .RSPLDA set will use : (XC) ',xcfun_fde_alda
253            write(unit,'(A,A20)') '                             (KE) ',kefun_fde_alda
254            write(unit,'(A)') ' '
255         end subroutine fde_print_nadd_functionals
256
257
258!   ----------------------------------------------------------------------------
259         subroutine get_fde_nadd_fun_derv(order,        &
260                                          block_length, &
261                                          w,            &
262                                          r_0,          &
263                                          r_0_t,        &
264                                          z_0,          &
265                                          z_0_t,        &
266                                          derv)
267!   ----------------------------------------------------------------------------
268         integer,          intent(in)  :: order
269         integer,          intent(in)  :: block_length
270         real(8),          intent(in)  :: w(*)
271         real(8),          intent(in)  :: r_0(*)
272         real(8),          intent(in)  :: z_0(*)
273         real(8),          intent(in)  :: r_0_t(*)
274         real(8),          intent(in)  :: z_0_t(*)
275         real(8),          intent(inout) :: derv(max_block_length, 0:derv_length)
276
277          derv      = 0.0d0
278
279         if (.not.skip_xc_contribution) &
280            call get_nadd_fun_derv(fde_xcf,       &
281                                   fde_xcf_alda,  &
282                                   fde_xcf_alda,  &
283                                   order,         &
284                                   block_length,  &
285                                   w,             &
286                                   r_0,           &
287                                   r_0_t,         &
288                                   z_0,           &
289                                   z_0_t,         &
290                                   derv)
291
292         if (.not.skip_kin_contribution) &
293            call get_nadd_fun_derv(fde_kef,       &
294                                   fde_kef_alda,  &
295                                   fde_kef_alda,  &
296                                   order,         &
297                                   block_length,  &
298                                   w,             &
299                                   r_0,           &
300                                   r_0_t,         &
301                                   z_0,           &
302                                   z_0_t,         &
303                                   derv)
304
305         end subroutine get_fde_nadd_fun_derv
306
307
308!   ----------------------------------------------------------------------------
309         subroutine get_xc_nadd_fun_derv(order,        &
310                                         block_length, &
311                                         w,            &
312                                         r_0,          &
313                                         r_0_t,        &
314                                         z_0,          &
315                                         z_0_t,        &
316                                         derv)
317!   ----------------------------------------------------------------------------
318            integer,          intent(in)  :: order
319            integer,          intent(in)  :: block_length
320            real(8),          intent(in)  :: w(*)
321            real(8),          intent(in)  :: r_0(*)
322            real(8),          intent(in)  :: z_0(*)
323            real(8),          intent(in)  :: r_0_t(*)
324            real(8),          intent(in)  :: z_0_t(*)
325            real(8),          intent(inout) :: derv(max_block_length, 0:derv_length)
326
327            call get_nadd_fun_derv(fde_xcf,       &
328                                fde_xcf_alda,  &
329                                fde_xcf_alda,  &
330                                order,         &
331                                block_length,  &
332                                w,             &
333                                r_0,           &
334                                r_0_t,         &
335                                z_0,           &
336                                z_0_t,         &
337                                derv)
338
339         end subroutine get_xc_nadd_fun_derv
340
341
342!   ----------------------------------------------------------------------------
343         subroutine get_ke_nadd_fun_derv(order,        &
344                                      block_length, &
345                                      w,            &
346                                      r_0,          &
347                                      r_0_t,        &
348                                      z_0,          &
349                                      z_0_t,        &
350                                      derv)
351!   ----------------------------------------------------------------------------
352            integer, intent(in)  :: order
353            integer, intent(in)  :: block_length
354            real(8), intent(in)  :: w(*)
355            real(8), intent(in)  :: r_0(*)
356            real(8), intent(in)  :: z_0(*)
357            real(8), intent(in)  :: r_0_t(*)
358            real(8), intent(in)  :: z_0_t(*)
359            real(8), intent(inout) :: derv(max_block_length, 0:derv_length)
360
361            call get_nadd_fun_derv(fde_kef,       &
362                                fde_kef_alda,  &
363                                fde_kef_alda,  &
364                                order,         &
365                                block_length,  &
366                                w,             &
367                                r_0,           &
368                                r_0_t,         &
369                                z_0,           &
370                                z_0_t,         &
371                                derv)
372
373      end subroutine get_ke_nadd_fun_derv
374
375
376
377!   ----------------------------------------------------------------------------
378      subroutine get_nadd_fun_derv(f,            &
379                                    f_alda,       &
380                                    f_xalda,      &
381                                    order,        &
382                                    block_length, &
383                                    w,            &
384                                    r_0,          &
385                                    r_0_t,        &
386                                    z_0,          &
387                                    z_0_t,        &
388                                    derv)
389!   ----------------------------------------------------------------------------
390         type(functional)              :: f
391         type(functional)              :: f_alda, f_xalda
392         integer,          intent(in)  :: order
393         integer,          intent(in)  :: block_length
394         real(8),          intent(in)  :: w(*)
395         real(8),          intent(in)  :: r_0(*)
396         real(8),          intent(in)  :: z_0(*)
397         real(8),          intent(in)  :: r_0_t(*)
398         real(8),          intent(in)  :: z_0_t(*)
399         real(8),          intent(inout) :: derv(max_block_length, 0:derv_length)
400         real(8)                       :: dervtmp(max_block_length, 0:derv_length)
401         integer :: i, j
402
403         dervtmp = 0.0d0
404
405         call get_functional_derv(f,          &
406                               f_alda,        &
407                               f_xalda,       &
408                               order,         &
409                               block_length,  &
410                               w,             &
411                               r_0_t,         &
412                               z_0_t,         &
413                               derv_length,   &
414                               dervtmp,       &
415                               alda_real=is_fdersp_alda, &
416                               alda_imag=is_fdersp_alda)
417
418         do i = 1, max_block_length
419            do j = 0, derv_length
420               derv(i,j) = derv(i,j) + dervtmp(i,j)
421            enddo
422         enddo
423
424         call get_functional_derv(f,          &
425                               f_alda,        &
426                               f_xalda,       &
427                               order,         &
428                               block_length,  &
429                               w,             &
430                               r_0,           &
431                               z_0,           &
432                               derv_length,   &
433                               dervtmp,       &
434                               alda_real=is_fdersp_alda, &
435                               alda_imag=is_fdersp_alda)
436
437         do i = 1, max_block_length
438            do j = 0, derv_length
439               derv(i,j) = derv(i,j) - dervtmp(i,j)
440            enddo
441         enddo
442
443      end subroutine get_nadd_fun_derv
444
445
446!   ----------------------------------------------------------------------------
447      subroutine get_xcke_fun_derv(order,        &
448                                   block_length, &
449                                   w,            &
450                                   r_0,          &
451                                   z_0,          &
452                                   derv)
453!   ----------------------------------------------------------------------------
454         integer,          intent(in)  :: order
455         integer,          intent(in)  :: block_length
456         real(8),          intent(in)  :: w(*)
457         real(8),          intent(in)  :: r_0(*)
458         real(8),          intent(in)  :: z_0(*)
459         real(8),        intent(inout) :: derv(max_block_length, 0:derv_length)
460         integer                       :: i, j
461         real(8)                       :: dervtmp(max_block_length, 0:derv_length)
462
463         derv = 0.0d0
464         dervtmp = 0.0d0
465
466         if (.not.skip_xc_contribution) then
467         call get_functional_derv(fde_xcf,       &
468                                  fde_xcf_alda,  &
469                                  fde_xcf_alda,  &
470                                  order,         &
471                                  block_length,  &
472                                  w,             &
473                                  r_0,           &
474                                  z_0,           &
475                                  derv_length,   &
476                                  dervtmp,       &
477                                  alda_real=is_fdersp_alda, &
478                                  alda_imag=is_fdersp_alda)
479
480         do i = 1, max_block_length
481            do j = 0, derv_length
482               derv(i,j) = derv(i,j) + dervtmp(i,j)
483            enddo
484         enddo
485         endif
486
487         if (.not.skip_kin_contribution) then
488         call get_functional_derv(fde_kef,       &
489                                  fde_kef_alda,  &
490                                  fde_kef_alda,  &
491                                  order,         &
492                                  block_length,  &
493                                  w,             &
494                                  r_0,           &
495                                  z_0,           &
496                                  derv_length,   &
497                                  dervtmp,       &
498                                  alda_real=is_fdersp_alda, &
499                                  alda_imag=is_fdersp_alda)
500
501
502         do i = 1, max_block_length
503            do j = 0, derv_length
504               derv(i,j) = derv(i,j) + dervtmp(i,j)
505            enddo
506         enddo
507         endif
508
509     end subroutine get_xcke_fun_derv
510
511end module fde_nadd_derv
512
513