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
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_api, ONLY: dbcsr_p_type
16 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
17 : dbcsr_deallocate_matrix_set
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_det,&
19 : cp_fm_invert,&
20 : cp_fm_transpose
21 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
22 : fm_pool_get_el_struct
23 : USE cp_fm_struct, ONLY: cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_get_info,&
26 : cp_fm_release,&
27 : cp_fm_type
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type,&
30 : cp_to_string
31 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32 : cp_print_key_unit_nr
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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
71 : LOGICAL :: is_spin_constraint
72 : REAL(KIND=dp) :: Sda, strength, Waa, Wbb, Wda
73 0 : 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 0 : 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 0 : TYPE(cp_fm_type), DIMENSION(2) :: rest_MO
80 : TYPE(cp_logger_type), POINTER :: logger
81 0 : 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 0 : NULLIFY (mo_mo_fmstruct, energy, matrix_s, dft_control, para_env)
88 :
89 0 : CALL timeset(routineN, handle)
90 :
91 0 : logger => cp_get_default_logger()
92 : et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
93 0 : "PROPERTIES%ET_COUPLING")
94 :
95 0 : 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 0 : extension=".log")
99 :
100 0 : is_spin_constraint = .FALSE.
101 0 : ALLOCATE (a(dft_control%nspins))
102 0 : ALLOCATE (b(dft_control%nspins))
103 0 : ALLOCATE (S_det(dft_control%nspins))
104 :
105 0 : CALL mpools_get(qs_env%mpools, mo_mo_fm_pools=mo_mo_fm_pools)
106 0 : mo_mo_fmstruct => fm_pool_get_el_struct(mo_mo_fm_pools(1)%pool)
107 0 : DO i = 1, dft_control%nspins
108 0 : 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 0 : nmo=nmo)
113 :
114 : CALL cp_fm_create(matrix=tmp2, &
115 : matrix_struct=qs_env%mos(i)%mo_coeff%matrix_struct, &
116 0 : 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 0 : name="INVERSE"//TRIM(ADJUSTL(cp_to_string(2)))//"MATRIX")
120 : CALL cp_fm_create(matrix=Tinverse, &
121 : matrix_struct=mo_mo_fmstruct, &
122 0 : name="T_INVERSE"//TRIM(ADJUSTL(cp_to_string(2)))//"MATRIX")
123 : CALL cp_fm_create(matrix=SMO, &
124 : matrix_struct=mo_mo_fmstruct, &
125 0 : name="ET_SMO"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
126 0 : DO j = 1, 2
127 : CALL cp_fm_create(matrix=rest_MO(j), &
128 : matrix_struct=mo_mo_fmstruct, &
129 0 : name="ET_rest_MO"//TRIM(ADJUSTL(cp_to_string(j)))//"MATRIX")
130 : END DO
131 :
132 : ! calculate MO-overlap
133 :
134 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : CALL cp_fm_det(SMO, S_det(i))
162 0 : 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 0 : row_indices=row_indices, col_indices=col_indices)
166 0 : b(i) = 0.0_dp
167 :
168 0 : DO j = 1, ncol_local
169 0 : DO k = 1, nrow_local
170 0 : 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 0 : CALL cp_fm_transpose(inverse_mat, Tinverse)
175 0 : a(i) = 0.0_dp
176 0 : DO j = 1, ncol_local
177 0 : DO k = 1, nrow_local
178 0 : 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 0 : CALL para_env%sum(a(i))
186 :
187 0 : CALL para_env%sum(b(i))
188 :
189 0 : CALL cp_fm_release(tmp2)
190 0 : DO j = 1, 2
191 0 : CALL cp_fm_release(rest_MO(j))
192 : END DO
193 0 : CALL cp_fm_release(SMO)
194 0 : CALL cp_fm_release(Tinverse)
195 0 : CALL cp_fm_release(inverse_mat)
196 : END DO
197 :
198 : ! solve eigenstates for the projector matrix
199 :
200 0 : IF (dft_control%nspins == 2) THEN
201 0 : Sda = S_det(1)*S_det(2)
202 0 : Wda = ((a(1) + a(2)) + (b(1) + b(2)))*0.5_dp*Sda
203 : ELSE
204 0 : Sda = S_det(1)**2
205 0 : Wda = (a(1) + b(1))*Sda
206 : END IF
207 :
208 0 : IF (dft_control%qs_control%ddapc_restraint) THEN
209 0 : Waa = qs_env%et_coupling%order_p
210 0 : Wbb = dft_control%qs_control%ddapc_restraint_control(1)%ddapc_order_p
211 0 : strength = dft_control%qs_control%ddapc_restraint_control(1)%strength
212 : END IF
213 :
214 : !! construct S and W !!!
215 0 : S_mat(1, 1) = 1.0_dp
216 0 : S_mat(2, 2) = 1.0_dp
217 0 : S_mat(2, 1) = Sda
218 0 : S_mat(1, 2) = Sda
219 :
220 0 : W_mat(1, 1) = Wbb
221 0 : W_mat(2, 2) = Waa
222 0 : W_mat(2, 1) = Wda
223 0 : W_mat(1, 2) = Wda
224 :
225 : !! solve WC=SCN
226 0 : CALL diamat_all(S_mat, eigenv, .TRUE.)
227 : ! U = S**(-1/2)
228 0 : U = 0.0_dp
229 0 : U(1, 1) = 1.0_dp/SQRT(eigenv(1))
230 0 : U(2, 2) = 1.0_dp/SQRT(eigenv(2))
231 0 : tmp_mat = MATMUL(U, TRANSPOSE(S_mat))
232 0 : U = MATMUL(S_mat, tmp_mat)
233 0 : tmp_mat = MATMUL(W_mat, U)
234 0 : W_mat = MATMUL(U, tmp_mat)
235 0 : CALL diamat_all(W_mat, eigenv, .TRUE.)
236 0 : tmp_mat = MATMUL(U, W_mat)
237 :
238 0 : CALL get_qs_env(qs_env, energy=energy)
239 0 : W_mat(1, 1) = energy%total
240 0 : W_mat(2, 2) = qs_env%et_coupling%energy
241 0 : a(1) = (energy%total + strength*Wbb)*Sda - strength*Wda
242 0 : a(2) = (qs_env%et_coupling%energy + qs_env%et_coupling%e1*Waa)*Sda - qs_env%et_coupling%e1*Wda
243 0 : W_mat(1, 2) = (a(1) + a(2))*0.5_dp
244 0 : W_mat(2, 1) = W_mat(1, 2)
245 :
246 0 : S_mat = MATMUL(W_mat, (tmp_mat))
247 0 : W_mat = MATMUL(TRANSPOSE(tmp_mat), S_mat)
248 :
249 0 : IF (iw > 0) THEN
250 0 : WRITE (iw, *)
251 0 : WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint A :', qs_env%et_coupling%e1
252 0 : WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint B :', strength
253 0 : WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint A:', Waa
254 0 : WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint B:', Wbb
255 0 : WRITE (iw, *)
256 : WRITE (iw, '(T3,A,T60,(3X,F12.6))') &
257 0 : 'Diabatic electronic coupling matrix element(mHartree):', ABS(W_mat(1, 2)*1000.0_dp)
258 :
259 : END IF
260 :
261 0 : 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 0 : "PROGRAM_RUN_INFO")
265 0 : CALL timestop(handle)
266 0 : END SUBROUTINE calc_et_coupling
267 :
268 : END MODULE et_coupling
269 :
|
