Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 20932 : 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 20932 : CALL timeset(routineN, handle)
55 :
56 20932 : IF (mixing_method >= gspace_mixing_nr) THEN
57 556 : CPASSERT(ASSOCIATED(mixing_store))
58 556 : mixing_store%ncall = mixing_store%ncall + 1
59 556 : ns = SIZE(charges, 2)
60 556 : ns = MIN(ns, mixing_store%max_shell)
61 556 : alpha = mixing_store%alpha
62 556 : nbuffer = mixing_store%nbuffer
63 556 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
64 556 : imin = inow - 1
65 556 : IF (imin == 0) imin = nbuffer
66 556 : IF (mixing_store%ncall > nbuffer) THEN
67 334 : nvec = nbuffer
68 : ELSE
69 222 : nvec = mixing_store%ncall - 1
70 : END IF
71 556 : IF (mixing_store%ncall > 1) THEN
72 : ! store in/out charge difference
73 2322 : DO ia = 1, mixing_store%nat_local
74 1802 : ii = mixing_store%atlist(ia)
75 11332 : 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 556 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
79 : ! skip mixing
80 36 : mixing_store%iter_method = "NoMix"
81 520 : ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
82 36 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
83 36 : mixing_store%iter_method = "Mixing"
84 484 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
85 0 : CPABORT("Kerker method not available for Charge Mixing")
86 484 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
87 0 : CPABORT("Pulay method not available for Charge Mixing")
88 484 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
89 484 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
90 484 : 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 2452 : DO ia = 1, mixing_store%nat_local
97 1896 : ii = mixing_store%atlist(ia)
98 11932 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
99 : END DO
100 :
101 : END IF
102 :
103 20932 : CALL timestop(handle)
104 :
105 20932 : 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 36 : 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 1156 : charges = 0.0_dp
127 :
128 130 : DO ia = 1, mixing_store%nat_local
129 94 : ii = mixing_store%atlist(ia)
130 600 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
131 : END DO
132 :
133 2276 : CALL para_env%sum(charges)
134 :
135 36 : 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 484 : 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 484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
157 484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
158 484 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
159 :
160 0 : CPASSERT(nvec > 1)
161 :
162 484 : nbuffer = mixing_store%nbuffer
163 484 : alpha = mixing_store%alpha
164 484 : imin = inow - 1
165 484 : IF (imin == 0) imin = nvec
166 484 : nv = nvec - 1
167 :
168 : ! charge vectors
169 484 : q_now => mixing_store%acharge(:, :, inow)
170 484 : q_last => mixing_store%acharge(:, :, imin)
171 484 : dq_now => mixing_store%dacharge(:, :, inow)
172 484 : dq_last => mixing_store%dacharge(:, :, imin)
173 :
174 484 : IF (nvec == nbuffer) THEN
175 : ! reshuffel Broyden storage n->n-1
176 1194 : DO i = 1, nv - 1
177 860 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
178 21080 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
179 21414 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
180 : END DO
181 1194 : DO i = 1, nv - 1
182 4450 : DO j = 1, nv - 1
183 4116 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
184 : END DO
185 : END DO
186 : END IF
187 :
188 484 : omega0 = 0.01_dp
189 484 : minw = 1.0_dp
190 484 : maxw = 100000.0_dp
191 484 : wfac = 0.01_dp
192 :
193 11444 : mixing_store%wbroy(nv) = SUM(dq_now(:, :)**2)
194 484 : CALL para_env%sum(mixing_store%wbroy(nv))
195 484 : mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
196 484 : IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
197 472 : mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
198 : ELSE
199 12 : mixing_store%wbroy(nv) = maxw
200 : END IF
201 484 : IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
202 :
203 : ! dfbroy
204 22404 : mixing_store%dfbroy(:, :, nv) = dq_now(:, :) - dq_last(:, :)
205 11444 : wdf = SUM(mixing_store%dfbroy(:, :, nv)**2)
206 484 : CALL para_env%sum(wdf)
207 484 : wdf = 1.0_dp/SQRT(wdf)
208 11444 : mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)
209 :
210 : ! abroy matrix
211 2350 : DO i = 1, nv
212 45846 : wfac = SUM(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
213 1866 : CALL para_env%sum(wfac)
214 1866 : mixing_store%abroy(i, nv) = wfac
215 2350 : mixing_store%abroy(nv, i) = wfac
216 : END DO
217 :
218 : ! broyden matrices
219 4356 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
220 2350 : DO i = 1, nv
221 45846 : wfac = SUM(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
222 1866 : CALL para_env%sum(wfac)
223 2350 : cvec(i) = mixing_store%wbroy(i)*wfac
224 : END DO
225 :
226 2350 : DO i = 1, nv
227 12508 : DO j = 1, nv
228 12508 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
229 : END DO
230 2350 : beta(i, i) = beta(i, i) + omega0*omega0
231 : END DO
232 :
233 484 : rskip = 1.e-12_dp
234 484 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
235 14858 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
236 :
237 : ! build ubroy
238 22404 : mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv) + wdf*(q_now(:, :) - q_last(:, :))
239 :
240 19984 : charges = 0.0_dp
241 2192 : DO ia = 1, mixing_store%nat_local
242 1708 : ii = mixing_store%atlist(ia)
243 10732 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
244 : END DO
245 2350 : DO i = 1, nv
246 9280 : DO ia = 1, mixing_store%nat_local
247 6930 : ii = mixing_store%atlist(ia)
248 43446 : 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 39484 : CALL para_env%sum(charges)
252 :
253 484 : DEALLOCATE (amat, beta, cvec, gammab)
254 :
255 484 : END SUBROUTINE broyden_mixing
256 :
257 : END MODULE qs_charge_mixing
|