(git:b279b6b)
et_coupling.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 calculates the electron transfer coupling elements
10 !> Wu, Van Voorhis, JCP 125, 164105 (2006)
11 !> \author fschiff (01.2007)
12 ! **************************************************************************************************
13 MODULE et_coupling
14  USE cp_control_types, ONLY: dft_control_type
15  USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
16  dbcsr_deallocate_matrix_set
17  USE cp_fm_basic_linalg, ONLY: cp_fm_det,&
18  cp_fm_invert,&
19  cp_fm_transpose
20  USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
21  fm_pool_get_el_struct
22  USE cp_fm_struct, ONLY: cp_fm_struct_type
23  USE cp_fm_types, ONLY: cp_fm_create,&
24  cp_fm_get_info,&
25  cp_fm_release,&
26  cp_fm_type
27  USE cp_log_handling, ONLY: cp_get_default_logger,&
28  cp_logger_type,&
29  cp_to_string
30  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31  cp_print_key_unit_nr
32  USE dbcsr_api, ONLY: dbcsr_p_type
33  USE input_section_types, ONLY: section_vals_get_subs_vals,&
34  section_vals_type
35  USE kinds, ONLY: dp
36  USE mathlib, ONLY: diamat_all
37  USE message_passing, ONLY: mp_para_env_type
38  USE parallel_gemm_api, ONLY: parallel_gemm
39  USE qs_energy_types, ONLY: qs_energy_type
40  USE qs_environment_types, ONLY: get_qs_env,&
41  qs_environment_type
42  USE qs_matrix_pools, ONLY: mpools_get
43  USE qs_mo_types, ONLY: get_mo_set
44 #include "./base/base_uses.f90"
45 
46  IMPLICIT NONE
47 
48  PRIVATE
49 
50  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'et_coupling'
51  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
52 
53 ! *** Public subroutines ***
54 
55  PUBLIC :: calc_et_coupling
56 
57 CONTAINS
58 ! **************************************************************************************************
59 !> \brief ...
60 !> \param qs_env ...
61 ! **************************************************************************************************
62  SUBROUTINE calc_et_coupling(qs_env)
63 
64  TYPE(qs_environment_type), POINTER :: qs_env
65 
66  CHARACTER(len=*), PARAMETER :: routinen = 'calc_et_coupling'
67 
68  INTEGER :: handle, i, iw, j, k, nao, ncol_local, &
69  nmo, nrow_local
70  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
71  LOGICAL :: is_spin_constraint
72  REAL(kind=dp) :: sda, strength, waa, wbb, wda
73  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: a, b, s_det
74  REAL(kind=dp), DIMENSION(2) :: eigenv
75  REAL(kind=dp), DIMENSION(2, 2) :: s_mat, tmp_mat, u, w_mat
76  TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: mo_mo_fm_pools
77  TYPE(cp_fm_struct_type), POINTER :: mo_mo_fmstruct
78  TYPE(cp_fm_type) :: inverse_mat, smo, tinverse, tmp2
79  TYPE(cp_fm_type), DIMENSION(2) :: rest_mo
80  TYPE(cp_logger_type), POINTER :: logger
81  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
82  TYPE(dft_control_type), POINTER :: dft_control
83  TYPE(mp_para_env_type), POINTER :: para_env
84  TYPE(qs_energy_type), POINTER :: energy
85  TYPE(section_vals_type), POINTER :: et_coupling_section
86 
87  NULLIFY (mo_mo_fmstruct, energy, matrix_s, dft_control, para_env)
88 
89  CALL timeset(routinen, handle)
90 
91  logger => cp_get_default_logger()
92  et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
93  "PROPERTIES%ET_COUPLING")
94 
95  CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
96 
97  iw = cp_print_key_unit_nr(logger, et_coupling_section, "PROGRAM_RUN_INFO", &
98  extension=".log")
99 
100  is_spin_constraint = .false.
101  ALLOCATE (a(dft_control%nspins))
102  ALLOCATE (b(dft_control%nspins))
103  ALLOCATE (s_det(dft_control%nspins))
104 
105  CALL mpools_get(qs_env%mpools, mo_mo_fm_pools=mo_mo_fm_pools)
106  mo_mo_fmstruct => fm_pool_get_el_struct(mo_mo_fm_pools(1)%pool)
107  DO i = 1, dft_control%nspins
108  mo_mo_fmstruct => fm_pool_get_el_struct(mo_mo_fm_pools(i)%pool)
109 
110  CALL get_mo_set(mo_set=qs_env%mos(i), &
111  nao=nao, &
112  nmo=nmo)
113 
114  CALL cp_fm_create(matrix=tmp2, &
115  matrix_struct=qs_env%mos(i)%mo_coeff%matrix_struct, &
116  name="ET_TMP"//trim(adjustl(cp_to_string(2)))//"MATRIX")
117  CALL cp_fm_create(matrix=inverse_mat, &
118  matrix_struct=mo_mo_fmstruct, &
119  name="INVERSE"//trim(adjustl(cp_to_string(2)))//"MATRIX")
120  CALL cp_fm_create(matrix=tinverse, &
121  matrix_struct=mo_mo_fmstruct, &
122  name="T_INVERSE"//trim(adjustl(cp_to_string(2)))//"MATRIX")
123  CALL cp_fm_create(matrix=smo, &
124  matrix_struct=mo_mo_fmstruct, &
125  name="ET_SMO"//trim(adjustl(cp_to_string(1)))//"MATRIX")
126  DO j = 1, 2
127  CALL cp_fm_create(matrix=rest_mo(j), &
128  matrix_struct=mo_mo_fmstruct, &
129  name="ET_rest_MO"//trim(adjustl(cp_to_string(j)))//"MATRIX")
130  END DO
131 
132 ! calculate MO-overlap
133 
134  CALL get_qs_env(qs_env, matrix_s=matrix_s)
135  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, qs_env%et_coupling%et_mo_coeff(i), &
136  tmp2, nmo, 1.0_dp, 0.0_dp)
137  CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
138  qs_env%mos(i)%mo_coeff, &
139  tmp2, 0.0_dp, smo)
140 
141 ! calculate the MO-representation of the restraint matrix A
142 
143  CALL cp_dbcsr_sm_fm_multiply(qs_env%et_coupling%rest_mat(1)%matrix, &
144  qs_env%et_coupling%et_mo_coeff(i), &
145  tmp2, nmo, 1.0_dp, 0.0_dp)
146 
147  CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
148  qs_env%mos(i)%mo_coeff, &
149  tmp2, 0.0_dp, rest_mo(1))
150 
151 ! calculate the MO-representation of the restraint matrix D
152 
153  CALL cp_dbcsr_sm_fm_multiply(qs_env%et_coupling%rest_mat(2)%matrix, &
154  qs_env%mos(i)%mo_coeff, &
155  tmp2, nmo, 1.0_dp, 0.0_dp)
156 
157  CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
158  qs_env%et_coupling%et_mo_coeff(i), &
159  tmp2, 0.0_dp, rest_mo(2))
160 ! TODO: could fix cp_fm_invert/LU determinant pivoting instead of calling cp_fm_det to save a pdgetrf call
161  CALL cp_fm_det(smo, s_det(i))
162  CALL cp_fm_invert(smo, inverse_mat)
163 
164  CALL cp_fm_get_info(inverse_mat, nrow_local=nrow_local, ncol_local=ncol_local, &
165  row_indices=row_indices, col_indices=col_indices)
166  b(i) = 0.0_dp
167 
168  DO j = 1, ncol_local
169  DO k = 1, nrow_local
170  b(i) = b(i) + rest_mo(2)%local_data(k, j)*inverse_mat%local_data(k, j)
171  END DO
172  END DO
173 
174  CALL cp_fm_transpose(inverse_mat, tinverse)
175  a(i) = 0.0_dp
176  DO j = 1, ncol_local
177  DO k = 1, nrow_local
178  a(i) = a(i) + rest_mo(1)%local_data(k, j)*tinverse%local_data(k, j)
179  END DO
180  END DO
181  IF (is_spin_constraint) THEN
182  a(i) = -a(i)
183  b(i) = -b(i)
184  END IF
185  CALL para_env%sum(a(i))
186 
187  CALL para_env%sum(b(i))
188 
189  CALL cp_fm_release(tmp2)
190  DO j = 1, 2
191  CALL cp_fm_release(rest_mo(j))
192  END DO
193  CALL cp_fm_release(smo)
194  CALL cp_fm_release(tinverse)
195  CALL cp_fm_release(inverse_mat)
196  END DO
197 
198 ! solve eigenstates for the projector matrix
199 
200  IF (dft_control%nspins == 2) THEN
201  sda = s_det(1)*s_det(2)
202  wda = ((a(1) + a(2)) + (b(1) + b(2)))*0.5_dp*sda
203  ELSE
204  sda = s_det(1)**2
205  wda = (a(1) + b(1))*sda
206  END IF
207 
208  IF (dft_control%qs_control%ddapc_restraint) THEN
209  waa = qs_env%et_coupling%order_p
210  wbb = dft_control%qs_control%ddapc_restraint_control(1)%ddapc_order_p
211  strength = dft_control%qs_control%ddapc_restraint_control(1)%strength
212  END IF
213 
214 !! construct S and W !!!
215  s_mat(1, 1) = 1.0_dp
216  s_mat(2, 2) = 1.0_dp
217  s_mat(2, 1) = sda
218  s_mat(1, 2) = sda
219 
220  w_mat(1, 1) = wbb
221  w_mat(2, 2) = waa
222  w_mat(2, 1) = wda
223  w_mat(1, 2) = wda
224 
225 !! solve WC=SCN
226  CALL diamat_all(s_mat, eigenv, .true.)
227  ! U = S**(-1/2)
228  u = 0.0_dp
229  u(1, 1) = 1.0_dp/sqrt(eigenv(1))
230  u(2, 2) = 1.0_dp/sqrt(eigenv(2))
231  tmp_mat = matmul(u, transpose(s_mat))
232  u = matmul(s_mat, tmp_mat)
233  tmp_mat = matmul(w_mat, u)
234  w_mat = matmul(u, tmp_mat)
235  CALL diamat_all(w_mat, eigenv, .true.)
236  tmp_mat = matmul(u, w_mat)
237 
238  CALL get_qs_env(qs_env, energy=energy)
239  w_mat(1, 1) = energy%total
240  w_mat(2, 2) = qs_env%et_coupling%energy
241  a(1) = (energy%total + strength*wbb)*sda - strength*wda
242  a(2) = (qs_env%et_coupling%energy + qs_env%et_coupling%e1*waa)*sda - qs_env%et_coupling%e1*wda
243  w_mat(1, 2) = (a(1) + a(2))*0.5_dp
244  w_mat(2, 1) = w_mat(1, 2)
245 
246  s_mat = matmul(w_mat, (tmp_mat))
247  w_mat = matmul(transpose(tmp_mat), s_mat)
248 
249  IF (iw > 0) THEN
250  WRITE (iw, *)
251  WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint A :', qs_env%et_coupling%e1
252  WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint B :', strength
253  WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint A:', waa
254  WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint B:', wbb
255  WRITE (iw, *)
256  WRITE (iw, '(T3,A,T60,(3X,F12.6))') &
257  'Diabatic electronic coupling matrix element(mHartree):', abs(w_mat(1, 2)*1000.0_dp)
258 
259  END IF
260 
261  CALL dbcsr_deallocate_matrix_set(qs_env%et_coupling%rest_mat)
262 
263  CALL cp_print_key_finished_output(iw, logger, et_coupling_section, &
264  "PROGRAM_RUN_INFO")
265  CALL timestop(handle)
266  END SUBROUTINE calc_et_coupling
267 
268 END MODULE et_coupling
269 
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_det
subroutine, public cp_fm_det(matrix_a, det_a)
Computes the determinant (with a correct sign even in parallel environment!) of a real square matrix.
Definition: cp_fm_basic_linalg.F:97
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition: cp_fm_basic_linalg.F:1710
cp_fm_basic_linalg::cp_fm_invert
subroutine, public cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix.
Definition: cp_fm_basic_linalg.F:2015
cp_fm_pool_types
pool for for elements that are retained and released
Definition: cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_get_el_struct
type(cp_fm_struct_type) function, pointer, public fm_pool_get_el_struct(pool)
returns the structure of the elements in this pool
Definition: cp_fm_pool_types.F:249
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
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
cp_fm_types::cp_fm_create
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
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)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
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,...
Definition: cp_output_handling.F:1013
et_coupling
calculates the electron transfer coupling elements Wu, Van Voorhis, JCP 125, 164105 (2006)
Definition: et_coupling.F:13
et_coupling::calc_et_coupling
subroutine, public calc_et_coupling(qs_env)
...
Definition: et_coupling.F:63
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:372
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition: qs_environment_types.F:502
qs_matrix_pools
wrapper for the pools of matrixes
Definition: qs_matrix_pools.F:14
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition: qs_matrix_pools.F:141
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397