Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate MAO's and analyze wavefunctions
10 : !> \par History
11 : !> 03.2016 created [JGH]
12 : !> 12.2016 split into four modules [JGH]
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE mao_optimizer
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, &
18 : dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_symmetric
19 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm,&
20 : dbcsr_reserve_diag_blocks
21 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
22 : dbcsr_deallocate_matrix_set
23 : USE kinds, ONLY: dp
24 : USE mao_methods, ONLY: mao_function,&
25 : mao_function_gradient,&
26 : mao_initialization,&
27 : mao_scalar_product,&
28 : mao_update_coef
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_optimizer'
35 :
36 : PUBLIC :: mao_optimize
37 :
38 : ! **************************************************************************************************
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief ...
44 : !> \param mao_coef ...
45 : !> \param matrix_q ...
46 : !> \param matrix_smm ...
47 : !> \param electra ...
48 : !> \param max_iter ...
49 : !> \param eps_grad ...
50 : !> \param eps1 ...
51 : !> \param iolevel ...
52 : !> \param iw ...
53 : ! **************************************************************************************************
54 14 : SUBROUTINE mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, &
55 : iolevel, iw)
56 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_q, matrix_smm
57 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: electra
58 : INTEGER, INTENT(IN) :: max_iter
59 : REAL(KIND=dp), INTENT(IN) :: eps_grad, eps1
60 : INTEGER, INTENT(IN) :: iolevel, iw
61 :
62 : CHARACTER(len=*), PARAMETER :: routineN = 'mao_optimize'
63 :
64 : INTEGER :: handle, i, ispin, iter, nspin
65 14 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_sizes_a, mao_sizes_b
66 : LOGICAL :: spin_warning
67 : REAL(KIND=dp) :: a1, a2, alpha, an, beta, eps_fun, fa1, &
68 : fa2, fnnew, fnold, fval, grad_norm
69 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
70 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_grad
71 : TYPE(dbcsr_type) :: amat, binv, cgmat
72 : TYPE(dbcsr_type), POINTER :: smat
73 :
74 14 : CALL timeset(routineN, handle)
75 :
76 14 : eps_fun = 10._dp*eps_grad
77 :
78 14 : nspin = SIZE(mao_coef, 1)
79 :
80 14 : IF (iw > 0) THEN
81 7 : WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
82 7 : WRITE (UNIT=iw, FMT="(T32,A)") "MAO BASIS GENERATION"
83 7 : WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
84 : END IF
85 :
86 : ! initialize MAO coeficients from diagonal blocks of the Q matrix
87 14 : smat => matrix_smm(1)%matrix
88 30 : DO ispin = 1, nspin
89 30 : CALL mao_initialization(mao_coef(ispin)%matrix, matrix_q(ispin)%matrix, smat, eps1, iolevel, iw)
90 : END DO
91 : ! check for mao sizes
92 14 : spin_warning = .FALSE.
93 14 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=mao_sizes_a, distribution=dbcsr_dist)
94 14 : IF (nspin == 2) THEN
95 2 : CALL dbcsr_get_info(mao_coef(2)%matrix, col_blk_size=mao_sizes_b, distribution=dbcsr_dist)
96 8 : DO i = 1, SIZE(mao_sizes_a)
97 8 : IF (mao_sizes_a(i) /= mao_sizes_b(i)) spin_warning = .TRUE.
98 : END DO
99 : END IF
100 14 : IF (iw > 0 .AND. spin_warning) THEN
101 0 : WRITE (UNIT=iw, FMT="(/,A)") ">>>> WARNING: Different number of MAOs for different spins"
102 : END IF
103 14 : IF (iw > 0) THEN
104 15 : DO ispin = 1, nspin
105 15 : IF (ispin == 1) THEN
106 43 : WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_a)
107 43 : WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_a(:)
108 1 : ELSE IF (ispin == 2) THEN
109 4 : WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_b)
110 4 : WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_b(:)
111 : END IF
112 : END DO
113 7 : WRITE (UNIT=iw, FMT="(A)")
114 : END IF
115 :
116 14 : IF (max_iter < 1) THEN
117 : ! projection only
118 12 : DO ispin = 1, nspin
119 6 : CALL dbcsr_get_info(mao_coef(ispin)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
120 : CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
121 6 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
122 : CALL mao_function(mao_coef(ispin)%matrix, fval, matrix_q(ispin)%matrix, &
123 6 : smat, binv, .FALSE.)
124 6 : IF (iw > 0) THEN
125 : WRITE (UNIT=iw, FMT="(T2,A,T31,A,I2,T56,A,F12.8)") &
126 3 : "MAO Projection", "Spin:", ispin, "Completeness:", fval/electra(ispin)
127 : END IF
128 18 : CALL dbcsr_release(binv)
129 : END DO
130 6 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)")
131 : ELSE
132 : ! optimize MAOs
133 8 : NULLIFY (mao_grad)
134 8 : CALL dbcsr_allocate_matrix_set(mao_grad, nspin)
135 18 : DO ispin = 1, nspin
136 10 : ALLOCATE (mao_grad(ispin)%matrix)
137 10 : CALL dbcsr_create(matrix=mao_grad(ispin)%matrix, name="MAO_GRAD", template=mao_coef(1)%matrix)
138 18 : CALL dbcsr_reserve_diag_blocks(matrix=mao_grad(ispin)%matrix)
139 : END DO
140 8 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
141 : CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
142 8 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
143 8 : CALL dbcsr_create(cgmat, template=mao_grad(1)%matrix)
144 8 : CALL dbcsr_create(amat, template=mao_coef(1)%matrix)
145 18 : DO ispin = 1, nspin
146 10 : alpha = 0.25_dp
147 : beta = 0.0_dp
148 : CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
149 10 : matrix_q(ispin)%matrix, smat, binv, .FALSE.)
150 10 : grad_norm = dbcsr_frobenius_norm(mao_grad(ispin)%matrix)
151 10 : CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
152 10 : fnold = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
153 10 : IF (iw > 0) THEN
154 5 : WRITE (UNIT=iw, FMT="(/,T2,A,T73,A,I2)") "MAO OPTIMIZATION", "Spin =", ispin
155 : WRITE (UNIT=iw, FMT="(T2,A,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
156 5 : "Initialization", "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
157 : END IF
158 216 : DO iter = 1, max_iter
159 202 : IF (grad_norm < eps_grad) EXIT
160 202 : IF ((1.0_dp - fval/electra(ispin)) < eps_fun) EXIT
161 198 : CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
162 198 : CALL mao_update_coef(mao_coef(ispin)%matrix, amat, cgmat, smat, alpha)
163 : CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
164 198 : matrix_q(ispin)%matrix, smat, binv, .TRUE.)
165 198 : grad_norm = dbcsr_frobenius_norm(mao_grad(ispin)%matrix)
166 198 : IF (iw > 0) THEN
167 : WRITE (UNIT=iw, FMT="(T2,A,i8,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
168 99 : "iter=", iter, "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
169 : END IF
170 198 : fnnew = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
171 198 : beta = fnnew/fnold
172 198 : CALL dbcsr_add(cgmat, mao_grad(ispin)%matrix, beta, 1.0_dp)
173 198 : fnold = fnnew
174 : ! line search, update alpha
175 198 : CALL mao_update_coef(amat, mao_coef(ispin)%matrix, cgmat, smat, 0.5_dp*alpha)
176 198 : CALL mao_function(amat, fa1, matrix_q(ispin)%matrix, smat, binv, .TRUE.)
177 198 : CALL mao_update_coef(amat, mao_coef(ispin)%matrix, cgmat, smat, alpha)
178 198 : CALL mao_function(amat, fa2, matrix_q(ispin)%matrix, smat, binv, .TRUE.)
179 198 : a2 = (4._dp*fa1 - fa2 - 3._dp*fval)/alpha
180 198 : a1 = (fa2 - fval - a2*alpha)/(alpha*alpha)
181 198 : IF (ABS(a1) > 1.e-14_dp) THEN
182 198 : an = -a2/(2._dp*a1)
183 198 : an = MIN(an, 2.0_dp*alpha)
184 : ELSE
185 0 : an = 2.0_dp*alpha
186 : END IF
187 208 : IF (an < 0.05_dp .OR. a1 > 0.0_dp) THEN
188 0 : CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
189 0 : alpha = 0.25_dp
190 : ELSE
191 198 : alpha = an
192 : END IF
193 : END DO
194 : END DO
195 8 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)")
196 8 : CALL dbcsr_release(binv)
197 8 : CALL dbcsr_release(cgmat)
198 8 : CALL dbcsr_release(amat)
199 8 : IF (ASSOCIATED(mao_grad)) CALL dbcsr_deallocate_matrix_set(mao_grad)
200 : END IF
201 :
202 14 : CALL timestop(handle)
203 :
204 14 : END SUBROUTINE mao_optimize
205 :
206 : END MODULE mao_optimizer
|
