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 : MODULE qs_charge_mixing
10 :
11 : USE kinds, ONLY: dp
12 : USE mathlib, ONLY: get_pseudo_inverse_svd
13 : USE message_passing, ONLY: mp_para_env_type
14 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
15 : gspace_mixing_nr,&
16 : mixing_storage_type,&
17 : multisecant_mixing_nr,&
18 : pulay_mixing_nr
19 : #include "./base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 :
23 : PRIVATE
24 :
25 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_charge_mixing'
26 :
27 : PUBLIC :: charge_mixing
28 :
29 : CONTAINS
30 :
31 : ! **************************************************************************************************
32 : !> \brief Driver for the charge mixing, calls the proper routine given the requested method
33 : !> \param mixing_method ...
34 : !> \param mixing_store ...
35 : !> \param charges ...
36 : !> \param para_env ...
37 : !> \param iter_count ...
38 : !> \par History
39 : !> \author JGH
40 : ! **************************************************************************************************
41 43498 : SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count)
42 : INTEGER, INTENT(IN) :: mixing_method
43 : TYPE(mixing_storage_type), POINTER :: mixing_store
44 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
45 : TYPE(mp_para_env_type), POINTER :: para_env
46 : INTEGER, INTENT(IN) :: iter_count
47 :
48 : CHARACTER(len=*), PARAMETER :: routineN = 'charge_mixing'
49 :
50 : INTEGER :: handle, ia, ii, imin, inow, nbuffer, ns, &
51 : nvec
52 : REAL(dp) :: alpha
53 :
54 43498 : CALL timeset(routineN, handle)
55 :
56 43498 : IF (mixing_method >= gspace_mixing_nr) THEN
57 492 : CPASSERT(ASSOCIATED(mixing_store))
58 492 : mixing_store%ncall = mixing_store%ncall + 1
59 492 : ns = SIZE(charges, 2)
60 492 : ns = MIN(ns, mixing_store%max_shell)
61 492 : alpha = mixing_store%alpha
62 492 : nbuffer = mixing_store%nbuffer
63 492 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
64 492 : imin = inow - 1
65 492 : IF (imin == 0) imin = nbuffer
66 492 : IF (mixing_store%ncall > nbuffer) THEN
67 108 : nvec = nbuffer
68 : ELSE
69 384 : nvec = mixing_store%ncall - 1
70 : END IF
71 492 : IF (mixing_store%ncall > 1) THEN
72 : ! store in/out charge difference
73 1796 : DO ia = 1, mixing_store%nat_local
74 1358 : ii = mixing_store%atlist(ia)
75 8586 : mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin) - charges(ii, 1:ns)
76 : END DO
77 : END IF
78 492 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
79 : ! skip mixing
80 54 : mixing_store%iter_method = "NoMix"
81 438 : ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
82 50 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
83 50 : mixing_store%iter_method = "Mixing"
84 388 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
85 0 : CPABORT("Kerker method not available for Charge Mixing")
86 388 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
87 0 : CPABORT("Pulay method not available for Charge Mixing")
88 388 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
89 388 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
90 388 : mixing_store%iter_method = "Broy."
91 0 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
92 0 : CPABORT("Multisecant_mixing method not available for Charge Mixing")
93 : END IF
94 :
95 : ! store new 'input' charges
96 1980 : DO ia = 1, mixing_store%nat_local
97 1488 : ii = mixing_store%atlist(ia)
98 9420 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
99 : END DO
100 :
101 : END IF
102 :
103 43498 : CALL timestop(handle)
104 :
105 43498 : END SUBROUTINE charge_mixing
106 :
107 : ! **************************************************************************************************
108 : !> \brief Simple charge mixing
109 : !> \param mixing_store ...
110 : !> \param charges ...
111 : !> \param alpha ...
112 : !> \param imin ...
113 : !> \param ns ...
114 : !> \param para_env ...
115 : !> \author JGH
116 : ! **************************************************************************************************
117 50 : SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
118 : TYPE(mixing_storage_type), POINTER :: mixing_store
119 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
120 : REAL(KIND=dp), INTENT(IN) :: alpha
121 : INTEGER, INTENT(IN) :: imin, ns
122 : TYPE(mp_para_env_type), POINTER :: para_env
123 :
124 : INTEGER :: ia, ii
125 :
126 1540 : charges = 0.0_dp
127 :
128 174 : DO ia = 1, mixing_store%nat_local
129 124 : ii = mixing_store%atlist(ia)
130 794 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
131 : END DO
132 :
133 3030 : CALL para_env%sum(charges)
134 :
135 50 : END SUBROUTINE mix_charges_only
136 :
137 : ! **************************************************************************************************
138 : !> \brief Broyden charge mixing
139 : !> \param mixing_store ...
140 : !> \param charges ...
141 : !> \param inow ...
142 : !> \param nvec ...
143 : !> \param ns ...
144 : !> \param para_env ...
145 : !> \author JGH
146 : ! **************************************************************************************************
147 388 : SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
148 : TYPE(mixing_storage_type), POINTER :: mixing_store
149 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
150 : INTEGER, INTENT(IN) :: inow, nvec, ns
151 : TYPE(mp_para_env_type), POINTER :: para_env
152 :
153 : INTEGER :: i, ia, ii, imin, j, nbuffer, nv
154 : REAL(KIND=dp) :: alpha, maxw, minw, omega0, rskip, wdf, &
155 : wfac
156 388 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
157 388 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
158 388 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
159 :
160 0 : CPASSERT(nvec > 1)
161 :
162 388 : nbuffer = mixing_store%nbuffer
163 388 : alpha = mixing_store%alpha
164 388 : imin = inow - 1
165 388 : IF (imin == 0) imin = nvec
166 388 : nv = nvec - 1
167 :
168 : ! charge vectors
169 388 : q_now => mixing_store%acharge(:, :, inow)
170 388 : q_last => mixing_store%acharge(:, :, imin)
171 388 : dq_now => mixing_store%dacharge(:, :, inow)
172 388 : dq_last => mixing_store%dacharge(:, :, imin)
173 :
174 388 : IF (nvec == nbuffer) THEN
175 : ! reshuffel Broyden storage n->n-1
176 956 : DO i = 1, nv - 1
177 848 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
178 21468 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
179 21576 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
180 : END DO
181 956 : DO i = 1, nv - 1
182 8556 : DO j = 1, nv - 1
183 8448 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
184 : END DO
185 : END DO
186 : END IF
187 :
188 388 : omega0 = 0.01_dp
189 388 : minw = 1.0_dp
190 388 : maxw = 100000.0_dp
191 388 : wfac = 0.01_dp
192 :
193 8498 : mixing_store%wbroy(nv) = SUM(dq_now(:, :)**2)
194 388 : CALL para_env%sum(mixing_store%wbroy(nv))
195 388 : mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
196 388 : IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
197 378 : mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
198 : ELSE
199 10 : mixing_store%wbroy(nv) = maxw
200 : END IF
201 388 : IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
202 :
203 : ! dfbroy
204 16608 : mixing_store%dfbroy(:, :, nv) = dq_now(:, :) - dq_last(:, :)
205 8498 : wdf = SUM(mixing_store%dfbroy(:, :, nv)**2)
206 388 : CALL para_env%sum(wdf)
207 388 : wdf = 1.0_dp/SQRT(wdf)
208 8498 : mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)
209 :
210 : ! abroy matrix
211 2624 : DO i = 1, nv
212 53636 : wfac = SUM(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
213 2236 : CALL para_env%sum(wfac)
214 2236 : mixing_store%abroy(i, nv) = wfac
215 2624 : mixing_store%abroy(nv, i) = wfac
216 : END DO
217 :
218 : ! broyden matrices
219 3492 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
220 2624 : DO i = 1, nv
221 53636 : wfac = SUM(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
222 2236 : CALL para_env%sum(wfac)
223 2624 : cvec(i) = mixing_store%wbroy(i)*wfac
224 : END DO
225 :
226 2624 : DO i = 1, nv
227 20648 : DO j = 1, nv
228 20648 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
229 : END DO
230 2624 : beta(i, i) = beta(i, i) + omega0*omega0
231 : END DO
232 :
233 388 : rskip = 1.e-12_dp
234 388 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
235 23272 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
236 :
237 : ! build ubroy
238 16608 : mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv) + wdf*(q_now(:, :) - q_last(:, :))
239 :
240 14668 : charges = 0.0_dp
241 1622 : DO ia = 1, mixing_store%nat_local
242 1234 : ii = mixing_store%atlist(ia)
243 7792 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
244 : END DO
245 2624 : DO i = 1, nv
246 10668 : DO ia = 1, mixing_store%nat_local
247 8044 : ii = mixing_store%atlist(ia)
248 50500 : charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
249 : END DO
250 : END DO
251 28948 : CALL para_env%sum(charges)
252 :
253 388 : DEALLOCATE (amat, beta, cvec, gammab)
254 :
255 388 : END SUBROUTINE broyden_mixing
256 :
257 : END MODULE qs_charge_mixing
|
