(git:ed6f26b)
Loading...
Searching...
No Matches
qs_charge_mixing.F
Go to the documentation of this file.
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! **************************************************************************************************
9MODULE 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
29CONTAINS
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 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 CALL timeset(routinen, handle)
55
56 IF (mixing_method >= gspace_mixing_nr) THEN
57 cpassert(ASSOCIATED(mixing_store))
58 mixing_store%ncall = mixing_store%ncall + 1
59 ns = SIZE(charges, 2)
60 ns = min(ns, mixing_store%max_shell)
61 alpha = mixing_store%alpha
62 nbuffer = mixing_store%nbuffer
63 inow = mod(mixing_store%ncall - 1, nbuffer) + 1
64 imin = inow - 1
65 IF (imin == 0) imin = nbuffer
66 IF (mixing_store%ncall > nbuffer) THEN
67 nvec = nbuffer
68 ELSE
69 nvec = mixing_store%ncall - 1
70 END IF
71 IF (mixing_store%ncall > 1) THEN
72 ! store in/out charge difference
73 DO ia = 1, mixing_store%nat_local
74 ii = mixing_store%atlist(ia)
75 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 IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
79 ! skip mixing
80 mixing_store%iter_method = "NoMix"
81 ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
82 CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
83 mixing_store%iter_method = "Mixing"
84 ELSEIF (mixing_method == gspace_mixing_nr) THEN
85 cpabort("Kerker method not available for Charge Mixing")
86 ELSEIF (mixing_method == pulay_mixing_nr) THEN
87 cpabort("Pulay method not available for Charge Mixing")
88 ELSEIF (mixing_method == broyden_mixing_nr) THEN
89 CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
90 mixing_store%iter_method = "Broy."
91 ELSEIF (mixing_method == multisecant_mixing_nr) THEN
92 cpabort("Multisecant_mixing method not available for Charge Mixing")
93 END IF
94
95 ! store new 'input' charges
96 DO ia = 1, mixing_store%nat_local
97 ii = mixing_store%atlist(ia)
98 mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
99 END DO
100
101 END IF
102
103 CALL timestop(handle)
104
105 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 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 charges = 0.0_dp
127
128 DO ia = 1, mixing_store%nat_local
129 ii = mixing_store%atlist(ia)
130 charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
131 END DO
132
133 CALL para_env%sum(charges)
134
135 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
157 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
158 REAL(kind=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
159
160 cpassert(nvec > 1)
161
162 nbuffer = mixing_store%nbuffer
163 alpha = mixing_store%alpha
164 imin = inow - 1
165 IF (imin == 0) imin = nvec
166 nv = nvec - 1
167
168 ! charge vectors
169 q_now => mixing_store%acharge(:, :, inow)
170 q_last => mixing_store%acharge(:, :, imin)
171 dq_now => mixing_store%dacharge(:, :, inow)
172 dq_last => mixing_store%dacharge(:, :, imin)
173
174 IF (nvec == nbuffer) THEN
175 ! reshuffel Broyden storage n->n-1
176 DO i = 1, nv - 1
177 mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
178 mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
179 mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
180 END DO
181 DO i = 1, nv - 1
182 DO j = 1, nv - 1
183 mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
184 END DO
185 END DO
186 END IF
187
188 omega0 = 0.01_dp
189 minw = 1.0_dp
190 maxw = 100000.0_dp
191 wfac = 0.01_dp
192
193 mixing_store%wbroy(nv) = sum(dq_now(:, :)**2)
194 CALL para_env%sum(mixing_store%wbroy(nv))
195 mixing_store%wbroy(nv) = sqrt(mixing_store%wbroy(nv))
196 IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
197 mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
198 ELSE
199 mixing_store%wbroy(nv) = maxw
200 END IF
201 IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
202
203 ! dfbroy
204 mixing_store%dfbroy(:, :, nv) = dq_now(:, :) - dq_last(:, :)
205 wdf = sum(mixing_store%dfbroy(:, :, nv)**2)
206 CALL para_env%sum(wdf)
207 wdf = 1.0_dp/sqrt(wdf)
208 mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)
209
210 ! abroy matrix
211 DO i = 1, nv
212 wfac = sum(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
213 CALL para_env%sum(wfac)
214 mixing_store%abroy(i, nv) = wfac
215 mixing_store%abroy(nv, i) = wfac
216 END DO
217
218 ! broyden matrices
219 ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
220 DO i = 1, nv
221 wfac = sum(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
222 CALL para_env%sum(wfac)
223 cvec(i) = mixing_store%wbroy(i)*wfac
224 END DO
225
226 DO i = 1, nv
227 DO j = 1, nv
228 beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
229 END DO
230 beta(i, i) = beta(i, i) + omega0*omega0
231 END DO
232
233 rskip = 1.e-12_dp
234 CALL get_pseudo_inverse_svd(beta, amat, rskip)
235 gammab(1:nv) = matmul(cvec(1:nv), amat(1:nv, 1:nv))
236
237 ! build ubroy
238 mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv) + wdf*(q_now(:, :) - q_last(:, :))
239
240 charges = 0.0_dp
241 DO ia = 1, mixing_store%nat_local
242 ii = mixing_store%atlist(ia)
243 charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
244 END DO
245 DO i = 1, nv
246 DO ia = 1, mixing_store%nat_local
247 ii = mixing_store%atlist(ia)
248 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 CALL para_env%sum(charges)
252
253 DEALLOCATE (amat, beta, cvec, gammab)
254
255 END SUBROUTINE broyden_mixing
256
257END MODULE qs_charge_mixing
imin
static int imin(int x, int y)
Returns the smaller of the two integers (missing from the C standard).
Definition dbm_miniapp.c:38
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::get_pseudo_inverse_svd
subroutine, public get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
returns the pseudoinverse of a real, square matrix using singular value decomposition
Definition mathlib.F:938
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
qs_charge_mixing
Definition qs_charge_mixing.F:9
qs_charge_mixing::charge_mixing
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count)
Driver for the charge mixing, calls the proper routine given the requested method.
Definition qs_charge_mixing.F:42
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::broyden_mixing_nr
integer, parameter, public broyden_mixing_nr
Definition qs_density_mixing_types.F:45
qs_density_mixing_types::multisecant_mixing_nr
integer, parameter, public multisecant_mixing_nr
Definition qs_density_mixing_types.F:45
qs_density_mixing_types::pulay_mixing_nr
integer, parameter, public pulay_mixing_nr
Definition qs_density_mixing_types.F:45
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:45
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_density_mixing_types::mixing_storage_type
Definition qs_density_mixing_types.F:55