13 USE dbcsr_api,
ONLY: dbcsr_get_block_p,&
18 #include "./base/base_uses.f90"
24 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pao_param_fock'
43 TYPE(pao_env_type),
POINTER :: pao
44 INTEGER,
INTENT(IN) :: iatom
45 REAL(
dp),
DIMENSION(:, :),
POINTER :: v, u
46 REAL(
dp),
INTENT(INOUT),
OPTIONAL :: penalty
47 REAL(
dp),
INTENT(OUT) :: gap
48 REAL(
dp),
DIMENSION(:),
INTENT(OUT),
OPTIONAL :: evals
49 REAL(
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: m1, g
51 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_calc_U_block_fock'
53 INTEGER :: handle, i, j, m, n
54 INTEGER,
DIMENSION(:),
POINTER :: blk_sizes_pao, blk_sizes_pri
56 REAL(
dp) :: alpha, beta, denom, diff
57 REAL(
dp),
DIMENSION(:),
POINTER :: h_evals
58 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_n, d1, d2, h, h0, h_evecs, m2, m3, &
61 CALL timeset(routinen, handle)
63 CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=h0, found=found)
64 cpassert(
ASSOCIATED(h0))
65 CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_n, found=found)
66 cpassert(
ASSOCIATED(block_n))
67 IF (maxval(abs(v - transpose(v))) > 1e-14_dp) cpabort(
"Expect symmetric matrix")
70 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
71 n = blk_sizes_pri(iatom)
72 m = blk_sizes_pao(iatom)
76 h = matmul(matmul(block_n, h0 + v), block_n)
79 ALLOCATE (h_evals(n), h_evecs(n, n))
87 IF (
PRESENT(evals))
THEN
88 cpassert(mod(
SIZE(evals), 2) == 0)
91 evals(1 + i - j:i) = h_evals(1 + m - j:m)
93 evals(i:i + j) = h_evals(m:m + j)
99 gap = h_evals(m + 1) - h_evals(m)
101 IF (
PRESENT(penalty))
THEN
103 alpha = pao%penalty_strength
104 beta = pao%penalty_dist
107 diff = h_evals(i) - h_evals(j)
108 penalty = penalty + alpha*exp(-(diff/beta)**2)
113 penalty = penalty + pao%regularization*sum(v**2)
118 cpassert(
PRESENT(m1))
121 ALLOCATE (d1(n, n), m2(n, n), m3(n, n), m4(n, n))
127 IF (i <= m .EQV. j <= m)
THEN
130 denom = h_evals(i) - h_evals(j)
131 IF (abs(denom) > 1e-9_dp)
THEN
132 d1(i, j) = 1.0_dp/denom
134 d1(i, j) = sign(1e+9_dp, denom)
139 IF (
ASSOCIATED(m1))
THEN
140 m2 = matmul(transpose(m1), h_evecs)
145 m4 = matmul(matmul(h_evecs, m3), transpose(h_evecs))
148 IF (
PRESENT(penalty))
THEN
153 IF (i <= m .EQV. j <= m) cycle
154 diff = h_evals(i) - h_evals(j)
155 d2(i, i) = d2(i, i) - 2.0_dp*alpha*diff/beta**2*exp(-(diff/beta)**2)
158 m4 = m4 + matmul(matmul(h_evecs, d2), transpose(h_evecs))
164 m5 = matmul(matmul(block_n, m4), block_n)
167 IF (
PRESENT(penalty)) &
168 m5 = m5 + 2.0_dp*pao%regularization*v
171 g = 0.5_dp*(m5 + transpose(m5))
173 DEALLOCATE (d1, m2, m3, m4, m5)
176 DEALLOCATE (h, h_evals, h_evecs)
178 CALL timestop(handle)
Defines the basic variable types.
integer, parameter, public dp
Collection of simple mathematical functions and subroutines.
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Common framework for using eigenvectors of a Fock matrix as PAO basis.
subroutine, public pao_calc_u_block_fock(pao, iatom, V, U, penalty, gap, evals, M1, G)
Calculate new matrix U and optinally its gradient G.
Types used by the PAO machinery.