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 : !> \brief Optimizer for the atomic code
10 : ! **************************************************************************************************
11 : MODULE atom_optimization
12 : USE atom_types, ONLY: atom_optimization_type
13 : USE kinds, ONLY: dp
14 : #include "./base/base_uses.f90"
15 :
16 : IMPLICIT NONE
17 :
18 : PRIVATE
19 :
20 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_optimization'
21 :
22 : TYPE hmat_type
23 : REAL(dp) :: energy = 0.0_dp
24 : REAL(dp) :: error = 0.0_dp
25 : REAL(dp), DIMENSION(:, :, :), POINTER :: emat => NULL()
26 : REAL(dp), DIMENSION(:, :, :), POINTER :: fmat => NULL()
27 : REAL(dp), DIMENSION(:, :, :), POINTER :: pmat => NULL()
28 : END TYPE hmat_type
29 :
30 : TYPE atom_history_type
31 : INTEGER :: max_history = 0
32 : INTEGER :: hlen = 0
33 : INTEGER :: hpos = 0
34 : REAL(dp) :: damping = 0.0_dp
35 : REAL(dp) :: eps_diis = 0.0_dp
36 : REAL(dp), DIMENSION(:, :), POINTER :: dmat => NULL()
37 : TYPE(hmat_type), DIMENSION(:), POINTER :: hmat => NULL()
38 : END TYPE atom_history_type
39 :
40 : PUBLIC :: atom_opt_fmat, &
41 : atom_history_type, atom_history_init, atom_history_update, atom_history_release
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief Initialise a circular buffer to keep Kohn-Sham and density matrices from previous iteration.
47 : !> \param history object to initialise
48 : !> \param optimization optimisation parameters
49 : !> \param matrix reference matrix. Historic matrices will have the same size as
50 : !> this reference matrix
51 : !> \par History
52 : !> * 08.2016 new structure element: density matrix [Juerg Hutter]
53 : !> * 08.2008 created [Juerg Hutter]
54 : ! **************************************************************************************************
55 11441 : PURE SUBROUTINE atom_history_init(history, optimization, matrix)
56 : TYPE(atom_history_type), INTENT(INOUT) :: history
57 : TYPE(atom_optimization_type), INTENT(IN) :: optimization
58 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: matrix
59 :
60 : INTEGER :: i, n1, n2, n3, ndiis
61 : REAL(dp) :: damp, eps
62 :
63 11441 : ndiis = optimization%n_diis
64 11441 : eps = optimization%eps_diis
65 11441 : damp = optimization%damping
66 :
67 11441 : CALL atom_history_release(history)
68 :
69 11441 : history%max_history = ndiis
70 11441 : history%hlen = 0
71 11441 : history%hpos = 0
72 11441 : history%damping = damp
73 11441 : history%eps_diis = eps
74 45764 : ALLOCATE (history%dmat(ndiis + 1, ndiis + 1))
75 :
76 91530 : ALLOCATE (history%hmat(ndiis))
77 11441 : n1 = SIZE(matrix, 1)
78 11441 : n2 = SIZE(matrix, 2)
79 11441 : n3 = SIZE(matrix, 3)
80 68648 : DO i = 1, ndiis
81 57207 : history%hmat(i)%energy = 0.0_dp
82 57207 : history%hmat(i)%error = 0.0_dp
83 285995 : ALLOCATE (history%hmat(i)%emat(n1, n2, n3))
84 228788 : ALLOCATE (history%hmat(i)%fmat(n1, n2, n3))
85 240229 : ALLOCATE (history%hmat(i)%pmat(n1, n2, n3))
86 : END DO
87 :
88 11441 : END SUBROUTINE atom_history_init
89 :
90 : ! **************************************************************************************************
91 : !> \brief Add matrices from the current iteration into the circular buffer.
92 : !> \param history object to keep historic matrices
93 : !> \param pmat density matrix
94 : !> \param fmat Kohn-Sham matrix
95 : !> \param emat error matrix
96 : !> \param energy total energy
97 : !> \param error convergence
98 : !> \par History
99 : !> * 08.2016 new formal argument: density matrix [Juerg Hutter]
100 : !> * 08.2008 created [Juerg Hutter]
101 : ! **************************************************************************************************
102 66856 : PURE SUBROUTINE atom_history_update(history, pmat, fmat, emat, energy, error)
103 : TYPE(atom_history_type), INTENT(INOUT) :: history
104 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: pmat, fmat, emat
105 : REAL(dp), INTENT(IN) :: energy, error
106 :
107 : INTEGER :: nlen, nmax, nnow
108 :
109 66856 : nmax = history%max_history
110 66856 : nlen = MIN(history%hlen + 1, nmax)
111 66856 : nnow = history%hpos + 1
112 66856 : IF (nnow > nmax) nnow = 1
113 :
114 66856 : history%hmat(nnow)%energy = energy
115 66856 : history%hmat(nnow)%error = error
116 34414108 : history%hmat(nnow)%pmat = pmat
117 34414108 : history%hmat(nnow)%fmat = fmat
118 34414108 : history%hmat(nnow)%emat = emat
119 :
120 66856 : history%hlen = nlen
121 66856 : history%hpos = nnow
122 :
123 66856 : END SUBROUTINE atom_history_update
124 :
125 : ! **************************************************************************************************
126 : !> \brief Release circular buffer to keep historic matrices.
127 : !> \param history object to release
128 : !> \par History
129 : !> * 08.2008 created [Juerg Hutter]
130 : ! **************************************************************************************************
131 22882 : PURE SUBROUTINE atom_history_release(history)
132 : TYPE(atom_history_type), INTENT(INOUT) :: history
133 :
134 : INTEGER :: i
135 :
136 22882 : history%max_history = 0
137 22882 : history%hlen = 0
138 22882 : history%hpos = 0
139 22882 : history%damping = 0._dp
140 22882 : history%eps_diis = 0._dp
141 22882 : IF (ASSOCIATED(history%dmat)) THEN
142 11441 : DEALLOCATE (history%dmat)
143 : END IF
144 22882 : IF (ASSOCIATED(history%hmat)) THEN
145 68648 : DO i = 1, SIZE(history%hmat)
146 57207 : IF (ASSOCIATED(history%hmat(i)%emat)) THEN
147 57207 : DEALLOCATE (history%hmat(i)%emat)
148 : END IF
149 57207 : IF (ASSOCIATED(history%hmat(i)%fmat)) THEN
150 57207 : DEALLOCATE (history%hmat(i)%fmat)
151 : END IF
152 68648 : IF (ASSOCIATED(history%hmat(i)%pmat)) THEN
153 57207 : DEALLOCATE (history%hmat(i)%pmat)
154 : END IF
155 : END DO
156 11441 : DEALLOCATE (history%hmat)
157 : END IF
158 :
159 22882 : END SUBROUTINE atom_history_release
160 :
161 : ! **************************************************************************************************
162 : !> \brief Construct a Kohn-Sham matrix for the next iteration based on the historic data.
163 : !> \param fmat new Kohn-Sham matrix
164 : !> \param history historic matrices
165 : !> \param err convergence
166 : !> \par History
167 : !> * 08.2016 renamed to atom_opt_fmat() [Juerg Hutter]
168 : !> * 08.2008 created as atom_opt() [Juerg Hutter]
169 : ! **************************************************************************************************
170 66856 : SUBROUTINE atom_opt_fmat(fmat, history, err)
171 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: fmat
172 : TYPE(atom_history_type), INTENT(INOUT) :: history
173 : REAL(dp), INTENT(IN) :: err
174 :
175 : INTEGER :: i, info, j, lwork, na, nb, nlen, nm, &
176 : nmax, nnow, rank
177 : REAL(dp) :: a, rcond, t
178 66856 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: s, work
179 66856 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vec
180 :
181 66856 : nmax = history%max_history
182 66856 : nnow = history%hpos
183 66856 : a = history%damping
184 66856 : IF (history%hlen > 1) THEN
185 57595 : IF (err < history%eps_diis) THEN
186 : ! DIIS
187 57595 : rcond = 1.e-10_dp
188 57595 : lwork = 25*nmax
189 460760 : ALLOCATE (vec(nmax + 1, 2), s(nmax + 1), work(lwork))
190 : nlen = history%hlen
191 863949 : vec = 0._dp
192 57595 : vec(nlen + 1, 1) = 1._dp
193 298876 : history%dmat(1:nlen, nlen + 1) = 1._dp
194 298876 : history%dmat(nlen + 1, 1:nlen) = 1._dp
195 57595 : history%dmat(nlen + 1, nlen + 1) = 0._dp
196 298876 : DO i = 1, nlen
197 241281 : na = nnow + 1 - i
198 241281 : IF (na < 1) na = nmax + na
199 963750 : DO j = i, nlen
200 664874 : nb = nnow + 1 - j
201 664874 : IF (nb < 1) nb = nmax + nb
202 374805326 : t = SUM(history%hmat(na)%emat*history%hmat(nb)%emat)
203 664874 : history%dmat(i, j) = t
204 906155 : history%dmat(j, i) = t
205 : END DO
206 : END DO
207 : CALL dgelss(nlen + 1, nlen + 1, 1, history%dmat, nmax + 1, vec, nmax + 1, s, &
208 57595 : rcond, rank, work, lwork, info)
209 57595 : CPASSERT(info == 0)
210 30911245 : fmat = 0._dp
211 298876 : DO i = 1, nlen
212 241281 : na = nnow + 1 - i
213 241281 : IF (na < 1) na = nmax + na
214 134061022 : fmat = fmat + vec(i, 1)*history%hmat(na)%fmat
215 : END DO
216 :
217 57595 : DEALLOCATE (vec, s, work)
218 : ELSE
219 : ! damping
220 0 : nm = nnow - 1
221 0 : IF (nm < 1) nm = history%max_history
222 0 : fmat = a*history%hmat(nnow)%fmat + (1._dp - a)*history%hmat(nm)%fmat
223 : END IF
224 9261 : ELSEIF (history%hlen == 1) THEN
225 3502863 : fmat = history%hmat(nnow)%fmat
226 : ELSE
227 0 : CPABORT("")
228 : END IF
229 :
230 66856 : END SUBROUTINE atom_opt_fmat
231 :
232 0 : END MODULE atom_optimization
|
