Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the calculation of wannier states
10 : !> \author Alin M Elena
11 : ! **************************************************************************************************
12 : MODULE wannier_states
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind,&
15 : get_atomic_kind_set
16 : USE basis_set_types, ONLY: get_gto_basis_set,&
17 : gto_basis_set_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_type
19 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
20 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
21 : cp_fm_struct_release,&
22 : cp_fm_struct_type
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_get_element,&
25 : cp_fm_get_info,&
26 : cp_fm_get_submatrix,&
27 : cp_fm_release,&
28 : cp_fm_to_fm,&
29 : cp_fm_type
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_get_default_io_unit,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34 : cp_print_key_unit_nr
35 : USE cp_units, ONLY: cp_unit_from_cp2k
36 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
37 : section_vals_type,&
38 : section_vals_val_get
39 : USE kinds, ONLY: default_string_length,&
40 : dp
41 : USE message_passing, ONLY: mp_para_env_type
42 : USE orbital_pointers, ONLY: indco,&
43 : nco,&
44 : nso
45 : USE orbital_symbols, ONLY: cgf_symbol,&
46 : sgf_symbol
47 : USE orbital_transformation_matrices, ONLY: orbtramat
48 : USE parallel_gemm_api, ONLY: parallel_gemm
49 : USE particle_types, ONLY: particle_type
50 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
51 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_kind_types, ONLY: get_qs_kind,&
55 : get_qs_kind_set,&
56 : qs_kind_type
57 : USE wannier_states_types, ONLY: wannier_centres_type
58 : !!!! this ones are needed to mapping
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : ! *** Global parameters ***
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wannier_states'
68 :
69 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
70 :
71 : ! *** Public subroutines ***
72 :
73 : PUBLIC :: construct_wannier_states
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief constructs wannier states. mo_localized should not be overwritten!
79 : !> \param mo_localized ...
80 : !> \param Hks ...
81 : !> \param qs_env ...
82 : !> \param loc_print_section ...
83 : !> \param WannierCentres ...
84 : !> \param ns ...
85 : !> \param states ...
86 : !> \par History
87 : !> - Printout Wannier states in AO basis (11.09.2025, H. Elgabarty)
88 : ! **************************************************************************************************
89 16 : SUBROUTINE construct_wannier_states(mo_localized, &
90 : Hks, qs_env, loc_print_section, WannierCentres, ns, states)
91 :
92 : TYPE(cp_fm_type), INTENT(in) :: mo_localized
93 : TYPE(dbcsr_type), POINTER :: Hks
94 : TYPE(qs_environment_type), POINTER :: qs_env
95 : TYPE(section_vals_type), POINTER :: loc_print_section
96 : TYPE(wannier_centres_type), INTENT(INOUT) :: WannierCentres
97 : INTEGER, INTENT(IN) :: ns
98 : INTEGER, INTENT(IN), POINTER :: states(:)
99 :
100 : CHARACTER(len=*), PARAMETER :: routineN = 'construct_wannier_states'
101 :
102 : CHARACTER(default_string_length) :: unit_str
103 : CHARACTER(LEN=12) :: symbol
104 16 : CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
105 : CHARACTER(LEN=2) :: element_symbol
106 : CHARACTER(LEN=40) :: fmtstr1, fmtstr2, fmtstr3
107 16 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
108 : INTEGER :: after, before, from, handle, i, iatom, icgf, ico, icol, ikind, iproc, irow, iset, &
109 : isgf, ishell, iso, jcol, left, lmax, lshell, natom, ncgf, ncol, ncol_global, nrow_global, &
110 : nset, nsgf, nstates(2), output_unit, right, to, unit_mat
111 16 : INTEGER, DIMENSION(:), POINTER :: nshell
112 16 : INTEGER, DIMENSION(:, :), POINTER :: l
113 : LOGICAL :: print_cartesian
114 : REAL(KIND=dp) :: unit_conv
115 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
116 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
117 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
118 : TYPE(cp_fm_type) :: b, c, d
119 : TYPE(cp_logger_type), POINTER :: logger
120 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
121 : TYPE(mp_para_env_type), POINTER :: para_env
122 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
124 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(section_vals_type), POINTER :: print_key
126 :
127 : !-----------------------------------------------------------------------
128 : !-----------------------------------------------------------------------
129 :
130 16 : CALL timeset(routineN, handle)
131 16 : NULLIFY (logger, para_env)
132 :
133 : CALL get_qs_env(qs_env, para_env=para_env, &
134 : atomic_kind_set=atomic_kind_set, &
135 : qs_kind_set=qs_kind_set, &
136 16 : particle_set=particle_set)
137 :
138 16 : logger => cp_get_default_logger()
139 :
140 16 : output_unit = cp_logger_get_default_io_unit(logger)
141 :
142 : CALL cp_fm_get_info(mo_localized, &
143 : ncol_global=ncol_global, &
144 16 : nrow_global=nrow_global)
145 :
146 16 : nstates(1) = ns
147 16 : nstates(2) = para_env%mepos
148 16 : iproc = nstates(2)
149 16 : NULLIFY (fm_struct_tmp, print_key)
150 :
151 16 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CENTERS")
152 16 : CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str)
153 16 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
154 :
155 16 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES")
156 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, &
157 : ncol_global=1, &
158 : para_env=mo_localized%matrix_struct%para_env, &
159 16 : context=mo_localized%matrix_struct%context)
160 :
161 16 : CALL cp_fm_create(b, fm_struct_tmp, name="b")
162 16 : CALL cp_fm_create(c, fm_struct_tmp, name="c")
163 :
164 16 : CALL cp_fm_struct_release(fm_struct_tmp)
165 :
166 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=1, ncol_global=1, &
167 : para_env=mo_localized%matrix_struct%para_env, &
168 16 : context=mo_localized%matrix_struct%context)
169 :
170 16 : CALL cp_fm_create(d, fm_struct_tmp, name="d")
171 16 : CALL cp_fm_struct_release(fm_struct_tmp)
172 :
173 132 : WannierCentres%WannierHamDiag = 0.0_dp
174 :
175 : unit_mat = cp_print_key_unit_nr(logger, loc_print_section, &
176 : "WANNIER_STATES", extension=".whks", &
177 16 : ignore_should_output=.FALSE.)
178 16 : IF (unit_mat > 0) THEN
179 8 : WRITE (unit_mat, '(a16,1(i0,1x))') "Wannier states: ", ns
180 8 : WRITE (unit_mat, '(a16)') "#No x y z energy "
181 : END IF
182 132 : DO i = 1, ns
183 116 : CALL cp_fm_to_fm(mo_localized, b, 1, states(i), 1)
184 116 : CALL cp_dbcsr_sm_fm_multiply(Hks, b, c, 1)
185 : CALL parallel_gemm('T', 'N', 1, 1, nrow_global, 1.0_dp, &
186 116 : b, c, 0.0_dp, d)
187 116 : CALL cp_fm_get_element(d, 1, 1, WannierCentres%WannierHamDiag(i))
188 : ! if (iproc==para_env%mepos) WRITE(unit_mat,'(f16.8,2x)', advance='no')WannierCentres%WannierHamDiag(i)
189 174 : IF (unit_mat > 0) WRITE (unit_mat, '(i0,1x,4(f16.8,2x))') states(i), &
190 306 : WannierCentres%centres(1:3, states(i))*unit_conv, WannierCentres%WannierHamDiag(states(i))
191 : END DO
192 :
193 16 : IF (unit_mat > 0) WRITE (unit_mat, *)
194 :
195 16 : IF (output_unit > 0) THEN
196 8 : WRITE (output_unit, *) ""
197 8 : WRITE (output_unit, *) "NUMBER OF Wannier STATES ", ns
198 8 : WRITE (output_unit, *) "ENERGY original MO-index"
199 66 : DO i = 1, ns
200 66 : WRITE (output_unit, '(f16.8,2x,i0)') WannierCentres%WannierHamDiag(i), states(i)
201 : END DO
202 : END IF
203 :
204 16 : CALL cp_fm_release(b)
205 16 : CALL cp_fm_release(c)
206 16 : CALL cp_fm_release(d)
207 :
208 : ! Print the states in AO basis
209 16 : CALL section_vals_val_get(print_key, "CARTESIAN", l_val=print_cartesian)
210 :
211 64 : ALLOCATE (smatrix(nrow_global, ncol_global))
212 16 : CALL cp_fm_get_submatrix(mo_localized, smatrix(1:nrow_global, 1:ncol_global))
213 :
214 16 : IF (unit_mat > 0) THEN
215 :
216 8 : NULLIFY (nshell)
217 8 : NULLIFY (bsgf_symbol)
218 8 : NULLIFY (l)
219 : NULLIFY (bsgf_symbol)
220 :
221 8 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
222 8 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
223 :
224 : ! Print header, define column widths and string templates
225 8 : after = 6
226 8 : before = 4
227 8 : ncol = INT(56/(before + after + 3))
228 :
229 8 : fmtstr1 = "(T2,A,21X, ( X,I5, X))"
230 8 : fmtstr2 = "(T2,A,9X, (1X,F . ))"
231 8 : fmtstr3 = "(T2,A,I5,1X,I5,1X,A,1X,A6, (1X,F . ))"
232 :
233 8 : right = MAX((after - 2), 1)
234 8 : left = (before + after + 3) - right - 5
235 :
236 8 : IF (print_cartesian) THEN
237 0 : WRITE (UNIT=unit_mat, FMT="(T2,A,16X,A)") "WS|", "Wannier states in the cartesian AO basis"
238 : ELSE
239 8 : WRITE (UNIT=unit_mat, FMT="(T2,A,16X,A)") "WS|", "Wannier states in the spherical AO basis"
240 : END IF
241 8 : WRITE (UNIT=fmtstr1(11:12), FMT="(I2)") ncol
242 8 : WRITE (UNIT=fmtstr1(14:15), FMT="(I2)") left
243 8 : WRITE (UNIT=fmtstr1(21:22), FMT="(I2)") right
244 :
245 8 : WRITE (UNIT=fmtstr2(10:11), FMT="(I2)") ncol
246 8 : WRITE (UNIT=fmtstr2(17:18), FMT="(I2)") before + after + 2
247 8 : WRITE (UNIT=fmtstr2(20:21), FMT="(I2)") after
248 :
249 8 : WRITE (UNIT=fmtstr3(27:28), FMT="(I2)") ncol
250 8 : WRITE (UNIT=fmtstr3(34:35), FMT="(I2)") before + after + 2
251 8 : WRITE (UNIT=fmtstr3(37:38), FMT="(I2)") after
252 :
253 : ! get MO coefficients in terms of contracted cartesian functions
254 8 : IF (print_cartesian) THEN
255 :
256 0 : ALLOCATE (cmatrix(ncgf, ncgf))
257 0 : cmatrix = 0.0_dp
258 :
259 : ! Transform spherical to Cartesian AO basis
260 0 : icgf = 1
261 0 : isgf = 1
262 0 : DO iatom = 1, natom
263 0 : NULLIFY (orb_basis_set, dftb_parameter)
264 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
265 : CALL get_qs_kind(qs_kind_set(ikind), &
266 : basis_set=orb_basis_set, &
267 0 : dftb_parameter=dftb_parameter)
268 0 : IF (ASSOCIATED(orb_basis_set)) THEN
269 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
270 : nset=nset, &
271 : nshell=nshell, &
272 0 : l=l)
273 0 : DO iset = 1, nset
274 0 : DO ishell = 1, nshell(iset)
275 0 : lshell = l(ishell, iset)
276 : CALL dgemm("T", "N", nco(lshell), ncol_global, nso(lshell), 1.0_dp, &
277 : orbtramat(lshell)%c2s, nso(lshell), &
278 : smatrix(isgf, 1), nsgf, 0.0_dp, &
279 0 : cmatrix(icgf, 1), ncgf)
280 0 : icgf = icgf + nco(lshell)
281 0 : isgf = isgf + nso(lshell)
282 : END DO
283 : END DO
284 0 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
285 0 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
286 0 : DO ishell = 1, lmax + 1
287 0 : lshell = ishell - 1
288 : CALL dgemm("T", "N", nco(lshell), nsgf, nso(lshell), 1.0_dp, &
289 : orbtramat(lshell)%c2s, nso(lshell), &
290 : smatrix(isgf, 1), nsgf, 0.0_dp, &
291 0 : cmatrix(icgf, 1), ncgf)
292 0 : icgf = icgf + nco(lshell)
293 0 : isgf = isgf + nso(lshell)
294 : END DO
295 : ELSE
296 : ! assume atom without basis set
297 : END IF
298 : END DO ! iatom
299 :
300 : END IF
301 :
302 : ! Print to file
303 23 : DO icol = 1, ncol_global, ncol
304 :
305 15 : from = icol
306 15 : to = MIN((from + ncol - 1), ncol_global)
307 :
308 15 : WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|"
309 73 : WRITE (UNIT=unit_mat, FMT=fmtstr1) "WS|", (jcol, jcol=from, to)
310 15 : WRITE (UNIT=unit_mat, FMT=fmtstr2) "WS| Energies", &
311 88 : (WannierCentres%WannierHamDiag(states(jcol)), jcol=from, to)
312 15 : WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|"
313 :
314 15 : irow = 1
315 :
316 110 : DO iatom = 1, natom
317 :
318 87 : IF (iatom /= 1) WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|"
319 :
320 87 : NULLIFY (orb_basis_set, dftb_parameter)
321 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
322 87 : element_symbol=element_symbol, kind_number=ikind)
323 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
324 87 : dftb_parameter=dftb_parameter)
325 :
326 189 : IF (print_cartesian) THEN
327 :
328 0 : NULLIFY (bcgf_symbol)
329 0 : IF (ASSOCIATED(orb_basis_set)) THEN
330 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
331 : nset=nset, &
332 : nshell=nshell, &
333 : l=l, &
334 0 : cgf_symbol=bcgf_symbol)
335 :
336 0 : icgf = 1
337 0 : DO iset = 1, nset
338 0 : DO ishell = 1, nshell(iset)
339 0 : lshell = l(ishell, iset)
340 0 : DO ico = 1, nco(lshell)
341 : WRITE (UNIT=unit_mat, FMT=fmtstr3) &
342 0 : "WS|", irow, iatom, ADJUSTR(element_symbol), bcgf_symbol(icgf), &
343 0 : (cmatrix(irow, jcol), jcol=from, to)
344 0 : icgf = icgf + 1
345 0 : irow = irow + 1
346 : END DO
347 : END DO
348 : END DO
349 0 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
350 0 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
351 0 : icgf = 1
352 0 : DO ishell = 1, lmax + 1
353 0 : lshell = ishell - 1
354 0 : DO ico = 1, nco(lshell)
355 0 : symbol = cgf_symbol(1, indco(1:3, icgf))
356 0 : symbol(1:2) = " "
357 : WRITE (UNIT=unit_mat, FMT=fmtstr3) &
358 0 : "WS|", irow, iatom, ADJUSTR(element_symbol), symbol, &
359 0 : (cmatrix(irow, jcol), jcol=from, to)
360 0 : icgf = icgf + 1
361 0 : irow = irow + 1
362 : END DO
363 : END DO
364 : ELSE
365 : ! assume atom without basis set
366 : END IF
367 :
368 : ELSE !print in spherical AO basis
369 :
370 87 : IF (ASSOCIATED(orb_basis_set)) THEN
371 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
372 : nset=nset, &
373 : nshell=nshell, &
374 : l=l, &
375 87 : sgf_symbol=bsgf_symbol)
376 87 : isgf = 1
377 255 : DO iset = 1, nset
378 525 : DO ishell = 1, nshell(iset)
379 270 : lshell = l(ishell, iset)
380 970 : DO iso = 1, nso(lshell)
381 : WRITE (UNIT=unit_mat, FMT=fmtstr3) &
382 532 : "WS|", irow, iatom, ADJUSTR(element_symbol), bsgf_symbol(isgf), &
383 1064 : (smatrix(irow, jcol), jcol=from, to)
384 532 : isgf = isgf + 1
385 802 : irow = irow + 1
386 : END DO
387 : END DO
388 : END DO
389 0 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
390 0 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
391 0 : isgf = 1
392 0 : DO ishell = 1, lmax + 1
393 0 : lshell = ishell - 1
394 0 : DO iso = 1, nso(lshell)
395 0 : symbol = sgf_symbol(1, lshell, -lshell + iso - 1)
396 0 : symbol(1:2) = " "
397 : WRITE (UNIT=unit_mat, FMT=fmtstr3) &
398 0 : "WS|", irow, iatom, ADJUSTR(element_symbol), symbol, &
399 0 : (smatrix(irow, jcol), jcol=from, to)
400 0 : isgf = isgf + 1
401 0 : irow = irow + 1
402 : END DO
403 : END DO
404 : ELSE
405 : ! assume atom without basis set
406 : END IF
407 :
408 : END IF ! print cartesian
409 :
410 : END DO ! iatom
411 :
412 : END DO ! icol
413 :
414 8 : WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|"
415 :
416 8 : IF (print_cartesian) THEN
417 0 : DEALLOCATE (cmatrix)
418 : END IF
419 :
420 : END IF ! output Wannier states in AO
421 16 : DEALLOCATE (smatrix)
422 :
423 : CALL cp_print_key_finished_output(unit_mat, logger, loc_print_section, &
424 16 : "WANNIER_STATES")
425 :
426 16 : CALL timestop(handle)
427 96 : END SUBROUTINE construct_wannier_states
428 :
429 : END MODULE wannier_states
430 :
|
