1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for the calculation of molecular states 8!> \author CJM 9! ************************************************************************************************** 10MODULE molecular_states 11 USE atomic_kind_types, ONLY: atomic_kind_type 12 USE bibliography, ONLY: Hunt2003,& 13 cite_reference 14 USE cell_types, ONLY: cell_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 17 USE cp_fm_diag, ONLY: choose_eigv_solver 18 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 19 cp_fm_struct_release,& 20 cp_fm_struct_type 21 USE cp_fm_types, ONLY: cp_fm_create,& 22 cp_fm_get_element,& 23 cp_fm_get_info,& 24 cp_fm_release,& 25 cp_fm_to_fm,& 26 cp_fm_type 27 USE cp_gemm_interface, ONLY: cp_gemm 28 USE cp_log_handling, ONLY: cp_get_default_logger,& 29 cp_logger_get_default_io_unit,& 30 cp_logger_type 31 USE cp_output_handling, ONLY: cp_p_file,& 32 cp_print_key_finished_output,& 33 cp_print_key_should_output,& 34 cp_print_key_unit_nr 35 USE cp_para_types, ONLY: cp_para_env_type 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_bcast,& 45 mp_maxloc 46 USE molecule_types, ONLY: molecule_type 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_p_type 51 USE qs_collocate_density, ONLY: calculate_wavefunction 52 USE qs_environment_types, ONLY: get_qs_env,& 53 qs_environment_type 54 USE qs_kind_types, ONLY: qs_kind_type 55#include "./base/base_uses.f90" 56 57 IMPLICIT NONE 58 59 PRIVATE 60 61! *** Global parameters *** 62 63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_states' 64 65 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 66 67! *** Public subroutines *** 68 69 PUBLIC :: construct_molecular_states 70 71CONTAINS 72 73! ************************************************************************************************** 74!> \brief constructs molecular states. mo_localized gets overwritten! 75!> \param molecule_set ... 76!> \param mo_localized ... 77!> \param mo_coeff ... 78!> \param mo_eigenvalues ... 79!> \param Hks ... 80!> \param matrix_S ... 81!> \param qs_env ... 82!> \param wf_r ... 83!> \param wf_g ... 84!> \param loc_print_section ... 85!> \param particles ... 86!> \param tag ... 87!> \param marked_states ... 88!> \param ispin ... 89! ************************************************************************************************** 90 SUBROUTINE construct_molecular_states(molecule_set, mo_localized, & 91 mo_coeff, mo_eigenvalues, Hks, matrix_S, qs_env, wf_r, wf_g, & 92 loc_print_section, particles, tag, marked_states, ispin) 93 94 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 95 TYPE(cp_fm_type), POINTER :: mo_localized, mo_coeff 96 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 97 TYPE(dbcsr_type), POINTER :: Hks, matrix_S 98 TYPE(qs_environment_type), POINTER :: qs_env 99 TYPE(pw_p_type), INTENT(INOUT) :: wf_r, wf_g 100 TYPE(section_vals_type), POINTER :: loc_print_section 101 TYPE(particle_list_type), POINTER :: particles 102 CHARACTER(LEN=*), INTENT(IN) :: tag 103 INTEGER, DIMENSION(:), POINTER :: marked_states 104 INTEGER, INTENT(IN) :: ispin 105 106 CHARACTER(len=*), PARAMETER :: routineN = 'construct_molecular_states', & 107 routineP = moduleN//':'//routineN 108 109 CHARACTER(LEN=default_path_length) :: filename 110 CHARACTER(LEN=default_string_length) :: title 111 INTEGER :: handle, i, imol, iproc, k, n_rep, & 112 ncol_global, nproc, nrow_global, ns, & 113 output_unit, unit_nr, unit_report 114 INTEGER, DIMENSION(:), POINTER :: ind, mark_list 115 INTEGER, DIMENSION(:, :), POINTER :: mark_states 116 INTEGER, POINTER :: nstates(:), states(:) 117 LOGICAL :: explicit, mpi_io 118 REAL(KIND=dp) :: tmp 119 REAL(KIND=dp), ALLOCATABLE :: evals(:) 120 REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range 121 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 122 TYPE(cell_type), POINTER :: cell 123 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 124 TYPE(cp_fm_type), POINTER :: b, c, d, D_igk, e_vectors, & 125 rot_e_vectors, smo, storage 126 TYPE(cp_logger_type), POINTER :: logger 127 TYPE(cp_para_env_type), POINTER :: para_env 128 TYPE(dft_control_type), POINTER :: dft_control 129 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 130 TYPE(pw_env_type), POINTER :: pw_env 131 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 132 133 CALL cite_reference(Hunt2003) 134 135 CALL timeset(routineN, handle) 136 137 NULLIFY (logger, mark_states, mark_list, para_env) 138 logger => cp_get_default_logger() 139 !----------------------------------------------------------------------------- 140 ! 1. 141 !----------------------------------------------------------------------------- 142 CALL get_qs_env(qs_env, para_env=para_env) 143 nproc = para_env%num_pe 144 output_unit = cp_logger_get_default_io_unit(logger) 145 CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", & 146 explicit=explicit) 147 IF (explicit) THEN 148 CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%CUBE_EVAL_RANGE", & 149 r_vals=eval_range) 150 ELSE 151 ALLOCATE (eval_range(2)) 152 eval_range(1) = -HUGE(0.0_dp) 153 eval_range(2) = +HUGE(0.0_dp) 154 ENDIF 155 CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", & 156 n_rep_val=n_rep) 157 IF (n_rep .GT. 0) THEN 158 ALLOCATE (mark_states(2, n_rep)) 159 IF (.NOT. ASSOCIATED(marked_states)) THEN 160 ALLOCATE (marked_states(n_rep)) 161 END IF 162 DO i = 1, n_rep 163 CALL section_vals_val_get(loc_print_section, "MOLECULAR_STATES%MARK_STATES", & 164 i_rep_val=i, i_vals=mark_list) 165 mark_states(:, i) = mark_list(:) 166 END DO 167 ELSE 168 NULLIFY (marked_states) 169 END IF 170 171 !----------------------------------------------------------------------------- 172 !----------------------------------------------------------------------------- 173 ! 2. 174 !----------------------------------------------------------------------------- 175 unit_report = cp_print_key_unit_nr(logger, loc_print_section, "MOLECULAR_STATES", & 176 extension=".data", middle_name="Molecular_DOS", log_filename=.FALSE.) 177 IF (unit_report > 0) THEN 178 WRITE (unit_report, *) SIZE(mo_eigenvalues), " number of states " 179 DO i = 1, SIZE(mo_eigenvalues) 180 WRITE (unit_report, *) mo_eigenvalues(i) 181 ENDDO 182 ENDIF 183 184 !----------------------------------------------------------------------------- 185 !----------------------------------------------------------------------------- 186 ! 3. 187 !----------------------------------------------------------------------------- 188 CALL cp_fm_get_info(mo_localized, & 189 ncol_global=ncol_global, & 190 nrow_global=nrow_global) 191 NULLIFY (smo) 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 ENDIF 209 nstates(2) = para_env%mepos 210 211 CALL mp_maxloc(nstates, para_env%group) 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 mp_bcast(states, iproc, para_env%group) 224 225 ns = nstates(1) 226 ind => states(:) 227 ALLOCATE (evals(ns)) 228 229 NULLIFY (b, c, d, fm_struct_tmp, e_vectors, rot_e_vectors) 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 cp_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 ENDDO 281 ENDIF 282 283 CALL cp_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 cp_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 ENDIF 306 ENDDO 307 ENDDO 308 CALL cp_fm_release(D_igk) 309 ENDIF 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%pw, 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 ENDIF 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 373END MODULE molecular_states 374