(git:374b731)
Loading...
Searching...
No Matches
qs_tddfpt2_assign.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
8MODULE qs_tddfpt2_assign
9 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
10 USE cp_fm_basic_linalg, ONLY: cp_fm_geadd,&
11 cp_fm_scale
12 USE cp_fm_diag, ONLY: cp_fm_power
13 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
14 cp_fm_struct_release,&
15 cp_fm_struct_type
16 USE cp_fm_types, ONLY: cp_fm_create,&
17 cp_fm_get_diag,&
18 cp_fm_get_info,&
19 cp_fm_release,&
20 cp_fm_to_fm,&
21 cp_fm_type
22 USE dbcsr_api, ONLY: dbcsr_p_type,&
23 dbcsr_type
24 USE exstates_types, ONLY: wfn_history_type
25 USE kinds, ONLY: dp
26 USE message_passing, ONLY: mp_para_env_type
27 USE parallel_gemm_api, ONLY: parallel_gemm
28 USE qs_environment_types, ONLY: get_qs_env,&
29 qs_environment_type
30#include "./base/base_uses.f90"
31
32 IMPLICIT NONE
33
34 PRIVATE
35
36 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_assign'
37
38 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
39
40 PUBLIC :: assign_state
41
42! **************************************************************************************************
43
44CONTAINS
45
46! **************************************************************************************************
47!> \brief ...
48!> \param qs_env ...
49!> \param matrix_s ...
50!> \param evects ...
51!> \param psi0 ...
52!> \param wfn_history ...
53!> \param my_state ...
54! **************************************************************************************************
55 SUBROUTINE assign_state(qs_env, matrix_s, evects, psi0, wfn_history, my_state)
56 TYPE(qs_environment_type), POINTER :: qs_env
57 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
58 TYPE(cp_fm_type), DIMENSION(:, :) :: evects
59 TYPE(cp_fm_type), DIMENSION(:) :: psi0
60 TYPE(wfn_history_type) :: wfn_history
61 INTEGER, INTENT(INOUT) :: my_state
62
63 CHARACTER(LEN=*), PARAMETER :: routinen = 'assign_state'
64
65 INTEGER :: handle, is, ispin, natom, ncol, nspins, &
66 nstate
67 REAL(kind=dp) :: xsum
68 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dv, rdiag
69 TYPE(dbcsr_type), POINTER :: smat
70 TYPE(mp_para_env_type), POINTER :: para_env
71
72 CALL timeset(routinen, handle)
73
74 CALL get_qs_env(qs_env, natom=natom, para_env=para_env)
75 nspins = SIZE(psi0)
76 nstate = SIZE(evects, 2)
77 !
78 smat => matrix_s(1)%matrix
79 !
80 IF (ASSOCIATED(wfn_history%evect)) THEN
81 ALLOCATE (dv(nstate))
82 !
83 wfn_history%gsval = 0.0_dp
84 wfn_history%gsmin = 1.0_dp
85 DO ispin = 1, nspins
86 CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
87 CALL lowdin_orthogonalization(wfn_history%cpmos(ispin), wfn_history%evect(ispin), &
88 ncol, smat)
89 ALLOCATE (rdiag(ncol))
90 CALL wfn_align(psi0(ispin), wfn_history%cpmos(ispin), wfn_history%evect(ispin), &
91 rdiag, smat)
92 wfn_history%gsval = wfn_history%gsval + sum(rdiag)/real(ncol*nspins, kind=dp)
93 wfn_history%gsmin = min(wfn_history%gsmin, minval(rdiag))
94 DEALLOCATE (rdiag)
95 END DO
96 DO is = 1, nstate
97 xsum = 0.0_dp
98 DO ispin = 1, nspins
99 CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
100 ALLOCATE (rdiag(ncol))
101 CALL xvec_ovlp(evects(ispin, is), wfn_history%evect(ispin), rdiag, smat)
102 xsum = xsum + sum(rdiag)
103 DEALLOCATE (rdiag)
104 END DO
105 dv(is) = abs(xsum)/sqrt(real(nspins, dp))
106 END DO
107 my_state = maxval(maxloc(dv))
108 wfn_history%xsval = dv(my_state)
109 IF (wfn_history%xsval < 0.75_dp) THEN
110 dv(my_state) = 0.0_dp
111 IF (wfn_history%xsval/maxval(dv) < 0.5_dp) THEN
112 CALL cp_warn(__location__, "Uncertain assignment for State following."// &
113 " Reduce trust radius in Geometry Optimization or timestep"// &
114 " in MD runs.")
115 END IF
116 END IF
117 DO ispin = 1, nspins
118 CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
119 CALL cp_fm_to_fm(evects(ispin, my_state), wfn_history%evect(ispin))
120 CALL cp_fm_to_fm(psi0(ispin), wfn_history%cpmos(ispin), ncol, 1, 1)
121 END DO
122 !
123 DEALLOCATE (dv)
124 ELSE
125 !
126 ALLOCATE (wfn_history%evect(nspins))
127 ALLOCATE (wfn_history%cpmos(nspins))
128 DO ispin = 1, nspins
129 CALL cp_fm_create(wfn_history%evect(ispin), evects(ispin, 1)%matrix_struct, "Xvec")
130 CALL cp_fm_create(wfn_history%cpmos(ispin), evects(ispin, 1)%matrix_struct, "Cvec")
131 END DO
132 DO ispin = 1, nspins
133 CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
134 CALL cp_fm_to_fm(evects(ispin, my_state), wfn_history%evect(ispin))
135 CALL cp_fm_to_fm(psi0(ispin), wfn_history%cpmos(ispin), ncol, 1, 1)
136 END DO
137 wfn_history%xsval = 1.0_dp
138 wfn_history%gsval = 1.0_dp
139 wfn_history%gsmin = 1.0_dp
140 END IF
141
142 CALL timestop(handle)
143
144 END SUBROUTINE assign_state
145
146! **************************************************************************************************
147!> \brief return a set of S orthonormal vectors (C^T S C == 1) where
148!> a Lodwin transformation is applied to keep the rotated vectors as close
149!> as possible to the original ones
150!> \param vmatrix ...
151!> \param xmatrix ...
152!> \param ncol ...
153!> \param matrix_s ...
154!> \param
155!> \par History
156!> 05.2009 created [MI]
157!> 06.2023 adapted to include a second set of vectors [JGH]
158!> \note
159! **************************************************************************************************
160 SUBROUTINE lowdin_orthogonalization(vmatrix, xmatrix, ncol, matrix_s)
161
162 TYPE(cp_fm_type), INTENT(IN) :: vmatrix, xmatrix
163 INTEGER, INTENT(IN) :: ncol
164 TYPE(dbcsr_type) :: matrix_s
165
166 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin_orthogonalization'
167 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
168
169 INTEGER :: handle, n, ncol_global, ndep
170 REAL(dp) :: threshold, xsum
171 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rdiag
172 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
173 TYPE(cp_fm_type) :: csc, sc, work
174
175 IF (ncol .EQ. 0) RETURN
176
177 CALL timeset(routinen, handle)
178
179 threshold = 1.0e-7_dp
180 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
181 IF (ncol .GT. ncol_global) cpabort("Wrong ncol value")
182
183 CALL cp_fm_create(sc, xmatrix%matrix_struct, "SC")
184 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
185
186 NULLIFY (fm_struct_tmp)
187 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
188 para_env=vmatrix%matrix_struct%para_env, &
189 context=vmatrix%matrix_struct%context)
190 CALL cp_fm_create(csc, fm_struct_tmp, "csc")
191 CALL cp_fm_create(work, fm_struct_tmp, "work")
192 CALL cp_fm_struct_release(fm_struct_tmp)
193
194 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
195 CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
196 CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
197 CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
198 !
199 CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, xmatrix, csc, rzero, sc)
200 CALL cp_fm_to_fm(sc, xmatrix, ncol, 1, 1)
201 ! projecton for xSv = 0
202 CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
203 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
204 CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
205 CALL cp_fm_geadd(-1.0_dp, 'N', sc, 1.0_dp, xmatrix)
206 ! normalisation
207 CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
208 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, xmatrix, sc, rzero, csc)
209 ALLOCATE (rdiag(ncol))
210 CALL cp_fm_get_diag(csc, rdiag)
211 xsum = sum(rdiag)
212 DEALLOCATE (rdiag)
213 xsum = 1._dp/sqrt(xsum)
214 CALL cp_fm_scale(xsum, xmatrix)
215
216 CALL cp_fm_release(csc)
217 CALL cp_fm_release(sc)
218 CALL cp_fm_release(work)
219
220 CALL timestop(handle)
221
222 END SUBROUTINE lowdin_orthogonalization
223
224! **************************************************************************************************
225!> \brief ...
226!> \param gmatrix ...
227!> \param vmatrix ...
228!> \param xmatrix ...
229!> \param rdiag ...
230!> \param matrix_s ...
231! **************************************************************************************************
232 SUBROUTINE wfn_align(gmatrix, vmatrix, xmatrix, rdiag, matrix_s)
233
234 TYPE(cp_fm_type), INTENT(IN) :: gmatrix, vmatrix, xmatrix
235 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: rdiag
236 TYPE(dbcsr_type) :: matrix_s
237
238 CHARACTER(LEN=*), PARAMETER :: routinen = 'wfn_align'
239 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
240
241 INTEGER :: handle, n, ncol, ncol_global
242 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
243 TYPE(cp_fm_type) :: csc, sc
244
245 CALL timeset(routinen, handle)
246
247 ncol = SIZE(rdiag)
248 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
249 IF (ncol .GT. ncol_global) cpabort("Wrong ncol value")
250
251 CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
252 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
253
254 NULLIFY (fm_struct_tmp)
255 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
256 para_env=vmatrix%matrix_struct%para_env, &
257 context=vmatrix%matrix_struct%context)
258 CALL cp_fm_create(csc, fm_struct_tmp, "csc")
259 CALL cp_fm_struct_release(fm_struct_tmp)
260
261 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, gmatrix, sc, rzero, csc)
262 CALL parallel_gemm('N', 'T', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
263 CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
264 CALL parallel_gemm('N', 'T', n, ncol, ncol, rone, xmatrix, csc, rzero, sc)
265 CALL cp_fm_to_fm(sc, xmatrix, ncol, 1, 1)
266 !
267 CALL lowdin_orthogonalization(vmatrix, xmatrix, ncol, matrix_s)
268 !
269 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
270 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, gmatrix, sc, rzero, csc)
271 CALL cp_fm_get_diag(csc, rdiag)
272
273 CALL cp_fm_release(csc)
274 CALL cp_fm_release(sc)
275
276 CALL timestop(handle)
277
278 END SUBROUTINE wfn_align
279
280! **************************************************************************************************
281!> \brief ...
282!> \param ematrix ...
283!> \param xmatrix ...
284!> \param rdiag ...
285!> \param matrix_s ...
286! **************************************************************************************************
287 SUBROUTINE xvec_ovlp(ematrix, xmatrix, rdiag, matrix_s)
288
289 TYPE(cp_fm_type), INTENT(IN) :: ematrix, xmatrix
290 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: rdiag
291 TYPE(dbcsr_type) :: matrix_s
292
293 CHARACTER(LEN=*), PARAMETER :: routinen = 'xvec_ovlp'
294 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
295
296 INTEGER :: handle, n, ncol, ncol_global
297 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
298 TYPE(cp_fm_type) :: csc, sc
299
300 CALL timeset(routinen, handle)
301
302 ncol = SIZE(rdiag)
303 CALL cp_fm_get_info(matrix=xmatrix, nrow_global=n, ncol_global=ncol_global)
304 IF (ncol .GT. ncol_global) cpabort("Wrong ncol value")
305
306 CALL cp_fm_create(sc, xmatrix%matrix_struct, "SC")
307 CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
308
309 NULLIFY (fm_struct_tmp)
310 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
311 para_env=xmatrix%matrix_struct%para_env, &
312 context=xmatrix%matrix_struct%context)
313 CALL cp_fm_create(csc, fm_struct_tmp, "csc")
314 CALL cp_fm_struct_release(fm_struct_tmp)
315
316 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, ematrix, sc, rzero, csc)
317 CALL cp_fm_get_diag(csc, rdiag)
318
319 CALL cp_fm_release(csc)
320 CALL cp_fm_release(sc)
321
322 CALL timestop(handle)
323
324 END SUBROUTINE xvec_ovlp
325
326END MODULE qs_tddfpt2_assign
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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_geadd
subroutine, public cp_fm_geadd(alpha, trans, matrix_a, beta, matrix_b)
interface to BLACS geadd: matrix_b = beta*matrix_b + alpha*opt(matrix_a) where opt(matrix_a) can be e...
Definition cp_fm_basic_linalg.F:274
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1680
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:896
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:570
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
exstates_types
Types for excited states potential energies.
Definition exstates_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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_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_tddfpt2_assign
Definition qs_tddfpt2_assign.F:8
qs_tddfpt2_assign::assign_state
subroutine, public assign_state(qs_env, matrix_s, evects, psi0, wfn_history, my_state)
...
Definition qs_tddfpt2_assign.F:56
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
exstates_types::wfn_history_type
Definition exstates_types.F:40
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215