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_operations, ONLY: copy_fm_to_dbcsr
16 : USE cp_files, ONLY: close_file,&
17 : open_file
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
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_release,&
24 : cp_fm_to_fm,&
25 : cp_fm_type
26 : USE dbcsr_api, ONLY: dbcsr_p_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 9365 : 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
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 9365 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_new, mos_orig_ext
100 : TYPE(section_vals_type), POINTER :: update_section, wfn_mix_section
101 :
102 9365 : CALL timeset(routineN, handle)
103 9365 : wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
104 9365 : CALL section_vals_get(wfn_mix_section, explicit=explicit)
105 :
106 : ! only perform action if explicitly required
107 9365 : 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 'b' 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_scal is the 'a' 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 : ! First get a copy of the proper orig
166 : ! Origin is in the MO matrix
167 30 : IF (orig_type == wfn_mix_orig_occ) THEN
168 : CALL cp_fm_to_fm(mos(orig_spin_index)%mo_coeff, matrix_x, 1, &
169 4 : mos(orig_spin_index)%nmo - orig_mo_index + 1, 1)
170 :
171 : ! Orgin is in the virtual matrix
172 26 : ELSE IF (orig_type == wfn_mix_orig_virtual) THEN
173 22 : IF (.NOT. ASSOCIATED(unoccupied_orbs)) &
174 : CALL cp_abort(__LOCATION__, &
175 : "If ORIG_TYPE is set to VIRTUAL, the array unoccupied_orbs must be associated! "// &
176 0 : "For instance, ask in the SCF section to compute virtual orbitals after the GS optimization.")
177 22 : CALL cp_fm_to_fm(unoccupied_orbs(orig_spin_index), matrix_x, 1, orig_mo_index, 1)
178 :
179 : ! Orgin is to be read from an external .wfn file
180 4 : ELSE IF (orig_type == wfn_mix_orig_external) THEN
181 : CALL section_vals_val_get(update_section, "ORIG_EXT_FILE_NAME", i_rep_section=i_rep, &
182 4 : c_val=read_file_name)
183 4 : IF (read_file_name == "EMPTY") &
184 : CALL cp_abort(__LOCATION__, &
185 : "If ORIG_TYPE is set to EXTERNAL, a file name should be set in ORIG_EXT_FILE_NAME "// &
186 0 : "so that it can be used as the orginal MO.")
187 :
188 18 : ALLOCATE (mos_orig_ext(SIZE(mos)))
189 10 : DO ispin = 1, SIZE(mos)
190 10 : CALL duplicate_mo_set(mos_orig_ext(ispin), mos(ispin))
191 : END DO
192 :
193 4 : IF (para_env%is_source()) THEN
194 2 : INQUIRE (FILE=TRIM(read_file_name), exist=is_file)
195 2 : IF (.NOT. is_file) &
196 : CALL cp_abort(__LOCATION__, &
197 0 : "Reference file not found! Name of the file CP2K looked for: "//TRIM(read_file_name))
198 :
199 : CALL open_file(file_name=read_file_name, &
200 : file_action="READ", &
201 : file_form="UNFORMATTED", &
202 : file_status="OLD", &
203 2 : unit_number=restart_unit)
204 : END IF
205 : CALL read_mos_restart_low(mos_orig_ext, para_env=para_env, qs_kind_set=qs_kind_set, &
206 : particle_set=particle_set, natom=SIZE(particle_set, 1), &
207 4 : rst_unit=restart_unit)
208 4 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
209 :
210 : CALL cp_fm_to_fm(mos_orig_ext(orig_spin_index)%mo_coeff, matrix_x, 1, &
211 4 : mos_orig_ext(orig_spin_index)%nmo - orig_mo_index + 1, 1)
212 :
213 10 : DO ispin = 1, SIZE(mos_orig_ext)
214 10 : CALL deallocate_mo_set(mos_orig_ext(ispin))
215 : END DO
216 4 : DEALLOCATE (mos_orig_ext)
217 : END IF
218 :
219 : ! Second, get a copy of the target
220 : CALL cp_fm_to_fm(mos_new(result_spin_index)%mo_coeff, matrix_y, &
221 30 : 1, mos_new(result_spin_index)%nmo - result_mo_index + 1, 1)
222 :
223 : ! Third, perform the mix
224 30 : CALL cp_fm_scale_and_add(result_scale, matrix_y, orig_scale, matrix_x)
225 : ! and copy back in the result mos
226 : CALL cp_fm_to_fm(matrix_y, mos_new(result_spin_index)%mo_coeff, &
227 144 : 1, 1, mos_new(result_spin_index)%nmo - result_mo_index + 1)
228 : END DO
229 :
230 24 : CALL cp_fm_release(matrix_x)
231 24 : CALL cp_fm_release(matrix_y)
232 :
233 24 : IF (my_for_rtp) THEN
234 12 : DO ispin = 1, SIZE(mos_new)
235 8 : CALL cp_fm_to_fm(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff)
236 8 : IF (mos_new(1)%use_mo_coeff_b) &
237 0 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos_new(ispin)%mo_coeff_b)
238 8 : IF (mos(1)%use_mo_coeff_b) &
239 4 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff_b)
240 : END DO
241 : ELSE
242 20 : IF (scf_env%method == special_diag_method_nr) THEN
243 0 : CALL calculate_orthonormality(orthonormality, mos)
244 : ELSE
245 20 : CALL calculate_orthonormality(orthonormality, mos, matrix_s(1)%matrix)
246 : END IF
247 :
248 20 : IF (output_unit > 0) THEN
249 10 : WRITE (output_unit, '()')
250 : WRITE (output_unit, '(T2,A,T61,E20.4)') &
251 10 : "Maximum deviation from MO S-orthonormality", orthonormality
252 10 : WRITE (output_unit, '(T2,A)') "Writing new MOs to file"
253 : END IF
254 :
255 : ! *** Write WaveFunction restart file ***
256 :
257 52 : DO ispin = 1, SIZE(mos_new)
258 32 : IF (overwrite_mos) THEN
259 12 : CALL cp_fm_to_fm(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff)
260 12 : IF (mos_new(1)%use_mo_coeff_b) &
261 4 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos_new(ispin)%mo_coeff_b)
262 : END IF
263 32 : IF (mos(1)%use_mo_coeff_b) &
264 24 : CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_coeff, mos(ispin)%mo_coeff_b)
265 : END DO
266 20 : CALL write_mo_set_to_restart(mos_new, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
267 : END IF
268 :
269 64 : DO ispin = 1, SIZE(mos_new)
270 64 : CALL deallocate_mo_set(mos_new(ispin))
271 : END DO
272 48 : DEALLOCATE (mos_new)
273 :
274 : END IF
275 :
276 9365 : CALL timestop(handle)
277 :
278 18730 : END SUBROUTINE wfn_mix
279 :
280 : END MODULE qs_scf_wfn_mix
|
