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