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 : 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 41260 : 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 41260 : CALL timeset(routineN, handle)
55 :
56 41260 : IF (mixing_method >= gspace_mixing_nr) THEN
57 638 : CPASSERT(ASSOCIATED(mixing_store))
58 638 : mixing_store%ncall = mixing_store%ncall + 1
59 638 : ns = SIZE(charges, 2)
60 638 : ns = MIN(ns, mixing_store%max_shell)
61 638 : alpha = mixing_store%alpha
62 638 : nbuffer = mixing_store%nbuffer
63 638 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
64 638 : imin = inow - 1
65 638 : IF (imin == 0) imin = nbuffer
66 638 : IF (mixing_store%ncall > nbuffer) THEN
67 350 : nvec = nbuffer
68 : ELSE
69 288 : nvec = mixing_store%ncall - 1
70 : END IF
71 638 : IF (mixing_store%ncall > 1) THEN
72 : ! store in/out charge difference
73 2514 : DO ia = 1, mixing_store%nat_local
74 1930 : ii = mixing_store%atlist(ia)
75 12164 : 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 638 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
79 : ! skip mixing
80 54 : mixing_store%iter_method = "NoMix"
81 584 : ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
82 52 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
83 52 : mixing_store%iter_method = "Mixing"
84 532 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
85 0 : CPABORT("Kerker method not available for Charge Mixing")
86 532 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
87 0 : CPABORT("Pulay method not available for Charge Mixing")
88 532 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
89 532 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
90 532 : 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 2698 : DO ia = 1, mixing_store%nat_local
97 2060 : ii = mixing_store%atlist(ia)
98 12998 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
99 : END DO
100 :
101 : END IF
102 :
103 41260 : CALL timestop(handle)
104 :
105 41260 : 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 52 : 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 1572 : charges = 0.0_dp
127 :
128 178 : DO ia = 1, mixing_store%nat_local
129 126 : ii = mixing_store%atlist(ia)
130 808 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
131 : END DO
132 :
133 3092 : CALL para_env%sum(charges)
134 :
135 52 : 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 532 : 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 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
157 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
158 532 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
159 :
160 0 : CPASSERT(nvec > 1)
161 :
162 532 : nbuffer = mixing_store%nbuffer
163 532 : alpha = mixing_store%alpha
164 532 : imin = inow - 1
165 532 : IF (imin == 0) imin = nvec
166 532 : nv = nvec - 1
167 :
168 : ! charge vectors
169 532 : q_now => mixing_store%acharge(:, :, inow)
170 532 : q_last => mixing_store%acharge(:, :, imin)
171 532 : dq_now => mixing_store%dacharge(:, :, inow)
172 532 : dq_last => mixing_store%dacharge(:, :, imin)
173 :
174 532 : IF (nvec == nbuffer) THEN
175 : ! reshuffel Broyden storage n->n-1
176 1242 : DO i = 1, nv - 1
177 892 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
178 21592 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
179 21942 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
180 : END DO
181 1242 : DO i = 1, nv - 1
182 4562 : DO j = 1, nv - 1
183 4212 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
184 : END DO
185 : END DO
186 : END IF
187 :
188 532 : omega0 = 0.01_dp
189 532 : minw = 1.0_dp
190 532 : maxw = 100000.0_dp
191 532 : wfac = 0.01_dp
192 :
193 12212 : mixing_store%wbroy(nv) = SUM(dq_now(:, :)**2)
194 532 : CALL para_env%sum(mixing_store%wbroy(nv))
195 532 : mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
196 532 : IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
197 520 : mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
198 : ELSE
199 12 : mixing_store%wbroy(nv) = maxw
200 : END IF
201 532 : IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
202 :
203 : ! dfbroy
204 23892 : mixing_store%dfbroy(:, :, nv) = dq_now(:, :) - dq_last(:, :)
205 12212 : wdf = SUM(mixing_store%dfbroy(:, :, nv)**2)
206 532 : CALL para_env%sum(wdf)
207 532 : wdf = 1.0_dp/SQRT(wdf)
208 12212 : mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)
209 :
210 : ! abroy matrix
211 2494 : DO i = 1, nv
212 47382 : wfac = SUM(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
213 1962 : CALL para_env%sum(wfac)
214 1962 : mixing_store%abroy(i, nv) = wfac
215 2494 : mixing_store%abroy(nv, i) = wfac
216 : END DO
217 :
218 : ! broyden matrices
219 4788 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
220 2494 : DO i = 1, nv
221 47382 : wfac = SUM(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
222 1962 : CALL para_env%sum(wfac)
223 2494 : cvec(i) = mixing_store%wbroy(i)*wfac
224 : END DO
225 :
226 2494 : DO i = 1, nv
227 12828 : DO j = 1, nv
228 12828 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
229 : END DO
230 2494 : beta(i, i) = beta(i, i) + omega0*omega0
231 : END DO
232 :
233 532 : rskip = 1.e-12_dp
234 532 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
235 15322 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
236 :
237 : ! build ubroy
238 23892 : mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv) + wdf*(q_now(:, :) - q_last(:, :))
239 :
240 21232 : charges = 0.0_dp
241 2336 : DO ia = 1, mixing_store%nat_local
242 1804 : ii = mixing_store%atlist(ia)
243 11356 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
244 : END DO
245 2494 : DO i = 1, nv
246 9616 : DO ia = 1, mixing_store%nat_local
247 7122 : ii = mixing_store%atlist(ia)
248 44694 : 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 41932 : CALL para_env%sum(charges)
252 :
253 532 : DEALLOCATE (amat, beta, cvec, gammab)
254 :
255 532 : END SUBROUTINE broyden_mixing
256 :
257 : END MODULE qs_charge_mixing
|
