Line data Source code
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 1 : 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 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: lamat
75 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: admm1_k, admm2_k, admm_xcmat, admm_xcmata, &
76 1 : 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 1 : IF (iw > 0) THEN
86 : WRITE (iw, '(/,T2,A)') &
87 1 : '!-----------------------------------------------------------------------------!'
88 1 : WRITE (iw, '(T30,A)') "Analysis of ADMM approximations"
89 : WRITE (iw, '(T2,A)') &
90 1 : '!-----------------------------------------------------------------------------!'
91 : END IF
92 :
93 : ! setup an xc section
94 1 : NULLIFY (xc_section, xc_pbex, xc_optx)
95 1 : CALL section_vals_duplicate(atom_info(1, 1)%atom%xc_section, xc_section)
96 1 : 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 1 : i_val=xc_funct_no_shortcut)
100 1 : ifun = 0
101 1 : DO
102 2 : ifun = ifun + 1
103 2 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
104 2 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
105 1 : CALL section_vals_remove_values(xc_fun)
106 : END DO
107 : ! PBEX
108 1 : CALL section_vals_duplicate(xc_section, xc_pbex)
109 1 : xc_fun_section => section_vals_get_subs_vals(xc_pbex, "XC_FUNCTIONAL")
110 1 : CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", l_val=.TRUE.)
111 1 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", r_val=1.0_dp)
112 1 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", r_val=0.0_dp)
113 : ! OPTX
114 1 : CALL section_vals_duplicate(xc_section, xc_optx)
115 1 : xc_fun_section => section_vals_get_subs_vals(xc_optx, "XC_FUNCTIONAL")
116 1 : CALL section_vals_val_set(xc_fun_section, "OPTX%_SECTION_PARAMETERS_", l_val=.TRUE.)
117 1 : CALL section_vals_val_set(xc_fun_section, "OPTX%SCALE_X", r_val=1.0_dp)
118 :
119 : ! ADMM basis set
120 1 : zval = atom_info(1, 1)%atom%z
121 1 : pp_calc = atom_info(1, 1)%atom%pp_calc
122 1 : btyp = "AE"
123 1 : IF (pp_calc) btyp = "PP"
124 19 : ALLOCATE (admm_basis)
125 1 : basis_section => section_vals_get_subs_vals(admm_section, "ADMM_BASIS")
126 1 : NULLIFY (admm_basis%grid)
127 1 : CALL init_atom_basis(admm_basis, basis_section, zval, btyp)
128 1 : IF (iw > 0) THEN
129 1 : CALL atom_print_basis(admm_basis, iw, " ADMM Basis")
130 : END IF
131 : ! reference basis set
132 1 : ref_basis => atom_info(1, 1)%atom%basis
133 : ! integrals
134 425 : ALLOCATE (ref_int, admm_int)
135 1 : CALL atom_int_setup(ref_int, ref_basis, eri_exchange=.TRUE.)
136 1 : CALL atom_int_setup(admm_int, admm_basis, eri_exchange=.TRUE.)
137 7 : DO l = 0, lmat
138 7 : IF (admm_int%n(l) /= admm_int%nne(l)) THEN
139 0 : IF (iw > 0) WRITE (iw, *) "ADMM Basis is linear dependent ", l, admm_int%n(l), admm_int%nne(l)
140 0 : CPABORT("ADMM basis is linear dependent")
141 : END IF
142 : END DO
143 : ! mixed overlap
144 7 : na = MAXVAL(admm_basis%nbas(:))
145 7 : nb = MAXVAL(ref_basis%nbas(:))
146 5 : ALLOCATE (ovlap(1:na, 1:nb, 0:lmat))
147 1 : CALL atom_basis_projection_overlap(ovlap, admm_basis, ref_basis)
148 : ! Inverse of ADMM overlap matrix
149 5 : ALLOCATE (sinv(1:na, 1:na, 0:lmat))
150 7 : DO l = 0, lmat
151 6 : n = admm_basis%nbas(l)
152 6 : IF (n < 1) CYCLE
153 20 : sinv(1:n, 1:n, l) = admm_int%ovlp(1:n, 1:n, l)
154 7 : CALL invmat_symm(sinv(1:n, 1:n, l))
155 : END DO
156 : ! ADMM transformation matrix
157 3 : ALLOCATE (tmat(1:na, 1:nb, 0:lmat))
158 7 : DO l = 0, lmat
159 6 : n = admm_basis%nbas(l)
160 6 : m = ref_basis%nbas(l)
161 6 : IF (n < 1 .OR. m < 1) CYCLE
162 132 : 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 5 : ALLOCATE (siref(1:nb, 1:nb, 0:lmat))
166 7 : DO l = 0, lmat
167 6 : n = ref_basis%nbas(l)
168 6 : IF (n < 1) CYCLE
169 72 : siref(1:n, 1:n, l) = ref_int%ovlp(1:n, 1:n, l)
170 7 : CALL invmat_symm(siref(1:n, 1:n, l))
171 : END DO
172 :
173 2 : DO i = 1, SIZE(atom_info, 1)
174 3 : DO j = 1, SIZE(atom_info, 2)
175 1 : atom => atom_info(i, j)%atom
176 2 : IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN
177 1 : ref_orbs => atom%orbitals
178 3 : ALLOCATE (ref_k(1:nb, 1:nb, 0:lmat))
179 1 : SELECT CASE (atom%method_type)
180 : CASE (do_rks_atom, do_rhf_atom)
181 : ! restricted
182 7 : rhf = .TRUE.
183 187 : ref_k = 0.0_dp
184 1 : CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmat, ref_basis%nbas)
185 1 : ref_energy = 0.5_dp*atom_trace(ref_k, ref_orbs%pmat)
186 : CASE (do_uks_atom, do_uhf_atom)
187 : ! unrestricted
188 0 : rhf = .FALSE.
189 0 : ref_k = 0.0_dp
190 0 : CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmata, ref_basis%nbas)
191 0 : ref_energy = atom_trace(ref_k, ref_orbs%pmata)
192 0 : ref_k = 0.0_dp
193 0 : CALL eeri_contract(ref_k, ref_int%eeri, ref_orbs%pmatb, ref_basis%nbas)
194 0 : ref_energy = ref_energy + atom_trace(ref_k, ref_orbs%pmatb)
195 : CASE (do_rohf_atom)
196 0 : CPABORT("ADMM not available")
197 : CASE DEFAULT
198 1 : CPABORT("ADMM not available")
199 : END SELECT
200 1 : DEALLOCATE (ref_k)
201 : ! reference number of electrons
202 1 : elref = atom_trace(ref_int%ovlp, ref_orbs%pmat)
203 : ! admm orbitals and density matrices
204 7 : mo = MAXVAL(atom%state%maxn_calc)
205 1 : NULLIFY (admm1_orbs, admm2_orbs, admmq_orbs)
206 1 : CALL create_atom_orbs(admm1_orbs, na, mo)
207 1 : CALL create_atom_orbs(admm2_orbs, na, mo)
208 1 : CALL create_atom_orbs(admmq_orbs, na, mo)
209 4 : ALLOCATE (lamat(1:mo, 1:mo))
210 5 : ALLOCATE (admm1_k(1:na, 1:na, 0:lmat))
211 3 : ALLOCATE (admm2_k(1:na, 1:na, 0:lmat))
212 3 : ALLOCATE (admmq_k(1:na, 1:na, 0:lmat))
213 1 : IF (rhf) THEN
214 7 : DO l = 0, lmat
215 6 : n = admm_basis%nbas(l)
216 6 : m = ref_basis%nbas(l)
217 6 : mo = atom%state%maxn_calc(l)
218 6 : IF (n < 1 .OR. m < 1 .OR. mo < 1) CYCLE
219 123 : admm2_orbs%wfn(1:n, 1:mo, l) = MATMUL(tmat(1:n, 1:m, l), ref_orbs%wfn(1:m, 1:mo, l))
220 2 : 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 61 : 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 1 : 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 1 : atom%state%maxl_occ, atom%state%maxn_occ)
227 1 : el1 = atom_trace(admm_int%ovlp, admm1_orbs%pmat)
228 1 : el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmat)
229 1 : xsi = elref/el2
230 157 : admmq_orbs%pmat = xsi*admm2_orbs%pmat
231 1 : elq = atom_trace(admm_int%ovlp, admmq_orbs%pmat)
232 109 : admmq_orbs%wfn = SQRT(xsi)*admm2_orbs%wfn
233 : !
234 79 : admm1_k = 0.0_dp
235 1 : CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmat, admm_basis%nbas)
236 1 : admm1_k_energy = 0.5_dp*atom_trace(admm1_k, admm1_orbs%pmat)
237 79 : admm2_k = 0.0_dp
238 1 : CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmat, admm_basis%nbas)
239 1 : admm2_k_energy = 0.5_dp*atom_trace(admm2_k, admm2_orbs%pmat)
240 79 : admmq_k = 0.0_dp
241 1 : CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmat, admm_basis%nbas)
242 1 : admmq_k_energy = 0.5_dp*atom_trace(admmq_k, admmq_orbs%pmat)
243 : ELSE
244 0 : DO l = 0, lmat
245 0 : n = admm_basis%nbas(l)
246 0 : m = ref_basis%nbas(l)
247 0 : mo = atom%state%maxn_calc(l)
248 0 : IF (n < 1 .OR. m < 1 .OR. mo < 1) CYCLE
249 0 : admm2_orbs%wfna(1:n, 1:mo, l) = MATMUL(tmat(1:n, 1:m, l), ref_orbs%wfna(1:m, 1:mo, l))
250 0 : CALL lowdin_matrix(admm2_orbs%wfna(1:n, 1:mo, l), lamat, admm_int%ovlp(1:n, 1:n, l))
251 0 : admm1_orbs%wfna(1:n, 1:mo, l) = MATMUL(admm2_orbs%wfna(1:n, 1:mo, l), lamat(1:mo, 1:mo))
252 0 : admm2_orbs%wfnb(1:n, 1:mo, l) = MATMUL(tmat(1:n, 1:m, l), ref_orbs%wfnb(1:m, 1:mo, l))
253 0 : CALL lowdin_matrix(admm2_orbs%wfnb(1:n, 1:mo, l), lamat, admm_int%ovlp(1:n, 1:n, l))
254 0 : 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 0 : 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 0 : atom%state%maxl_occ, atom%state%maxn_occ)
260 0 : 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 0 : 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 0 : atom%state%maxl_occ, atom%state%maxn_occ)
265 0 : admm2_orbs%pmat = admm2_orbs%pmata + admm2_orbs%pmatb
266 0 : elref = atom_trace(ref_int%ovlp, ref_orbs%pmata)
267 0 : el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmata)
268 0 : xsi = elref/el2
269 0 : admmq_orbs%pmata = xsi*admm2_orbs%pmata
270 0 : admmq_orbs%wfna = SQRT(xsi)*admm2_orbs%wfna
271 0 : elref = atom_trace(ref_int%ovlp, ref_orbs%pmatb)
272 0 : el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmatb)
273 0 : xsi = elref/el2
274 0 : admmq_orbs%pmatb = xsi*admm2_orbs%pmatb
275 0 : admmq_orbs%wfnb = SQRT(xsi)*admm2_orbs%wfnb
276 0 : admmq_orbs%pmat = admmq_orbs%pmata + admmq_orbs%pmatb
277 0 : el1 = atom_trace(admm_int%ovlp, admm1_orbs%pmat)
278 0 : el2 = atom_trace(admm_int%ovlp, admm2_orbs%pmat)
279 : elq = atom_trace(admm_int%ovlp, admmq_orbs%pmat)
280 0 : elref = atom_trace(ref_int%ovlp, ref_orbs%pmat)
281 : !
282 0 : admm1_k = 0.0_dp
283 0 : CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmata, admm_basis%nbas)
284 0 : admm1_k_energy = atom_trace(admm1_k, admm1_orbs%pmata)
285 0 : admm1_k = 0.0_dp
286 0 : CALL eeri_contract(admm1_k, admm_int%eeri, admm1_orbs%pmatb, admm_basis%nbas)
287 0 : admm1_k_energy = admm1_k_energy + atom_trace(admm1_k, admm1_orbs%pmatb)
288 0 : admm2_k = 0.0_dp
289 0 : CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmata, admm_basis%nbas)
290 0 : admm2_k_energy = atom_trace(admm2_k, admm2_orbs%pmata)
291 0 : admm2_k = 0.0_dp
292 0 : CALL eeri_contract(admm2_k, admm_int%eeri, admm2_orbs%pmatb, admm_basis%nbas)
293 0 : admm2_k_energy = admm2_k_energy + atom_trace(admm2_k, admm2_orbs%pmatb)
294 0 : admmq_k = 0.0_dp
295 0 : CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmata, admm_basis%nbas)
296 0 : admmq_k_energy = atom_trace(admmq_k, admmq_orbs%pmata)
297 0 : admmq_k = 0.0_dp
298 0 : CALL eeri_contract(admmq_k, admm_int%eeri, admmq_orbs%pmatb, admm_basis%nbas)
299 0 : admmq_k_energy = admmq_k_energy + atom_trace(admmq_k, admmq_orbs%pmatb)
300 : END IF
301 1 : DEALLOCATE (lamat)
302 : !
303 : ! ADMM correction terms
304 1 : maxl = atom%state%maxl_occ
305 1 : IF (rhf) THEN
306 3 : ALLOCATE (ref_xcmat(1:nb, 1:nb, 0:lmat))
307 5 : ALLOCATE (admm_xcmat(1:na, 1:na, 0:lmat))
308 : ! PBEX
309 1 : CALL atom_vxc_lda(ref_basis, ref_orbs%pmat, maxl, xc_pbex, fexc_pbex_ref, ref_xcmat)
310 1 : CALL atom_vxc_lda(admm_basis, admm1_orbs%pmat, maxl, xc_pbex, fexc_pbex_admm1, admm_xcmat)
311 1 : CALL atom_vxc_lda(admm_basis, admm2_orbs%pmat, maxl, xc_pbex, fexc_pbex_admm2, admm_xcmat)
312 1 : CALL atom_vxc_lda(admm_basis, admmq_orbs%pmat, maxl, xc_pbex, fexc_pbex_admmq, admm_xcmat)
313 : ! OPTX
314 1 : CALL atom_vxc_lda(ref_basis, ref_orbs%pmat, maxl, xc_optx, fexc_optx_ref, ref_xcmat)
315 1 : CALL atom_vxc_lda(admm_basis, admm1_orbs%pmat, maxl, xc_optx, fexc_optx_admm1, admm_xcmat)
316 1 : CALL atom_vxc_lda(admm_basis, admm2_orbs%pmat, maxl, xc_optx, fexc_optx_admm2, admm_xcmat)
317 1 : 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 1 : maxl, "LINX", dfexc_admm1)
321 : CALL atom_dpot_lda(ref_basis, ref_orbs%pmat, admm_basis, admm2_orbs%pmat, &
322 1 : maxl, "LINX", dfexc_admm2)
323 : CALL atom_dpot_lda(ref_basis, ref_orbs%pmat, admm_basis, admmq_orbs%pmat, &
324 1 : maxl, "LINX", dfexc_admmq)
325 1 : DEALLOCATE (ref_xcmat, admm_xcmat)
326 : ELSE
327 0 : ALLOCATE (ref_xcmata(1:nb, 1:nb, 0:lmat))
328 0 : ALLOCATE (ref_xcmatb(1:nb, 1:nb, 0:lmat))
329 0 : ALLOCATE (admm_xcmata(1:na, 1:na, 0:lmat))
330 0 : 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 0 : ref_xcmata, ref_xcmatb)
334 : CALL atom_vxc_lsd(admm_basis, admm1_orbs%pmata, admm1_orbs%pmatb, maxl, xc_pbex, &
335 0 : fexc_pbex_admm1, admm_xcmata, admm_xcmatb)
336 : CALL atom_vxc_lsd(admm_basis, admm2_orbs%pmata, admm2_orbs%pmatb, maxl, xc_pbex, &
337 0 : fexc_pbex_admm2, admm_xcmata, admm_xcmatb)
338 : CALL atom_vxc_lsd(admm_basis, admmq_orbs%pmata, admmq_orbs%pmatb, maxl, xc_pbex, &
339 0 : 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 0 : ref_xcmata, ref_xcmatb)
342 : CALL atom_vxc_lsd(admm_basis, admm1_orbs%pmata, admm1_orbs%pmatb, maxl, xc_optx, &
343 0 : fexc_optx_admm1, admm_xcmata, admm_xcmatb)
344 : CALL atom_vxc_lsd(admm_basis, admm2_orbs%pmata, admm2_orbs%pmatb, maxl, xc_optx, &
345 0 : fexc_optx_admm2, admm_xcmata, admm_xcmatb)
346 : CALL atom_vxc_lsd(admm_basis, admmq_orbs%pmata, admmq_orbs%pmatb, maxl, xc_optx, &
347 0 : fexc_optx_admmq, admm_xcmata, admm_xcmatb)
348 0 : DEALLOCATE (ref_xcmata, ref_xcmatb, admm_xcmata, admm_xcmatb)
349 : END IF
350 : !
351 :
352 1 : IF (iw > 0) THEN
353 1 : WRITE (iw, "(/,A,I3,T48,A,I3)") " Electronic Structure Setting: ", i, &
354 2 : " Electronic Structure Method: ", j
355 1 : WRITE (iw, "(' Norm of ADMM Basis projection ',T61,F20.10)") el2/elref
356 1 : WRITE (iw, "(' Reference Exchange Energy [Hartree]',T61,F20.10)") ref_energy
357 : ! ADMM1
358 1 : dxk = ref_energy - admm1_k_energy
359 1 : WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMM1 METHOD: Energy ", admm1_k_energy, &
360 2 : " Error: ", dxk
361 1 : dxc = fexc_pbex_ref - fexc_pbex_admm1
362 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
363 2 : " Error: ", dxk - dxc
364 1 : dxc = fexc_optx_ref - fexc_optx_admm1
365 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
366 2 : " Error: ", dxk - dxc
367 1 : dxc = dfexc_admm1
368 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
369 2 : " Error: ", dxk - dxc
370 : ! ADMM2
371 1 : dxk = ref_energy - admm2_k_energy
372 1 : WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMM2 METHOD: Energy ", admm2_k_energy, &
373 2 : " Error: ", dxk
374 1 : dxc = fexc_pbex_ref - fexc_pbex_admm2
375 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
376 2 : " Error: ", dxk - dxc
377 1 : dxc = fexc_optx_ref - fexc_optx_admm2
378 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
379 2 : " Error: ", dxk - dxc
380 1 : dxc = dfexc_admm2
381 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
382 2 : " Error: ", dxk - dxc
383 : ! ADMMQ
384 1 : dxk = ref_energy - admmq_k_energy
385 1 : WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMMQ METHOD: Energy ", admmq_k_energy, &
386 2 : " Error: ", dxk
387 1 : dxc = fexc_pbex_ref - fexc_pbex_admmq
388 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
389 2 : " Error: ", dxk - dxc
390 1 : dxc = fexc_optx_ref - fexc_optx_admmq
391 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
392 2 : " Error: ", dxk - dxc
393 1 : dxc = dfexc_admmq
394 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "LINX Correction ", dxc, dxc/dxk*100._dp, &
395 2 : " Error: ", dxk - dxc
396 : ! ADMMS
397 : dxk = ref_energy - admmq_k_energy
398 1 : WRITE (iw, "(A,F20.10,T60,A,F13.10)") " ADMMS METHOD: Energy ", admmq_k_energy, &
399 2 : " Error: ", dxk
400 1 : dxc = fexc_pbex_ref - fexc_pbex_admmq*xsi**(2._dp/3._dp)
401 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "PBEX Correction ", dxc, dxc/dxk*100._dp, &
402 2 : " Error: ", dxk - dxc
403 1 : dxc = fexc_optx_ref - fexc_optx_admmq*xsi**(2._dp/3._dp)
404 1 : WRITE (iw, "(T10,A,F12.6,F12.3,'%',T60,A,F13.10)") "OPTX Correction ", dxc, dxc/dxk*100._dp, &
405 2 : " Error: ", dxk - dxc
406 : END IF
407 : !
408 1 : DEALLOCATE (admm1_k, admm2_k, admmq_k)
409 : !
410 1 : CALL release_atom_orbs(admm1_orbs)
411 1 : CALL release_atom_orbs(admm2_orbs)
412 1 : CALL release_atom_orbs(admmq_orbs)
413 : END IF
414 : END DO
415 : END DO
416 :
417 : ! clean up
418 1 : CALL atom_int_release(ref_int)
419 1 : CALL atom_int_release(admm_int)
420 1 : CALL release_atom_basis(admm_basis)
421 1 : DEALLOCATE (ref_int, admm_int, admm_basis)
422 1 : DEALLOCATE (ovlap, sinv, tmat, siref)
423 :
424 1 : CALL section_vals_release(xc_pbex)
425 1 : CALL section_vals_release(xc_optx)
426 1 : CALL section_vals_release(xc_section)
427 :
428 1 : IF (iw > 0) THEN
429 : WRITE (iw, '(/,T2,A)') &
430 1 : '!------------------------------End of ADMM analysis---------------------------!'
431 : END IF
432 :
433 1 : END SUBROUTINE atom_admm
434 :
435 : ! **************************************************************************************************
436 : !> \brief ...
437 : !> \param wfn ...
438 : !> \param lamat ...
439 : !> \param ovlp ...
440 : ! **************************************************************************************************
441 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w
448 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vmat
449 :
450 2 : n = SIZE(wfn, 2)
451 2 : IF (n > 0) THEN
452 96 : lamat = MATMUL(TRANSPOSE(wfn), MATMUL(ovlp, wfn))
453 12 : ALLOCATE (w(n), vmat(n, n))
454 2 : CALL jacobi(lamat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
455 5 : w(1:n) = 1.0_dp/SQRT(w(1:n))
456 5 : DO i = 1, n
457 10 : DO j = 1, n
458 5 : lamat(i, j) = 0.0_dp
459 17 : DO k = 1, n
460 14 : lamat(i, j) = lamat(i, j) + vmat(i, k)*w(k)*vmat(j, k)
461 : END DO
462 : END DO
463 : END DO
464 2 : DEALLOCATE (vmat, w)
465 : END IF
466 :
467 2 : END SUBROUTINE lowdin_matrix
468 :
469 2 : END MODULE atom_admm_methods
|
