(git:1f285aa)
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-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  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 
257 END MODULE qs_charge_mixing
imin
static int imin(int x, int y)
Returns the smaller of two given integer (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:949
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