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 :
8 : ! **************************************************************************************************
9 : !> \brief Does all kind of post scf calculations for GPW/GAPW
10 : !> \par History
11 : !> Started as a copy from the relevant part of qs_scf
12 : !> \author Joost VandeVondele (10.2003)
13 : ! **************************************************************************************************
14 : MODULE qs_scf_wfn_mix
15 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
16 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
17 : USE cp_files, ONLY: close_file,&
18 : open_file
19 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
20 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
21 : cp_fm_struct_release,&
22 : cp_fm_struct_type
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_release,&
25 : cp_fm_to_fm,&
26 : cp_fm_type
27 : USE input_constants, ONLY: wfn_mix_orig_external,&
28 : wfn_mix_orig_occ,&
29 : wfn_mix_orig_virtual
30 : USE input_section_types, ONLY: section_vals_get,&
31 : section_vals_get_subs_vals,&
32 : section_vals_type,&
33 : section_vals_val_get
34 : USE kinds, ONLY: default_path_length,&
35 : dp
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE particle_types, ONLY: particle_type
38 : USE qs_kind_types, ONLY: qs_kind_type
39 : USE qs_mo_io, ONLY: read_mos_restart_low,&
40 : write_mo_set_to_restart
41 : USE qs_mo_methods, ONLY: calculate_orthonormality
42 : USE qs_mo_types, ONLY: deallocate_mo_set,&
43 : duplicate_mo_set,&
44 : get_mo_set,&
45 : mo_set_type
46 : USE qs_scf_types, ONLY: qs_scf_env_type,&
47 : special_diag_method_nr
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 : PRIVATE
52 :
53 : ! Global parameters
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_wfn_mix'
55 : PUBLIC :: wfn_mix
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief writes a new 'mixed' set of mos to restart file, without touching the current MOs
61 : !> \param mos ...
62 : !> \param particle_set ...
63 : !> \param dft_section ...
64 : !> \param qs_kind_set ...
65 : !> \param para_env ...
66 : !> \param output_unit ...
67 : !> \param unoccupied_orbs ...
68 : !> \param scf_env ...
69 : !> \param matrix_s ...
70 : !> \param marked_states ...
71 : !> \param for_rtp ...
72 : ! **************************************************************************************************
73 9601 : SUBROUTINE wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, output_unit, &
74 : unoccupied_orbs, scf_env, matrix_s, marked_states, for_rtp)
75 :
76 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
77 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
78 : TYPE(section_vals_type), POINTER :: dft_section
79 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
80 : TYPE(mp_para_env_type), POINTER :: para_env
81 : INTEGER :: output_unit
82 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
83 : OPTIONAL, POINTER :: unoccupied_orbs
84 : TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env
85 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
86 : POINTER :: matrix_s
87 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: marked_states
88 : LOGICAL, OPTIONAL :: for_rtp
89 :
90 : CHARACTER(len=*), PARAMETER :: routineN = 'wfn_mix'
91 :
92 : CHARACTER(LEN=default_path_length) :: read_file_name
93 : INTEGER :: handle, homo, i_rep, ispin, mark_ind, mark_number, n_rep, nmo, orig_mo_index, &
94 : orig_spin_index, orig_type, restart_unit, result_mo_index, result_spin_index
95 : LOGICAL :: explicit, is_file, my_for_rtp, &
96 : overwrite_mos, reverse_mo_index
97 : REAL(KIND=dp) :: orig_scale, orthonormality, result_scale
98 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_vector
99 : TYPE(cp_fm_type) :: matrix_x, matrix_y
100 9601 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_new, mos_orig_ext
101 : TYPE(section_vals_type), POINTER :: update_section, wfn_mix_section
102 :
103 9601 : CALL timeset(routineN, handle)
104 9601 : wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
105 9601 : CALL section_vals_get(wfn_mix_section, explicit=explicit)
106 :
107 : ! only perform action if explicitly required
108 9601 : IF (explicit) THEN
109 :
110 24 : IF (PRESENT(for_rtp)) THEN
111 4 : my_for_rtp = for_rtp
112 : ELSE
113 : my_for_rtp = .FALSE.
114 : END IF
115 :
116 24 : IF (output_unit > 0) THEN
117 12 : WRITE (output_unit, '()')
118 12 : WRITE (output_unit, '(T2,A)') "Performing wfn mixing"
119 12 : WRITE (output_unit, '(T2,A)') "====================="
120 : END IF
121 :
122 112 : ALLOCATE (mos_new(SIZE(mos)))
123 64 : DO ispin = 1, SIZE(mos)
124 64 : CALL duplicate_mo_set(mos_new(ispin), mos(ispin))
125 : END DO
126 :
127 : ! a single vector matrix structure
128 24 : NULLIFY (fm_struct_vector)
129 : CALL cp_fm_struct_create(fm_struct_vector, template_fmstruct=mos(1)%mo_coeff%matrix_struct, &
130 24 : ncol_global=1)
131 24 : CALL cp_fm_create(matrix_x, fm_struct_vector, name="x")
132 24 : CALL cp_fm_create(matrix_y, fm_struct_vector, name="y")
133 24 : CALL cp_fm_struct_release(fm_struct_vector)
134 :
135 24 : update_section => section_vals_get_subs_vals(wfn_mix_section, "UPDATE")
136 24 : CALL section_vals_get(update_section, n_repetition=n_rep)
137 24 : CALL section_vals_get(update_section, explicit=explicit)
138 24 : IF (.NOT. explicit) n_rep = 0
139 :
140 : ! Mix the MOs as : y = ay + bx
141 54 : DO i_rep = 1, n_rep
142 : ! The occupied MO that will be modified or saved, 'y'
143 30 : CALL section_vals_val_get(update_section, "RESULT_MO_INDEX", i_rep_section=i_rep, i_val=result_mo_index)
144 30 : CALL section_vals_val_get(update_section, "RESULT_MARKED_STATE", i_rep_section=i_rep, i_val=mark_number)
145 30 : CALL section_vals_val_get(update_section, "RESULT_SPIN_INDEX", i_rep_section=i_rep, i_val=result_spin_index)
146 : ! result_scale is the 'a' coefficient
147 30 : CALL section_vals_val_get(update_section, "RESULT_SCALE", i_rep_section=i_rep, r_val=result_scale)
148 :
149 30 : mark_ind = 1
150 30 : IF (mark_number .GT. 0) result_mo_index = marked_states(mark_number, result_spin_index, mark_ind)
151 :
152 : ! The MO that will be added to the previous one, 'x'
153 : CALL section_vals_val_get(update_section, "ORIG_TYPE", i_rep_section=i_rep, &
154 30 : i_val=orig_type)
155 30 : CALL section_vals_val_get(update_section, "ORIG_MO_INDEX", i_rep_section=i_rep, i_val=orig_mo_index)
156 30 : CALL section_vals_val_get(update_section, "ORIG_MARKED_STATE", i_rep_section=i_rep, i_val=mark_number)
157 30 : CALL section_vals_val_get(update_section, "ORIG_SPIN_INDEX", i_rep_section=i_rep, i_val=orig_spin_index)
158 : ! orig_scale is the 'b' coefficient
159 30 : CALL section_vals_val_get(update_section, "ORIG_SCALE", i_rep_section=i_rep, r_val=orig_scale)
160 :
161 30 : IF (orig_type == wfn_mix_orig_virtual) mark_ind = 2
162 30 : IF (mark_number .GT. 0) orig_mo_index = marked_states(mark_number, orig_spin_index, mark_ind)
163 :
164 30 : CALL section_vals_val_get(wfn_mix_section, "OVERWRITE_MOS", l_val=overwrite_mos)
165 :
166 30 : CALL section_vals_val_get(update_section, "REVERSE_MO_INDEX", i_rep_section=i_rep, l_val=reverse_mo_index)
167 :
168 : ! First get a copy of the proper orig
169 : ! Origin is in an occupied MO
170 30 : IF (orig_type == wfn_mix_orig_occ) THEN
171 4 : IF (reverse_mo_index) THEN
172 : CALL cp_fm_to_fm(mos(orig_spin_index)%mo_coeff, matrix_x, 1, &
173 0 : orig_mo_index, 1)
174 : ELSE
175 : CALL cp_fm_to_fm(mos(orig_spin_index)%mo_coeff, matrix_x, 1, &
176 4 : mos(orig_spin_index)%nmo - orig_mo_index + 1, 1)
177 : END IF
178 : ! Origin is an unoccupied MO
179 26 : ELSE IF (orig_type == wfn_mix_orig_virtual) THEN
180 22 : CALL get_mo_set(mos(orig_spin_index), homo=homo, nmo=nmo)
181 : ! check whether the MO is already available in the original MO set
182 22 : IF (orig_mo_index + homo <= nmo) THEN
183 0 : CALL cp_fm_to_fm(mos(orig_spin_index)%mo_coeff, matrix_x, 1, orig_mo_index + homo, 1)
184 22 : ELSE IF (.NOT. ASSOCIATED(unoccupied_orbs)) THEN
185 : CALL cp_abort(__LOCATION__, &
186 : "If ORIG_TYPE is set to VIRTUAL, the array unoccupied_orbs must be associated! "// &
187 0 : "For instance, ask in the SCF section to compute virtual orbitals after the GS optimization.")
188 : ELSE
189 22 : CALL cp_fm_to_fm(unoccupied_orbs(orig_spin_index), matrix_x, 1, orig_mo_index, 1)
190 : END IF
191 :
192 : ! Origin is to be read from an external .wfn file
193 4 : ELSE IF (orig_type == wfn_mix_orig_external) THEN
194 : CALL section_vals_val_get(update_section, "ORIG_EXT_FILE_NAME", i_rep_section=i_rep, &
195 4 : c_val=read_file_name)
196 4 : IF (read_file_name == "EMPTY") &
197 : CALL cp_abort(__LOCATION__, &
198 : "If ORIG_TYPE is set to EXTERNAL, a file name should be set in ORIG_EXT_FILE_NAME "// &
199 0 : "so that it can be used as the orginal MO.")
200 :
201 18 : ALLOCATE (mos_orig_ext(SIZE(mos)))
202 10 : DO ispin = 1, SIZE(mos)
203 10 : CALL duplicate_mo_set(mos_orig_ext(ispin), mos(ispin))
204 : END DO
205 :
206 4 : IF (para_env%is_source()) THEN
207 2 : INQUIRE (FILE=TRIM(read_file_name), exist=is_file)
208 2 : IF (.NOT. is_file) &
209 : CALL cp_abort(__LOCATION__, &
210 0 : "Reference file not found! Name of the file CP2K looked for: "//TRIM(read_file_name))
211 :
212 : CALL open_file(file_name=read_file_name, &
213 : file_action="READ", &
214 : file_form="UNFORMATTED", &
215 : file_status="OLD", &
216 2 : unit_number=restart_unit)
217 : END IF
218 : CALL read_mos_restart_low(mos_orig_ext, para_env=para_env, qs_kind_set=qs_kind_set, &
219 : particle_set=particle_set, natom=SIZE(particle_set, 1), &
220 4 : rst_unit=restart_unit)
221 4 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
222 :
223 4 : IF (reverse_mo_index) THEN
224 : CALL cp_fm_to_fm(mos_orig_ext(orig_spin_index)%mo_coeff, matrix_x, 1, &
225 0 : orig_mo_index, 1)
226 : ELSE
227 : CALL cp_fm_to_fm(mos_orig_ext(orig_spin_index)%mo_coeff, matrix_x, 1, &
228 4 : mos_orig_ext(orig_spin_index)%nmo - orig_mo_index + 1, 1)
229 : END IF
230 10 : DO ispin = 1, SIZE(mos_orig_ext)
231 10 : CALL deallocate_mo_set(mos_orig_ext(ispin))
232 : END DO
233 4 : DEALLOCATE (mos_orig_ext)
234 : END IF
235 :
236 : ! Second, get a copy of the target
237 30 : IF (reverse_mo_index) THEN
238 : CALL cp_fm_to_fm(mos_new(result_spin_index)%mo_coeff, matrix_y, &
239 0 : 1, result_mo_index, 1)
240 : ELSE
241 : CALL cp_fm_to_fm(mos_new(result_spin_index)%mo_coeff, matrix_y, &
242 30 : 1, mos_new(result_spin_index)%nmo - result_mo_index + 1, 1)
243 : END IF
244 :
245 : ! Third, perform the mix
246 30 : CALL cp_fm_scale_and_add(result_scale, matrix_y, orig_scale, matrix_x)
247 :
248 : ! and copy back in the result mos
249 144 : IF (reverse_mo_index) THEN
250 : CALL cp_fm_to_fm(matrix_y, mos_new(result_spin_index)%mo_coeff, &
251 0 : 1, 1, result_mo_index)
252 : ELSE
253 : CALL cp_fm_to_fm(matrix_y, mos_new(result_spin_index)%mo_coeff, &
254 30 : 1, 1, mos_new(result_spin_index)%nmo - result_mo_index + 1)
255 : END IF
256 : END DO
257 :
258 24 : CALL cp_fm_release(matrix_x)
259 24 : CALL cp_fm_release(matrix_y)
260 :
261 24 : IF (my_for_rtp) THEN
262 12 : DO ispin = 1, SIZE(mos_new)
263 8 : CALL cp_fm_to_fm(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff)
264 8 : IF (mos_new(1)%use_mo_coeff_b) &
265 0 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos_new(ispin)%mo_coeff_b)
266 8 : IF (mos(1)%use_mo_coeff_b) &
267 4 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff_b)
268 : END DO
269 : ELSE
270 20 : IF (scf_env%method == special_diag_method_nr) THEN
271 0 : CALL calculate_orthonormality(orthonormality, mos)
272 : ELSE
273 20 : CALL calculate_orthonormality(orthonormality, mos, matrix_s(1)%matrix)
274 : END IF
275 :
276 20 : IF (output_unit > 0) THEN
277 10 : WRITE (output_unit, '()')
278 : WRITE (output_unit, '(T2,A,T61,E20.4)') &
279 10 : "Maximum deviation from MO S-orthonormality", orthonormality
280 10 : WRITE (output_unit, '(T2,A)') "Writing new MOs to file"
281 : END IF
282 :
283 : ! *** Write WaveFunction restart file ***
284 :
285 52 : DO ispin = 1, SIZE(mos_new)
286 32 : IF (overwrite_mos) THEN
287 12 : CALL cp_fm_to_fm(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff)
288 12 : IF (mos_new(1)%use_mo_coeff_b) &
289 4 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos_new(ispin)%mo_coeff_b)
290 : END IF
291 32 : IF (mos(1)%use_mo_coeff_b) &
292 24 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff_b)
293 : END DO
294 20 : CALL write_mo_set_to_restart(mos_new, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
295 : END IF
296 :
297 64 : DO ispin = 1, SIZE(mos_new)
298 64 : CALL deallocate_mo_set(mos_new(ispin))
299 : END DO
300 72 : DEALLOCATE (mos_new)
301 :
302 : END IF
303 :
304 9601 : CALL timestop(handle)
305 :
306 19202 : END SUBROUTINE wfn_mix
307 :
308 : END MODULE qs_scf_wfn_mix
|
