(git:58e3e09)
molecular_states.F
Go to the documentation of this file.
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 ! **************************************************************************************************
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
22  cp_fm_struct_type
23  USE cp_fm_types, ONLY: cp_fm_create,&
26  cp_fm_release,&
27  cp_fm_to_fm,&
28  cp_fm_type
31  cp_logger_type
32  USE cp_output_handling, ONLY: cp_p_file,&
37  USE dbcsr_api, ONLY: dbcsr_type
39  section_vals_type,&
41  USE kinds, ONLY: default_path_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
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 
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  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  INTEGER, DIMENSION(:), POINTER :: ind, mark_list
116  INTEGER, DIMENSION(:, :), POINTER :: mark_states
117  INTEGER, POINTER :: nstates(:), states(:)
118  LOGICAL :: explicit, mpi_io
119  REAL(kind=dp) :: tmp
120  REAL(kind=dp), ALLOCATABLE :: evals(:)
121  REAL(kind=dp), DIMENSION(:), POINTER :: eval_range
122  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
131  TYPE(pw_env_type), POINTER :: pw_env
132  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
133 
134  CALL cite_reference(hunt2003)
135 
136  CALL timeset(routinen, handle)
137 
138  NULLIFY (logger, mark_states, mark_list, para_env)
139  logger => cp_get_default_logger()
140  !-----------------------------------------------------------------------------
141  ! 1.
142  !-----------------------------------------------------------------------------
143  CALL get_qs_env(qs_env, para_env=para_env)
144  nproc = para_env%num_pe
145  output_unit = cp_logger_get_default_io_unit(logger)
146  CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", &
147  explicit=explicit)
148  IF (explicit) THEN
149  CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", &
150  r_vals=eval_range)
151  ELSE
152  ALLOCATE (eval_range(2))
153  eval_range(1) = -huge(0.0_dp)
154  eval_range(2) = +huge(0.0_dp)
155  END IF
156  CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", &
157  n_rep_val=n_rep)
158  IF (n_rep .GT. 0) THEN
159  ALLOCATE (mark_states(2, n_rep))
160  IF (.NOT. ASSOCIATED(marked_states)) THEN
161  ALLOCATE (marked_states(n_rep))
162  END IF
163  DO i = 1, n_rep
164  CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", &
165  i_rep_val=i, i_vals=mark_list)
166  mark_states(:, i) = mark_list(:)
167  END DO
168  ELSE
169  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  extension=".data", middle_name="Molecular_DOS", log_filename=.false.)
178  IF (unit_report > 0) THEN
179  WRITE (unit_report, *) SIZE(mo_eigenvalues), " number of states "
180  DO i = 1, SIZE(mo_eigenvalues)
181  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  nrow_global=nrow_global)
192  CALL cp_fm_create(smo, mo_coeff%matrix_struct)
193  CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_coeff, smo, ncol_global)
194 
195  !-----------------------------------------------------------------------------
196  !-----------------------------------------------------------------------------
197  ! 4.
198  !-----------------------------------------------------------------------------
199  ALLOCATE (nstates(2))
200 
201  CALL cp_fm_create(storage, mo_localized%matrix_struct, name='storage')
202 
203  DO imol = 1, SIZE(molecule_set)
204  IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN
205  nstates(1) = molecule_set(imol)%lmi(ispin)%nstates
206  ELSE
207  nstates(1) = 0
208  END IF
209  nstates(2) = para_env%mepos
210 
211  CALL para_env%maxloc(nstates)
212 
213  IF (nstates(1) == 0) cycle
214  NULLIFY (states)
215  ALLOCATE (states(nstates(1)))
216  states(:) = 0
217 
218  iproc = nstates(2)
219  IF (iproc == para_env%mepos) THEN
220  states(:) = molecule_set(imol)%lmi(ispin)%states(:)
221  END IF
222  !!BCAST from here root = iproc
223  CALL para_env%bcast(states, iproc)
224 
225  ns = nstates(1)
226  ind => states(:)
227  ALLOCATE (evals(ns))
228 
229  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  context=mo_localized%matrix_struct%context)
235 
236  CALL cp_fm_create(b, fm_struct_tmp, name="b")
237  CALL cp_fm_create(c, fm_struct_tmp, name="c")
238  CALL cp_fm_create(rot_e_vectors, fm_struct_tmp, name="rot_e_vectors")
239 
240  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  context=mo_localized%matrix_struct%context)
245 
246  CALL cp_fm_create(d, fm_struct_tmp, name="d")
247  CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e_vectors")
248  CALL cp_fm_struct_release(fm_struct_tmp)
249 
250  DO i = 1, ns
251  CALL cp_fm_to_fm(mo_localized, b, 1, ind(i), i)
252  END DO
253 
254  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  b, c, 0.0_dp, d)
258 
259  CALL choose_eigv_solver(d, e_vectors, evals)
260 
261  IF (output_unit > 0) WRITE (output_unit, *) ""
262  IF (output_unit > 0) WRITE (output_unit, *) "MOLECULE ", imol
263  IF (output_unit > 0) WRITE (output_unit, *) "NUMBER OF STATES ", ns
264  IF (output_unit > 0) WRITE (output_unit, *) "EIGENVALUES"
265  IF (output_unit > 0) WRITE (output_unit, *) ""
266  IF (output_unit > 0) WRITE (output_unit, *) "ENERGY original MO-index"
267 
268  DO k = 1, ns
269  IF (ASSOCIATED(mark_states)) THEN
270  DO i = 1, n_rep
271  IF (imol == mark_states(1, i) .AND. k == mark_states(2, i)) marked_states(i) = ind(k)
272  END DO
273  END IF
274  IF (output_unit > 0) WRITE (output_unit, *) evals(k), ind(k)
275  END DO
276  IF (unit_report > 0) THEN
277  WRITE (unit_report, *) imol, ns, " imol, number of states"
278  DO k = 1, ns
279  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  b, e_vectors, 0.0_dp, rot_e_vectors)
285 
286  DO i = 1, ns
287  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  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  pw_env=pw_env)
321 
322  DO i = 1, ns
323  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  pw_env)
328 
329  WRITE (filename, '(a9,I4.4,a1,I5.5,a4)') "MOLECULE_", imol, "_", i, tag
330  WRITE (title, '(A,I0,A,I0,A,F14.6,A,I0)') "Mol. Eigenstate ", i, " of ", ns, " E [a.u.] = ", &
331  evals(i), " Orig. index ", ind(i)
332  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  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  "MOLECULAR_STATES%CUBES%STRIDE"), mpi_io=mpi_io)
339  CALL cp_print_key_finished_output(unit_nr, logger, loc_print_section, &
340  "MOLECULAR_STATES%CUBES", mpi_io=mpi_io)
341  END DO
342  END IF
343 
344  DEALLOCATE (evals)
345  CALL cp_fm_release(b)
346  CALL cp_fm_release(c)
347  CALL cp_fm_release(d)
348  CALL cp_fm_release(e_vectors)
349  CALL cp_fm_release(rot_e_vectors)
350 
351  DEALLOCATE (states)
352 
353  END DO
354  CALL cp_fm_release(smo)
355  CALL cp_fm_to_fm(storage, mo_localized)
356  CALL cp_fm_release(storage)
357  IF (ASSOCIATED(mark_states)) THEN
358  DEALLOCATE (mark_states)
359  END IF
360  DEALLOCATE (nstates)
361  CALL cp_print_key_finished_output(unit_report, logger, loc_print_section, &
362  "MOLECULAR_STATES")
363  !------------------------------------------------------------------------------
364 
365  IF (.NOT. explicit) THEN
366  DEALLOCATE (eval_range)
367  END IF
368 
369  CALL timestop(handle)
370 
371  END SUBROUTINE construct_molecular_states
372 
373 END MODULE molecular_states
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public hunt2003
Definition: bibliography.F:43
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition: cp_fm_types.F:643
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
Interface to the message passing library MPI.
Routines for the calculation of molecular states.
subroutine, public construct_molecular_states(molecule_set, mo_localized, mo_coeff, mo_eigenvalues, Hks, matrix_S, qs_env, wf_r, wf_g, loc_print_section, particles, tag, marked_states, ispin)
constructs molecular states. mo_localized gets overwritten!
Define the data structure for the molecule information.
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23