Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_maxabs,&
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_orthogonalization,&
28 : mao_scalar_product
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 :
73 14 : CALL timeset(routineN, handle)
74 :
75 14 : eps_fun = 10._dp*eps_grad
76 :
77 14 : nspin = SIZE(mao_coef, 1)
78 :
79 14 : IF (iw > 0) THEN
80 7 : WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
81 7 : WRITE (UNIT=iw, FMT="(T32,A)") "MAO BASIS GENERATION"
82 7 : WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
83 : END IF
84 :
85 : ! initialize MAO coeficients from diagonal blocks of the Q matrix
86 30 : DO ispin = 1, nspin
87 : CALL mao_initialization(mao_coef(ispin)%matrix, &
88 30 : matrix_q(ispin)%matrix, matrix_smm(1)%matrix, eps1, iolevel, iw)
89 : END DO
90 : ! check for mao sizes
91 14 : spin_warning = .FALSE.
92 14 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=mao_sizes_a, distribution=dbcsr_dist)
93 14 : IF (nspin == 2) THEN
94 2 : CALL dbcsr_get_info(mao_coef(2)%matrix, col_blk_size=mao_sizes_b, distribution=dbcsr_dist)
95 8 : DO i = 1, SIZE(mao_sizes_a)
96 8 : IF (mao_sizes_a(i) /= mao_sizes_b(i)) spin_warning = .TRUE.
97 : END DO
98 : END IF
99 14 : IF (iw > 0 .AND. spin_warning) THEN
100 0 : WRITE (UNIT=iw, FMT="(/,A)") ">>>> WARNING: Different number of MAOs for different spins"
101 : END IF
102 14 : IF (iw > 0) THEN
103 15 : DO ispin = 1, nspin
104 15 : IF (ispin == 1) THEN
105 43 : WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_a)
106 7 : WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_a(:)
107 1 : ELSE IF (ispin == 2) THEN
108 4 : WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_b)
109 1 : WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_b(:)
110 : END IF
111 : END DO
112 7 : WRITE (UNIT=iw, FMT="(A)")
113 : END IF
114 :
115 14 : IF (max_iter < 1) THEN
116 : ! projection only
117 12 : DO ispin = 1, nspin
118 6 : CALL dbcsr_get_info(mao_coef(ispin)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
119 : CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
120 6 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
121 : CALL mao_function(mao_coef(ispin)%matrix, fval, matrix_q(ispin)%matrix, &
122 6 : matrix_smm(1)%matrix, binv, .FALSE.)
123 6 : IF (iw > 0) THEN
124 : WRITE (UNIT=iw, FMT="(T2,A,T31,A,I2,T56,A,F12.8)") &
125 3 : "MAO Projection", "Spin:", ispin, "Completeness:", fval/electra(ispin)
126 : END IF
127 18 : CALL dbcsr_release(binv)
128 : END DO
129 6 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)")
130 : ELSE
131 : ! optimize MAOs
132 8 : NULLIFY (mao_grad)
133 8 : CALL dbcsr_allocate_matrix_set(mao_grad, nspin)
134 18 : DO ispin = 1, nspin
135 10 : ALLOCATE (mao_grad(ispin)%matrix)
136 10 : CALL dbcsr_create(matrix=mao_grad(ispin)%matrix, name="MAO_GRAD", template=mao_coef(1)%matrix)
137 18 : CALL dbcsr_reserve_diag_blocks(matrix=mao_grad(ispin)%matrix)
138 : END DO
139 8 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
140 : CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
141 8 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
142 8 : CALL dbcsr_create(cgmat, template=mao_grad(1)%matrix)
143 8 : CALL dbcsr_create(amat, template=mao_coef(1)%matrix)
144 18 : DO ispin = 1, nspin
145 10 : alpha = 0.25_dp
146 10 : beta = 0.0_dp
147 : CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
148 10 : matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .FALSE.)
149 10 : CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
150 10 : grad_norm = dbcsr_maxabs(mao_grad(ispin)%matrix)
151 10 : fnold = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
152 10 : IF (iw > 0) THEN
153 5 : WRITE (UNIT=iw, FMT="(/,T2,A,T73,A,I2)") "MAO OPTIMIZATION", "Spin =", ispin
154 : WRITE (UNIT=iw, FMT="(T2,A,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
155 5 : "Initialization", "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
156 : END IF
157 174 : DO iter = 1, max_iter
158 160 : IF (grad_norm < eps_grad) EXIT
159 160 : IF ((1.0_dp - fval/electra(ispin)) < eps_fun) EXIT
160 156 : CALL dbcsr_add(mao_coef(ispin)%matrix, cgmat, 1.0_dp, alpha)
161 156 : CALL mao_orthogonalization(mao_coef(ispin)%matrix, matrix_smm(1)%matrix)
162 : CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
163 156 : matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
164 156 : grad_norm = dbcsr_maxabs(mao_grad(ispin)%matrix)
165 156 : IF (iw > 0) THEN
166 : WRITE (UNIT=iw, FMT="(T2,A,i8,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
167 78 : "iter=", iter, "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
168 : END IF
169 156 : fnnew = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
170 156 : beta = fnnew/fnold
171 156 : CALL dbcsr_add(cgmat, mao_grad(ispin)%matrix, beta, 1.0_dp)
172 156 : fnold = fnnew
173 : ! line search, update alpha
174 156 : CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
175 156 : CALL dbcsr_add(amat, cgmat, 1.0_dp, 0.5_dp*alpha)
176 156 : CALL mao_orthogonalization(amat, matrix_smm(1)%matrix)
177 156 : CALL mao_function(amat, fa1, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
178 156 : CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
179 156 : CALL dbcsr_add(amat, cgmat, 1.0_dp, alpha)
180 156 : CALL mao_orthogonalization(amat, matrix_smm(1)%matrix)
181 156 : CALL mao_function(amat, fa2, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
182 156 : a2 = (4._dp*fa1 - fa2 - 3._dp*fval)/alpha
183 156 : a1 = (fa2 - fval - a2*alpha)/(alpha*alpha)
184 156 : IF (ABS(a1) > 1.e-14_dp) THEN
185 156 : an = -a2/(2._dp*a1)
186 156 : an = MIN(an, 2.0_dp*alpha)
187 : ELSE
188 0 : an = 2.0_dp*alpha
189 : END IF
190 166 : IF (an < 0.05_dp .OR. a1 > 0.0_dp) THEN
191 0 : CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
192 0 : alpha = 0.25_dp
193 : ELSE
194 156 : alpha = an
195 : END IF
196 : END DO
197 : END DO
198 8 : CALL dbcsr_release(binv)
199 8 : CALL dbcsr_release(cgmat)
200 8 : CALL dbcsr_release(amat)
201 8 : IF (ASSOCIATED(mao_grad)) CALL dbcsr_deallocate_matrix_set(mao_grad)
202 : END IF
203 :
204 14 : CALL timestop(handle)
205 :
206 14 : END SUBROUTINE mao_optimize
207 :
208 : END MODULE mao_optimizer
|
