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 : MODULE optbas_opt_utils
8 : USE atomic_kind_types, ONLY: atomic_kind_type,&
9 : get_atomic_kind
10 : USE basis_set_types, ONLY: get_gto_basis_set,&
11 : gto_basis_set_type
12 : USE cell_types, ONLY: cell_type
13 : USE cp_blacs_env, ONLY: cp_blacs_env_type
14 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
15 : dbcsr_distribution_type,&
16 : dbcsr_get_info,&
17 : dbcsr_p_type,&
18 : dbcsr_release,&
19 : dbcsr_transposed,&
20 : dbcsr_type,&
21 : dbcsr_type_no_symmetry
22 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 : cp_dbcsr_sm_fm_multiply
24 : USE cp_fm_basic_linalg, ONLY: cp_fm_invert,&
25 : cp_fm_trace
26 : USE cp_fm_diag, ONLY: cp_fm_power
27 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_get_info,&
32 : cp_fm_release,&
33 : cp_fm_type
34 : USE distribution_1d_types, ONLY: distribution_1d_type
35 : USE distribution_2d_types, ONLY: distribution_2d_type
36 : USE input_section_types, ONLY: section_vals_val_get
37 : USE kinds, ONLY: dp
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE molecule_types, ONLY: molecule_type
40 : USE parallel_gemm_api, ONLY: parallel_gemm
41 : USE particle_types, ONLY: particle_type
42 : USE qs_condnum, ONLY: overlap_condnum
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kind_types, ONLY: get_qs_kind,&
46 : qs_kind_type
47 : USE qs_ks_types, ONLY: qs_ks_env_type
48 : USE qs_mo_types, ONLY: get_mo_set,&
49 : mo_set_type
50 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
51 : USE qs_neighbor_lists, ONLY: atom2d_build,&
52 : atom2d_cleanup,&
53 : build_neighbor_lists,&
54 : local_atoms_type,&
55 : pair_radius_setup
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 : PRIVATE
60 :
61 : PUBLIC :: evaluate_optvals, fit_mo_coeffs, optbas_build_neighborlist
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optbas_opt_utils'
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param mos ...
70 : !> \param mos_aux_fit ...
71 : !> \param matrix_ks ...
72 : !> \param Q ...
73 : !> \param Snew ...
74 : !> \param S_inv_orb ...
75 : !> \param fval ...
76 : !> \param energy ...
77 : !> \param S_cond_number ...
78 : ! **************************************************************************************************
79 234 : SUBROUTINE evaluate_optvals(mos, mos_aux_fit, matrix_ks, Q, Snew, S_inv_orb, &
80 : fval, energy, S_cond_number)
81 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
83 : TYPE(dbcsr_type), POINTER :: Q, Snew
84 : TYPE(cp_fm_type), INTENT(IN) :: S_inv_orb
85 : REAL(KIND=dp) :: fval, energy, S_cond_number
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'evaluate_optvals'
88 :
89 : INTEGER :: handle, ispin, iunit, naux, nmo, norb, &
90 : nspins
91 234 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
92 : REAL(KIND=dp) :: tmp_energy, trace
93 : REAL(KIND=dp), DIMENSION(2) :: condnum
94 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
95 : TYPE(cp_fm_type) :: tmp1, tmp2
96 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
97 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
98 234 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
99 : TYPE(dbcsr_type) :: Qt
100 :
101 234 : CALL timeset(routineN, handle)
102 :
103 234 : nspins = SIZE(mos)
104 :
105 234 : NULLIFY (col_blk_sizes, row_blk_sizes)
106 : CALL dbcsr_get_info(Q, distribution=dbcsr_dist, &
107 : nfullrows_total=naux, nfullcols_total=norb, &
108 234 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
109 : CALL dbcsr_create(matrix=Qt, name="Qt", &
110 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
111 234 : row_blk_size=col_blk_sizes, col_blk_size=row_blk_sizes)
112 234 : CALL dbcsr_transposed(Qt, Q)
113 : !
114 234 : fval = 0.0_dp
115 234 : energy = 0.0_dp
116 468 : DO ispin = 1, nspins
117 234 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
118 234 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
119 234 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
120 234 : CALL cp_fm_create(tmp1, matrix_struct=mo_coeff%matrix_struct)
121 234 : CALL cp_dbcsr_sm_fm_multiply(Qt, mo_coeff_aux_fit, tmp1, nmo)
122 234 : CALL cp_fm_trace(tmp1, mo_coeff, trace)
123 234 : fval = fval - 2.0_dp*trace + 2.0_dp*nmo
124 : !
125 234 : CALL cp_fm_create(tmp2, matrix_struct=mo_coeff%matrix_struct)
126 234 : CALL parallel_gemm('N', 'N', norb, nmo, norb, 1.0_dp, S_inv_orb, tmp1, 0.0_dp, tmp2)
127 234 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, tmp2, tmp1, nmo)
128 234 : CALL cp_fm_trace(tmp2, tmp1, tmp_energy)
129 234 : energy = energy + tmp_energy*(3.0_dp - REAL(nspins, KIND=dp))
130 234 : CALL cp_fm_release(tmp1)
131 1170 : CALL cp_fm_release(tmp2)
132 : END DO
133 234 : CALL dbcsr_release(Qt)
134 :
135 702 : ALLOCATE (smat(1, 1))
136 234 : smat(1, 1)%matrix => Snew
137 234 : iunit = -1
138 234 : CALL cp_fm_get_info(S_inv_orb, context=blacs_env)
139 234 : CALL overlap_condnum(smat, condnum, iunit, .FALSE., .TRUE., .FALSE., blacs_env)
140 234 : S_cond_number = condnum(2)
141 234 : DEALLOCATE (smat)
142 :
143 234 : CALL timestop(handle)
144 :
145 234 : END SUBROUTINE evaluate_optvals
146 :
147 : ! **************************************************************************************************
148 : !> \brief ...
149 : !> \param saux ...
150 : !> \param sauxorb ...
151 : !> \param mos ...
152 : !> \param mosaux ...
153 : ! **************************************************************************************************
154 234 : SUBROUTINE fit_mo_coeffs(saux, sauxorb, mos, mosaux)
155 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: saux, sauxorb
156 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mosaux
157 :
158 : CHARACTER(len=*), PARAMETER :: routineN = 'fit_mo_coeffs'
159 : REAL(KIND=dp), PARAMETER :: threshold = 1.E-12_dp
160 :
161 : INTEGER :: handle, ispin, naux, ndep, nmo, norb, &
162 : nspins
163 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
164 : TYPE(cp_fm_type) :: fm_s, fm_sinv, tmat, tmp1, tmp2, work
165 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux
166 :
167 234 : CALL timeset(routineN, handle)
168 :
169 234 : CALL dbcsr_get_info(saux(1)%matrix, nfullrows_total=naux)
170 234 : CALL dbcsr_get_info(sauxorb(1)%matrix, nfullcols_total=norb)
171 234 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
172 :
173 : CALL cp_fm_struct_create(fm_struct, nrow_global=naux, ncol_global=naux, &
174 : context=mo_coeff%matrix_struct%context, &
175 234 : para_env=mo_coeff%matrix_struct%para_env)
176 234 : CALL cp_fm_create(fm_s, fm_struct, name="s_aux")
177 234 : CALL cp_fm_create(fm_sinv, fm_struct, name="s_aux_inv")
178 234 : CALL copy_dbcsr_to_fm(saux(1)%matrix, fm_s)
179 234 : CALL cp_fm_invert(fm_s, fm_sinv)
180 234 : CALL cp_fm_release(fm_s)
181 234 : CALL cp_fm_struct_release(fm_struct)
182 234 : nspins = SIZE(mos)
183 468 : DO ispin = 1, nspins
184 234 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
185 234 : CALL get_mo_set(mosaux(ispin), mo_coeff=mo_coeff_aux)
186 234 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
187 234 : CALL cp_fm_create(tmp1, matrix_struct=mo_coeff_aux%matrix_struct)
188 234 : CALL cp_fm_create(tmp2, matrix_struct=mo_coeff_aux%matrix_struct)
189 : CALL cp_fm_struct_create(fm_struct, nrow_global=nmo, ncol_global=nmo, &
190 : context=mo_coeff%matrix_struct%context, &
191 234 : para_env=mo_coeff%matrix_struct%para_env)
192 234 : CALL cp_fm_create(tmat, fm_struct, name="tmat")
193 234 : CALL cp_fm_create(work, fm_struct, name="work")
194 234 : CALL cp_fm_struct_release(fm_struct)
195 : !
196 234 : CALL cp_dbcsr_sm_fm_multiply(sauxorb(1)%matrix, mo_coeff, tmp1, nmo)
197 234 : CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, fm_sinv, tmp1, 0.0_dp, tmp2)
198 234 : CALL parallel_gemm('T', 'N', nmo, nmo, naux, 1.0_dp, tmp1, tmp2, 0.0_dp, tmat)
199 234 : CALL cp_fm_power(tmat, work, -0.5_dp, threshold, ndep)
200 234 : CALL parallel_gemm('N', 'N', naux, nmo, nmo, 1.0_dp, tmp2, tmat, 0.0_dp, mo_coeff_aux)
201 : !
202 234 : CALL cp_fm_release(work)
203 234 : CALL cp_fm_release(tmat)
204 234 : CALL cp_fm_release(tmp1)
205 936 : CALL cp_fm_release(tmp2)
206 : END DO
207 234 : CALL cp_fm_release(fm_sinv)
208 :
209 234 : CALL timestop(handle)
210 :
211 234 : END SUBROUTINE fit_mo_coeffs
212 :
213 : ! **************************************************************************************************
214 : !> \brief rebuilds neighborlist for absis sets
215 : !> \param qs_env ...
216 : !> \param sab_aux ...
217 : !> \param sab_aux_orb ...
218 : !> \param basis_type ...
219 : !> \par History
220 : !> adapted from kg_build_neighborlist
221 : ! **************************************************************************************************
222 234 : SUBROUTINE optbas_build_neighborlist(qs_env, sab_aux, sab_aux_orb, basis_type)
223 : TYPE(qs_environment_type), POINTER :: qs_env
224 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
225 : POINTER :: sab_aux, sab_aux_orb
226 : CHARACTER(*) :: basis_type
227 :
228 : CHARACTER(LEN=*), PARAMETER :: routineN = 'optbas_build_neighborlist'
229 :
230 : INTEGER :: handle, ikind, nkind
231 : LOGICAL :: mic, molecule_only
232 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_fit_present, orb_present
233 : REAL(dp) :: subcells
234 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_fit_radius, orb_radius
235 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
236 234 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
237 : TYPE(cell_type), POINTER :: cell
238 : TYPE(distribution_1d_type), POINTER :: distribution_1d
239 : TYPE(distribution_2d_type), POINTER :: distribution_2d
240 : TYPE(gto_basis_set_type), POINTER :: aux_fit_basis_set, orb_basis_set
241 234 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
242 234 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
243 : TYPE(mp_para_env_type), POINTER :: para_env
244 234 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
245 234 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
246 : TYPE(qs_ks_env_type), POINTER :: ks_env
247 :
248 234 : CALL timeset(routineN, handle)
249 234 : NULLIFY (para_env)
250 :
251 : ! restrict lists to molecular subgroups
252 234 : molecule_only = .FALSE.
253 234 : mic = molecule_only
254 :
255 : CALL get_qs_env(qs_env=qs_env, &
256 : ks_env=ks_env, &
257 : atomic_kind_set=atomic_kind_set, &
258 : qs_kind_set=qs_kind_set, &
259 : cell=cell, &
260 : distribution_2d=distribution_2d, &
261 : molecule_set=molecule_set, &
262 : local_particles=distribution_1d, &
263 : particle_set=particle_set, &
264 234 : para_env=para_env)
265 :
266 234 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
267 :
268 : ! Allocate work storage
269 234 : nkind = SIZE(atomic_kind_set)
270 936 : ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind))
271 639 : orb_radius(:) = 0.0_dp
272 639 : aux_fit_radius(:) = 0.0_dp
273 936 : ALLOCATE (orb_present(nkind), aux_fit_present(nkind))
274 936 : ALLOCATE (pair_radius(nkind, nkind))
275 1107 : ALLOCATE (atom2d(nkind))
276 :
277 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
278 234 : molecule_set, molecule_only, particle_set=particle_set)
279 :
280 639 : DO ikind = 1, nkind
281 405 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
282 405 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
283 405 : IF (ASSOCIATED(orb_basis_set)) THEN
284 405 : orb_present(ikind) = .TRUE.
285 405 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
286 : ELSE
287 0 : orb_present(ikind) = .FALSE.
288 0 : orb_radius(ikind) = 0.0_dp
289 : END IF
290 405 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type=basis_type)
291 639 : IF (ASSOCIATED(aux_fit_basis_set)) THEN
292 405 : aux_fit_present(ikind) = .TRUE.
293 405 : CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
294 : ELSE
295 0 : aux_fit_present(ikind) = .FALSE.
296 0 : aux_fit_radius(ikind) = 0.0_dp
297 : END IF
298 : END DO
299 : !
300 234 : CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, pair_radius)
301 : CALL build_neighbor_lists(sab_aux, particle_set, atom2d, cell, pair_radius, &
302 234 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="sab_aux")
303 234 : CALL pair_radius_setup(aux_fit_present, orb_present, aux_fit_radius, orb_radius, pair_radius)
304 : CALL build_neighbor_lists(sab_aux_orb, particle_set, atom2d, cell, pair_radius, &
305 : mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, &
306 234 : nlname="sab_aux_orb")
307 :
308 : ! Release work storage
309 234 : CALL atom2d_cleanup(atom2d)
310 234 : DEALLOCATE (atom2d)
311 234 : DEALLOCATE (orb_present, aux_fit_present)
312 234 : DEALLOCATE (orb_radius, aux_fit_radius)
313 234 : DEALLOCATE (pair_radius)
314 :
315 234 : CALL timestop(handle)
316 :
317 702 : END SUBROUTINE optbas_build_neighborlist
318 :
319 : END MODULE optbas_opt_utils
|
