(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
9MODULE 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
49CONTAINS
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
469END 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
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_integrals
Definition atom_types.F:213
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:234
atom_types::atom_p_type
Definition atom_types.F:320
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:290
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126