1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE atom_admm_methods 8 USE atom_operators, ONLY: atom_basis_projection_overlap,& 9 atom_int_release,& 10 atom_int_setup 11 USE atom_output, ONLY: atom_print_basis 12 USE atom_types, ONLY: & 13 atom_basis_type, atom_integrals, atom_orbitals, atom_p_type, atom_type, create_atom_orbs, & 14 init_atom_basis, lmat, release_atom_basis, release_atom_orbs 15 USE atom_utils, ONLY: atom_consistent_method,& 16 atom_denmat,& 17 atom_trace,& 18 eeri_contract 19 USE atom_xc, ONLY: atom_dpot_lda,& 20 atom_vxc_lda,& 21 atom_vxc_lsd 22 USE input_constants, ONLY: do_rhf_atom,& 23 do_rks_atom,& 24 do_rohf_atom,& 25 do_uhf_atom,& 26 do_uks_atom,& 27 xc_funct_no_shortcut 28 USE input_section_types, ONLY: section_vals_duplicate,& 29 section_vals_get_subs_vals,& 30 section_vals_get_subs_vals2,& 31 section_vals_release,& 32 section_vals_remove_values,& 33 section_vals_type,& 34 section_vals_val_set 35 USE kinds, ONLY: dp 36 USE mathlib, ONLY: invmat_symm,& 37 jacobi 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 42 PRIVATE 43 PUBLIC :: atom_admm 44 45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_admm_methods' 46 47CONTAINS 48 49! ************************************************************************************************** 50!> \brief Analysis of ADMM approximation to exact exchange. 51!> \param atom_info information about the atomic kind. Two-dimensional array of size 52!> (electronic-configuration, electronic-structure-method) 53!> \param admm_section ADMM print section 54!> \param iw output file unit 55!> \par History 56!> * 07.2016 created [Juerg Hutter] 57! ************************************************************************************************** 58 SUBROUTINE atom_admm(atom_info, admm_section, iw) 59 60 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info 61 TYPE(section_vals_type), POINTER :: admm_section 62 INTEGER, INTENT(IN) :: iw 63 64 CHARACTER(len=*), PARAMETER :: routineN = 'atom_admm', routineP = moduleN//':'//routineN 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 470