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 Routines for the calculation of molecular states
10 : !> \author CJM
11 : ! **************************************************************************************************
12 : MODULE molecular_states
13 : USE atomic_kind_types, ONLY: atomic_kind_type
14 : USE bibliography, ONLY: Hunt2003,&
15 : cite_reference
16 : USE cell_types, ONLY: cell_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
19 : USE cp_fm_diag, ONLY: choose_eigv_solver
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_get_element,&
25 : cp_fm_get_info,&
26 : cp_fm_release,&
27 : cp_fm_to_fm,&
28 : cp_fm_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_get_default_io_unit,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_p_file,&
33 : cp_print_key_finished_output,&
34 : cp_print_key_should_output,&
35 : cp_print_key_unit_nr
36 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
37 : USE dbcsr_api, ONLY: dbcsr_type
38 : USE input_section_types, ONLY: section_get_ivals,&
39 : section_vals_type,&
40 : section_vals_val_get
41 : USE kinds, ONLY: default_path_length,&
42 : default_string_length,&
43 : dp
44 : USE message_passing, ONLY: mp_para_env_type
45 : USE molecule_types, ONLY: molecule_type
46 : USE parallel_gemm_api, ONLY: parallel_gemm
47 : USE particle_list_types, ONLY: particle_list_type
48 : USE particle_types, ONLY: particle_type
49 : USE pw_env_types, ONLY: pw_env_type
50 : USE pw_types, ONLY: pw_c1d_gs_type,&
51 : pw_r3d_rs_type
52 : USE qs_collocate_density, ONLY: calculate_wavefunction
53 : USE qs_environment_types, ONLY: get_qs_env,&
54 : qs_environment_type
55 : USE qs_kind_types, ONLY: qs_kind_type
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : ! *** Global parameters ***
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_states'
65 :
66 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
67 :
68 : ! *** Public subroutines ***
69 :
70 : PUBLIC :: construct_molecular_states
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief constructs molecular states. mo_localized gets overwritten!
76 : !> \param molecule_set ...
77 : !> \param mo_localized ...
78 : !> \param mo_coeff ...
79 : !> \param mo_eigenvalues ...
80 : !> \param Hks ...
81 : !> \param matrix_S ...
82 : !> \param qs_env ...
83 : !> \param wf_r ...
84 : !> \param wf_g ...
85 : !> \param loc_print_section ...
86 : !> \param particles ...
87 : !> \param tag ...
88 : !> \param marked_states ...
89 : !> \param ispin ...
90 : ! **************************************************************************************************
91 28 : SUBROUTINE construct_molecular_states(molecule_set, mo_localized, &
92 : mo_coeff, mo_eigenvalues, Hks, matrix_S, qs_env, wf_r, wf_g, &
93 : loc_print_section, particles, tag, marked_states, ispin)
94 :
95 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
96 : TYPE(cp_fm_type), INTENT(IN) :: mo_localized, mo_coeff
97 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
98 : TYPE(dbcsr_type), POINTER :: Hks, matrix_S
99 : TYPE(qs_environment_type), POINTER :: qs_env
100 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: wf_r
101 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: wf_g
102 : TYPE(section_vals_type), POINTER :: loc_print_section
103 : TYPE(particle_list_type), POINTER :: particles
104 : CHARACTER(LEN=*), INTENT(IN) :: tag
105 : INTEGER, DIMENSION(:), POINTER :: marked_states
106 : INTEGER, INTENT(IN) :: ispin
107 :
108 : CHARACTER(len=*), PARAMETER :: routineN = 'construct_molecular_states'
109 :
110 : CHARACTER(LEN=default_path_length) :: filename
111 : CHARACTER(LEN=default_string_length) :: title
112 : INTEGER :: handle, i, imol, iproc, k, n_rep, &
113 : ncol_global, nproc, nrow_global, ns, &
114 : output_unit, unit_nr, unit_report
115 28 : INTEGER, DIMENSION(:), POINTER :: ind, mark_list
116 28 : INTEGER, DIMENSION(:, :), POINTER :: mark_states
117 28 : INTEGER, POINTER :: nstates(:), states(:)
118 : LOGICAL :: explicit, mpi_io
119 : REAL(KIND=dp) :: tmp
120 28 : REAL(KIND=dp), ALLOCATABLE :: evals(:)
121 28 : REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range
122 28 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
123 : TYPE(cell_type), POINTER :: cell
124 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
125 : TYPE(cp_fm_type) :: b, c, d, D_igk, e_vectors, &
126 : rot_e_vectors, smo, storage
127 : TYPE(cp_logger_type), POINTER :: logger
128 : TYPE(dft_control_type), POINTER :: dft_control
129 : TYPE(mp_para_env_type), POINTER :: para_env
130 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
131 : TYPE(pw_env_type), POINTER :: pw_env
132 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
133 :
134 28 : CALL cite_reference(Hunt2003)
135 :
136 28 : CALL timeset(routineN, handle)
137 :
138 28 : NULLIFY (logger, mark_states, mark_list, para_env)
139 28 : logger => cp_get_default_logger()
140 : !-----------------------------------------------------------------------------
141 : ! 1.
142 : !-----------------------------------------------------------------------------
143 28 : CALL get_qs_env(qs_env, para_env=para_env)
144 28 : nproc = para_env%num_pe
145 28 : output_unit = cp_logger_get_default_io_unit(logger)
146 : CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", &
147 28 : explicit=explicit)
148 28 : IF (explicit) THEN
149 : CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", &
150 0 : r_vals=eval_range)
151 : ELSE
152 28 : ALLOCATE (eval_range(2))
153 28 : eval_range(1) = -HUGE(0.0_dp)
154 28 : eval_range(2) = +HUGE(0.0_dp)
155 : END IF
156 : CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", &
157 28 : n_rep_val=n_rep)
158 28 : IF (n_rep .GT. 0) THEN
159 84 : ALLOCATE (mark_states(2, n_rep))
160 28 : IF (.NOT. ASSOCIATED(marked_states)) THEN
161 84 : ALLOCATE (marked_states(n_rep))
162 : END IF
163 64 : DO i = 1, n_rep
164 : CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", &
165 36 : i_rep_val=i, i_vals=mark_list)
166 208 : mark_states(:, i) = mark_list(:)
167 : END DO
168 : ELSE
169 0 : NULLIFY (marked_states)
170 : END IF
171 :
172 : !-----------------------------------------------------------------------------
173 : !-----------------------------------------------------------------------------
174 : ! 2.
175 : !-----------------------------------------------------------------------------
176 : unit_report = cp_print_key_unit_nr(logger, loc_print_section, "MOLECULAR_STATES", &
177 28 : extension=".data", middle_name="Molecular_DOS", log_filename=.FALSE.)
178 28 : IF (unit_report > 0) THEN
179 14 : WRITE (unit_report, *) SIZE(mo_eigenvalues), " number of states "
180 104 : DO i = 1, SIZE(mo_eigenvalues)
181 104 : WRITE (unit_report, *) mo_eigenvalues(i)
182 : END DO
183 : END IF
184 :
185 : !-----------------------------------------------------------------------------
186 : !-----------------------------------------------------------------------------
187 : ! 3.
188 : !-----------------------------------------------------------------------------
189 : CALL cp_fm_get_info(mo_localized, &
190 : ncol_global=ncol_global, &
191 28 : nrow_global=nrow_global)
192 28 : CALL cp_fm_create(smo, mo_coeff%matrix_struct)
193 28 : CALL cp_dbcsr_sm_fm_multiply(matrix_S, mo_coeff, smo, ncol_global)
194 :
195 : !-----------------------------------------------------------------------------
196 : !-----------------------------------------------------------------------------
197 : ! 4.
198 : !-----------------------------------------------------------------------------
199 28 : ALLOCATE (nstates(2))
200 :
201 28 : CALL cp_fm_create(storage, mo_localized%matrix_struct, name='storage')
202 :
203 72 : DO imol = 1, SIZE(molecule_set)
204 44 : IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN
205 22 : nstates(1) = molecule_set(imol)%lmi(ispin)%nstates
206 : ELSE
207 22 : nstates(1) = 0
208 : END IF
209 44 : nstates(2) = para_env%mepos
210 :
211 220 : CALL para_env%maxloc(nstates)
212 :
213 44 : IF (nstates(1) == 0) CYCLE
214 44 : NULLIFY (states)
215 132 : ALLOCATE (states(nstates(1)))
216 192 : states(:) = 0
217 :
218 44 : iproc = nstates(2)
219 44 : IF (iproc == para_env%mepos) THEN
220 96 : states(:) = molecule_set(imol)%lmi(ispin)%states(:)
221 : END IF
222 : !!BCAST from here root = iproc
223 340 : CALL para_env%bcast(states, iproc)
224 :
225 44 : ns = nstates(1)
226 44 : ind => states(:)
227 132 : ALLOCATE (evals(ns))
228 :
229 44 : NULLIFY (fm_struct_tmp)
230 :
231 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, &
232 : ncol_global=ns, &
233 : para_env=mo_localized%matrix_struct%para_env, &
234 44 : context=mo_localized%matrix_struct%context)
235 :
236 44 : CALL cp_fm_create(b, fm_struct_tmp, name="b")
237 44 : CALL cp_fm_create(c, fm_struct_tmp, name="c")
238 44 : CALL cp_fm_create(rot_e_vectors, fm_struct_tmp, name="rot_e_vectors")
239 :
240 44 : CALL cp_fm_struct_release(fm_struct_tmp)
241 :
242 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ns, ncol_global=ns, &
243 : para_env=mo_localized%matrix_struct%para_env, &
244 44 : context=mo_localized%matrix_struct%context)
245 :
246 44 : CALL cp_fm_create(d, fm_struct_tmp, name="d")
247 44 : CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e_vectors")
248 44 : CALL cp_fm_struct_release(fm_struct_tmp)
249 :
250 192 : DO i = 1, ns
251 192 : CALL cp_fm_to_fm(mo_localized, b, 1, ind(i), i)
252 : END DO
253 :
254 44 : CALL cp_dbcsr_sm_fm_multiply(Hks, b, c, ns)
255 :
256 : CALL parallel_gemm('T', 'N', ns, ns, nrow_global, 1.0_dp, &
257 44 : b, c, 0.0_dp, d)
258 :
259 44 : CALL choose_eigv_solver(d, e_vectors, evals)
260 :
261 44 : IF (output_unit > 0) WRITE (output_unit, *) ""
262 44 : IF (output_unit > 0) WRITE (output_unit, *) "MOLECULE ", imol
263 44 : IF (output_unit > 0) WRITE (output_unit, *) "NUMBER OF STATES ", ns
264 44 : IF (output_unit > 0) WRITE (output_unit, *) "EIGENVALUES"
265 44 : IF (output_unit > 0) WRITE (output_unit, *) ""
266 44 : IF (output_unit > 0) WRITE (output_unit, *) "ENERGY original MO-index"
267 :
268 192 : DO k = 1, ns
269 148 : IF (ASSOCIATED(mark_states)) THEN
270 328 : DO i = 1, n_rep
271 328 : IF (imol == mark_states(1, i) .AND. k == mark_states(2, i)) marked_states(i) = ind(k)
272 : END DO
273 : END IF
274 192 : IF (output_unit > 0) WRITE (output_unit, *) evals(k), ind(k)
275 : END DO
276 44 : IF (unit_report > 0) THEN
277 22 : WRITE (unit_report, *) imol, ns, " imol, number of states"
278 96 : DO k = 1, ns
279 96 : WRITE (unit_report, *) evals(k)
280 : END DO
281 : END IF
282 :
283 : CALL parallel_gemm('N', 'N', nrow_global, ns, ns, 1.0_dp, &
284 44 : b, e_vectors, 0.0_dp, rot_e_vectors)
285 :
286 192 : DO i = 1, ns
287 192 : CALL cp_fm_to_fm(rot_e_vectors, storage, 1, i, ind(i))
288 : END DO
289 :
290 : IF (.FALSE.) THEN ! this is too much data for large systems
291 : ! compute Eq. 6 from P. Hunt et al. (CPL 376, p. 68-74)
292 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ns, &
293 : ncol_global=ncol_global, &
294 : para_env=mo_localized%matrix_struct%para_env, &
295 : context=mo_localized%matrix_struct%context)
296 : CALL cp_fm_create(D_igk, fm_struct_tmp)
297 : CALL cp_fm_struct_release(fm_struct_tmp)
298 : CALL parallel_gemm('T', 'N', ns, ncol_global, nrow_global, 1.0_dp, &
299 : rot_e_vectors, smo, 0.0_dp, D_igk)
300 : DO i = 1, ns
301 : DO k = 1, ncol_global
302 : CALL cp_fm_get_element(D_igk, i, k, tmp)
303 : IF (unit_report > 0) THEN
304 : WRITE (unit_report, *) tmp**2
305 : END IF
306 : END DO
307 : END DO
308 : CALL cp_fm_release(D_igk)
309 : END IF
310 :
311 44 : IF (BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
312 : "MOLECULAR_STATES%CUBES"), cp_p_file)) THEN
313 :
314 : CALL get_qs_env(qs_env=qs_env, &
315 : atomic_kind_set=atomic_kind_set, &
316 : qs_kind_set=qs_kind_set, &
317 : cell=cell, &
318 : dft_control=dft_control, &
319 : particle_set=particle_set, &
320 16 : pw_env=pw_env)
321 :
322 48 : DO i = 1, ns
323 32 : IF (evals(i) < eval_range(1) .OR. evals(i) > eval_range(2)) CYCLE
324 :
325 : CALL calculate_wavefunction(rot_e_vectors, i, wf_r, &
326 : wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
327 32 : pw_env)
328 :
329 32 : WRITE (filename, '(a9,I4.4,a1,I5.5,a4)') "MOLECULE_", imol, "_", i, tag
330 32 : WRITE (title, '(A,I0,A,I0,A,F14.6,A,I0)') "Mol. Eigenstate ", i, " of ", ns, " E [a.u.] = ", &
331 64 : evals(i), " Orig. index ", ind(i)
332 32 : mpi_io = .TRUE.
333 : unit_nr = cp_print_key_unit_nr(logger, loc_print_section, "MOLECULAR_STATES%CUBES", &
334 : extension=".cube", middle_name=TRIM(filename), log_filename=.FALSE., &
335 32 : mpi_io=mpi_io)
336 : CALL cp_pw_to_cube(wf_r, unit_nr, particles=particles, title=title, &
337 : stride=section_get_ivals(loc_print_section, &
338 32 : "MOLECULAR_STATES%CUBES%STRIDE"), mpi_io=mpi_io)
339 : CALL cp_print_key_finished_output(unit_nr, logger, loc_print_section, &
340 48 : "MOLECULAR_STATES%CUBES", mpi_io=mpi_io)
341 : END DO
342 : END IF
343 :
344 44 : DEALLOCATE (evals)
345 44 : CALL cp_fm_release(b)
346 44 : CALL cp_fm_release(c)
347 44 : CALL cp_fm_release(d)
348 44 : CALL cp_fm_release(e_vectors)
349 44 : CALL cp_fm_release(rot_e_vectors)
350 :
351 72 : DEALLOCATE (states)
352 :
353 : END DO
354 28 : CALL cp_fm_release(smo)
355 28 : CALL cp_fm_to_fm(storage, mo_localized)
356 28 : CALL cp_fm_release(storage)
357 28 : IF (ASSOCIATED(mark_states)) THEN
358 28 : DEALLOCATE (mark_states)
359 : END IF
360 28 : DEALLOCATE (nstates)
361 : CALL cp_print_key_finished_output(unit_report, logger, loc_print_section, &
362 28 : "MOLECULAR_STATES")
363 : !------------------------------------------------------------------------------
364 :
365 28 : IF (.NOT. explicit) THEN
366 28 : DEALLOCATE (eval_range)
367 : END IF
368 :
369 28 : CALL timestop(handle)
370 :
371 140 : END SUBROUTINE construct_molecular_states
372 :
373 : END MODULE molecular_states
|
