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