(git:1f285aa)
atom_admm_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 MODULE atom_admm_methods
10  USE atom_operators, ONLY: atom_basis_projection_overlap,&
11  atom_int_release,&
12  atom_int_setup
13  USE atom_output, ONLY: atom_print_basis
14  USE atom_types, ONLY: &
15  atom_basis_type, atom_integrals, atom_orbitals, atom_p_type, atom_type, create_atom_orbs, &
16  init_atom_basis, lmat, release_atom_basis, release_atom_orbs
17  USE atom_utils, ONLY: atom_consistent_method,&
18  atom_denmat,&
19  atom_trace,&
20  eeri_contract
21  USE atom_xc, ONLY: atom_dpot_lda,&
22  atom_vxc_lda,&
23  atom_vxc_lsd
24  USE input_constants, ONLY: do_rhf_atom,&
25  do_rks_atom,&
26  do_rohf_atom,&
27  do_uhf_atom,&
28  do_uks_atom,&
29  xc_funct_no_shortcut
30  USE input_section_types, ONLY: section_vals_duplicate,&
31  section_vals_get_subs_vals,&
32  section_vals_get_subs_vals2,&
33  section_vals_release,&
34  section_vals_remove_values,&
35  section_vals_type,&
36  section_vals_val_set
37  USE kinds, ONLY: dp
38  USE mathlib, ONLY: invmat_symm,&
39  jacobi
40 #include "./base/base_uses.f90"
41 
42  IMPLICIT NONE
43 
44  PRIVATE
45  PUBLIC :: atom_admm
46 
47  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_admm_methods'
48 
49 CONTAINS
50 
51 ! **************************************************************************************************
52 !> \brief Analysis of ADMM approximation to exact exchange.
53 !> \param atom_info information about the atomic kind. Two-dimensional array of size
54 !> (electronic-configuration, electronic-structure-method)
55 !> \param admm_section ADMM print section
56 !> \param iw output file unit
57 !> \par History
58 !> * 07.2016 created [Juerg Hutter]
59 ! **************************************************************************************************
60  SUBROUTINE atom_admm(atom_info, admm_section, iw)
61 
62  TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
63  TYPE(section_vals_type), POINTER :: admm_section
64  INTEGER, INTENT(IN) :: iw
65 
66  CHARACTER(LEN=2) :: btyp
67  INTEGER :: i, ifun, j, l, m, maxl, mo, n, na, nb, &
68  zval
69  LOGICAL :: pp_calc, rhf
70  REAL(kind=dp) :: admm1_k_energy, admm2_k_energy, admmq_k_energy, dfexc_admm1, dfexc_admm2, &
71  dfexc_admmq, dxc, dxk, el1, el2, elq, elref, fexc_optx_admm1, fexc_optx_admm2, &
72  fexc_optx_admmq, fexc_optx_ref, fexc_pbex_admm1, fexc_pbex_admm2, fexc_pbex_admmq, &
73  fexc_pbex_ref, ref_energy, xsi
74  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lamat
75  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: admm1_k, admm2_k, admm_xcmat, admm_xcmata, &
76  admm_xcmatb, admmq_k, ovlap, ref_k, ref_xcmat, ref_xcmata, ref_xcmatb, sinv, siref, tmat
77  TYPE(atom_basis_type), POINTER :: admm_basis, ref_basis
78  TYPE(atom_integrals), POINTER :: admm_int, ref_int
79  TYPE(atom_orbitals), POINTER :: admm1_orbs, admm2_orbs, admmq_orbs, &
80  ref_orbs
81  TYPE(atom_type), POINTER :: atom
82  TYPE(section_vals_type), POINTER :: basis_section, xc_fun, xc_fun_section, &
83  xc_optx, xc_pbex, xc_section
84 
85  IF (iw > 0) THEN
86  WRITE (iw, '(/,T2,A)') &
87  '!-----------------------------------------------------------------------------!'
88  WRITE (iw, '(T30,A)') "Analysis of ADMM approximations"
89  WRITE (iw, '(T2,A)') &
90  '!-----------------------------------------------------------------------------!'
91  END IF
92 
93  ! setup an xc section
94  NULLIFY (xc_section, xc_pbex, xc_optx)
95  CALL section_vals_duplicate(atom_info(1, 1)%atom%xc_section, xc_section)
96  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
97  ! Overwrite possible shortcut
98  CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
99  i_val=xc_funct_no_shortcut)
100  ifun = 0
101  DO
102  ifun = ifun + 1
103  xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
104  IF (.NOT. ASSOCIATED(xc_fun)) EXIT
105  CALL section_vals_remove_values(xc_fun)
106  END DO
107  ! PBEX
108  CALL section_vals_duplicate(xc_section, xc_pbex)
109  xc_fun_section => section_vals_get_subs_vals(xc_pbex, "XC_FUNCTIONAL")
110  CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", l_val=.true.)
111  CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", r_val=1.0_dp)
112  CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", r_val=0.0_dp)
113  ! OPTX
114  CALL section_vals_duplicate(xc_section, xc_optx)
115  xc_fun_section => section_vals_get_subs_vals(xc_optx, "XC_FUNCTIONAL")
116  CALL section_vals_val_set(xc_fun_section, "OPTX%_SECTION_PARAMETERS_", l_val=.true.)
117  CALL section_vals_val_set(xc_fun_section, "OPTX%SCALE_X", r_val=1.0_dp)
118 
119  ! ADMM basis set
120  zval = atom_info(1, 1)%atom%z
121  pp_calc = atom_info(1, 1)%atom%pp_calc
122  btyp = "AE"
123  IF (pp_calc) btyp = "PP"
124  ALLOCATE (admm_basis)
125  basis_section => section_vals_get_subs_vals(admm_section, "ADMM_BASIS")
126  NULLIFY (admm_basis%grid)
127  CALL init_atom_basis(admm_basis, basis_section, zval, btyp)
128  IF (iw > 0) THEN
129  CALL atom_print_basis(admm_basis, iw, " ADMM Basis")
130  END IF
131  ! reference basis set
132  ref_basis => atom_info(1, 1)%atom%basis
133  ! integrals
134  ALLOCATE (ref_int, admm_int)
135  CALL atom_int_setup(ref_int, ref_basis, eri_exchange=.true.)
136  CALL atom_int_setup(admm_int, admm_basis, eri_exchange=.true.)
137  DO l = 0, lmat
138  IF (admm_int%n(l) /= admm_int%nne(l)) THEN
139  IF (iw > 0) WRITE (iw, *) "ADMM Basis is linear dependent ", l, admm_int%n(l), admm_int%nne(l)
140  cpabort("ADMM basis is linear dependent")
141  END IF
142  END DO
143  ! mixed overlap
144  na = maxval(admm_basis%nbas(:))
145  nb = maxval(ref_basis%nbas(:))
146  ALLOCATE (ovlap(1:na, 1:nb, 0:lmat))
147  CALL atom_basis_projection_overlap(ovlap, admm_basis, ref_basis)
148  ! Inverse of ADMM overlap matrix
149  ALLOCATE (sinv(1:na, 1:na, 0:lmat))
150  DO l = 0, lmat
151  n = admm_basis%nbas(l)
152  IF (n < 1) cycle
153  sinv(1:n, 1:n, l) = admm_int%ovlp(1:n, 1:n, l)
154  CALL invmat_symm(sinv(1:n, 1:n, l))
155  END DO
156  ! ADMM transformation matrix
157  ALLOCATE (tmat(1:na, 1:nb, 0:lmat))
158  DO l = 0, lmat
159  n = admm_basis%nbas(l)
160  m = ref_basis%nbas(l)
161  IF (n < 1 .OR. m < 1) cycle
162  tmat(1:n, 1:m, l) = matmul(sinv(1:n, 1:n, l), ovlap(1:n, 1:m, l))
163  END DO
164  ! Inverse of REF overlap matrix
165  ALLOCATE (siref(1:nb, 1:nb, 0:lmat))
166  DO l = 0, lmat
167  n = ref_basis%nbas(l)
168  IF (n < 1) cycle
169  siref(1:n, 1:n, l) = ref_int%ovlp(1:n, 1:n, l)
170  CALL invmat_symm(siref(1:n, 1:n, l))
171  END DO
172 
173  DO i = 1, SIZE(atom_info, 1)
174  DO j = 1, SIZE(atom_info, 2)
175  atom => atom_info(i, j)%atom
176  IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN
177  ref_orbs => atom%orbitals
178  ALLOCATE (ref_k(1:nb, 1:nb, 0:lmat))
179  SELECT CASE (atom%method_type)
180  CASE (do_rks_atom, do_rhf_atom)
181  ! restricted
182  rhf = .true.
183  ref_k = 0.0_dp
184  CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmat, ref_basis%nbas)
185  ref_energy = 0.5_dp*atom_trace(ref_k, ref_orbs%pmat)
186  CASE (do_uks_atom, do_uhf_atom)
187  ! unrestricted
188  rhf = .false.
189  ref_k = 0.0_dp
190  CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmata, ref_basis%nbas)
191  ref_energy = atom_trace(ref_k, ref_orbs%pmata)
192  ref_k = 0.0_dp
193  CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmatb, ref_basis%nbas)
194  ref_energy = ref_energy + atom_trace(ref_k, ref_orbs%pmatb)
195  CASE (do_rohf_atom)
196  cpabort("ADMM not available")
197  CASE DEFAULT
198  cpabort("ADMM not available")
199  END SELECT
200  DEALLOCATE (ref_k)
201  ! reference number of electrons
202  elref = atom_trace(ref_int%ovlp, ref_orbs%pmat)
203  ! admm orbitals and density matrices
204  mo = maxval(atom%state%maxn_calc)
205  NULLIFY (admm1_orbs, admm2_orbs, admmq_orbs)
206  CALL create_atom_orbs(admm1_orbs, na, mo)
207  CALL create_atom_orbs(admm2_orbs, na, mo)
208  CALL create_atom_orbs(admmq_orbs, na, mo)
209  ALLOCATE (lamat(1:mo, 1:mo))
210  ALLOCATE (admm1_k(1:na, 1:na, 0:lmat))
211  ALLOCATE (admm2_k(1:na, 1:na, 0:lmat))
212  ALLOCATE (admmq_k(1:na, 1:na, 0:lmat))
213  IF (rhf) THEN
214  DO l = 0, lmat
215  n = admm_basis%nbas(l)
216  m = ref_basis%nbas(l)
217  mo = atom%state%maxn_calc(l)
218  IF (n < 1 .OR. m < 1 .OR. mo < 1) cycle
219  admm2_orbs%wfn(1:n, 1:mo, l) = matmul(tmat(1:n, 1:m, l), ref_orbs%wfn(1:m, 1:mo, l))
220  CALL lowdin_matrix(admm2_orbs%wfn(1:n, 1:mo, l), lamat(1:mo, 1:mo), admm_int%ovlp(1:n, 1:n, l))
221  admm1_orbs%wfn(1:n, 1:mo, l) = matmul(admm2_orbs%wfn(1:n, 1:mo, l), lamat(1:mo, 1:mo))
222  END DO
223  CALL atom_denmat(admm1_orbs%pmat, admm1_orbs%wfn, admm_basis%nbas, atom%state%occupation, &
224  atom%state%maxl_occ, atom%state%maxn_occ)
225  CALL atom_denmat(admm2_orbs%pmat, admm2_orbs%wfn, admm_basis%nbas, atom%state%occupation, &
226  atom%state%maxl_occ, atom%state%maxn_occ)
227  el1 = atom_trace(admm_int%ovlp, admm1_orbs%pmat)
228  el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmat)
229  xsi = elref/el2
230  admmq_orbs%pmat = xsi*admm2_orbs%pmat
231  elq = atom_trace(admm_int%ovlp, admmq_orbs%pmat)
232  admmq_orbs%wfn = sqrt(xsi)*admm2_orbs%wfn
233  !
234  admm1_k = 0.0_dp
235  CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmat, admm_basis%nbas)
236  admm1_k_energy = 0.5_dp*atom_trace(admm1_k, admm1_orbs%pmat)
237  admm2_k = 0.0_dp
238  CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmat, admm_basis%nbas)
239  admm2_k_energy = 0.5_dp*atom_trace(admm2_k, admm2_orbs%pmat)
240  admmq_k = 0.0_dp
241  CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmat, admm_basis%nbas)
242  admmq_k_energy = 0.5_dp*atom_trace(admmq_k, admmq_orbs%pmat)
243  ELSE
244  DO l = 0, lmat
245  n = admm_basis%nbas(l)
246  m = ref_basis%nbas(l)
247  mo = atom%state%maxn_calc(l)
248  IF (n < 1 .OR. m < 1 .OR. mo < 1) cycle
249  admm2_orbs%wfna(1:n, 1:mo, l) = matmul(tmat(1:n, 1:m, l), ref_orbs%wfna(1:m, 1:mo, l))
250  CALL lowdin_matrix(admm2_orbs%wfna(1:n, 1:mo, l), lamat, admm_int%ovlp(1:n, 1:n, l))
251  admm1_orbs%wfna(1:n, 1:mo, l) = matmul(admm2_orbs%wfna(1:n, 1:mo, l), lamat(1:mo, 1:mo))
252  admm2_orbs%wfnb(1:n, 1:mo, l) = matmul(tmat(1:n, 1:m, l), ref_orbs%wfnb(1:m, 1:mo, l))
253  CALL lowdin_matrix(admm2_orbs%wfnb(1:n, 1:mo, l), lamat, admm_int%ovlp(1:n, 1:n, l))
254  admm1_orbs%wfnb(1:n, 1:mo, l) = matmul(admm2_orbs%wfnb(1:n, 1:mo, l), lamat(1:mo, 1:mo))
255  END DO
256  CALL atom_denmat(admm1_orbs%pmata, admm1_orbs%wfna, admm_basis%nbas, atom%state%occa, &
257  atom%state%maxl_occ, atom%state%maxn_occ)
258  CALL atom_denmat(admm1_orbs%pmatb, admm1_orbs%wfnb, admm_basis%nbas, atom%state%occb, &
259  atom%state%maxl_occ, atom%state%maxn_occ)
260  admm1_orbs%pmat = admm1_orbs%pmata + admm1_orbs%pmatb
261  CALL atom_denmat(admm2_orbs%pmata, admm2_orbs%wfna, admm_basis%nbas, atom%state%occa, &
262  atom%state%maxl_occ, atom%state%maxn_occ)
263  CALL atom_denmat(admm2_orbs%pmatb, admm2_orbs%wfnb, admm_basis%nbas, atom%state%occb, &
264  atom%state%maxl_occ, atom%state%maxn_occ)
265  admm2_orbs%pmat = admm2_orbs%pmata + admm2_orbs%pmatb
266  elref = atom_trace(ref_int%ovlp, ref_orbs%pmata)
267  el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmata)
268  xsi = elref/el2
269  admmq_orbs%pmata = xsi*admm2_orbs%pmata
270  admmq_orbs%wfna = sqrt(xsi)*admm2_orbs%wfna
271  elref = atom_trace(ref_int%ovlp, ref_orbs%pmatb)
272  el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmatb)
273  xsi = elref/el2
274  admmq_orbs%pmatb = xsi*admm2_orbs%pmatb
275  admmq_orbs%wfnb = sqrt(xsi)*admm2_orbs%wfnb
276  admmq_orbs%pmat = admmq_orbs%pmata + admmq_orbs%pmatb
277  el1 = atom_trace(admm_int%ovlp, admm1_orbs%pmat)
278  el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmat)
279  elq = atom_trace(admm_int%ovlp, admmq_orbs%pmat)
280  elref = atom_trace(ref_int%ovlp, ref_orbs%pmat)
281  !
282  admm1_k = 0.0_dp
283  CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmata, admm_basis%nbas)
284  admm1_k_energy = atom_trace(admm1_k, admm1_orbs%pmata)
285  admm1_k = 0.0_dp
286  CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmatb, admm_basis%nbas)
287  admm1_k_energy = admm1_k_energy + atom_trace(admm1_k, admm1_orbs%pmatb)
288  admm2_k = 0.0_dp
289  CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmata, admm_basis%nbas)
290  admm2_k_energy = atom_trace(admm2_k, admm2_orbs%pmata)
291  admm2_k = 0.0_dp
292  CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmatb, admm_basis%nbas)
293  admm2_k_energy = admm2_k_energy + atom_trace(admm2_k, admm2_orbs%pmatb)
294  admmq_k = 0.0_dp
295  CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmata, admm_basis%nbas)
296  admmq_k_energy = atom_trace(admmq_k, admmq_orbs%pmata)
297  admmq_k = 0.0_dp
298  CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmatb, admm_basis%nbas)
299  admmq_k_energy = admmq_k_energy + atom_trace(admmq_k, admmq_orbs%pmatb)
300  END IF
301  DEALLOCATE (lamat)
302  !
303  ! ADMM correction terms
304  maxl = atom%state%maxl_occ
305  IF (rhf) THEN
306  ALLOCATE (ref_xcmat(1:nb, 1:nb, 0:lmat))
307  ALLOCATE (admm_xcmat(1:na, 1:na, 0:lmat))
308  ! PBEX
309  CALL atom_vxc_lda(ref_basis, ref_orbs%pmat, maxl, xc_pbex, fexc_pbex_ref, ref_xcmat)
310  CALL atom_vxc_lda(admm_basis, admm1_orbs%pmat, maxl, xc_pbex, fexc_pbex_admm1, admm_xcmat)
311  CALL atom_vxc_lda(admm_basis, admm2_orbs%pmat, maxl, xc_pbex, fexc_pbex_admm2, admm_xcmat)
312  CALL atom_vxc_lda(admm_basis, admmq_orbs%pmat, maxl, xc_pbex, fexc_pbex_admmq, admm_xcmat)
313  ! OPTX
314  CALL atom_vxc_lda(ref_basis, ref_orbs%pmat, maxl, xc_optx, fexc_optx_ref, ref_xcmat)
315  CALL atom_vxc_lda(admm_basis, admm1_orbs%pmat, maxl, xc_optx, fexc_optx_admm1, admm_xcmat)
316  CALL atom_vxc_lda(admm_basis, admm2_orbs%pmat, maxl, xc_optx, fexc_optx_admm2, admm_xcmat)
317  CALL atom_vxc_lda(admm_basis, admmq_orbs%pmat, maxl, xc_optx, fexc_optx_admmq, admm_xcmat)
318  ! LINX
319  CALL atom_dpot_lda(ref_basis, ref_orbs%pmat, admm_basis, admm1_orbs%pmat, &
320  maxl, "LINX", dfexc_admm1)
321  CALL atom_dpot_lda(ref_basis, ref_orbs%pmat, admm_basis, admm2_orbs%pmat, &
322  maxl, "LINX", dfexc_admm2)
323  CALL atom_dpot_lda(ref_basis, ref_orbs%pmat, admm_basis, admmq_orbs%pmat, &
324  maxl, "LINX", dfexc_admmq)
325  DEALLOCATE (ref_xcmat, admm_xcmat)
326  ELSE
327  ALLOCATE (ref_xcmata(1:nb, 1:nb, 0:lmat))
328  ALLOCATE (ref_xcmatb(1:nb, 1:nb, 0:lmat))
329  ALLOCATE (admm_xcmata(1:na, 1:na, 0:lmat))
330  ALLOCATE (admm_xcmatb(1:na, 1:na, 0:lmat))
331  ! PBEX
332  CALL atom_vxc_lsd(ref_basis, ref_orbs%pmata, ref_orbs%pmatb, maxl, xc_pbex, fexc_pbex_ref, &
333  ref_xcmata, ref_xcmatb)
334  CALL atom_vxc_lsd(admm_basis, admm1_orbs%pmata, admm1_orbs%pmatb, maxl, xc_pbex, &
335  fexc_pbex_admm1, admm_xcmata, admm_xcmatb)
336  CALL atom_vxc_lsd(admm_basis, admm2_orbs%pmata, admm2_orbs%pmatb, maxl, xc_pbex, &
337  fexc_pbex_admm2, admm_xcmata, admm_xcmatb)
338  CALL atom_vxc_lsd(admm_basis, admmq_orbs%pmata, admmq_orbs%pmatb, maxl, xc_pbex, &
339  fexc_pbex_admmq, admm_xcmata, admm_xcmatb)
340  CALL atom_vxc_lsd(ref_basis, ref_orbs%pmata, ref_orbs%pmatb, maxl, xc_optx, fexc_optx_ref, &
341  ref_xcmata, ref_xcmatb)
342  CALL atom_vxc_lsd(admm_basis, admm1_orbs%pmata, admm1_orbs%pmatb, maxl, xc_optx, &
343  fexc_optx_admm1, admm_xcmata, admm_xcmatb)
344  CALL atom_vxc_lsd(admm_basis, admm2_orbs%pmata, admm2_orbs%pmatb, maxl, xc_optx, &
345  fexc_optx_admm2, admm_xcmata, admm_xcmatb)
346  CALL atom_vxc_lsd(admm_basis, admmq_orbs%pmata, admmq_orbs%pmatb, maxl, xc_optx, &
347  fexc_optx_admmq, admm_xcmata, admm_xcmatb)
348  DEALLOCATE (ref_xcmata, ref_xcmatb, admm_xcmata, admm_xcmatb)
349  END IF
350  !
351 
352  IF (iw > 0) THEN
353  WRITE (iw, "(/,A,I3,T48,A,I3)") " Electronic Structure Setting: ", i, &
354  " Electronic Structure Method: ", j
355  WRITE (iw, "(' Norm of ADMM Basis projection ',T61,F20.10)") el2/elref
356  WRITE (iw, "(' Reference Exchange Energy [Hartree]',T61,F20.10)") ref_energy
357  ! ADMM1
358  dxk = ref_energy - admm1_k_energy
359  WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMM1 METHOD: Energy ", admm1_k_energy, &
360  " Error: ", dxk
361  dxc = fexc_pbex_ref - fexc_pbex_admm1
362  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
363  " Error: ", dxk - dxc
364  dxc = fexc_optx_ref - fexc_optx_admm1
365  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
366  " Error: ", dxk - dxc
367  dxc = dfexc_admm1
368  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
369  " Error: ", dxk - dxc
370  ! ADMM2
371  dxk = ref_energy - admm2_k_energy
372  WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMM2 METHOD: Energy ", admm2_k_energy, &
373  " Error: ", dxk
374  dxc = fexc_pbex_ref - fexc_pbex_admm2
375  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
376  " Error: ", dxk - dxc
377  dxc = fexc_optx_ref - fexc_optx_admm2
378  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
379  " Error: ", dxk - dxc
380  dxc = dfexc_admm2
381  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
382  " Error: ", dxk - dxc
383  ! ADMMQ
384  dxk = ref_energy - admmq_k_energy
385  WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMMQ METHOD: Energy ", admmq_k_energy, &
386  " Error: ", dxk
387  dxc = fexc_pbex_ref - fexc_pbex_admmq
388  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
389  " Error: ", dxk - dxc
390  dxc = fexc_optx_ref - fexc_optx_admmq
391  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
392  " Error: ", dxk - dxc
393  dxc = dfexc_admmq
394  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
395  " Error: ", dxk - dxc
396  ! ADMMS
397  dxk = ref_energy - admmq_k_energy
398  WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMMS METHOD: Energy ", admmq_k_energy, &
399  " Error: ", dxk
400  dxc = fexc_pbex_ref - fexc_pbex_admmq*xsi**(2._dp/3._dp)
401  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
402  " Error: ", dxk - dxc
403  dxc = fexc_optx_ref - fexc_optx_admmq*xsi**(2._dp/3._dp)
404  WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
405  " Error: ", dxk - dxc
406  END IF
407  !
408  DEALLOCATE (admm1_k, admm2_k, admmq_k)
409  !
410  CALL release_atom_orbs(admm1_orbs)
411  CALL release_atom_orbs(admm2_orbs)
412  CALL release_atom_orbs(admmq_orbs)
413  END IF
414  END DO
415  END DO
416 
417  ! clean up
418  CALL atom_int_release(ref_int)
419  CALL atom_int_release(admm_int)
420  CALL release_atom_basis(admm_basis)
421  DEALLOCATE (ref_int, admm_int, admm_basis)
422  DEALLOCATE (ovlap, sinv, tmat, siref)
423 
424  CALL section_vals_release(xc_pbex)
425  CALL section_vals_release(xc_optx)
426  CALL section_vals_release(xc_section)
427 
428  IF (iw > 0) THEN
429  WRITE (iw, '(/,T2,A)') &
430  '!------------------------------End of ADMM analysis---------------------------!'
431  END IF
432 
433  END SUBROUTINE atom_admm
434 
435 ! **************************************************************************************************
436 !> \brief ...
437 !> \param wfn ...
438 !> \param lamat ...
439 !> \param ovlp ...
440 ! **************************************************************************************************
441  SUBROUTINE lowdin_matrix(wfn, lamat, ovlp)
442  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: wfn
443  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: lamat
444  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ovlp
445 
446  INTEGER :: i, j, k, n
447  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w
448  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vmat
449 
450  n = SIZE(wfn, 2)
451  IF (n > 0) THEN
452  lamat = matmul(transpose(wfn), matmul(ovlp, wfn))
453  ALLOCATE (w(n), vmat(n, n))
454  CALL jacobi(lamat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
455  w(1:n) = 1.0_dp/sqrt(w(1:n))
456  DO i = 1, n
457  DO j = 1, n
458  lamat(i, j) = 0.0_dp
459  DO k = 1, n
460  lamat(i, j) = lamat(i, j) + vmat(i, k)*w(k)*vmat(j, k)
461  END DO
462  END DO
463  END DO
464  DEALLOCATE (vmat, w)
465  END IF
466 
467  END SUBROUTINE lowdin_matrix
468 
469 END MODULE atom_admm_methods
atom_admm_methods
Definition: atom_admm_methods.F:9
atom_admm_methods::atom_admm
subroutine, public atom_admm(atom_info, admm_section, iw)
Analysis of ADMM approximation to exact exchange.
Definition: atom_admm_methods.F:61
atom_operators
Calculate the atomic operator matrices.
Definition: atom_operators.F:15
atom_operators::atom_int_setup
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
Definition: atom_operators.F:64
atom_operators::atom_basis_projection_overlap
subroutine, public atom_basis_projection_overlap(ovlap, basis_a, basis_b)
Calculate overlap matrix between two atomic basis sets.
Definition: atom_operators.F:1036
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition: atom_operators.F:1108
atom_output
Routines that print various information about an atomic kind.
Definition: atom_output.F:11
atom_output::atom_print_basis
subroutine, public atom_print_basis(atom_basis, iw, title)
Print atomic basis set.
Definition: atom_output.F:293
atom_types
Define the atom type and its sub types.
Definition: atom_types.F:15
atom_types::init_atom_basis
subroutine, public init_atom_basis(basis, basis_section, zval, btyp)
Initialize the basis for the atomic code.
Definition: atom_types.F:375
atom_types::release_atom_orbs
subroutine, public release_atom_orbs(orbs)
...
Definition: atom_types.F:1104
atom_types::lmat
integer, parameter, public lmat
Definition: atom_types.F:67
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition: atom_types.F:910
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition: atom_types.F:1055
atom_utils
Some basic routines for atomic calculations.
Definition: atom_utils.F:15
atom_utils::atom_denmat
pure subroutine, public atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
Calculate a density matrix using atomic orbitals.
Definition: atom_utils.F:328
atom_utils::eeri_contract
pure subroutine, public eeri_contract(kmat, erint, pmat, nsize)
Contract exchange Electron Repulsion Integrals.
Definition: atom_utils.F:1974
atom_utils::atom_consistent_method
pure logical function, public atom_consistent_method(method, multiplicity)
Check that the atomic multiplicity is consistent with the electronic structure method.
Definition: atom_utils.F:2189
atom_utils::atom_trace
pure real(kind=dp) function, public atom_trace(opmat, pmat)
Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
Definition: atom_utils.F:1823
atom_xc
routines that build the integrals of the Vxc potential calculated for the atomic code
Definition: atom_xc.F:12
atom_xc::atom_dpot_lda
subroutine, public atom_dpot_lda(basis0, pmat0, basis1, pmat1, maxl, functional, dfexc, linxpar)
Estimate the ADMM exchange energy correction term using a model exchange functional.
Definition: atom_xc.F:979
atom_xc::atom_vxc_lsd
subroutine, public atom_vxc_lsd(basis, pmata, pmatb, maxl, xc_section, fexc, xcmata, xcmatb)
Alternative subroutine that calculates atomic exchange-correlation potential in spin-polarized case.
Definition: atom_xc.F:800
atom_xc::atom_vxc_lda
subroutine, public atom_vxc_lda(basis, pmat, maxl, xc_section, fexc, xcmat)
Alternative subroutine that calculates atomic exchange-correlation potential in spin-restricted case.
Definition: atom_xc.F:649
atom
Definition: atom.F:9
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_rhf_atom
integer, parameter, public do_rhf_atom
Definition: input_constants.F:159
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition: input_constants.F:159
input_constants::xc_funct_no_shortcut
integer, parameter, public xc_funct_no_shortcut
Definition: input_constants.F:551
input_constants::do_uhf_atom
integer, parameter, public do_uhf_atom
Definition: input_constants.F:159
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition: input_constants.F:159
input_constants::do_rohf_atom
integer, parameter, public do_rohf_atom
Definition: input_constants.F:159
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
input_section_types::section_vals_get_subs_vals2
type(section_vals_type) function, pointer, public section_vals_get_subs_vals2(section_vals, i_section, i_rep_section)
returns the values of the n-th non default subsection (null if no such section exists (not so many no...
Definition: input_section_types.F:779
input_section_types::section_vals_remove_values
subroutine, public section_vals_remove_values(section_vals)
removes the values of a repetition of the section
Definition: input_section_types.F:889
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_duplicate
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
Definition: input_section_types.F:1732
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition: input_section_types.F:646
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition: mathlib.F:1479
mathlib::invmat_symm
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:574
xc_optx
calculate optx
Definition: xc_optx.F:14