18 USE dbcsr_api,
ONLY: &
19 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, dbcsr_norm, &
20 dbcsr_norm_maxabsnorm, dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type, &
28 #include "./base/base_uses.f90"
33 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mao_optimizer'
53 SUBROUTINE mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, &
55 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef, matrix_q, matrix_smm
56 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: electra
57 INTEGER,
INTENT(IN) :: max_iter
58 REAL(kind=
dp),
INTENT(IN) :: eps_grad, eps1
59 INTEGER,
INTENT(IN) :: iolevel, iw
61 CHARACTER(len=*),
PARAMETER :: routinen =
'mao_optimize'
63 INTEGER :: handle, i, ispin, iter, nspin
64 INTEGER,
DIMENSION(:),
POINTER :: col_blk_sizes, mao_sizes_a, mao_sizes_b
65 LOGICAL :: spin_warning
66 REAL(kind=
dp) :: a1, a2, alpha, an, beta, eps_fun, fa1, &
67 fa2, fnnew, fnold, fval, grad_norm
68 TYPE(dbcsr_distribution_type) :: dbcsr_dist
69 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_grad
70 TYPE(dbcsr_type) :: amat, binv, cgmat
72 CALL timeset(routinen, handle)
74 eps_fun = 10._dp*eps_grad
76 nspin =
SIZE(mao_coef, 1)
79 WRITE (iw,
'(/,T2,A)')
'!-----------------------------------------------------------------------------!'
80 WRITE (unit=iw, fmt=
"(T32,A)")
"MAO BASIS GENERATION"
81 WRITE (iw,
'(T2,A)')
'!-----------------------------------------------------------------------------!'
87 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, eps1, iolevel, iw)
90 spin_warning = .false.
91 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=mao_sizes_a, distribution=dbcsr_dist)
93 CALL dbcsr_get_info(mao_coef(2)%matrix, col_blk_size=mao_sizes_b, distribution=dbcsr_dist)
94 DO i = 1,
SIZE(mao_sizes_a)
95 IF (mao_sizes_a(i) /= mao_sizes_b(i)) spin_warning = .true.
98 IF (iw > 0 .AND. spin_warning)
THEN
99 WRITE (unit=iw, fmt=
"(/,A)")
">>>> WARNING: Different number of MAOs for different spins"
104 WRITE (unit=iw, fmt=
"(/,A,I2,T70,I11)")
" MAO's for SPIN ", ispin, sum(mao_sizes_a)
105 WRITE (unit=iw, fmt=
"(20I4)") mao_sizes_a(:)
106 ELSE IF (ispin == 2)
THEN
107 WRITE (unit=iw, fmt=
"(/,A,I2,T70,I11)")
" MAO's for SPIN ", ispin, sum(mao_sizes_b)
108 WRITE (unit=iw, fmt=
"(20I4)") mao_sizes_b(:)
111 WRITE (unit=iw, fmt=
"(A)")
114 IF (max_iter < 1)
THEN
117 CALL dbcsr_get_info(mao_coef(ispin)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
118 CALL dbcsr_create(binv, name=
"Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
119 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
120 CALL mao_function(mao_coef(ispin)%matrix, fval, matrix_q(ispin)%matrix, &
121 matrix_smm(1)%matrix, binv, .false.)
123 WRITE (unit=iw, fmt=
"(T2,A,T31,A,I2,T56,A,F12.8)") &
124 "MAO Projection",
"Spin:", ispin,
"Completeness:", fval/electra(ispin)
126 CALL dbcsr_release(binv)
128 IF (iw > 0)
WRITE (unit=iw, fmt=
"(A)")
134 ALLOCATE (mao_grad(ispin)%matrix)
135 CALL dbcsr_create(matrix=mao_grad(ispin)%matrix, name=
"MAO_GRAD", template=mao_coef(1)%matrix)
136 CALL dbcsr_reserve_diag_blocks(matrix=mao_grad(ispin)%matrix)
138 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
139 CALL dbcsr_create(binv, name=
"Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
140 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
141 CALL dbcsr_create(cgmat, template=mao_grad(1)%matrix)
142 CALL dbcsr_create(amat, template=mao_coef(1)%matrix)
147 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .false.)
148 CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
149 CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm)
152 WRITE (unit=iw, fmt=
"(/,T2,A,T73,A,I2)")
"MAO OPTIMIZATION",
"Spin =", ispin
153 WRITE (unit=iw, fmt=
"(T2,A,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
154 "Initialization",
"fval =", fval/electra(ispin),
"grad =", grad_norm,
"step =", alpha
156 DO iter = 1, max_iter
157 IF (grad_norm < eps_grad)
EXIT
158 IF ((1.0_dp - fval/electra(ispin)) < eps_fun)
EXIT
159 CALL dbcsr_add(mao_coef(ispin)%matrix, cgmat, 1.0_dp, alpha)
162 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .true.)
163 CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm)
165 WRITE (unit=iw, fmt=
"(T2,A,i8,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
166 "iter=", iter,
"fval =", fval/electra(ispin),
"grad =", grad_norm,
"step =", alpha
170 CALL dbcsr_add(cgmat, mao_grad(ispin)%matrix, beta, 1.0_dp)
173 CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
174 CALL dbcsr_add(amat, cgmat, 1.0_dp, 0.5_dp*alpha)
176 CALL mao_function(amat, fa1, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .true.)
177 CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
178 CALL dbcsr_add(amat, cgmat, 1.0_dp, alpha)
180 CALL mao_function(amat, fa2, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .true.)
181 a2 = (4._dp*fa1 - fa2 - 3._dp*fval)/alpha
182 a1 = (fa2 - fval - a2*alpha)/(alpha*alpha)
183 IF (abs(a1) > 1.e-14_dp)
THEN
185 an = min(an, 2.0_dp*alpha)
189 IF (an < 0.05_dp .OR. a1 > 0.0_dp)
THEN
190 CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
197 CALL dbcsr_release(binv)
198 CALL dbcsr_release(cgmat)
199 CALL dbcsr_release(amat)
203 CALL timestop(handle)
DBCSR operations in CP2K.
Defines the basic variable types.
integer, parameter, public dp
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
...
subroutine, public mao_orthogonalization(mao_coef, smat)
...
subroutine, public mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
...
real(kind=dp) function, public mao_scalar_product(fmat1, fmat2)
...
subroutine, public mao_function(mao_coef, fval, qmat, smat, binv, reuse)
...
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iolevel, iw)
...