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