Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods for preparing and committing k-point orbital states to the QS environment.
10 : ! **************************************************************************************************
11 : MODULE qs_kpoint_state
12 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,&
13 : cp_cfm_scale
14 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
15 : cp_cfm_heevd
16 : USE cp_cfm_types, ONLY: cp_cfm_create,&
17 : cp_cfm_release,&
18 : cp_cfm_to_cfm,&
19 : cp_cfm_to_fm,&
20 : cp_cfm_type,&
21 : cp_fm_to_cfm
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_info,&
28 : cp_fm_release,&
29 : cp_fm_set_all,&
30 : cp_fm_to_fm,&
31 : cp_fm_type,&
32 : cp_fm_vectorsnorm
33 : USE kinds, ONLY: dp
34 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
35 : kpoint_density_transform,&
36 : kpoint_set_mo_occupation
37 : USE kpoint_types, ONLY: get_kpoint_info,&
38 : kpoint_env_type,&
39 : kpoint_type
40 : USE mathconstants, ONLY: z_one,&
41 : z_zero
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE qs_energy_types, ONLY: qs_energy_type
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_ks_types, ONLY: qs_ks_did_change
47 : USE qs_mo_types, ONLY: get_mo_set,&
48 : mo_set_type
49 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
50 : USE qs_rho_methods, ONLY: qs_rho_update_rho
51 : USE qs_rho_types, ONLY: qs_rho_get,&
52 : qs_rho_type
53 : USE qs_scf_types, ONLY: qs_scf_env_type
54 : USE scf_control_types, ONLY: scf_control_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpoint_state'
61 :
62 : PUBLIC :: qs_kpoint_copy_spin_mos, qs_kpoint_mos_initialized, qs_kpoint_set_fixed_occupations, &
63 : qs_kpoint_state_canonicalize_fixed, qs_kpoint_state_commit, &
64 : qs_kpoint_state_prepare_fixed_density
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief Reconstructs fixed-rank complex k-point orbitals from a density kernel.
70 : !> The density and overlap must first be transferred to kp%ot_hmat and kp%ot_smat. The
71 : !> largest natural-orbital subspace gives a metric-defined fixed-rank projector without
72 : !> replacing the guess by the ground-state projector of its initial Hamiltonian.
73 : !> \param kpoints K-point environment containing the transferred density and overlap matrices.
74 : ! **************************************************************************************************
75 46 : SUBROUTINE qs_kpoint_state_prepare_fixed_density(kpoints)
76 : TYPE(kpoint_type), POINTER :: kpoints
77 :
78 : INTEGER :: ikpoint_local, ispin, nao, nmo, nocc
79 : LOGICAL :: use_real_wfn
80 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: density_eigenvalues
81 : TYPE(cp_cfm_type) :: density_metric, density_work, matrix_p, &
82 : matrix_s, matrix_s_work, projected
83 : TYPE(cp_fm_struct_type), POINTER :: active_struct, square_struct
84 : TYPE(cp_fm_type) :: active_im, active_re
85 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
86 : TYPE(kpoint_env_type), POINTER :: kp
87 :
88 46 : NULLIFY (active_struct, kp, mo_coeff_im, mo_coeff_re, square_struct)
89 :
90 0 : CPASSERT(ASSOCIATED(kpoints))
91 46 : CPASSERT(ASSOCIATED(kpoints%kp_env))
92 46 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
93 46 : IF (use_real_wfn) THEN
94 0 : CPABORT("Fixed-density k-point OT initialization requires complex wavefunctions.")
95 : END IF
96 :
97 153 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
98 107 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
99 107 : CPASSERT(ASSOCIATED(kp))
100 107 : CPASSERT(ASSOCIATED(kp%mos))
101 107 : CPASSERT(ASSOCIATED(kp%ot_hmat))
102 107 : CPASSERT(ASSOCIATED(kp%ot_smat))
103 107 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
104 107 : CPASSERT(SIZE(kp%ot_hmat, 2) >= SIZE(kp%mos, 2))
105 107 : CPASSERT(SIZE(kp%ot_smat) >= 2)
106 :
107 107 : CALL cp_fm_get_info(kp%ot_hmat(1, 1), matrix_struct=square_struct, nrow_global=nao)
108 107 : CALL cp_cfm_create(matrix_p, square_struct, name="natural-orbital density P(k)")
109 107 : CALL cp_cfm_create(matrix_s, square_struct, name="natural-orbital overlap S(k)")
110 107 : CALL cp_fm_to_cfm(kp%ot_smat(1), kp%ot_smat(2), matrix_s)
111 :
112 233 : DO ispin = 1, SIZE(kp%mos, 2)
113 126 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff_re, nmo=nmo)
114 126 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
115 126 : CPASSERT(ASSOCIATED(mo_coeff_re))
116 126 : CPASSERT(ASSOCIATED(mo_coeff_im))
117 126 : CPASSERT(nocc >= 1 .AND. nocc <= nmo)
118 :
119 : CALL cp_fm_struct_create(active_struct, template_fmstruct=mo_coeff_re%matrix_struct, &
120 126 : ncol_global=nocc)
121 126 : CALL cp_fm_to_cfm(kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin), matrix_p)
122 126 : CALL cp_cfm_create(projected, active_struct, name="fixed-density natural orbitals")
123 126 : CALL cp_cfm_create(density_metric, square_struct, name="natural-orbital metric")
124 126 : CALL cp_cfm_create(density_work, square_struct, name="natural-orbital work")
125 126 : CALL cp_cfm_create(matrix_s_work, square_struct, name="natural-orbital overlap work")
126 :
127 : ! The largest eigenvectors of S*P*S*C = S*C*n are the natural orbitals of
128 : ! S**(1/2)*P*S**(1/2). Negating the left-hand side lets the generalized
129 : ! eigensolver request only its lowest nocc eigenpairs when supported.
130 : CALL cp_cfm_gemm('N', 'N', nao, nao, nao, z_one, matrix_p, matrix_s, &
131 126 : z_zero, density_work)
132 : CALL cp_cfm_gemm('N', 'N', nao, nao, nao, z_one, matrix_s, density_work, &
133 126 : z_zero, density_metric)
134 126 : CALL cp_cfm_scale(-z_one, density_metric)
135 126 : CALL cp_cfm_to_cfm(matrix_s, matrix_s_work)
136 378 : ALLOCATE (density_eigenvalues(nocc))
137 : CALL cp_cfm_geeig(density_metric, matrix_s_work, projected, &
138 126 : density_eigenvalues, density_work, lowest_subset=.TRUE.)
139 126 : DEALLOCATE (density_eigenvalues)
140 :
141 126 : CALL cp_fm_create(active_re, active_struct, name="fixed-density orbitals real")
142 126 : CALL cp_fm_create(active_im, active_struct, name="fixed-density orbitals imaginary")
143 126 : CALL cp_cfm_to_fm(projected, active_re, active_im)
144 126 : CALL cp_fm_set_all(mo_coeff_re, 0.0_dp)
145 126 : CALL cp_fm_set_all(mo_coeff_im, 0.0_dp)
146 126 : CALL cp_fm_to_fm(active_re, mo_coeff_re, nocc)
147 126 : CALL cp_fm_to_fm(active_im, mo_coeff_im, nocc)
148 :
149 126 : CALL cp_fm_release(active_im)
150 126 : CALL cp_fm_release(active_re)
151 126 : CALL cp_cfm_release(matrix_s_work)
152 126 : CALL cp_cfm_release(density_work)
153 126 : CALL cp_cfm_release(density_metric)
154 126 : CALL cp_cfm_release(projected)
155 611 : CALL cp_fm_struct_release(active_struct)
156 : END DO
157 :
158 107 : CALL cp_cfm_release(matrix_s)
159 260 : CALL cp_cfm_release(matrix_p)
160 : END DO
161 :
162 92 : END SUBROUTINE qs_kpoint_state_prepare_fixed_density
163 :
164 : ! **************************************************************************************************
165 : !> \brief Rotates a uniformly occupied k-point subspace into Ritz states of the current
166 : !> Hamiltonian and stores their energy labels. The Hamiltonian must first be transferred
167 : !> to kp%ot_hmat. Since the rotation is internal to a uniformly occupied subspace, it does
168 : !> not change the density represented by the orbitals.
169 : !> \param kpoints K-point environment containing the orbitals and transferred Hamiltonians.
170 : ! **************************************************************************************************
171 46 : SUBROUTINE qs_kpoint_state_canonicalize_fixed(kpoints)
172 : TYPE(kpoint_type), POINTER :: kpoints
173 :
174 : INTEGER :: ikpoint_local, ispin, nao, nmo, nocc
175 : LOGICAL :: use_real_wfn
176 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ritz_values
177 46 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_im, eigenvalues_re
178 : TYPE(cp_cfm_type) :: coeff, eigenvectors, hblock, hc, &
179 : matrix_h, rotated
180 : TYPE(cp_fm_struct_type), POINTER :: active_struct, square_struct
181 : TYPE(cp_fm_type) :: active_im, active_re
182 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
183 : TYPE(kpoint_env_type), POINTER :: kp
184 :
185 46 : NULLIFY (active_struct, eigenvalues_im, eigenvalues_re, kp, mo_coeff_im, mo_coeff_re, &
186 46 : square_struct)
187 :
188 0 : CPASSERT(ASSOCIATED(kpoints))
189 46 : CPASSERT(ASSOCIATED(kpoints%kp_env))
190 46 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
191 46 : IF (use_real_wfn) THEN
192 0 : CPABORT("Fixed-density k-point OT canonicalization requires complex wavefunctions.")
193 : END IF
194 :
195 153 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
196 107 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
197 107 : CPASSERT(ASSOCIATED(kp))
198 107 : CPASSERT(ASSOCIATED(kp%mos))
199 107 : CPASSERT(ASSOCIATED(kp%ot_hmat))
200 107 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
201 107 : CPASSERT(SIZE(kp%ot_hmat, 2) >= SIZE(kp%mos, 2))
202 :
203 107 : CALL cp_fm_get_info(kp%ot_hmat(1, 1), matrix_struct=square_struct, nrow_global=nao)
204 107 : CALL cp_cfm_create(matrix_h, square_struct, name="fixed-density Ritz Hamiltonian")
205 :
206 233 : DO ispin = 1, SIZE(kp%mos, 2)
207 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff_re, nmo=nmo, &
208 126 : eigenvalues=eigenvalues_re)
209 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im, &
210 126 : eigenvalues=eigenvalues_im)
211 126 : CPASSERT(ASSOCIATED(mo_coeff_re))
212 126 : CPASSERT(ASSOCIATED(mo_coeff_im))
213 126 : CPASSERT(ASSOCIATED(eigenvalues_re))
214 126 : CPASSERT(ASSOCIATED(eigenvalues_im))
215 126 : CPASSERT(nocc >= 1 .AND. nocc <= nmo)
216 126 : CPASSERT(SIZE(eigenvalues_re) >= nocc)
217 126 : CPASSERT(SIZE(eigenvalues_im) >= nocc)
218 :
219 : CALL cp_fm_struct_create(active_struct, template_fmstruct=mo_coeff_re%matrix_struct, &
220 126 : ncol_global=nocc)
221 126 : CALL cp_fm_create(active_re, active_struct, name="fixed-density Ritz orbitals real")
222 126 : CALL cp_fm_create(active_im, active_struct, name="fixed-density Ritz orbitals imaginary")
223 126 : CALL cp_fm_to_fm(mo_coeff_re, active_re, nocc)
224 126 : CALL cp_fm_to_fm(mo_coeff_im, active_im, nocc)
225 :
226 126 : CALL cp_cfm_create(coeff, active_struct, name="fixed-density Ritz orbitals")
227 126 : CALL cp_cfm_create(hc, active_struct, name="fixed-density Ritz HC")
228 126 : CALL cp_cfm_create(rotated, active_struct, name="fixed-density Ritz rotated orbitals")
229 : CALL cp_cfm_create(hblock, active_struct, nrow=nocc, ncol=nocc, &
230 126 : name="fixed-density projected Hamiltonian")
231 : CALL cp_cfm_create(eigenvectors, active_struct, nrow=nocc, ncol=nocc, &
232 126 : name="fixed-density Ritz eigenvectors")
233 126 : CALL cp_fm_to_cfm(kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin), matrix_h)
234 126 : CALL cp_fm_to_cfm(active_re, active_im, coeff)
235 :
236 126 : CALL cp_cfm_gemm('N', 'N', nao, nocc, nao, z_one, matrix_h, coeff, z_zero, hc)
237 126 : CALL cp_cfm_gemm('C', 'N', nocc, nocc, nao, z_one, coeff, hc, z_zero, hblock)
238 378 : ALLOCATE (ritz_values(nocc))
239 126 : CALL cp_cfm_heevd(hblock, eigenvectors, ritz_values)
240 : CALL cp_cfm_gemm('N', 'N', nao, nocc, nocc, z_one, coeff, eigenvectors, &
241 126 : z_zero, rotated)
242 126 : CALL cp_cfm_to_fm(rotated, active_re, active_im)
243 126 : CALL cp_fm_to_fm(active_re, mo_coeff_re, nocc)
244 126 : CALL cp_fm_to_fm(active_im, mo_coeff_im, nocc)
245 483 : eigenvalues_re(1:nocc) = ritz_values
246 483 : eigenvalues_im(1:nocc) = ritz_values
247 126 : DEALLOCATE (ritz_values)
248 :
249 126 : CALL cp_cfm_release(eigenvectors)
250 126 : CALL cp_cfm_release(hblock)
251 126 : CALL cp_cfm_release(rotated)
252 126 : CALL cp_cfm_release(hc)
253 126 : CALL cp_cfm_release(coeff)
254 126 : CALL cp_fm_release(active_im)
255 126 : CALL cp_fm_release(active_re)
256 611 : CALL cp_fm_struct_release(active_struct)
257 : END DO
258 :
259 260 : CALL cp_cfm_release(matrix_h)
260 : END DO
261 :
262 92 : END SUBROUTINE qs_kpoint_state_canonicalize_fixed
263 :
264 : ! **************************************************************************************************
265 : !> \brief Rebuilds the physical density from the current k-point orbitals and occupations.
266 : !> \param qs_env QS environment containing the orbital state to commit.
267 : !> \param update_occupations Recompute occupations from the current MO energies and SCF settings.
268 : !> \param separate_spin_occupations Fill restricted energy channels independently.
269 : !> \param fixed_occupations Restore the fixed occupied rank required by OT.
270 : ! **************************************************************************************************
271 608 : SUBROUTINE qs_kpoint_state_commit(qs_env, update_occupations, &
272 : separate_spin_occupations, fixed_occupations)
273 : TYPE(qs_environment_type), POINTER :: qs_env
274 : LOGICAL, INTENT(IN) :: update_occupations
275 : LOGICAL, INTENT(IN), OPTIONAL :: separate_spin_occupations, &
276 : fixed_occupations
277 :
278 : INTEGER :: ispin
279 : LOGICAL :: fixed, separate_spins
280 608 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
281 : TYPE(kpoint_type), POINTER :: kpoints
282 608 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
283 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
284 608 : POINTER :: sab_nl
285 : TYPE(qs_energy_type), POINTER :: energy
286 : TYPE(qs_rho_type), POINTER :: rho
287 : TYPE(qs_scf_env_type), POINTER :: scf_env
288 : TYPE(scf_control_type), POINTER :: scf_control
289 :
290 608 : NULLIFY (energy, kpoints, matrix_s_kp, mos, rho_ao_kp, sab_nl, rho, scf_env, scf_control)
291 :
292 608 : fixed = .FALSE.
293 608 : separate_spins = .FALSE.
294 584 : IF (PRESENT(fixed_occupations)) fixed = fixed_occupations
295 608 : IF (PRESENT(separate_spin_occupations)) separate_spins = separate_spin_occupations
296 608 : IF (fixed .AND. .NOT. update_occupations) THEN
297 0 : CPABORT("Fixed occupations can only be requested while updating occupations.")
298 : END IF
299 :
300 : CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp, &
301 608 : energy=energy, rho=rho, scf_env=scf_env, scf_control=scf_control)
302 608 : CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
303 :
304 608 : IF (update_occupations) THEN
305 : CALL kpoint_set_mo_occupation( &
306 584 : kpoints, scf_control%smear, separate_spin_occupations=separate_spins)
307 584 : IF (fixed) CALL qs_kpoint_set_fixed_occupations(kpoints)
308 : END IF
309 608 : CALL kpoint_density_matrices(kpoints)
310 608 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
311 : CALL kpoint_density_transform(kpoints, rho_ao_kp, .FALSE., &
312 : matrix_s_kp(1, 1)%matrix, sab_nl, scf_env%scf_work1, &
313 608 : overlap_rs=matrix_s_kp)
314 608 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
315 608 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
316 :
317 608 : energy%kTS = 0.0_dp
318 608 : energy%efermi = 0.0_dp
319 608 : mos => kpoints%kp_env(1)%kpoint_env%mos
320 608 : CPASSERT(ASSOCIATED(mos))
321 1268 : DO ispin = 1, SIZE(mos, 2)
322 660 : energy%kTS = energy%kTS + mos(1, ispin)%kTS
323 1268 : energy%efermi = energy%efermi + mos(1, ispin)%mu
324 : END DO
325 608 : energy%efermi = energy%efermi/REAL(SIZE(mos, 2), KIND=dp)
326 :
327 608 : END SUBROUTINE qs_kpoint_state_commit
328 :
329 : ! **************************************************************************************************
330 : !> \brief Checks whether every local k-point channel contains nonzero occupied orbitals.
331 : !> \param kpoints K-point environment.
332 : !> \param require_full_space Require every available orbital instead of the occupied block.
333 : !> \return True if every channel contains a nonzero real or imaginary coefficient.
334 : ! **************************************************************************************************
335 154 : FUNCTION qs_kpoint_mos_initialized(kpoints, require_full_space) RESULT(initialized)
336 : TYPE(kpoint_type), POINTER :: kpoints
337 : LOGICAL, INTENT(IN), OPTIONAL :: require_full_space
338 : LOGICAL :: initialized
339 :
340 : INTEGER :: icomponent, ikpoint_local, ispin, &
341 : nactive, nmo, nocc
342 : LOGICAL :: full_space
343 : REAL(KIND=dp) :: coefficient_norm
344 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: column_norm, component_norm
345 : TYPE(cp_fm_type), POINTER :: mo_coeff
346 : TYPE(kpoint_env_type), POINTER :: kp
347 : TYPE(mp_para_env_type), POINTER :: para_env
348 :
349 154 : NULLIFY (kp, mo_coeff, para_env)
350 154 : CPASSERT(ASSOCIATED(kpoints))
351 154 : CPASSERT(ASSOCIATED(kpoints%kp_env))
352 :
353 154 : full_space = .FALSE.
354 154 : IF (PRESENT(require_full_space)) full_space = require_full_space
355 154 : coefficient_norm = HUGE(1.0_dp)
356 403 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
357 249 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
358 249 : CPASSERT(ASSOCIATED(kp))
359 249 : CPASSERT(ASSOCIATED(kp%mos))
360 693 : DO ispin = 1, SIZE(kp%mos, 2)
361 290 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff, nmo=nmo)
362 290 : CPASSERT(ASSOCIATED(mo_coeff))
363 290 : nactive = MERGE(nmo, nocc, full_space)
364 290 : IF (nactive < 1 .OR. nactive > nmo) THEN
365 0 : coefficient_norm = 0.0_dp
366 0 : CYCLE
367 : END IF
368 1160 : ALLOCATE (column_norm(nmo), component_norm(nmo))
369 290 : column_norm(:) = 0.0_dp
370 870 : DO icomponent = 1, SIZE(kp%mos, 1)
371 580 : CALL get_mo_set(kp%mos(icomponent, ispin), mo_coeff=mo_coeff)
372 580 : CPASSERT(ASSOCIATED(mo_coeff))
373 580 : CALL cp_fm_vectorsnorm(mo_coeff, component_norm)
374 4282 : column_norm(:) = column_norm(:) + component_norm(:)**2
375 : END DO
376 1973 : coefficient_norm = MIN(coefficient_norm, SQRT(MINVAL(column_norm(1:nactive))))
377 829 : DEALLOCATE (column_norm, component_norm)
378 : END DO
379 : END DO
380 :
381 154 : CALL get_kpoint_info(kpoints, para_env=para_env)
382 154 : CPASSERT(ASSOCIATED(para_env))
383 154 : CALL para_env%min(coefficient_norm)
384 154 : initialized = coefficient_norm > 100.0_dp*EPSILON(1.0_dp)
385 :
386 154 : END FUNCTION qs_kpoint_mos_initialized
387 :
388 : ! **************************************************************************************************
389 : !> \brief Copies the first k-point spin channel into the remaining channels.
390 : !> \param kpoints K-point environment.
391 : !> \param nspin Number of spin channels to populate.
392 : ! **************************************************************************************************
393 10 : SUBROUTINE qs_kpoint_copy_spin_mos(kpoints, nspin)
394 : TYPE(kpoint_type), POINTER :: kpoints
395 : INTEGER, INTENT(IN) :: nspin
396 :
397 : INTEGER :: icomponent, ikpoint_local, ispin
398 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_source, eigenvalues_target
399 : TYPE(cp_fm_type), POINTER :: mo_coeff_source, mo_coeff_target
400 : TYPE(kpoint_env_type), POINTER :: kp
401 :
402 10 : NULLIFY (eigenvalues_source, eigenvalues_target, kp, mo_coeff_source, mo_coeff_target)
403 0 : CPASSERT(ASSOCIATED(kpoints))
404 10 : CPASSERT(ASSOCIATED(kpoints%kp_env))
405 10 : CPASSERT(nspin >= 1)
406 :
407 21 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
408 11 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
409 11 : CPASSERT(ASSOCIATED(kp))
410 11 : CPASSERT(ASSOCIATED(kp%mos))
411 11 : CPASSERT(SIZE(kp%mos, 2) >= nspin)
412 43 : DO icomponent = 1, SIZE(kp%mos, 1)
413 : CALL get_mo_set(kp%mos(icomponent, 1), mo_coeff=mo_coeff_source, &
414 22 : eigenvalues=eigenvalues_source)
415 22 : CPASSERT(ASSOCIATED(mo_coeff_source))
416 22 : CPASSERT(ASSOCIATED(eigenvalues_source))
417 55 : DO ispin = 2, nspin
418 : CALL get_mo_set(kp%mos(icomponent, ispin), mo_coeff=mo_coeff_target, &
419 22 : eigenvalues=eigenvalues_target)
420 22 : CPASSERT(ASSOCIATED(mo_coeff_target))
421 22 : CPASSERT(ASSOCIATED(eigenvalues_target))
422 22 : CALL cp_fm_to_fm(mo_coeff_source, mo_coeff_target)
423 346 : eigenvalues_target(:) = eigenvalues_source(:)
424 : END DO
425 : END DO
426 : END DO
427 :
428 10 : END SUBROUTINE qs_kpoint_copy_spin_mos
429 :
430 : ! **************************************************************************************************
431 : !> \brief Restores the fixed occupied rank required by k-point OT.
432 : !> \param kpoints K-point environment.
433 : ! **************************************************************************************************
434 120 : SUBROUTINE qs_kpoint_set_fixed_occupations(kpoints)
435 : TYPE(kpoint_type), POINTER :: kpoints
436 :
437 : INTEGER :: icomponent, ikpoint_local, ispin, &
438 : nelectron, nmo, nocc
439 : REAL(KIND=dp) :: maxocc
440 120 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
441 : TYPE(kpoint_env_type), POINTER :: kp
442 :
443 120 : NULLIFY (kp, occupation_numbers)
444 0 : CPASSERT(ASSOCIATED(kpoints))
445 120 : CPASSERT(ASSOCIATED(kpoints%kp_env))
446 :
447 374 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
448 254 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
449 254 : CPASSERT(ASSOCIATED(kp))
450 254 : CPASSERT(ASSOCIATED(kp%mos))
451 678 : DO ispin = 1, SIZE(kp%mos, 2)
452 1166 : DO icomponent = 1, SIZE(kp%mos, 1)
453 : CALL get_mo_set(kp%mos(icomponent, ispin), maxocc=maxocc, &
454 608 : nelectron=nelectron, nmo=nmo, occupation_numbers=occupation_numbers)
455 608 : CPASSERT(ASSOCIATED(occupation_numbers))
456 : ! The fixed rank comes from the spin electron count, not a previous k-point filling.
457 608 : nocc = NINT(REAL(nelectron, KIND=dp)/maxocc)
458 608 : IF (nocc < 1 .OR. nocc > nmo) THEN
459 0 : CPABORT("Fixed-occupation k-point OT requires a nonempty occupied subspace.")
460 : END IF
461 608 : IF (ABS(maxocc*REAL(nocc, KIND=dp) - REAL(nelectron, KIND=dp)) > &
462 : 100.0_dp*EPSILON(1.0_dp)) THEN
463 0 : CPABORT("Fixed-occupation k-point OT requires uniform integer occupations.")
464 : END IF
465 3104 : occupation_numbers(:) = 0.0_dp
466 3004 : occupation_numbers(1:nocc) = maxocc
467 608 : kp%mos(icomponent, ispin)%homo = nocc
468 1520 : kp%mos(icomponent, ispin)%kTS = 0.0_dp
469 : END DO
470 : END DO
471 : END DO
472 :
473 120 : END SUBROUTINE qs_kpoint_set_fixed_occupations
474 :
475 : END MODULE qs_kpoint_state
|
