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 calculates the electron transfer coupling elements by projection-operator approach
10 : !> Kondov et al. J.Phys.Chem.C 2007, 111, 11970-11981
11 : !> \author Z. Futera (02.2017)
12 : ! **************************************************************************************************
13 : MODULE et_coupling_proj
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE bibliography, ONLY: Futera2017,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type
22 : USE cp_blacs_env, ONLY: cp_blacs_env_type
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : cp_dbcsr_sm_fm_multiply
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
28 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
29 : cp_fm_power
30 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
31 : cp_fm_struct_equivalent,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE cp_fm_types, ONLY: &
35 : cp_fm_create, cp_fm_get_element, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
36 : cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_vectorssum
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_type,&
39 : cp_to_string
40 : USE cp_output_handling, ONLY: cp_p_file,&
41 : cp_print_key_finished_output,&
42 : cp_print_key_should_output,&
43 : cp_print_key_unit_nr
44 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
45 : USE input_section_types, ONLY: section_get_ivals,&
46 : section_get_lval,&
47 : section_vals_get,&
48 : section_vals_get_subs_vals,&
49 : section_vals_type,&
50 : section_vals_val_get
51 : USE kinds, ONLY: default_path_length,&
52 : default_string_length,&
53 : dp
54 : USE kpoint_types, ONLY: kpoint_type
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE orbital_pointers, ONLY: nso
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE particle_list_types, ONLY: particle_list_type
59 : USE particle_types, ONLY: particle_type
60 : USE physcon, ONLY: evolt
61 : USE pw_env_types, ONLY: pw_env_get,&
62 : pw_env_type
63 : USE pw_pool_types, ONLY: pw_pool_p_type,&
64 : pw_pool_type
65 : USE pw_types, ONLY: pw_c1d_gs_type,&
66 : pw_r3d_rs_type
67 : USE qs_collocate_density, ONLY: calculate_wavefunction
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
70 : USE qs_kind_types, ONLY: get_qs_kind,&
71 : get_qs_kind_set,&
72 : qs_kind_type
73 : USE qs_mo_methods, ONLY: make_mo_eig
74 : USE qs_mo_occupation, ONLY: set_mo_occupation
75 : USE qs_mo_types, ONLY: allocate_mo_set,&
76 : deallocate_mo_set,&
77 : mo_set_type
78 : USE qs_subsys_types, ONLY: qs_subsys_get,&
79 : qs_subsys_type
80 : USE scf_control_types, ONLY: scf_control_type
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'et_coupling_proj'
88 :
89 : ! Electronic-coupling calculation data structure
90 : !
91 : ! n_atoms - number of atoms in the blocks
92 : ! n_blocks - number of atomic blocks (donor,acceptor,bridge,...)
93 : ! fermi - system Fermi level (alpha/beta spin component)
94 : ! m_transf - transformation matrix for basis-set orthogonalization (S^{-1/2})
95 : ! m_transf_inv - inversion transformation matrix
96 : ! block - atomic data blocks
97 : TYPE et_cpl
98 : INTEGER :: n_atoms = 0
99 : INTEGER :: n_blocks = 0
100 : REAL(KIND=dp), DIMENSION(:), POINTER :: fermi => NULL()
101 : TYPE(cp_fm_type), POINTER :: m_transf => NULL()
102 : TYPE(cp_fm_type), POINTER :: m_transf_inv => NULL()
103 : TYPE(et_cpl_block), DIMENSION(:), POINTER :: block => NULL()
104 : END TYPE et_cpl
105 :
106 : ! Electronic-coupling data block
107 : !
108 : ! n_atoms - number of atoms
109 : ! n_electrons - number of electrons
110 : ! n_ao - number of AO basis functions
111 : ! atom - list of atoms
112 : ! mo - electronic states
113 : ! hab - electronic-coupling elements
114 : TYPE et_cpl_block
115 : INTEGER :: n_atoms = 0
116 : INTEGER :: n_electrons = 0
117 : INTEGER :: n_ao = 0
118 : TYPE(et_cpl_atom), DIMENSION(:), POINTER :: atom => NULL()
119 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo => NULL()
120 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: hab => NULL()
121 : END TYPE et_cpl_block
122 :
123 : ! Electronic-coupling block-atom data
124 : ! id - atom ID
125 : ! n_ao - number of AO basis functions
126 : ! ao_pos - position of atom in array of AO functions
127 : TYPE et_cpl_atom
128 : INTEGER :: id = 0
129 : INTEGER :: n_ao = 0
130 : INTEGER :: ao_pos = 0
131 : END TYPE et_cpl_atom
132 :
133 : PUBLIC :: calc_et_coupling_proj
134 :
135 : CONTAINS
136 :
137 : ! **************************************************************************************************
138 : !> \brief Release memory allocate for electronic coupling data structures
139 : !> \param ec electronic coupling data structure
140 : !> \author Z. Futera (02.2017)
141 : ! **************************************************************************************************
142 10 : SUBROUTINE release_ec_data(ec)
143 :
144 : ! Routine arguments
145 : TYPE(et_cpl), POINTER :: ec
146 :
147 : INTEGER :: i, j
148 :
149 : ! Routine name for debug purposes
150 :
151 10 : IF (ASSOCIATED(ec)) THEN
152 :
153 10 : IF (ASSOCIATED(ec%fermi)) &
154 10 : DEALLOCATE (ec%fermi)
155 10 : IF (ASSOCIATED(ec%m_transf)) THEN
156 10 : CALL cp_fm_release(matrix=ec%m_transf)
157 10 : DEALLOCATE (ec%m_transf)
158 10 : NULLIFY (ec%m_transf)
159 : END IF
160 10 : IF (ASSOCIATED(ec%m_transf_inv)) THEN
161 10 : CALL cp_fm_release(matrix=ec%m_transf_inv)
162 10 : DEALLOCATE (ec%m_transf_inv)
163 10 : NULLIFY (ec%m_transf_inv)
164 : END IF
165 :
166 10 : IF (ASSOCIATED(ec%block)) THEN
167 :
168 30 : DO i = 1, SIZE(ec%block)
169 20 : IF (ASSOCIATED(ec%block(i)%atom)) &
170 20 : DEALLOCATE (ec%block(i)%atom)
171 20 : IF (ASSOCIATED(ec%block(i)%mo)) THEN
172 60 : DO j = 1, SIZE(ec%block(i)%mo)
173 60 : CALL deallocate_mo_set(ec%block(i)%mo(j))
174 : END DO
175 20 : DEALLOCATE (ec%block(i)%mo)
176 : END IF
177 30 : CALL cp_fm_release(ec%block(i)%hab)
178 : END DO
179 :
180 10 : DEALLOCATE (ec%block)
181 :
182 : END IF
183 :
184 10 : DEALLOCATE (ec)
185 :
186 : END IF
187 :
188 10 : END SUBROUTINE release_ec_data
189 :
190 : ! **************************************************************************************************
191 : !> \brief check the electronic-coupling input section and set the atomic block data
192 : !> \param qs_env QuickStep environment containing all system data
193 : !> \param et_proj_sec the electronic-coupling input section
194 : !> \param ec electronic coupling data structure
195 : !> \author Z. Futera (02.2017)
196 : ! **************************************************************************************************
197 10 : SUBROUTINE set_block_data(qs_env, et_proj_sec, ec)
198 :
199 : ! Routine arguments
200 : TYPE(qs_environment_type), POINTER :: qs_env
201 : TYPE(section_vals_type), POINTER :: et_proj_sec
202 : TYPE(et_cpl), POINTER :: ec
203 :
204 : INTEGER :: i, j, k, l, n, n_ao, n_atoms, n_set
205 10 : INTEGER, DIMENSION(:), POINTER :: atom_id, atom_nf, atom_ps, n_shell, t
206 10 : INTEGER, DIMENSION(:, :), POINTER :: ang_mom_id
207 : LOGICAL :: found
208 : TYPE(gto_basis_set_type), POINTER :: ao_basis_set
209 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
211 : TYPE(section_vals_type), POINTER :: block_sec
212 :
213 : ! Routine name for debug purposes
214 :
215 10 : NULLIFY (ao_basis_set)
216 10 : NULLIFY (particle_set)
217 10 : NULLIFY (qs_kind_set)
218 10 : NULLIFY (n_shell)
219 10 : NULLIFY (ang_mom_id)
220 10 : NULLIFY (atom_nf)
221 10 : NULLIFY (atom_id)
222 10 : NULLIFY (block_sec)
223 :
224 : ! Initialization
225 10 : ec%n_atoms = 0
226 10 : ec%n_blocks = 0
227 10 : NULLIFY (ec%fermi)
228 10 : NULLIFY (ec%m_transf)
229 10 : NULLIFY (ec%m_transf_inv)
230 10 : NULLIFY (ec%block)
231 :
232 : ! Number of atoms / atom types
233 10 : CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, natom=n_atoms)
234 : ! Number of AO basis functions
235 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
236 :
237 : ! Number of AO functions per atom
238 30 : ALLOCATE (atom_nf(n_atoms))
239 10 : CPASSERT(ASSOCIATED(atom_nf))
240 :
241 82 : atom_nf = 0
242 82 : DO i = 1, n_atoms
243 72 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=j)
244 72 : CALL get_qs_kind(qs_kind_set(j), basis_set=ao_basis_set)
245 72 : IF (.NOT. ASSOCIATED(ao_basis_set)) &
246 0 : CPABORT('Unsupported basis set type. ')
247 : CALL get_gto_basis_set(gto_basis_set=ao_basis_set, &
248 72 : nset=n_set, nshell=n_shell, l=ang_mom_id)
249 238 : DO j = 1, n_set
250 280 : DO k = 1, n_shell(j)
251 208 : atom_nf(i) = atom_nf(i) + nso(ang_mom_id(k, j))
252 : END DO
253 : END DO
254 : END DO
255 :
256 : ! Sanity check
257 10 : n = 0
258 82 : DO i = 1, n_atoms
259 82 : n = n + atom_nf(i)
260 : END DO
261 10 : CPASSERT(n == n_ao)
262 :
263 : ! Atom position in AO array
264 30 : ALLOCATE (atom_ps(n_atoms))
265 10 : CPASSERT(ASSOCIATED(atom_ps))
266 82 : atom_ps = 1
267 72 : DO i = 1, n_atoms - 1
268 72 : atom_ps(i + 1) = atom_ps(i) + atom_nf(i)
269 : END DO
270 :
271 : ! Number of blocks
272 10 : block_sec => section_vals_get_subs_vals(et_proj_sec, 'BLOCK')
273 10 : CALL section_vals_get(block_sec, n_repetition=ec%n_blocks)
274 50 : ALLOCATE (ec%block(ec%n_blocks))
275 10 : CPASSERT(ASSOCIATED(ec%block))
276 :
277 : ! Block data
278 30 : ALLOCATE (t(n_atoms))
279 10 : CPASSERT(ASSOCIATED(t))
280 :
281 10 : ec%n_atoms = 0
282 30 : DO i = 1, ec%n_blocks
283 :
284 : ! Initialization
285 20 : ec%block(i)%n_atoms = 0
286 20 : ec%block(i)%n_electrons = 0
287 20 : ec%block(i)%n_ao = 0
288 20 : NULLIFY (ec%block(i)%atom)
289 20 : NULLIFY (ec%block(i)%mo)
290 20 : NULLIFY (ec%block(i)%hab)
291 :
292 : ! Number of electrons
293 : CALL section_vals_val_get(block_sec, i_rep_section=i, &
294 20 : keyword_name='NELECTRON', i_val=ec%block(i)%n_electrons)
295 :
296 : ! User-defined atom array
297 : CALL section_vals_val_get(block_sec, i_rep_section=i, &
298 20 : keyword_name='ATOMS', i_vals=atom_id)
299 :
300 : ! Count unique atoms
301 92 : DO j = 1, SIZE(atom_id)
302 : ! Check atom ID validity
303 72 : IF (atom_id(j) < 1 .OR. atom_id(j) > n_atoms) &
304 0 : CPABORT('invalid fragment atom ID ('//TRIM(ADJUSTL(cp_to_string(atom_id(j))))//')')
305 : ! Check if the atom is not in previously-defined blocks
306 : found = .FALSE.
307 108 : DO k = 1, i - 1
308 348 : DO l = 1, ec%block(k)%n_atoms
309 276 : IF (ec%block(k)%atom(l)%id == atom_id(j)) THEN
310 0 : CPWARN('multiple definition of atom'//TRIM(ADJUSTL(cp_to_string(atom_id(j)))))
311 0 : found = .TRUE.
312 0 : EXIT
313 : END IF
314 : END DO
315 : END DO
316 : ! Check if the atom is not in already defined in the present block
317 72 : IF (.NOT. found) THEN
318 276 : DO k = 1, ec%block(i)%n_atoms
319 276 : IF (t(k) == atom_id(j)) THEN
320 0 : CPWARN('multiple definition of atom'//TRIM(ADJUSTL(cp_to_string(atom_id(j)))))
321 : found = .TRUE.
322 : EXIT
323 : END IF
324 : END DO
325 : END IF
326 : ! Save the atom
327 20 : IF (.NOT. found) THEN
328 72 : ec%block(i)%n_atoms = ec%block(i)%n_atoms + 1
329 72 : t(ec%block(i)%n_atoms) = atom_id(j)
330 : END IF
331 : END DO
332 :
333 : ! Memory allocation
334 132 : ALLOCATE (ec%block(i)%atom(ec%block(i)%n_atoms))
335 20 : CPASSERT(ASSOCIATED(ec%block(i)%atom))
336 :
337 : ! Save atom IDs and number of AOs
338 92 : DO j = 1, ec%block(i)%n_atoms
339 72 : ec%block(i)%atom(j)%id = t(j)
340 72 : ec%block(i)%atom(j)%n_ao = atom_nf(ec%block(i)%atom(j)%id)
341 72 : ec%block(i)%atom(j)%ao_pos = atom_ps(ec%block(i)%atom(j)%id)
342 92 : ec%block(i)%n_ao = ec%block(i)%n_ao + ec%block(i)%atom(j)%n_ao
343 : END DO
344 :
345 30 : ec%n_atoms = ec%n_atoms + ec%block(i)%n_atoms
346 : END DO
347 :
348 : ! Clean memory
349 10 : IF (ASSOCIATED(atom_nf)) &
350 10 : DEALLOCATE (atom_nf)
351 10 : IF (ASSOCIATED(atom_ps)) &
352 10 : DEALLOCATE (atom_ps)
353 10 : IF (ASSOCIATED(t)) &
354 10 : DEALLOCATE (t)
355 :
356 10 : END SUBROUTINE set_block_data
357 :
358 : ! **************************************************************************************************
359 : !> \brief check the electronic-coupling input section and set the atomic block data
360 : !> \param ec electronic coupling data structure
361 : !> \param fa system Fermi level (alpha spin)
362 : !> \param fb system Fermi level (beta spin)
363 : !> \author Z. Futera (02.2017)
364 : ! **************************************************************************************************
365 10 : SUBROUTINE set_fermi(ec, fa, fb)
366 :
367 : ! Routine arguments
368 : TYPE(et_cpl), POINTER :: ec
369 : REAL(KIND=dp) :: fa
370 : REAL(KIND=dp), OPTIONAL :: fb
371 :
372 : ! Routine name for debug purposes
373 :
374 10 : NULLIFY (ec%fermi)
375 :
376 10 : IF (PRESENT(fb)) THEN
377 :
378 10 : ALLOCATE (ec%fermi(2))
379 10 : CPASSERT(ASSOCIATED(ec%fermi))
380 10 : ec%fermi(1) = fa
381 10 : ec%fermi(2) = fb
382 :
383 : ELSE
384 :
385 0 : ALLOCATE (ec%fermi(1))
386 0 : CPASSERT(ASSOCIATED(ec%fermi))
387 0 : ec%fermi(1) = fa
388 :
389 : END IF
390 :
391 10 : END SUBROUTINE set_fermi
392 :
393 : ! **************************************************************************************************
394 : !> \brief reorder Hamiltonian matrix according to defined atomic blocks
395 : !> \param ec electronic coupling data structure
396 : !> \param mat_h the Hamiltonian matrix
397 : !> \param mat_w working matrix of the same dimension
398 : !> \author Z. Futera (02.2017)
399 : ! **************************************************************************************************
400 20 : SUBROUTINE reorder_hamiltonian_matrix(ec, mat_h, mat_w)
401 :
402 : ! Routine arguments
403 : TYPE(et_cpl), POINTER :: ec
404 : TYPE(cp_fm_type), INTENT(IN) :: mat_h, mat_w
405 :
406 : INTEGER :: ic, ir, jc, jr, kc, kr, mc, mr, nc, nr
407 : REAL(KIND=dp) :: xh
408 :
409 : ! Routine name for debug purposes
410 : ! Local variables
411 :
412 20 : IF (.NOT. cp_fm_struct_equivalent(mat_h%matrix_struct, mat_w%matrix_struct)) &
413 0 : CPABORT('cannot reorder Hamiltonian, working-matrix structure is not equivalent')
414 :
415 : ! Matrix-element reordering
416 20 : nr = 1
417 : ! Rows
418 60 : DO ir = 1, ec%n_blocks
419 204 : DO jr = 1, ec%block(ir)%n_atoms
420 592 : DO kr = 1, ec%block(ir)%atom(jr)%n_ao
421 : ! Columns
422 408 : nc = 1
423 1224 : DO ic = 1, ec%n_blocks
424 6072 : DO jc = 1, ec%block(ic)%n_atoms
425 18384 : DO kc = 1, ec%block(ic)%atom(jc)%n_ao
426 12720 : mr = ec%block(ir)%atom(jr)%ao_pos + kr - 1
427 12720 : mc = ec%block(ic)%atom(jc)%ao_pos + kc - 1
428 12720 : CALL cp_fm_get_element(mat_h, nr, nc, xh)
429 12720 : CALL cp_fm_set_element(mat_w, nr, nc, xh)
430 30288 : nc = nc + 1
431 : END DO
432 : END DO
433 : END DO
434 552 : nr = nr + 1
435 : END DO
436 : END DO
437 : END DO
438 :
439 : ! Copy the reordered matrix to original data array
440 20 : CALL cp_fm_to_fm(mat_w, mat_h)
441 :
442 20 : END SUBROUTINE reorder_hamiltonian_matrix
443 :
444 : ! **************************************************************************************************
445 : !> \brief calculated transformation matrix for basis-set orthogonalization (S^{-1/2})
446 : !> \param qs_env QuickStep environment containing all system data
447 : !> \param mat_t storage for the transformation matrix
448 : !> \param mat_i storage for the inversion transformation matrix
449 : !> \param mat_w working matrix of the same dimension
450 : !> \author Z. Futera (02.2017)
451 : ! **************************************************************************************************
452 10 : SUBROUTINE get_s_half_inv_matrix(qs_env, mat_t, mat_i, mat_w)
453 :
454 : ! Routine arguments
455 : TYPE(qs_environment_type), POINTER :: qs_env
456 : TYPE(cp_fm_type), INTENT(IN) :: mat_t, mat_i, mat_w
457 :
458 : INTEGER :: n_deps
459 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
460 : TYPE(scf_control_type), POINTER :: scf_cntrl
461 :
462 : ! Routine name for debug purposes
463 :
464 10 : NULLIFY (mat_s)
465 10 : NULLIFY (scf_cntrl)
466 :
467 10 : CALL get_qs_env(qs_env, matrix_s=mat_s)
468 10 : CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_t)
469 10 : CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_i)
470 :
471 : ! Transformation S -> S^{-1/2}
472 10 : CALL get_qs_env(qs_env, scf_control=scf_cntrl)
473 10 : CALL cp_fm_power(mat_t, mat_w, -0.5_dp, scf_cntrl%eps_eigval, n_deps)
474 10 : CALL cp_fm_power(mat_i, mat_w, +0.5_dp, scf_cntrl%eps_eigval, n_deps)
475 : ! Sanity check
476 10 : IF (n_deps /= 0) THEN
477 : CALL cp_warn(__LOCATION__, &
478 : "Overlap matrix exhibits linear dependencies. At least some "// &
479 0 : "eigenvalues have been quenched.")
480 : END IF
481 :
482 10 : END SUBROUTINE get_s_half_inv_matrix
483 :
484 : ! **************************************************************************************************
485 : !> \brief transform KS hamiltonian to orthogonalized block-separated basis set
486 : !> \param qs_env QuickStep environment containing all system data
487 : !> \param ec electronic coupling data structure
488 : !> \param fm_s full-matrix structure used for allocation of KS matrices
489 : !> \param mat_t storage for pointers to the transformed KS matrices
490 : !> \param mat_w working matrix of the same dimension
491 : !> \param n_ao total number of AO basis functions
492 : !> \param n_spins number of spin components
493 : !> \author Z. Futera (02.2017)
494 : ! **************************************************************************************************
495 10 : SUBROUTINE get_block_hamiltonian(qs_env, ec, fm_s, mat_t, mat_w, n_ao, n_spins)
496 :
497 : ! Routine arguments
498 : TYPE(qs_environment_type), POINTER :: qs_env
499 : TYPE(et_cpl), POINTER :: ec
500 : TYPE(cp_fm_struct_type), POINTER :: fm_s
501 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
502 : INTENT(OUT) :: mat_t
503 : TYPE(cp_fm_type), INTENT(IN) :: mat_w
504 : INTEGER :: n_ao, n_spins
505 :
506 : INTEGER :: i
507 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_h
508 :
509 : ! Routine name for debug purposes
510 :
511 10 : NULLIFY (mat_h)
512 :
513 : ! Memory allocation
514 50 : ALLOCATE (mat_t(n_spins))
515 :
516 : ! KS Hamiltonian
517 10 : CALL get_qs_env(qs_env, matrix_ks=mat_h)
518 : ! Transformation matrix
519 10 : ALLOCATE (ec%m_transf, ec%m_transf_inv)
520 : CALL cp_fm_create(matrix=ec%m_transf, matrix_struct=fm_s, &
521 10 : name='S^(-1/2) TRANSFORMATION MATRIX')
522 : CALL cp_fm_create(matrix=ec%m_transf_inv, matrix_struct=fm_s, &
523 10 : name='S^(+1/2) TRANSFORMATION MATRIX')
524 10 : CALL get_s_half_inv_matrix(qs_env, ec%m_transf, ec%m_transf_inv, mat_w)
525 :
526 30 : DO i = 1, n_spins
527 :
528 : ! Full-matrix format
529 : CALL cp_fm_create(matrix=mat_t(i), matrix_struct=fm_s, &
530 20 : name='KS HAMILTONIAN IN SEPARATED ORTHOGONALIZED BASIS SET')
531 20 : CALL copy_dbcsr_to_fm(mat_h(i)%matrix, mat_t(i))
532 :
533 : ! Transform KS Hamiltonian to the orthogonalized AO basis set
534 20 : CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, 1.0_dp, ec%m_transf, mat_t(i), 0.0_dp, mat_w)
535 20 : CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, 1.0_dp, mat_w, ec%m_transf, 0.0_dp, mat_t(i))
536 :
537 : ! Reorder KS Hamiltonain elements to defined block structure
538 30 : CALL reorder_hamiltonian_matrix(ec, mat_t(i), mat_w)
539 :
540 : END DO
541 :
542 10 : END SUBROUTINE get_block_hamiltonian
543 :
544 : ! **************************************************************************************************
545 : !> \brief Diagonalize diagonal blocks of the KS hamiltonian in separated orthogonalized basis set
546 : !> \param qs_env QuickStep environment containing all system data
547 : !> \param ec electronic coupling data structure
548 : !> \param mat_h Hamiltonian in separated orthogonalized basis set
549 : !> \author Z. Futera (02.2017)
550 : ! **************************************************************************************************
551 10 : SUBROUTINE hamiltonian_block_diag(qs_env, ec, mat_h)
552 :
553 : ! Routine arguments
554 : TYPE(qs_environment_type), POINTER :: qs_env
555 : TYPE(et_cpl), POINTER :: ec
556 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mat_h
557 :
558 : INTEGER :: i, j, k, l, n_spins, spin
559 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
560 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
561 : TYPE(cp_fm_struct_type), POINTER :: fm_s
562 : TYPE(cp_fm_type) :: mat_u
563 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: dat
564 : TYPE(mp_para_env_type), POINTER :: para_env
565 :
566 : ! Routine name for debug purposes
567 :
568 10 : NULLIFY (vec_e)
569 10 : NULLIFY (blacs_env)
570 10 : NULLIFY (para_env)
571 10 : NULLIFY (fm_s)
572 :
573 : ! Parallel environment
574 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
575 :
576 : ! Storage for block sub-matrices
577 50 : ALLOCATE (dat(ec%n_blocks))
578 10 : CPASSERT(ALLOCATED(dat))
579 :
580 : ! Storage for electronic states and couplings
581 10 : n_spins = SIZE(mat_h)
582 30 : DO i = 1, ec%n_blocks
583 100 : ALLOCATE (ec%block(i)%mo(n_spins))
584 20 : CPASSERT(ASSOCIATED(ec%block(i)%mo))
585 200 : ALLOCATE (ec%block(i)%hab(n_spins, ec%n_blocks))
586 30 : CPASSERT(ASSOCIATED(ec%block(i)%hab))
587 : END DO
588 :
589 : ! Spin components
590 30 : DO spin = 1, n_spins
591 :
592 : ! Diagonal blocks
593 20 : j = 1
594 60 : DO i = 1, ec%n_blocks
595 :
596 : ! Memory allocation
597 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
598 40 : nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(i)%n_ao)
599 : CALL cp_fm_create(matrix=dat(i), matrix_struct=fm_s, &
600 40 : name='H_KS DIAGONAL BLOCK')
601 :
602 120 : ALLOCATE (vec_e(ec%block(i)%n_ao))
603 40 : CPASSERT(ASSOCIATED(vec_e))
604 :
605 : ! Copy block data
606 : CALL cp_fm_to_fm_submat(mat_h(spin), &
607 : dat(i), ec%block(i)%n_ao, &
608 40 : ec%block(i)%n_ao, j, j, 1, 1)
609 :
610 : ! Diagonalization
611 40 : CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='UNITARY MATRIX')
612 40 : CALL choose_eigv_solver(dat(i), mat_u, vec_e)
613 40 : CALL cp_fm_to_fm(mat_u, dat(i))
614 :
615 : ! Save state energies / vectors
616 40 : CALL create_block_mo_set(qs_env, ec, i, spin, mat_u, vec_e)
617 :
618 : ! Clean memory
619 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
620 40 : CALL cp_fm_release(matrix=mat_u)
621 40 : DEALLOCATE (vec_e)
622 :
623 : ! Off-set for next block
624 100 : j = j + ec%block(i)%n_ao
625 :
626 : END DO
627 :
628 : ! Off-diagonal blocks
629 20 : k = 1
630 60 : DO i = 1, ec%n_blocks
631 40 : l = 1
632 120 : DO j = 1, ec%n_blocks
633 80 : IF (i /= j) THEN
634 :
635 : ! Memory allocation
636 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
637 40 : nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(j)%n_ao)
638 : CALL cp_fm_create(matrix=ec%block(i)%hab(spin, j), matrix_struct=fm_s, &
639 40 : name='H_KS OFF-DIAGONAL BLOCK')
640 :
641 : ! Copy block data
642 : CALL cp_fm_to_fm_submat(mat_h(spin), &
643 : ec%block(i)%hab(spin, j), ec%block(i)%n_ao, &
644 40 : ec%block(j)%n_ao, k, l, 1, 1)
645 :
646 : ! Transformation
647 40 : CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='FULL WORK MATRIX')
648 : CALL parallel_gemm("T", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(i)%n_ao, &
649 40 : 1.0_dp, dat(i), ec%block(i)%hab(spin, j), 0.0_dp, mat_u)
650 : CALL parallel_gemm("N", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(j)%n_ao, &
651 40 : 1.0_dp, mat_u, dat(j), 0.0_dp, ec%block(i)%hab(spin, j))
652 :
653 : ! Clean memory
654 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
655 40 : CALL cp_fm_release(matrix=mat_u)
656 :
657 : END IF
658 : ! Off-set for next block
659 120 : l = l + ec%block(j)%n_ao
660 : END DO
661 : ! Off-set for next block
662 60 : k = k + ec%block(i)%n_ao
663 : END DO
664 :
665 : ! Clean memory
666 30 : IF (ALLOCATED(dat)) THEN
667 60 : DO i = 1, SIZE(dat)
668 60 : CALL cp_fm_release(dat(i))
669 : END DO
670 : END IF
671 : END DO
672 :
673 : ! Clean memory
674 10 : IF (ALLOCATED(dat)) &
675 10 : DEALLOCATE (dat)
676 :
677 20 : END SUBROUTINE hamiltonian_block_diag
678 :
679 : ! **************************************************************************************************
680 : !> \brief Return sum of selected squared MO coefficients
681 : !> \param blk_at list of atoms in the block
682 : !> \param mo array of MO sets
683 : !> \param id state index
684 : !> \param atom list of atoms for MO coefficient summing
685 : !> \return ...
686 : !> \author Z. Futera (02.2017)
687 : ! **************************************************************************************************
688 0 : FUNCTION get_mo_c2_sum(blk_at, mo, id, atom) RESULT(c2)
689 :
690 : ! Routine arguments
691 : TYPE(et_cpl_atom), DIMENSION(:), POINTER :: blk_at
692 : TYPE(cp_fm_type), INTENT(IN) :: mo
693 : INTEGER, INTENT(IN) :: id
694 : INTEGER, DIMENSION(:), POINTER :: atom
695 : REAL(KIND=dp) :: c2
696 :
697 : INTEGER :: i, ir, j, k
698 : LOGICAL :: found
699 : REAL(KIND=dp) :: c
700 :
701 : ! Returning value
702 : ! Routine name for debug purposes
703 : ! Local variables
704 :
705 : ! initialization
706 0 : c2 = 0.0d0
707 :
708 : ! selected atoms
709 0 : DO i = 1, SIZE(atom)
710 :
711 : ! find atomic function offset
712 0 : found = .FALSE.
713 0 : DO j = 1, SIZE(blk_at)
714 0 : IF (blk_at(j)%id == atom(i)) THEN
715 : found = .TRUE.
716 : EXIT
717 : END IF
718 : END DO
719 :
720 0 : IF (.NOT. found) &
721 0 : CPABORT('MO-fraction atom ID not defined in the block')
722 :
723 : ! sum MO coefficients from the atom
724 0 : DO k = 1, blk_at(j)%n_ao
725 0 : ir = blk_at(j)%ao_pos + k - 1
726 0 : CALL cp_fm_get_element(mo, ir, id, c)
727 0 : c2 = c2 + c*c
728 : END DO
729 :
730 : END DO
731 :
732 0 : END FUNCTION get_mo_c2_sum
733 :
734 : ! **************************************************************************************************
735 : !> \brief Print out specific MO coefficients
736 : !> \param output_unit unit number of the open output stream
737 : !> \param qs_env QuickStep environment containing all system data
738 : !> \param ec electronic coupling data structure
739 : !> \param blk atomic-block ID
740 : !> \param n_spins number of spin components
741 : !> \author Z. Futera (02.2017)
742 : ! **************************************************************************************************
743 20 : SUBROUTINE print_mo_coeff(output_unit, qs_env, ec, blk, n_spins)
744 :
745 : ! Routine arguments
746 : INTEGER, INTENT(IN) :: output_unit
747 : TYPE(qs_environment_type), POINTER :: qs_env
748 : TYPE(et_cpl), POINTER :: ec
749 : INTEGER, INTENT(IN) :: blk, n_spins
750 :
751 : INTEGER :: j, k, l, m, n, n_ao, n_mo
752 20 : INTEGER, DIMENSION(:), POINTER :: list_at, list_mo
753 : REAL(KIND=dp) :: c1, c2
754 20 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_w
755 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
756 : TYPE(section_vals_type), POINTER :: block_sec, print_sec
757 :
758 : ! Routine name for debug purposes
759 :
760 20 : NULLIFY (block_sec)
761 20 : NULLIFY (print_sec)
762 20 : NULLIFY (qs_kind_set)
763 :
764 : ! Atomic block data
765 : block_sec => section_vals_get_subs_vals(qs_env%input, &
766 40 : 'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
767 :
768 20 : print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=blk)
769 :
770 : ! List of atoms
771 20 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', n_rep_val=n)
772 :
773 20 : IF (n > 0) THEN
774 :
775 0 : IF (output_unit > 0) &
776 0 : WRITE (output_unit, '(/,T3,A/)') 'Block state fractions:'
777 :
778 : ! Number of AO functions
779 0 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
780 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
781 :
782 : ! MOs in orthonormal basis set
783 0 : ALLOCATE (mat_w(n_spins))
784 0 : DO j = 1, n_spins
785 0 : n_mo = ec%block(blk)%n_ao
786 : CALL cp_fm_create(matrix=mat_w(j), &
787 : matrix_struct=ec%block(blk)%mo(j)%mo_coeff%matrix_struct, &
788 0 : name='BLOCK MOs IN ORTHONORMAL BASIS SET')
789 : CALL parallel_gemm("N", "N", n_ao, n_mo, n_ao, 1.0_dp, ec%m_transf_inv, &
790 0 : ec%block(blk)%mo(j)%mo_coeff, 0.0_dp, mat_w(j))
791 : END DO
792 :
793 0 : DO j = 1, n
794 0 : NULLIFY (list_at)
795 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', &
796 0 : i_rep_val=j, i_vals=list_at)
797 0 : IF (ASSOCIATED(list_at)) THEN
798 :
799 : ! List of states
800 0 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', n_rep_val=m)
801 :
802 0 : IF (m > 0) THEN
803 :
804 0 : DO k = 1, m
805 0 : NULLIFY (list_mo)
806 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', &
807 0 : i_rep_val=k, i_vals=list_mo)
808 0 : IF (ASSOCIATED(list_mo)) THEN
809 :
810 0 : IF (j > 1) THEN
811 0 : IF (output_unit > 0) &
812 0 : WRITE (output_unit, *)
813 : END IF
814 :
815 0 : DO l = 1, SIZE(list_mo)
816 :
817 0 : IF (n_spins > 1) THEN
818 : c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
819 0 : list_mo(l), list_at)
820 : c2 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(2), &
821 0 : list_mo(l), list_at)
822 0 : IF (output_unit > 0) &
823 0 : WRITE (output_unit, '(I5,A,I5,2F20.10)') j, ' /', list_mo(l), c1, c2
824 : ELSE
825 : c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
826 0 : list_mo(l), list_at)
827 0 : IF (output_unit > 0) &
828 0 : WRITE (output_unit, '(I5,A,I5,F20.10)') j, ' /', list_mo(l), c1
829 : END IF
830 :
831 : END DO
832 :
833 : END IF
834 : END DO
835 :
836 : END IF
837 :
838 : END IF
839 : END DO
840 :
841 : ! Clean memory
842 0 : CALL cp_fm_release(mat_w)
843 :
844 : END IF
845 :
846 40 : END SUBROUTINE print_mo_coeff
847 :
848 : ! **************************************************************************************************
849 : !> \brief Print out electronic states (MOs)
850 : !> \param output_unit unit number of the open output stream
851 : !> \param mo array of MO sets
852 : !> \param n_spins number of spin components
853 : !> \param label output label
854 : !> \param mx_mo_a maximum number of alpha states to print out
855 : !> \param mx_mo_b maximum number of beta states to print out
856 : !> \param fermi print out Fermi level and number of electrons
857 : !> \author Z. Futera (02.2017)
858 : ! **************************************************************************************************
859 40 : SUBROUTINE print_states(output_unit, mo, n_spins, label, mx_mo_a, mx_mo_b, fermi)
860 :
861 : ! Routine arguments
862 : INTEGER, INTENT(IN) :: output_unit
863 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
864 : INTEGER, INTENT(IN) :: n_spins
865 : CHARACTER(LEN=*), INTENT(IN) :: label
866 : INTEGER, INTENT(IN), OPTIONAL :: mx_mo_a, mx_mo_b
867 : LOGICAL, INTENT(IN), OPTIONAL :: fermi
868 :
869 : INTEGER :: i, mx_a, mx_b, n
870 : LOGICAL :: prnt_fm
871 :
872 : ! Routine name for debug purposes
873 :
874 20 : prnt_fm = .FALSE.
875 20 : IF (PRESENT(fermi)) &
876 20 : prnt_fm = fermi
877 :
878 20 : IF (output_unit > 0) THEN
879 :
880 15 : WRITE (output_unit, '(/,T3,A/)') 'State energies ('//TRIM(ADJUSTL(label))//'):'
881 :
882 : ! Spin-polarized calculation
883 15 : IF (n_spins > 1) THEN
884 :
885 15 : mx_a = mo(1)%nmo
886 15 : IF (PRESENT(mx_mo_a)) &
887 10 : mx_a = MIN(mo(1)%nmo, mx_mo_a)
888 15 : mx_b = mo(2)%nmo
889 15 : IF (PRESENT(mx_mo_b)) &
890 10 : mx_b = MIN(mo(2)%nmo, mx_mo_b)
891 15 : n = MAX(mx_a, mx_b)
892 :
893 181 : DO i = 1, n
894 166 : WRITE (output_unit, '(T3,I10)', ADVANCE='no') i
895 166 : IF (i <= mx_a) THEN
896 : WRITE (output_unit, '(2F12.4)', ADVANCE='no') &
897 166 : mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
898 : ELSE
899 0 : WRITE (output_unit, '(A)', ADVANCE='no') ' '
900 : END IF
901 166 : WRITE (output_unit, '(A)', ADVANCE='no') ' '
902 181 : IF (i <= mx_b) THEN
903 : WRITE (output_unit, '(2F12.4)') &
904 161 : mo(2)%occupation_numbers(i), mo(2)%eigenvalues(i)
905 : ELSE
906 5 : WRITE (output_unit, *)
907 : END IF
908 : END DO
909 :
910 15 : IF (prnt_fm) THEN
911 : WRITE (output_unit, '(/,T3,I10,F24.4,I10,F19.4)') &
912 15 : mo(1)%nelectron, mo(1)%mu, &
913 30 : mo(2)%nelectron, mo(2)%mu
914 : END IF
915 :
916 : ! Spin-restricted calculation
917 : ELSE
918 :
919 0 : mx_a = mo(1)%nmo
920 0 : IF (PRESENT(mx_mo_a)) &
921 0 : mx_a = MIN(mo(1)%nmo, mx_mo_a)
922 :
923 0 : DO i = 1, mx_a
924 : WRITE (output_unit, '(T3,I10,2F12.4)') &
925 0 : i, mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
926 : END DO
927 :
928 0 : IF (prnt_fm) THEN
929 : WRITE (output_unit, '(/,T3,I10,F24.4)') &
930 0 : mo(1)%nelectron, mo(1)%mu
931 : END IF
932 :
933 : END IF
934 :
935 : END IF
936 :
937 20 : END SUBROUTINE print_states
938 :
939 : ! **************************************************************************************************
940 : !> \brief Print out donor-acceptor state couplings
941 : !> \param ec_sec ...
942 : !> \param output_unit unit number of the open output stream
943 : !> \param logger ...
944 : !> \param ec electronic coupling data structure
945 : !> \param mo ...
946 : !> \author Z. Futera (02.2017)
947 : ! **************************************************************************************************
948 10 : SUBROUTINE print_couplings(ec_sec, output_unit, logger, ec, mo)
949 :
950 : ! Routine arguments
951 : TYPE(section_vals_type), POINTER :: ec_sec
952 : INTEGER, INTENT(IN) :: output_unit
953 : TYPE(cp_logger_type), POINTER :: logger
954 : TYPE(et_cpl), POINTER :: ec
955 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
956 :
957 : CHARACTER(LEN=default_path_length) :: filename, my_pos, title
958 : INTEGER :: i, j, k, l, n_states(2), nc, nr, nspins, &
959 : unit_nr
960 : LOGICAL :: append
961 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: w1, w2
962 : TYPE(section_vals_type), POINTER :: print_key
963 :
964 : ! Routine name for debug purposes
965 : ! Local variables
966 :
967 10 : n_states = 0
968 30 : DO i = 1, SIZE(mo)
969 30 : n_states(i) = mo(i)%nmo
970 : END DO
971 10 : nspins = 1
972 10 : IF (n_states(2) > 0) nspins = 2
973 :
974 : print_key => section_vals_get_subs_vals(section_vals=ec_sec, &
975 10 : subsection_name="PRINT%COUPLINGS")
976 :
977 10 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
978 : cp_p_file)) THEN
979 :
980 10 : my_pos = "REWIND"
981 10 : append = section_get_lval(print_key, "APPEND")
982 10 : IF (append) THEN
983 0 : my_pos = "APPEND"
984 : END IF
985 :
986 10 : IF (output_unit > 0) &
987 5 : WRITE (output_unit, '(/,T3,A/)') 'Printing coupling elements to output files'
988 :
989 30 : DO i = 1, ec%n_blocks
990 40 : DO j = i + 1, ec%n_blocks
991 :
992 10 : nr = ec%block(i)%hab(1, j)%matrix_struct%nrow_global
993 10 : nc = ec%block(i)%hab(1, j)%matrix_struct%ncol_global
994 :
995 40 : ALLOCATE (w1(nr, nc))
996 10 : CPASSERT(ASSOCIATED(w1))
997 10 : CALL cp_fm_get_submatrix(ec%block(i)%hab(1, j), w1)
998 10 : IF (nspins > 1) THEN
999 30 : ALLOCATE (w2(nr, nc))
1000 10 : CPASSERT(ASSOCIATED(w2))
1001 10 : CALL cp_fm_get_submatrix(ec%block(i)%hab(2, j), w2)
1002 : END IF
1003 :
1004 10 : IF (output_unit > 0) THEN
1005 :
1006 5 : WRITE (filename, '(a5,I1.1,a1,I1.1)') "ET_BL_", i, "-", j
1007 : unit_nr = cp_print_key_unit_nr(logger, ec_sec, "PRINT%COUPLINGS", extension=".elcoup", &
1008 5 : middle_name=TRIM(filename), file_position=my_pos, log_filename=.FALSE.)
1009 :
1010 5 : WRITE (title, *) 'Coupling elements [meV] between blocks:', i, j
1011 :
1012 5 : WRITE (unit_nr, *) TRIM(title)
1013 5 : IF (nspins > 1) THEN
1014 5 : WRITE (unit_nr, '(T3,A8,T13,A8,T28,A,A)') 'State A', 'State B', 'Coupling spin 1', ' Coupling spin 2'
1015 : ELSE
1016 0 : WRITE (unit_nr, '(T3,A8,T13,A8,T28,A)') 'State A', 'State B', 'Coupling'
1017 : END IF
1018 :
1019 55 : DO k = 1, MIN(ec%block(i)%n_ao, n_states(1))
1020 836 : DO l = 1, MIN(ec%block(j)%n_ao, n_states(1))
1021 :
1022 836 : IF (nspins > 1) THEN
1023 :
1024 : WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)', ADVANCE='no') &
1025 786 : k, l, w1(k, l)*evolt*1000.0_dp
1026 786 : IF ((k <= n_states(2)) .AND. (l <= n_states(2))) THEN
1027 : WRITE (unit_nr, '(E20.6)') &
1028 779 : w2(k, l)*evolt*1000.0_dp
1029 : ELSE
1030 7 : WRITE (unit_nr, *)
1031 : END IF
1032 :
1033 : ELSE
1034 :
1035 : WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)') &
1036 0 : k, l, w1(k, l)*evolt*1000.0_dp
1037 : END IF
1038 :
1039 : END DO
1040 55 : WRITE (unit_nr, *)
1041 : END DO
1042 5 : CALL cp_print_key_finished_output(unit_nr, logger, ec_sec, "PRINT%COUPLINGS")
1043 :
1044 : END IF
1045 :
1046 10 : IF (ASSOCIATED(w1)) DEALLOCATE (w1)
1047 30 : IF (ASSOCIATED(w2)) DEALLOCATE (w2)
1048 :
1049 : END DO
1050 : END DO
1051 :
1052 : END IF
1053 10 : END SUBROUTINE print_couplings
1054 :
1055 : ! **************************************************************************************************
1056 : !> \brief Normalize set of MO vectors
1057 : !> \param qs_env QuickStep environment containing all system data
1058 : !> \param mo storage for the MO data set
1059 : !> \param n_ao number of AO basis functions
1060 : !> \param n_mo number of block states
1061 : !> \author Z. Futera (02.2017)
1062 : ! **************************************************************************************************
1063 40 : SUBROUTINE normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1064 :
1065 : ! Routine arguments
1066 : TYPE(qs_environment_type), POINTER :: qs_env
1067 : TYPE(mo_set_type), POINTER :: mo
1068 : INTEGER, INTENT(IN) :: n_ao, n_mo
1069 :
1070 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_t
1071 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1072 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1073 : TYPE(cp_fm_type) :: mat_sc, mat_t
1074 40 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
1075 : TYPE(mp_para_env_type), POINTER :: para_env
1076 :
1077 : ! Routine name for debug purposes
1078 :
1079 : ! Initialization
1080 40 : NULLIFY (blacs_env)
1081 40 : NULLIFY (para_env)
1082 40 : NULLIFY (fm_s)
1083 40 : NULLIFY (mat_s)
1084 : NULLIFY (vec_t)
1085 :
1086 : ! Overlap matrix
1087 40 : CALL get_qs_env(qs_env, matrix_s=mat_s)
1088 :
1089 : ! Calculate S*C product
1090 : CALL cp_fm_create(matrix=mat_sc, matrix_struct=mo%mo_coeff%matrix_struct, &
1091 40 : name='S*C PRODUCT MATRIX')
1092 40 : CALL cp_dbcsr_sm_fm_multiply(mat_s(1)%matrix, mo%mo_coeff, mat_sc, n_mo)
1093 :
1094 : ! Calculate C^T*S*C
1095 40 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1096 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1097 40 : nrow_global=n_mo, ncol_global=n_mo)
1098 : CALL cp_fm_create(matrix=mat_t, matrix_struct=fm_s, &
1099 40 : name='C^T*S*C OVERLAP PRODUCT MATRIX')
1100 40 : CALL parallel_gemm('T', 'N', n_mo, n_mo, n_ao, 1.0_dp, mo%mo_coeff, mat_sc, 0.0_dp, mat_t)
1101 :
1102 : ! Normalization
1103 120 : ALLOCATE (vec_t(n_mo))
1104 40 : CPASSERT(ASSOCIATED(vec_t))
1105 40 : CALL cp_fm_vectorssum(mat_t, vec_t)
1106 448 : vec_t = 1.0_dp/DSQRT(vec_t)
1107 40 : CALL cp_fm_column_scale(mo%mo_coeff, vec_t)
1108 :
1109 : ! Clean memory
1110 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
1111 40 : CALL cp_fm_release(matrix=mat_sc)
1112 40 : CALL cp_fm_release(matrix=mat_t)
1113 40 : IF (ASSOCIATED(vec_t)) &
1114 40 : DEALLOCATE (vec_t)
1115 :
1116 80 : END SUBROUTINE normalize_mo_vectors
1117 :
1118 : ! **************************************************************************************************
1119 : !> \brief Transform block MO coefficients to original non-orthogonal basis set and save them
1120 : !> \param qs_env QuickStep environment containing all system data
1121 : !> \param ec electronic coupling data structure
1122 : !> \param id block ID
1123 : !> \param mo storage for the MO data set
1124 : !> \param mat_u matrix of the block states
1125 : !> \param n_ao number of AO basis functions
1126 : !> \param n_mo number of block states
1127 : !> \author Z. Futera (02.2017)
1128 : ! **************************************************************************************************
1129 40 : SUBROUTINE set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1130 :
1131 : ! Routine arguments
1132 : TYPE(qs_environment_type), POINTER :: qs_env
1133 : TYPE(et_cpl), POINTER :: ec
1134 : INTEGER, INTENT(IN) :: id
1135 : TYPE(mo_set_type), POINTER :: mo
1136 : TYPE(cp_fm_type), INTENT(IN) :: mat_u
1137 : INTEGER, INTENT(IN) :: n_ao, n_mo
1138 :
1139 : INTEGER :: ic, ir, jc, jr, mr, nc, nr
1140 : REAL(KIND=dp) :: xu
1141 : TYPE(cp_fm_type) :: mat_w
1142 :
1143 : ! Routine name for debug purposes
1144 : ! Local variables
1145 :
1146 : ! Working matrix
1147 : CALL cp_fm_create(matrix=mat_w, matrix_struct=mo%mo_coeff%matrix_struct, &
1148 40 : name='BLOCK MO-TRANSFORMATION WORKING MATRIX')
1149 40 : CALL cp_fm_set_all(mat_w, 0.0_dp)
1150 :
1151 : ! Matrix-element reordering
1152 40 : nr = 1
1153 : ! Rows
1154 184 : DO ir = 1, ec%block(id)%n_atoms
1155 592 : DO jr = 1, ec%block(id)%atom(ir)%n_ao
1156 : ! Columns
1157 408 : nc = 1
1158 2832 : DO ic = 1, ec%block(id)%n_atoms
1159 9192 : DO jc = 1, ec%block(id)%atom(ic)%n_ao
1160 6360 : mr = ec%block(id)%atom(ir)%ao_pos + jr - 1
1161 6360 : CALL cp_fm_get_element(mat_u, nr, nc, xu)
1162 6360 : CALL cp_fm_set_element(mat_w, mr, nc, xu)
1163 15144 : nc = nc + 1
1164 : END DO
1165 : END DO
1166 552 : nr = nr + 1
1167 : END DO
1168 : END DO
1169 :
1170 : ! Transformation to original non-orthogonal basis set
1171 40 : CALL parallel_gemm("N", "N", n_ao, n_mo, n_ao, 1.0_dp, ec%m_transf, mat_w, 0.0_dp, mo%mo_coeff)
1172 40 : CALL normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1173 :
1174 : ! Clean memory
1175 40 : CALL cp_fm_release(matrix=mat_w)
1176 :
1177 40 : END SUBROUTINE set_mo_coefficients
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief Creates MO set corresponding to one atomic data block
1181 : !> \param qs_env QuickStep environment containing all system data
1182 : !> \param ec electronic coupling data structure
1183 : !> \param id block ID
1184 : !> \param spin spin component
1185 : !> \param mat_u matrix of the block states
1186 : !> \param vec_e array of the block eigenvalues
1187 : !> \author Z. Futera (02.2017)
1188 : ! **************************************************************************************************
1189 40 : SUBROUTINE create_block_mo_set(qs_env, ec, id, spin, mat_u, vec_e)
1190 :
1191 : ! Routine arguments
1192 : TYPE(qs_environment_type), POINTER :: qs_env
1193 : TYPE(et_cpl), POINTER :: ec
1194 : INTEGER, INTENT(IN) :: id, spin
1195 : TYPE(cp_fm_type), INTENT(IN) :: mat_u
1196 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
1197 :
1198 : INTEGER :: n_ao, n_el, n_mo
1199 : REAL(KIND=dp) :: mx_occ
1200 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1201 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1202 : TYPE(dft_control_type), POINTER :: dft_cntrl
1203 : TYPE(mo_set_type), POINTER :: mo
1204 : TYPE(mp_para_env_type), POINTER :: para_env
1205 40 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1206 : TYPE(scf_control_type), POINTER :: scf_cntrl
1207 :
1208 : ! Routine name for debug purposes
1209 :
1210 40 : NULLIFY (blacs_env)
1211 40 : NULLIFY (dft_cntrl)
1212 40 : NULLIFY (para_env)
1213 40 : NULLIFY (qs_kind_set)
1214 40 : NULLIFY (fm_s)
1215 40 : NULLIFY (scf_cntrl)
1216 40 : NULLIFY (mo)
1217 :
1218 : ! Number of basis functions
1219 40 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1220 40 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1221 :
1222 : ! Number of states
1223 40 : n_mo = mat_u%matrix_struct%nrow_global
1224 40 : IF (n_mo /= mat_u%matrix_struct%ncol_global) &
1225 0 : CPABORT('block state matrix is not square')
1226 40 : IF (n_mo /= SIZE(vec_e)) &
1227 0 : CPABORT('inconsistent number of states / energies')
1228 :
1229 : ! Maximal occupancy
1230 40 : CALL get_qs_env(qs_env, dft_control=dft_cntrl)
1231 40 : mx_occ = 2.0_dp
1232 40 : IF (dft_cntrl%nspins > 1) &
1233 40 : mx_occ = 1.0_dp
1234 :
1235 : ! Number of electrons
1236 40 : n_el = ec%block(id)%n_electrons
1237 40 : IF (dft_cntrl%nspins > 1) THEN
1238 40 : n_el = n_el/2
1239 40 : IF (MOD(ec%block(id)%n_electrons, 2) == 1) THEN
1240 28 : IF (spin == 1) &
1241 14 : n_el = n_el + 1
1242 : END IF
1243 : END IF
1244 :
1245 : ! Memory allocation (Use deallocate_mo_set to prevent accidental memory leaks)
1246 40 : CALL deallocate_mo_set(ec%block(id)%mo(spin))
1247 40 : CALL allocate_mo_set(ec%block(id)%mo(spin), n_ao, n_mo, n_el, REAL(n_el, dp), mx_occ, 0.0_dp)
1248 40 : mo => ec%block(id)%mo(spin)
1249 :
1250 : ! State energies
1251 120 : ALLOCATE (mo%eigenvalues(n_mo))
1252 40 : CPASSERT(ASSOCIATED(mo%eigenvalues))
1253 896 : mo%eigenvalues = vec_e
1254 :
1255 : ! States coefficients
1256 40 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1257 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1258 40 : nrow_global=n_ao, ncol_global=n_mo)
1259 40 : ALLOCATE (mo%mo_coeff)
1260 40 : CALL cp_fm_create(matrix=mo%mo_coeff, matrix_struct=fm_s, name='BLOCK STATES')
1261 :
1262 : ! Transform MO coefficients to original non-orthogonal basis set
1263 40 : CALL set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1264 :
1265 : ! Occupancies
1266 120 : ALLOCATE (mo%occupation_numbers(n_mo))
1267 40 : CPASSERT(ASSOCIATED(mo%occupation_numbers))
1268 448 : mo%occupation_numbers = 0.0_dp
1269 :
1270 40 : IF (n_el > 0) THEN
1271 28 : CALL get_qs_env(qs_env, scf_control=scf_cntrl)
1272 28 : CALL set_mo_occupation(mo_set=mo, smear=scf_cntrl%smear)
1273 : END IF
1274 :
1275 : ! Clean memory
1276 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
1277 :
1278 40 : END SUBROUTINE create_block_mo_set
1279 :
1280 : ! **************************************************************************************************
1281 : !> \brief save given electronic state to cube files
1282 : !> \param qs_env QuickStep environment containing all system data
1283 : !> \param logger output logger
1284 : !> \param input input-file block print setting section
1285 : !> \param mo electronic states data
1286 : !> \param ib block ID
1287 : !> \param im state ID
1288 : !> \param is spin ID
1289 : !> \author Z. Futera (02.2017)
1290 : ! **************************************************************************************************
1291 0 : SUBROUTINE save_mo_cube(qs_env, logger, input, mo, ib, im, is)
1292 :
1293 : ! Routine arguments
1294 : TYPE(qs_environment_type), POINTER :: qs_env
1295 : TYPE(cp_logger_type), POINTER :: logger
1296 : TYPE(section_vals_type), POINTER :: input
1297 : TYPE(mo_set_type), POINTER :: mo
1298 : INTEGER, INTENT(IN) :: ib, im, is
1299 :
1300 : CHARACTER(LEN=default_path_length) :: filename
1301 : CHARACTER(LEN=default_string_length) :: title
1302 : INTEGER :: unit_nr
1303 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1304 : TYPE(cell_type), POINTER :: cell
1305 : TYPE(dft_control_type), POINTER :: dft_control
1306 : TYPE(particle_list_type), POINTER :: particles
1307 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1308 : TYPE(pw_c1d_gs_type) :: wf_g
1309 : TYPE(pw_env_type), POINTER :: pw_env
1310 0 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1311 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1312 : TYPE(pw_r3d_rs_type) :: wf_r
1313 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1314 : TYPE(qs_subsys_type), POINTER :: subsys
1315 :
1316 : ! Routine name for debug purposes
1317 :
1318 : ! Initialization
1319 0 : NULLIFY (particles)
1320 0 : NULLIFY (subsys)
1321 :
1322 0 : NULLIFY (pw_env)
1323 0 : NULLIFY (pw_pools)
1324 0 : NULLIFY (auxbas_pw_pool)
1325 :
1326 0 : NULLIFY (atomic_kind_set)
1327 0 : NULLIFY (cell)
1328 0 : NULLIFY (dft_control)
1329 0 : NULLIFY (particle_set)
1330 0 : NULLIFY (qs_kind_set)
1331 :
1332 : ! Name of the cube file
1333 0 : WRITE (filename, '(A4,I1.1,A1,I5.5,A1,I1.1)') 'BWF_', ib, '_', im, '_', is
1334 : ! Open the file
1335 : unit_nr = cp_print_key_unit_nr(logger, input, 'MO_CUBES', extension='.cube', &
1336 0 : middle_name=TRIM(filename), file_position='REWIND', log_filename=.FALSE.)
1337 : ! Title of the file
1338 0 : WRITE (title, *) 'WAVEFUNCTION ', im, ' block ', ib, ' spin ', is
1339 :
1340 : ! List of all atoms
1341 0 : CALL get_qs_env(qs_env, subsys=subsys)
1342 0 : CALL qs_subsys_get(subsys, particles=particles)
1343 :
1344 : ! Grids for wavefunction
1345 0 : CALL get_qs_env(qs_env, pw_env=pw_env)
1346 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1347 0 : CALL auxbas_pw_pool%create_pw(wf_r)
1348 0 : CALL auxbas_pw_pool%create_pw(wf_g)
1349 :
1350 : ! Calculate the grid values
1351 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1352 0 : cell=cell, dft_control=dft_control, particle_set=particle_set)
1353 : CALL calculate_wavefunction(mo%mo_coeff, im, wf_r, wf_g, atomic_kind_set, &
1354 0 : qs_kind_set, cell, dft_control, particle_set, pw_env)
1355 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, &
1356 0 : stride=section_get_ivals(input, 'MO_CUBES%STRIDE'))
1357 :
1358 : ! Close file
1359 0 : CALL cp_print_key_finished_output(unit_nr, logger, input, 'MO_CUBES')
1360 :
1361 : ! Clean memory
1362 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1363 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1364 :
1365 0 : END SUBROUTINE save_mo_cube
1366 :
1367 : ! **************************************************************************************************
1368 : !> \brief save specified electronic states to cube files
1369 : !> \param qs_env QuickStep environment containing all system data
1370 : !> \param ec electronic coupling data structure
1371 : !> \param n_spins number of spin states
1372 : !> \author Z. Futera (02.2017)
1373 : ! **************************************************************************************************
1374 10 : SUBROUTINE save_el_states(qs_env, ec, n_spins)
1375 :
1376 : ! Routine arguments
1377 : TYPE(qs_environment_type), POINTER :: qs_env
1378 : TYPE(et_cpl), POINTER :: ec
1379 : INTEGER, INTENT(IN) :: n_spins
1380 :
1381 : INTEGER :: i, j, k, l, n
1382 10 : INTEGER, DIMENSION(:), POINTER :: list
1383 : TYPE(cp_logger_type), POINTER :: logger
1384 : TYPE(mo_set_type), POINTER :: mo
1385 : TYPE(section_vals_type), POINTER :: block_sec, mo_sec, print_sec
1386 :
1387 : ! Routine name for debug purposes
1388 :
1389 10 : NULLIFY (logger)
1390 10 : NULLIFY (block_sec)
1391 10 : NULLIFY (print_sec)
1392 10 : NULLIFY (mo_sec)
1393 :
1394 : ! Output logger
1395 20 : logger => cp_get_default_logger()
1396 : block_sec => section_vals_get_subs_vals(qs_env%input, &
1397 10 : 'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
1398 :
1399 : ! Print states of all blocks
1400 30 : DO i = 1, ec%n_blocks
1401 :
1402 20 : print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=i)
1403 :
1404 : ! Check if the print input section is active
1405 20 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1406 10 : print_sec, 'MO_CUBES'), cp_p_file)) THEN
1407 :
1408 0 : mo_sec => section_vals_get_subs_vals(print_sec, 'MO_CUBES')
1409 :
1410 : ! Spin states
1411 0 : DO j = 1, n_spins
1412 :
1413 0 : mo => ec%block(i)%mo(j)
1414 :
1415 0 : CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', n_rep_val=n)
1416 :
1417 : ! List of specific MOs
1418 0 : IF (n > 0) THEN
1419 :
1420 0 : DO k = 1, n
1421 0 : NULLIFY (list)
1422 : CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', &
1423 0 : i_rep_val=k, i_vals=list)
1424 0 : IF (ASSOCIATED(list)) THEN
1425 0 : DO l = 1, SIZE(list)
1426 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, list(l), j)
1427 : END DO
1428 : END IF
1429 : END DO
1430 :
1431 : ! Frontier MOs
1432 : ELSE
1433 :
1434 : ! Occupied states
1435 0 : CALL section_vals_val_get(mo_sec, keyword_name='NHOMO', i_val=n)
1436 :
1437 0 : IF (n > 0) THEN
1438 0 : DO k = MAX(1, mo%homo - n + 1), mo%homo
1439 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1440 : END DO
1441 : END IF
1442 :
1443 : ! Unoccupied states
1444 0 : CALL section_vals_val_get(mo_sec, keyword_name='NLUMO', i_val=n)
1445 :
1446 0 : IF (n > 0) THEN
1447 0 : DO k = mo%lfomo, MIN(mo%lfomo + n - 1, mo%nmo)
1448 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1449 : END DO
1450 : END IF
1451 :
1452 : END IF
1453 :
1454 : END DO
1455 :
1456 : END IF
1457 :
1458 : END DO
1459 :
1460 10 : END SUBROUTINE save_el_states
1461 :
1462 : ! **************************************************************************************************
1463 : !> \brief calculates the electron transfer coupling elements by projection-operator approach
1464 : !> Kondov et al. J.Phys.Chem.C 2007, 111, 11970-11981
1465 : !> \param qs_env QuickStep environment containing all system data
1466 : !> \author Z. Futera (02.2017)
1467 : ! **************************************************************************************************
1468 10 : SUBROUTINE calc_et_coupling_proj(qs_env)
1469 :
1470 : ! Routine arguments
1471 : TYPE(qs_environment_type), POINTER :: qs_env
1472 :
1473 : INTEGER :: i, j, k, n_ao, n_atoms, output_unit
1474 : LOGICAL :: do_kp, master
1475 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1476 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1477 : TYPE(cp_fm_type) :: mat_w
1478 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_h
1479 : TYPE(cp_logger_type), POINTER :: logger
1480 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks, mo_der
1481 : TYPE(dft_control_type), POINTER :: dft_cntrl
1482 : TYPE(et_cpl), POINTER :: ec
1483 : TYPE(kpoint_type), POINTER :: kpoints
1484 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo
1485 : TYPE(mp_para_env_type), POINTER :: para_env
1486 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1487 : TYPE(scf_control_type), POINTER :: scf_control
1488 : TYPE(section_vals_type), POINTER :: et_proj_sec
1489 :
1490 : ! Routine name for debug purposes
1491 :
1492 : ! Pointer initialization
1493 10 : NULLIFY (logger)
1494 :
1495 10 : NULLIFY (blacs_env)
1496 10 : NULLIFY (para_env)
1497 10 : NULLIFY (dft_cntrl)
1498 10 : NULLIFY (kpoints)
1499 10 : NULLIFY (qs_kind_set)
1500 : NULLIFY (et_proj_sec)
1501 :
1502 10 : NULLIFY (fm_s)
1503 10 : NULLIFY (ks, mo_der)
1504 :
1505 10 : NULLIFY (ec)
1506 :
1507 : ! Reference
1508 10 : CALL cite_reference(Futera2017)
1509 :
1510 : ! Stream for output to LOG file
1511 10 : logger => cp_get_default_logger()
1512 :
1513 10 : et_proj_sec => section_vals_get_subs_vals(qs_env%input, 'PROPERTIES%ET_COUPLING%PROJECTION')
1514 :
1515 : output_unit = cp_print_key_unit_nr(logger, et_proj_sec, &
1516 10 : 'PROGRAM_RUN_INFO', extension='.log')
1517 :
1518 : ! Parallel calculation - master thread
1519 10 : master = .FALSE.
1520 10 : IF (output_unit > 0) &
1521 : master = .TRUE.
1522 :
1523 : ! Header
1524 : IF (master) THEN
1525 : WRITE (output_unit, '(/,T2,A)') &
1526 5 : '!-----------------------------------------------------------------------------!'
1527 : WRITE (output_unit, '(T17,A)') &
1528 5 : 'Electronic coupling - Projection-operator method'
1529 : END IF
1530 :
1531 : ! Main data structure
1532 10 : ALLOCATE (ec)
1533 10 : CPASSERT(ASSOCIATED(ec))
1534 10 : CALL set_block_data(qs_env, et_proj_sec, ec)
1535 :
1536 : ! Number of atoms and AO functions
1537 10 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=n_atoms)
1538 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1539 :
1540 : ! Print out info about system partitioning
1541 10 : IF (master) THEN
1542 :
1543 : WRITE (output_unit, '(/,T3,A,I10)') &
1544 5 : 'Number of atoms = ', n_atoms
1545 : WRITE (output_unit, '(T3,A,I10)') &
1546 5 : 'Number of fragments = ', ec%n_blocks
1547 : WRITE (output_unit, '(T3,A,I10)') &
1548 5 : 'Number of fragment atoms = ', ec%n_atoms
1549 : WRITE (output_unit, '(T3,A,I10)') &
1550 5 : 'Number of unassigned atoms = ', n_atoms - ec%n_atoms
1551 : WRITE (output_unit, '(T3,A,I10)') &
1552 5 : 'Number of AO basis functions = ', n_ao
1553 :
1554 15 : DO i = 1, ec%n_blocks
1555 :
1556 : WRITE (output_unit, '(/,T3,A,I0,A)') &
1557 10 : 'Block ', i, ':'
1558 : WRITE (output_unit, '(T3,A,I10)') &
1559 10 : 'Number of block atoms = ', ec%block(i)%n_atoms
1560 : WRITE (output_unit, '(T3,A,I10)') &
1561 10 : 'Number of block electrons = ', ec%block(i)%n_electrons
1562 : WRITE (output_unit, '(T3,A,I10)') &
1563 10 : 'Number of block AO functions = ', ec%block(i)%n_ao
1564 :
1565 15 : IF (ec%block(i)%n_atoms < 10) THEN
1566 :
1567 : WRITE (output_unit, '(T3,A,10I6)') &
1568 10 : 'Block atom IDs = ', &
1569 56 : (ec%block(i)%atom(j)%id, j=1, ec%block(i)%n_atoms)
1570 :
1571 : ELSE
1572 :
1573 0 : WRITE (output_unit, '(T3,A)') 'Block atom IDs ='
1574 0 : DO j = 1, ec%block(i)%n_atoms/10
1575 0 : WRITE (output_unit, '(T3,A,10I6)') ' ', &
1576 0 : (ec%block(i)%atom((j - 1)*10 + k)%id, k=1, 10)
1577 : END DO
1578 0 : IF (MOD(ec%block(i)%n_atoms, 10) /= 0) THEN
1579 0 : WRITE (output_unit, '(T3,A,10I6)') ' ', &
1580 0 : (ec%block(i)%atom(k + 10*(ec%block(i)%n_atoms/10))%id, &
1581 0 : k=1, MOD(ec%block(i)%n_atoms, 10))
1582 : END IF
1583 :
1584 : END IF
1585 :
1586 : END DO
1587 :
1588 : END IF
1589 :
1590 : ! Full matrix data structure
1591 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1592 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1593 10 : nrow_global=n_ao, ncol_global=n_ao)
1594 10 : CALL cp_fm_create(matrix=mat_w, matrix_struct=fm_s, name='FULL WORK MATRIX')
1595 :
1596 : ! Spin polarization / K-point sampling
1597 10 : CALL get_qs_env(qs_env, dft_control=dft_cntrl, do_kpoints=do_kp)
1598 10 : CALL get_qs_env(qs_env, mos=mo, matrix_ks=ks, mo_derivs=mo_der, scf_control=scf_control)
1599 10 : CALL make_mo_eig(mo, dft_cntrl%nspins, ks, scf_control, mo_der)
1600 :
1601 10 : IF (do_kp) THEN
1602 0 : CPABORT('ET_COUPLING not implemented with kpoints')
1603 : ELSE
1604 : ! no K-points
1605 10 : IF (master) &
1606 5 : WRITE (output_unit, '(T3,A)') 'No K-point sampling (Gamma point only)'
1607 : END IF
1608 :
1609 10 : IF (dft_cntrl%nspins == 2) THEN
1610 :
1611 10 : IF (master) &
1612 5 : WRITE (output_unit, '(/,T3,A)') 'Spin-polarized calculation'
1613 :
1614 : !<--- Open shell / No K-points ------------------------------------------------>!
1615 :
1616 : ! State eneries of the whole system
1617 10 : IF (mo(1)%nao /= mo(2)%nao) &
1618 0 : CPABORT('different number of alpha/beta AO basis functions')
1619 10 : IF (master) THEN
1620 : WRITE (output_unit, '(/,T3,A,I10)') &
1621 5 : 'Number of AO basis functions = ', mo(1)%nao
1622 : WRITE (output_unit, '(T3,A,I10)') &
1623 5 : 'Number of alpha states = ', mo(1)%nmo
1624 : WRITE (output_unit, '(T3,A,I10)') &
1625 5 : 'Number of beta states = ', mo(2)%nmo
1626 : END IF
1627 10 : CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.TRUE.)
1628 10 : CALL set_fermi(ec, mo(1)%mu, mo(2)%mu)
1629 :
1630 : ! KS Hamiltonian
1631 10 : CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1632 :
1633 : ! Block diagonization
1634 10 : CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1635 :
1636 : ! Print out energies and couplings
1637 30 : DO i = 1, ec%n_blocks
1638 20 : IF (output_unit > 0) THEN
1639 : CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1640 : 'block '//TRIM(ADJUSTL(cp_to_string(i)))//' states', &
1641 10 : mx_mo_a=mo(1)%nmo, mx_mo_b=mo(2)%nmo, fermi=.TRUE.)
1642 : END IF
1643 30 : CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1644 : END DO
1645 :
1646 10 : CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1647 :
1648 : ELSE
1649 :
1650 0 : IF (master) &
1651 0 : WRITE (output_unit, '(/,T3,A)') 'Spin-restricted calculation'
1652 :
1653 : !<--- Close shell / No K-points ----------------------------------------------->!
1654 :
1655 : ! State eneries of the whole system
1656 : IF (master) THEN
1657 : WRITE (output_unit, '(/,T3,A,I10)') &
1658 0 : 'Number of AO basis functions = ', mo(1)%nao
1659 : WRITE (output_unit, '(T3,A,I10)') &
1660 0 : 'Number of states = ', mo(1)%nmo
1661 : END IF
1662 0 : CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.TRUE.)
1663 0 : CALL set_fermi(ec, mo(1)%mu)
1664 :
1665 : ! KS Hamiltonian
1666 0 : CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1667 :
1668 : ! Block diagonization
1669 0 : CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1670 :
1671 : ! Print out energies and couplings
1672 0 : DO i = 1, ec%n_blocks
1673 0 : IF (output_unit > 0) THEN
1674 : CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1675 : 'block '//TRIM(ADJUSTL(cp_to_string(i)))//' states', &
1676 0 : mx_mo_a=mo(1)%nmo, fermi=.TRUE.)
1677 : END IF
1678 0 : CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1679 : END DO
1680 :
1681 0 : CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1682 :
1683 : END IF
1684 :
1685 : ! Save electronic states
1686 10 : CALL save_el_states(qs_env, ec, dft_cntrl%nspins)
1687 :
1688 : ! Footer
1689 10 : IF (master) WRITE (output_unit, '(/,T2,A)') &
1690 5 : '!-----------------------------------------------------------------------------!'
1691 :
1692 : ! Clean memory
1693 10 : CALL cp_fm_struct_release(fmstruct=fm_s)
1694 10 : CALL cp_fm_release(matrix=mat_w)
1695 10 : IF (ALLOCATED(mat_h)) THEN
1696 30 : DO i = 1, SIZE(mat_h)
1697 30 : CALL cp_fm_release(matrix=mat_h(i))
1698 : END DO
1699 10 : DEALLOCATE (mat_h)
1700 : END IF
1701 10 : CALL release_ec_data(ec)
1702 :
1703 : ! Close output stream
1704 10 : CALL cp_print_key_finished_output(output_unit, logger, et_proj_sec, 'PROGRAM_RUN_INFO')
1705 :
1706 20 : END SUBROUTINE calc_et_coupling_proj
1707 :
1708 0 : END MODULE et_coupling_proj
|
