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 deltaSCF calculations
10 : ! **************************************************************************************************
11 : MODULE qs_mom_methods
12 : USE bibliography, ONLY: Barca2018,&
13 : Gilbert2008,&
14 : cite_reference
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
17 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
19 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
20 : cp_fm_struct_release,&
21 : cp_fm_struct_type
22 : USE cp_fm_types, ONLY: cp_fm_create,&
23 : cp_fm_get_info,&
24 : cp_fm_maxabsval,&
25 : cp_fm_to_fm,&
26 : cp_fm_type,&
27 : cp_fm_vectorsnorm,&
28 : cp_fm_vectorssum
29 : USE input_constants, ONLY: momproj_norm,&
30 : momproj_sum,&
31 : momtype_imom,&
32 : momtype_mom
33 : USE input_section_types, ONLY: section_vals_type
34 : USE kinds, ONLY: dp
35 : USE parallel_gemm_api, ONLY: parallel_gemm
36 : USE qs_density_matrices, ONLY: calculate_density_matrix
37 : USE qs_mo_types, ONLY: get_mo_set,&
38 : mo_set_type,&
39 : set_mo_set
40 : USE qs_scf_diagonalization, ONLY: general_eigenproblem
41 : USE qs_scf_types, ONLY: qs_scf_env_type
42 : USE scf_control_types, ONLY: scf_control_type
43 : USE string_utilities, ONLY: integer_to_string
44 : USE util, ONLY: sort,&
45 : sort_unique
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mom_methods'
53 :
54 : PUBLIC :: do_mom_guess, do_mom_diag
55 : PRIVATE :: mom_is_unique_orbital_indices, mom_reoccupy_orbitals
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief check that every molecular orbital index appears only once in each
61 : !> (de-)occupation list supplied by user. Check that all the indices
62 : !> are positive integers and abort if it is not the case.
63 : !> \param iarr list of molecular orbital indices to be checked
64 : !> \return .true. if all the elements are unique or the list contains
65 : !> exactly one 0 element (meaning no excitation)
66 : !> \par History
67 : !> 01.2016 created [Sergey Chulkov]
68 : ! **************************************************************************************************
69 80 : FUNCTION mom_is_unique_orbital_indices(iarr) RESULT(is_unique)
70 : INTEGER, DIMENSION(:), POINTER :: iarr
71 : LOGICAL :: is_unique
72 :
73 : CHARACTER(len=*), PARAMETER :: routineN = 'mom_is_unique_orbital_indices'
74 :
75 : INTEGER :: handle, norbs
76 80 : INTEGER, DIMENSION(:), POINTER :: tmp_iarr
77 :
78 80 : CALL timeset(routineN, handle)
79 :
80 80 : CPASSERT(ASSOCIATED(iarr))
81 80 : norbs = SIZE(iarr)
82 :
83 80 : IF (norbs > 0) THEN
84 240 : ALLOCATE (tmp_iarr(norbs))
85 :
86 320 : tmp_iarr(:) = iarr(:)
87 80 : CALL sort_unique(tmp_iarr, is_unique)
88 :
89 : ! Ensure that all orbital indices are positive integers.
90 : ! A special value '0' means 'disabled keyword',
91 : ! it must appear once to be interpreted in such a way
92 80 : IF (tmp_iarr(1) < 0 .OR. (tmp_iarr(1) == 0 .AND. norbs > 1)) THEN
93 0 : CPABORT("MOM: all molecular orbital indices must be positive integer numbers")
94 : END IF
95 :
96 160 : DEALLOCATE (tmp_iarr)
97 : END IF
98 :
99 : is_unique = .TRUE.
100 :
101 80 : CALL timestop(handle)
102 :
103 80 : END FUNCTION mom_is_unique_orbital_indices
104 :
105 : ! **************************************************************************************************
106 : !> \brief swap occupation numbers between molecular orbitals
107 : !> from occupation and de-occupation lists
108 : !> \param mo_set set of molecular orbitals
109 : !> \param deocc_orb_set list of de-occupied orbital indices
110 : !> \param occ_orb_set list of newly occupied orbital indices
111 : !> \param spin spin component of the molecular orbitals;
112 : !> to be used for diagnostic messages
113 : !> \par History
114 : !> 01.2016 created [Sergey Chulkov]
115 : ! **************************************************************************************************
116 40 : SUBROUTINE mom_reoccupy_orbitals(mo_set, deocc_orb_set, occ_orb_set, spin)
117 : TYPE(mo_set_type), INTENT(INOUT) :: mo_set
118 : INTEGER, DIMENSION(:), POINTER :: deocc_orb_set, occ_orb_set
119 : CHARACTER(len=*), INTENT(in) :: spin
120 :
121 : CHARACTER(len=*), PARAMETER :: routineN = 'mom_reoccupy_orbitals'
122 :
123 : CHARACTER(len=10) :: str_iorb, str_norbs
124 : CHARACTER(len=3) :: str_prefix
125 : INTEGER :: handle, homo, iorb, lfomo, nao, nmo, &
126 : norbs
127 : REAL(kind=dp) :: maxocc
128 40 : REAL(kind=dp), DIMENSION(:), POINTER :: occ_nums
129 :
130 40 : CALL timeset(routineN, handle)
131 :
132 : ! MOM electron excitation should preserve both the number of electrons and
133 : ! multiplicity of the electronic system thus ensuring the following constraint :
134 : ! norbs = SIZE(deocc_orb_set) == SIZE(occ_orb_set)
135 40 : norbs = SIZE(deocc_orb_set)
136 :
137 : ! the following assertion should never raise an exception
138 40 : CPASSERT(SIZE(deocc_orb_set) == SIZE(occ_orb_set))
139 :
140 : ! MOM does not follow aufbau principle producing non-uniformly occupied orbitals
141 40 : CALL set_mo_set(mo_set=mo_set, uniform_occupation=.FALSE.)
142 :
143 40 : IF (deocc_orb_set(1) /= 0 .AND. occ_orb_set(1) /= 0) THEN
144 : CALL get_mo_set(mo_set=mo_set, maxocc=maxocc, &
145 20 : nao=nao, nmo=nmo, occupation_numbers=occ_nums)
146 :
147 20 : IF (deocc_orb_set(norbs) > nao .OR. occ_orb_set(norbs) > nao) THEN
148 : ! STOP: one of the molecular orbital index exceeds the number of atomic basis functions available
149 0 : CALL integer_to_string(nao, str_norbs)
150 :
151 0 : IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
152 0 : iorb = deocc_orb_set(norbs)
153 0 : str_prefix = 'de-'
154 : ELSE
155 0 : iorb = occ_orb_set(norbs)
156 0 : str_prefix = ''
157 : END IF
158 0 : CALL integer_to_string(iorb, str_iorb)
159 :
160 : CALL cp_abort(__LOCATION__, "Unable to "//TRIM(str_prefix)//"occupy "// &
161 : TRIM(spin)//" orbital No. "//TRIM(str_iorb)// &
162 : " since its index exceeds the number of atomic orbital functions available ("// &
163 0 : TRIM(str_norbs)//"). Please consider using a larger basis set.")
164 : END IF
165 :
166 20 : IF (deocc_orb_set(norbs) > nmo .OR. occ_orb_set(norbs) > nmo) THEN
167 : ! STOP: one of the molecular orbital index exceeds the number of constructed molecular orbitals
168 0 : IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
169 0 : iorb = deocc_orb_set(norbs)
170 : ELSE
171 0 : iorb = occ_orb_set(norbs)
172 : END IF
173 :
174 0 : IF (iorb - nmo > 1) THEN
175 0 : CALL integer_to_string(iorb - nmo, str_iorb)
176 0 : str_prefix = 's'
177 : ELSE
178 0 : str_iorb = 'an'
179 0 : str_prefix = ''
180 : END IF
181 :
182 0 : CALL integer_to_string(nmo, str_norbs)
183 :
184 : CALL cp_abort(__LOCATION__, "The number of molecular orbitals ("//TRIM(str_norbs)// &
185 : ") is not enough to perform MOM calculation. Please add "// &
186 : TRIM(str_iorb)//" extra orbital"//TRIM(str_prefix)// &
187 0 : " using the ADDED_MOS keyword in the SCF section of your input file.")
188 : END IF
189 :
190 40 : DO iorb = 1, norbs
191 : ! swap occupation numbers between two adjoint molecular orbitals
192 20 : IF (occ_nums(deocc_orb_set(iorb)) <= 0.0_dp) THEN
193 0 : CALL integer_to_string(deocc_orb_set(iorb), str_iorb)
194 :
195 : CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
196 0 : TRIM(str_iorb)//" is not occupied thus it cannot be deoccupied.")
197 : END IF
198 :
199 20 : IF (occ_nums(occ_orb_set(iorb)) > 0.0_dp) THEN
200 0 : CALL integer_to_string(occ_orb_set(iorb), str_iorb)
201 :
202 : CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
203 0 : TRIM(str_iorb)//" is already occupied thus it cannot be reoccupied.")
204 : END IF
205 :
206 20 : occ_nums(occ_orb_set(iorb)) = occ_nums(deocc_orb_set(iorb))
207 40 : occ_nums(deocc_orb_set(iorb)) = 0.0_dp
208 : END DO
209 :
210 : ! locate the lowest non-maxocc occupied orbital
211 78 : DO lfomo = 1, nmo
212 78 : IF (occ_nums(lfomo) /= maxocc) EXIT
213 : END DO
214 :
215 : ! locate the highest occupied orbital
216 90 : DO homo = nmo, 1, -1
217 90 : IF (occ_nums(homo) > 0.0_dp) EXIT
218 : END DO
219 :
220 20 : CALL set_mo_set(mo_set=mo_set, homo=homo, lfomo=lfomo)
221 :
222 20 : ELSE IF (deocc_orb_set(1) /= 0 .OR. occ_orb_set(1) /= 0) THEN
223 : CALL cp_abort(__LOCATION__, &
224 0 : "Incorrect multiplicity of the MOM reference electronic state")
225 : END IF
226 :
227 40 : CALL timestop(handle)
228 :
229 40 : END SUBROUTINE mom_reoccupy_orbitals
230 :
231 : ! **************************************************************************************************
232 : !> \brief initial guess for the maximum overlap method
233 : !> \param nspins number of spin components
234 : !> \param mos array of molecular orbitals
235 : !> \param scf_control SCF control variables
236 : !> \param p_rmpv density matrix to be computed
237 : !> \par History
238 : !> * 01.2016 created [Sergey Chulkov]
239 : ! **************************************************************************************************
240 20 : SUBROUTINE do_mom_guess(nspins, mos, scf_control, p_rmpv)
241 : INTEGER, INTENT(in) :: nspins
242 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
243 : TYPE(scf_control_type), POINTER :: scf_control
244 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
245 :
246 : CHARACTER(len=*), PARAMETER :: routineN = 'do_mom_guess'
247 :
248 : CHARACTER(len=10) :: str_iter
249 : INTEGER :: handle, ispin, scf_iter
250 : LOGICAL :: is_mo
251 : REAL(kind=dp) :: maxa
252 : TYPE(cp_fm_type), POINTER :: mo_coeff
253 :
254 20 : CALL timeset(routineN, handle)
255 :
256 : ! we are about to initialise the maximum overlap method,
257 : ! so cite the relevant reference first
258 20 : IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
259 20 : CALL cite_reference(Gilbert2008)
260 0 : ELSE IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
261 0 : CALL cite_reference(Barca2018)
262 : END IF
263 :
264 : ! ensure we do not have duplicated orbital indices
265 20 : IF (.NOT. &
266 : (mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccA) .AND. &
267 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccB) .AND. &
268 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occA) .AND. &
269 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occB))) THEN
270 : CALL cp_abort(__LOCATION__, &
271 0 : "Duplicate orbital indices were found in the MOM section")
272 : END IF
273 :
274 : ! ignore beta orbitals for spin-unpolarized calculations
275 20 : IF (nspins == 1 .AND. (scf_control%diagonalization%mom_deoccB(1) /= 0 &
276 : .OR. scf_control%diagonalization%mom_occB(1) /= 0)) THEN
277 :
278 : CALL cp_warn(__LOCATION__, "Maximum overlap method will"// &
279 0 : " ignore beta orbitals since neither UKS nor ROKS calculation is performed")
280 : END IF
281 :
282 : ! compute the change in multiplicity and number of electrons
283 : IF (SIZE(scf_control%diagonalization%mom_deoccA) /= &
284 20 : SIZE(scf_control%diagonalization%mom_occA) .OR. &
285 : (nspins > 1 .AND. &
286 : SIZE(scf_control%diagonalization%mom_deoccB) /= &
287 : SIZE(scf_control%diagonalization%mom_occB))) THEN
288 :
289 : CALL cp_abort(__LOCATION__, "Incorrect multiplicity of the MOM reference"// &
290 0 : " electronic state or inconsistent number of electrons")
291 : END IF
292 :
293 20 : is_mo = .FALSE.
294 : ! by default activate MOM at the second SCF iteration as the
295 : ! 'old' molecular orbitals are unavailable from the very beginning
296 20 : scf_iter = 2
297 : ! check if the molecular orbitals are actually there
298 : ! by finding at least one MO coefficient > 0
299 28 : DO ispin = 1, nspins
300 24 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
301 24 : CALL cp_fm_maxabsval(mo_coeff, maxa)
302 : ! is_mo |= maxa > 0.0_dp
303 28 : IF (maxa > 0.0_dp) THEN
304 16 : is_mo = .TRUE.
305 : ! we already have the molecular orbitals (e.g. from a restart file);
306 : ! activate MOM immediately if the input keyword START_ITER is not given
307 16 : scf_iter = 1
308 16 : EXIT
309 : END IF
310 : END DO
311 :
312 : ! proceed alpha orbitals
313 20 : IF (nspins >= 1) THEN
314 : CALL mom_reoccupy_orbitals(mos(1), &
315 : scf_control%diagonalization%mom_deoccA, &
316 20 : scf_control%diagonalization%mom_occA, 'alpha')
317 : END IF
318 :
319 : ! proceed beta orbitals (if any)
320 20 : IF (nspins >= 2) THEN
321 : CALL mom_reoccupy_orbitals(mos(2), &
322 : scf_control%diagonalization%mom_deoccB, &
323 20 : scf_control%diagonalization%mom_occB, 'beta')
324 : END IF
325 :
326 : ! recompute the density matrix if the molecular orbitals are here;
327 : ! otherwise do nothing to prevent zeroing out the density matrix
328 : ! obtained from atomic guess
329 20 : IF (is_mo) THEN
330 48 : DO ispin = 1, nspins
331 48 : CALL calculate_density_matrix(mos(ispin), p_rmpv(ispin)%matrix)
332 : END DO
333 : END IF
334 :
335 : ! adjust the start SCF iteration number if needed
336 20 : IF (scf_control%diagonalization%mom_start < scf_iter) THEN
337 18 : IF (scf_control%diagonalization%mom_start > 0) THEN
338 : ! inappropriate iteration number has been provided through the input file;
339 : ! fix it and issue a warning message
340 0 : CALL integer_to_string(scf_iter, str_iter)
341 : CALL cp_warn(__LOCATION__, &
342 : "The maximum overlap method will be activated at the SCF iteration No. "// &
343 0 : TRIM(str_iter)//" due to the SCF guess method used.")
344 : END IF
345 18 : scf_control%diagonalization%mom_start = scf_iter
346 2 : ELSE IF (scf_control%diagonalization%mom_start > scf_iter .AND. &
347 : (scf_control%diagonalization%mom_occA(1) > 0 .OR. scf_control%diagonalization%mom_occB(1) > 0)) THEN
348 : ! the keyword START_ITER has been provided for an excited state calculation, ignore it
349 2 : CALL integer_to_string(scf_iter, str_iter)
350 : CALL cp_warn(__LOCATION__, &
351 : "The maximum overlap method will be activated at the SCF iteration No. "// &
352 2 : TRIM(str_iter)//" because an excited state calculation has been requested")
353 2 : scf_control%diagonalization%mom_start = scf_iter
354 : END IF
355 :
356 : ! MOM is now initialised properly
357 20 : scf_control%diagonalization%mom_didguess = .TRUE.
358 :
359 20 : CALL timestop(handle)
360 :
361 20 : END SUBROUTINE do_mom_guess
362 :
363 : ! **************************************************************************************************
364 : !> \brief do an SCF iteration, then compute occupation numbers of the new
365 : !> molecular orbitals according to their overlap with the previous ones
366 : !> \param scf_env SCF environment information
367 : !> \param mos array of molecular orbitals
368 : !> \param matrix_ks sparse Kohn-Sham matrix
369 : !> \param matrix_s sparse overlap matrix
370 : !> \param scf_control SCF control variables
371 : !> \param scf_section SCF input section
372 : !> \param diis_step have we done a DIIS step
373 : !> \par History
374 : !> * 07.2014 created [Matt Watkins]
375 : !> * 01.2016 release version [Sergey Chulkov]
376 : !> * 03.2018 initial maximum overlap method [Sergey Chulkov]
377 : ! **************************************************************************************************
378 324 : SUBROUTINE do_mom_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
379 : TYPE(qs_scf_env_type), POINTER :: scf_env
380 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
381 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
382 : TYPE(scf_control_type), POINTER :: scf_control
383 : TYPE(section_vals_type), POINTER :: scf_section
384 : LOGICAL, INTENT(INOUT) :: diis_step
385 :
386 : CHARACTER(len=*), PARAMETER :: routineN = 'do_mom_diag'
387 :
388 : INTEGER :: handle, homo, iproj, ispin, lfomo, nao, &
389 : nmo, nspins
390 324 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
391 : REAL(kind=dp) :: maxocc
392 324 : REAL(kind=dp), DIMENSION(:), POINTER :: occ_nums, proj, tmp_occ_nums
393 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
394 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_fmstruct, mo_mo_fmstruct
395 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_ref, overlap, svec
396 :
397 324 : CALL timeset(routineN, handle)
398 :
399 324 : IF (.NOT. scf_control%diagonalization%mom_didguess) THEN
400 : CALL cp_abort(__LOCATION__, &
401 0 : "The current implementation of the maximum overlap method is incompatible with the initial SCF guess")
402 : END IF
403 :
404 : ! number of spins == dft_control%nspins
405 324 : nspins = SIZE(matrix_ks)
406 :
407 : ! copy old molecular orbitals
408 324 : IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
409 318 : IF (.NOT. ASSOCIATED(scf_env%mom_ref_mo_coeff)) THEN
410 110 : ALLOCATE (scf_env%mom_ref_mo_coeff(nspins))
411 66 : DO ispin = 1, nspins
412 44 : NULLIFY (ao_mo_fmstruct)
413 44 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
414 44 : CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
415 44 : CALL cp_fm_create(scf_env%mom_ref_mo_coeff(ispin), ao_mo_fmstruct)
416 :
417 : ! Initial Maximum Overlap Method: keep initial molecular orbitals
418 66 : IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
419 0 : CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
420 0 : CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
421 : END IF
422 : END DO
423 : END IF
424 :
425 : ! allocate the molecular orbitals overlap matrix
426 318 : IF (.NOT. ASSOCIATED(scf_env%mom_overlap)) THEN
427 100 : ALLOCATE (scf_env%mom_overlap(nspins))
428 60 : DO ispin = 1, nspins
429 40 : NULLIFY (blacs_env, mo_mo_fmstruct)
430 40 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, mo_coeff=mo_coeff)
431 40 : CALL cp_fm_get_info(mo_coeff, context=blacs_env)
432 40 : CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, context=blacs_env)
433 40 : CALL cp_fm_create(scf_env%mom_overlap(ispin), mo_mo_fmstruct)
434 100 : CALL cp_fm_struct_release(mo_mo_fmstruct)
435 : END DO
436 : END IF
437 :
438 : ! allocate a matrix to store the product S * mo_coeff
439 318 : IF (.NOT. ASSOCIATED(scf_env%mom_s_mo_coeff)) THEN
440 100 : ALLOCATE (scf_env%mom_s_mo_coeff(nspins))
441 60 : DO ispin = 1, nspins
442 40 : NULLIFY (ao_mo_fmstruct)
443 40 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
444 40 : CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
445 60 : CALL cp_fm_create(scf_env%mom_s_mo_coeff(ispin), ao_mo_fmstruct)
446 : END DO
447 : END IF
448 :
449 : ! Original Maximum Overlap Method: keep orbitals from the previous SCF iteration
450 318 : IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
451 954 : DO ispin = 1, nspins
452 636 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
453 636 : CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
454 954 : CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
455 : END DO
456 : END IF
457 : END IF
458 :
459 : ! solve the eigenproblem
460 324 : CALL general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
461 :
462 324 : IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
463 954 : DO ispin = 1, nspins
464 :
465 : ! TO DO: sparse-matrix variant; check if use_mo_coeff_b is set, and if yes use mo_coeff_b instead
466 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc, mo_coeff=mo_coeff, &
467 636 : nao=nao, nmo=nmo, occupation_numbers=occ_nums)
468 :
469 636 : mo_coeff_ref => scf_env%mom_ref_mo_coeff(ispin)
470 636 : overlap => scf_env%mom_overlap(ispin)
471 636 : svec => scf_env%mom_s_mo_coeff(ispin)
472 :
473 : ! svec = S * C(new)
474 636 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, svec, nmo)
475 :
476 : ! overlap = C(reference occupied)^T * S * C(new)
477 636 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff_ref, svec, 0.0_dp, overlap)
478 :
479 1908 : ALLOCATE (proj(nmo))
480 1908 : ALLOCATE (inds(nmo))
481 1272 : ALLOCATE (tmp_occ_nums(nmo))
482 :
483 : ! project the new molecular orbitals into the space of the reference occupied orbitals
484 636 : SELECT CASE (scf_control%diagonalization%mom_proj_formula)
485 : CASE (momproj_sum)
486 : ! proj_j = abs( \sum_i overlap(i, j) )
487 0 : CALL cp_fm_vectorssum(overlap, proj)
488 :
489 0 : DO iproj = 1, nmo
490 0 : proj(iproj) = ABS(proj(iproj))
491 : END DO
492 :
493 : CASE (momproj_norm)
494 : ! proj_j = (\sum_i overlap(i, j)**2) ** 0.5
495 636 : CALL cp_fm_vectorsnorm(overlap, proj)
496 :
497 : CASE DEFAULT
498 636 : CPABORT("Unimplemented projection formula")
499 : END SELECT
500 :
501 11688 : tmp_occ_nums(:) = occ_nums(:)
502 : ! sort occupation numbers in ascending order
503 636 : CALL sort(tmp_occ_nums, nmo, inds)
504 : ! sort overlap projection in ascending order
505 636 : CALL sort(proj, nmo, inds)
506 :
507 : ! reorder occupation numbers according to overlap projections
508 5844 : DO iproj = 1, nmo
509 5844 : occ_nums(inds(iproj)) = tmp_occ_nums(iproj)
510 : END DO
511 :
512 636 : DEALLOCATE (tmp_occ_nums)
513 636 : DEALLOCATE (inds)
514 636 : DEALLOCATE (proj)
515 :
516 : ! locate the lowest non-fully occupied orbital
517 2882 : DO lfomo = 1, nmo
518 2882 : IF (occ_nums(lfomo) /= maxocc) EXIT
519 : END DO
520 :
521 : ! locate the highest occupied orbital
522 2918 : DO homo = nmo, 1, -1
523 2918 : IF (occ_nums(homo) > 0.0_dp) EXIT
524 : END DO
525 :
526 1590 : CALL set_mo_set(mo_set=mos(ispin), homo=homo, lfomo=lfomo)
527 : END DO
528 : END IF
529 :
530 : ! recompute density matrix
531 972 : DO ispin = 1, nspins
532 972 : CALL calculate_density_matrix(mos(ispin), scf_env%p_mix_new(ispin, 1)%matrix)
533 : END DO
534 :
535 324 : CALL timestop(handle)
536 :
537 648 : END SUBROUTINE do_mom_diag
538 :
539 : END MODULE qs_mom_methods
|
