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 Calculation of band structures
10 : !> \par History
11 : !> 2015.06 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_band_structure
15 : USE cell_types, ONLY: cell_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_files, ONLY: close_file,&
19 : open_file
20 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
21 : USE cp_parser_methods, ONLY: read_float_object
22 : USE dbcsr_api, ONLY: dbcsr_p_type
23 : USE input_section_types, ONLY: section_vals_get,&
24 : section_vals_get_subs_vals,&
25 : section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: default_string_length,&
28 : dp,&
29 : max_line_length
30 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
31 : kpoint_init_cell_index,&
32 : kpoint_initialize_mo_set,&
33 : kpoint_initialize_mos
34 : USE kpoint_types, ONLY: get_kpoint_info,&
35 : kpoint_create,&
36 : kpoint_env_type,&
37 : kpoint_release,&
38 : kpoint_sym_create,&
39 : kpoint_type
40 : USE machine, ONLY: m_walltime
41 : USE mathconstants, ONLY: twopi
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE physcon, ONLY: angstrom,&
44 : evolt
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_env_release,&
47 : qs_environment_type
48 : USE qs_gamma2kp, ONLY: create_kp_from_gamma
49 : USE qs_mo_types, ONLY: get_mo_set
50 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
51 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
52 : USE qs_scf_types, ONLY: qs_scf_env_type
53 : USE scf_control_types, ONLY: scf_control_type
54 : USE string_utilities, ONLY: uppercase
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_band_structure'
62 :
63 : PUBLIC :: calculate_band_structure, calculate_kp_orbitals, calculate_kpoints_for_bs
64 :
65 : ! **************************************************************************************************
66 :
67 : CONTAINS
68 :
69 : ! **************************************************************************************************
70 : !> \brief Main routine for band structure calculation
71 : !> \param qs_env ...
72 : !> \author JGH
73 : ! **************************************************************************************************
74 35838 : SUBROUTINE calculate_band_structure(qs_env)
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 :
77 : LOGICAL :: do_kpoints, explicit
78 : TYPE(section_vals_type), POINTER :: bs_input
79 :
80 17919 : bs_input => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%BAND_STRUCTURE")
81 17919 : CALL section_vals_get(bs_input, explicit=explicit)
82 17919 : IF (explicit) THEN
83 8 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
84 8 : IF (do_kpoints) THEN
85 6 : CALL do_calculate_band_structure(qs_env)
86 : ELSE
87 : BLOCK
88 : TYPE(qs_environment_type), POINTER :: qs_env_kp
89 2 : CALL create_kp_from_gamma(qs_env, qs_env_kp)
90 2 : CALL do_calculate_band_structure(qs_env_kp)
91 2 : CALL qs_env_release(qs_env_kp)
92 2 : DEALLOCATE (qs_env_kp)
93 : END BLOCK
94 : END IF
95 : END IF
96 :
97 17919 : END SUBROUTINE calculate_band_structure
98 :
99 : ! **************************************************************************************************
100 : !> \brief band structure calculation
101 : !> \param qs_env ...
102 : !> \author JGH
103 : ! **************************************************************************************************
104 16 : SUBROUTINE do_calculate_band_structure(qs_env)
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 :
107 : CHARACTER(LEN=default_string_length) :: filename, ustr
108 : CHARACTER(LEN=default_string_length), &
109 8 : DIMENSION(:), POINTER :: spname, strptr
110 : CHARACTER(LEN=max_line_length) :: error_message
111 : INTEGER :: bs_data_unit, i, i_rep, ik, ikk, ikpgr, &
112 : imo, ip, ispin, n_ptr, n_rep, nadd, &
113 : nkp, nmo, npline, npoints, nspins, &
114 : unit_nr
115 : INTEGER, DIMENSION(2) :: kp_range
116 : LOGICAL :: explicit, io_default, my_kpgrp
117 : REAL(KIND=dp) :: t1, t2
118 : REAL(KIND=dp), DIMENSION(3) :: kpptr
119 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, eigval, occnum, &
120 8 : occupation_numbers, wkp
121 8 : REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral, kspecial, xkp
122 : TYPE(cell_type), POINTER :: cell
123 : TYPE(dft_control_type), POINTER :: dft_control
124 : TYPE(kpoint_env_type), POINTER :: kp
125 : TYPE(kpoint_type), POINTER :: kpoint
126 : TYPE(mp_para_env_type), POINTER :: para_env
127 : TYPE(section_vals_type), POINTER :: bs_input, kpset
128 :
129 16 : bs_input => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%BAND_STRUCTURE")
130 8 : CALL section_vals_get(bs_input, explicit=explicit)
131 8 : IF (explicit) THEN
132 8 : CALL section_vals_val_get(bs_input, "FILE_NAME", c_val=filename)
133 8 : CALL section_vals_val_get(bs_input, "ADDED_MOS", i_val=nadd)
134 8 : unit_nr = cp_logger_get_default_io_unit()
135 8 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
136 8 : CALL get_qs_env(qs_env, cell=cell)
137 8 : kpset => section_vals_get_subs_vals(bs_input, "KPOINT_SET")
138 8 : CALL section_vals_get(kpset, n_repetition=n_rep)
139 8 : IF (unit_nr > 0) THEN
140 4 : WRITE (unit_nr, FMT="(/,T2,A)") "KPOINTS| Band Structure Calculation"
141 4 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of k-point sets", n_rep
142 4 : IF (nadd /= 0) THEN
143 1 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of added MOs/bands", nadd
144 : END IF
145 : END IF
146 8 : IF (filename == "") THEN
147 : ! use standard output file
148 2 : bs_data_unit = unit_nr
149 2 : io_default = .TRUE.
150 : ELSE
151 6 : io_default = .FALSE.
152 6 : IF (para_env%is_source()) THEN
153 : CALL open_file(filename, unit_number=bs_data_unit, file_status="UNKNOWN", file_action="WRITE", &
154 3 : file_position="APPEND")
155 : ELSE
156 3 : bs_data_unit = -1
157 : END IF
158 : END IF
159 18 : DO i_rep = 1, n_rep
160 10 : t1 = m_walltime()
161 10 : CALL section_vals_val_get(kpset, "NPOINTS", i_rep_section=i_rep, i_val=npline)
162 10 : CALL section_vals_val_get(kpset, "UNITS", i_rep_section=i_rep, c_val=ustr)
163 10 : CALL uppercase(ustr)
164 10 : CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_section=i_rep, n_rep_val=n_ptr)
165 30 : ALLOCATE (kspecial(3, n_ptr))
166 30 : ALLOCATE (spname(n_ptr))
167 30 : DO ip = 1, n_ptr
168 20 : CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_section=i_rep, i_rep_val=ip, c_vals=strptr)
169 20 : IF (SIZE(strptr(:), 1) == 4) THEN
170 6 : spname(ip) = strptr(1)
171 24 : DO i = 1, 3
172 18 : CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
173 24 : IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
174 : END DO
175 14 : ELSE IF (SIZE(strptr(:), 1) == 3) THEN
176 14 : spname(ip) = "not specified"
177 56 : DO i = 1, 3
178 42 : CALL read_float_object(strptr(i), kpptr(i), error_message)
179 56 : IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
180 : END DO
181 : ELSE
182 0 : CPABORT("Input SPECIAL_POINT invalid")
183 : END IF
184 10 : SELECT CASE (ustr)
185 : CASE ("B_VECTOR")
186 48 : kspecial(1:3, ip) = kpptr(1:3)
187 : CASE ("CART_ANGSTROM")
188 : kspecial(1:3, ip) = (kpptr(1)*cell%hmat(1, 1:3) + &
189 : kpptr(2)*cell%hmat(2, 1:3) + &
190 0 : kpptr(3)*cell%hmat(3, 1:3))/twopi*angstrom
191 : CASE ("CART_BOHR")
192 : kspecial(1:3, ip) = (kpptr(1)*cell%hmat(1, 1:3) + &
193 : kpptr(2)*cell%hmat(2, 1:3) + &
194 56 : kpptr(3)*cell%hmat(3, 1:3))/twopi
195 : CASE DEFAULT
196 20 : CPABORT("Unknown unit <"//TRIM(ustr)//"> specified for k-point definition")
197 : END SELECT
198 : END DO
199 10 : npoints = (n_ptr - 1)*npline + 1
200 10 : CPASSERT(npoints >= 1)
201 :
202 : ! Initialize environment and calculate MOs
203 30 : ALLOCATE (kpgeneral(3, npoints))
204 70 : kpgeneral(1:3, 1) = kspecial(1:3, 1)
205 10 : ikk = 1
206 20 : DO ik = 2, n_ptr
207 100 : DO ip = 1, npline
208 80 : ikk = ikk + 1
209 : kpgeneral(1:3, ikk) = kspecial(1:3, ik - 1) + &
210 : REAL(ip, KIND=dp)/REAL(npline, KIND=dp)* &
211 570 : (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
212 : END DO
213 : END DO
214 10 : NULLIFY (kpoint)
215 10 : CALL calculate_kp_orbitals(qs_env, kpoint, "GENERAL", nadd, kpgeneral=kpgeneral)
216 10 : DEALLOCATE (kpgeneral)
217 :
218 10 : CALL get_qs_env(qs_env, dft_control=dft_control)
219 10 : nspins = dft_control%nspins
220 10 : kp => kpoint%kp_env(1)%kpoint_env
221 10 : CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
222 50 : ALLOCATE (eigval(nmo), occnum(nmo))
223 10 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp, wkp=wkp)
224 :
225 10 : IF (unit_nr > 0) THEN
226 : WRITE (UNIT=unit_nr, FMT="(T2,A,I4,T71,I10)") &
227 5 : "KPOINTS| Number of k-points in set ", i_rep, npoints
228 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
229 5 : "KPOINTS| In units of b-vector [2pi/Bohr]"
230 15 : DO ip = 1, n_ptr
231 : WRITE (UNIT=unit_nr, FMT="(T2,A,I5,1X,A11,3(1X,F12.6))") &
232 15 : "KPOINTS| Special point ", ip, ADJUSTL(TRIM(spname(ip))), kspecial(1:3, ip)
233 : END DO
234 : END IF
235 10 : IF (bs_data_unit > 0 .AND. (bs_data_unit /= unit_nr)) THEN
236 : WRITE (UNIT=bs_data_unit, FMT="(4(A,I0),A)") &
237 4 : "# Set ", i_rep, ": ", n_ptr, " special points, ", npoints, " k-points, ", nmo, " bands"
238 13 : DO ip = 1, n_ptr
239 : WRITE (UNIT=bs_data_unit, FMT="(A,I0,T20,T24,3(1X,F14.8),2X,A)") &
240 13 : "# Special point ", ip, kspecial(1:3, ip), ADJUSTL(TRIM(spname(ip)))
241 : END DO
242 : END IF
243 :
244 100 : DO ik = 1, nkp
245 90 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
246 190 : DO ispin = 1, nspins
247 90 : IF (my_kpgrp) THEN
248 90 : ikpgr = ik - kp_range(1) + 1
249 90 : kp => kpoint%kp_env(ikpgr)%kpoint_env
250 90 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues, occupation_numbers=occupation_numbers)
251 3370 : eigval(1:nmo) = eigenvalues(1:nmo)
252 3370 : occnum(1:nmo) = occupation_numbers(1:nmo)
253 : ELSE
254 0 : eigval(1:nmo) = 0.0_dp
255 0 : occnum(1:nmo) = 0.0_dp
256 : END IF
257 3370 : CALL kpoint%para_env_inter_kp%sum(eigval)
258 3370 : CALL kpoint%para_env_inter_kp%sum(occnum)
259 180 : IF (bs_data_unit > 0) THEN
260 : WRITE (UNIT=bs_data_unit, FMT="(A,I0,T15,A,I0,A,T24,3(1X,F14.8),3X,F14.8)") &
261 45 : "# Point ", ik, " Spin ", ispin, ":", xkp(1:3, ik), wkp(ik)
262 : WRITE (UNIT=bs_data_unit, FMT="(A)") &
263 45 : "# Band Energy [eV] Occupation"
264 865 : DO imo = 1, nmo
265 : WRITE (UNIT=bs_data_unit, FMT="(T2,I7,2(1X,F14.8))") &
266 865 : imo, eigval(imo)*evolt, occnum(imo)
267 : END DO
268 : END IF
269 : END DO
270 : END DO
271 :
272 10 : DEALLOCATE (kspecial, spname)
273 10 : DEALLOCATE (eigval, occnum)
274 10 : CALL kpoint_release(kpoint)
275 10 : t2 = m_walltime()
276 48 : IF (unit_nr > 0) THEN
277 5 : WRITE (UNIT=unit_nr, FMT="(T2,A,T67,F14.3)") "KPOINTS| Time for k-point line ", t2 - t1
278 : END IF
279 :
280 : END DO
281 :
282 : ! Close output files
283 8 : IF (.NOT. io_default) THEN
284 6 : IF (para_env%is_source()) CALL close_file(bs_data_unit)
285 : END IF
286 :
287 : END IF
288 :
289 8 : END SUBROUTINE do_calculate_band_structure
290 :
291 : ! **************************************************************************************************
292 : !> \brief diagonalize KS matrices at a set of kpoints
293 : !> \param qs_env ...
294 : !> \param kpoint ...
295 : !> \param scheme ...
296 : !> \param nadd ...
297 : !> \param mp_grid ...
298 : !> \param kpgeneral ...
299 : !> \param group_size_ext ...
300 : !> \author JGH
301 : ! **************************************************************************************************
302 10 : SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgeneral, group_size_ext)
303 : TYPE(qs_environment_type), POINTER :: qs_env
304 : TYPE(kpoint_type), POINTER :: kpoint
305 : CHARACTER(LEN=*), INTENT(IN) :: scheme
306 : INTEGER, INTENT(IN) :: nadd
307 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid
308 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
309 : OPTIONAL :: kpgeneral
310 : INTEGER, INTENT(IN), OPTIONAL :: group_size_ext
311 :
312 : LOGICAL :: diis_step
313 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
314 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
315 : TYPE(dft_control_type), POINTER :: dft_control
316 : TYPE(mp_para_env_type), POINTER :: para_env
317 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
318 10 : POINTER :: sab_nl
319 : TYPE(qs_scf_env_type), POINTER :: scf_env
320 : TYPE(scf_control_type), POINTER :: scf_control
321 :
322 10 : CALL calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral)
323 :
324 10 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
325 10 : CALL kpoint_env_initialize(kpoint, para_env, blacs_env)
326 :
327 10 : CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd)
328 10 : CALL kpoint_initialize_mo_set(kpoint)
329 :
330 10 : CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control)
331 10 : CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
332 :
333 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
334 10 : scf_env=scf_env, scf_control=scf_control)
335 10 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
336 :
337 10 : END SUBROUTINE calculate_kp_orbitals
338 :
339 : ! **************************************************************************************************
340 : !> \brief ...
341 : !> \param kpoint ...
342 : !> \param scheme ...
343 : !> \param group_size_ext ...
344 : !> \param mp_grid ...
345 : !> \param kpgeneral ...
346 : ! **************************************************************************************************
347 46 : SUBROUTINE calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral)
348 :
349 : TYPE(kpoint_type), POINTER :: kpoint
350 : CHARACTER(LEN=*), INTENT(IN) :: scheme
351 : INTEGER, INTENT(IN), OPTIONAL :: group_size_ext
352 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid
353 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
354 : OPTIONAL :: kpgeneral
355 :
356 : INTEGER :: i, ix, iy, iz, npoints
357 :
358 64 : CPASSERT(.NOT. ASSOCIATED(kpoint))
359 :
360 64 : CALL kpoint_create(kpoint)
361 :
362 64 : kpoint%kp_scheme = scheme
363 64 : kpoint%symmetry = .FALSE.
364 64 : kpoint%verbose = .FALSE.
365 64 : kpoint%full_grid = .FALSE.
366 64 : kpoint%use_real_wfn = .FALSE.
367 64 : kpoint%eps_geo = 1.e-6_dp
368 64 : IF (PRESENT(group_size_ext)) THEN
369 54 : kpoint%parallel_group_size = group_size_ext
370 : ELSE
371 10 : kpoint%parallel_group_size = -1
372 : END IF
373 0 : SELECT CASE (scheme)
374 : CASE ("GAMMA")
375 0 : kpoint%nkp = 1
376 0 : ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
377 0 : kpoint%xkp(1:3, 1) = 0.0_dp
378 0 : kpoint%wkp(1) = 1.0_dp
379 0 : kpoint%symmetry = .TRUE.
380 0 : ALLOCATE (kpoint%kp_sym(1))
381 0 : NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
382 0 : CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
383 : CASE ("MONKHORST-PACK")
384 18 : CPASSERT(PRESENT(mp_grid))
385 18 : npoints = mp_grid(1)*mp_grid(2)*mp_grid(3)
386 72 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
387 18 : kpoint%full_grid = .TRUE.
388 18 : kpoint%nkp = npoints
389 90 : ALLOCATE (kpoint%xkp(3, npoints), kpoint%wkp(npoints))
390 36 : kpoint%wkp(:) = 1._dp/REAL(npoints, KIND=dp)
391 : i = 0
392 36 : DO ix = 1, mp_grid(1)
393 54 : DO iy = 1, mp_grid(2)
394 54 : DO iz = 1, mp_grid(3)
395 18 : i = i + 1
396 18 : kpoint%xkp(1, i) = REAL(2*ix - mp_grid(1) - 1, KIND=dp)/(2._dp*REAL(mp_grid(1), KIND=dp))
397 18 : kpoint%xkp(2, i) = REAL(2*iy - mp_grid(2) - 1, KIND=dp)/(2._dp*REAL(mp_grid(2), KIND=dp))
398 36 : kpoint%xkp(3, i) = REAL(2*iz - mp_grid(3) - 1, KIND=dp)/(2._dp*REAL(mp_grid(3), KIND=dp))
399 : END DO
400 : END DO
401 : END DO
402 : ! default: no symmetry settings
403 54 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
404 36 : DO i = 1, kpoint%nkp
405 18 : NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
406 36 : CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
407 : END DO
408 : CASE ("MACDONALD")
409 0 : CPABORT("MACDONALD not implemented")
410 : CASE ("GENERAL")
411 46 : CPASSERT(PRESENT(kpgeneral))
412 46 : npoints = SIZE(kpgeneral, 2)
413 46 : kpoint%nkp = npoints
414 230 : ALLOCATE (kpoint%xkp(3, npoints), kpoint%wkp(npoints))
415 408 : kpoint%wkp(:) = 1._dp/REAL(npoints, KIND=dp)
416 1494 : kpoint%xkp(1:3, 1:npoints) = kpgeneral(1:3, 1:npoints)
417 : ! default: no symmetry settings
418 138 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
419 408 : DO i = 1, kpoint%nkp
420 362 : NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
421 408 : CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
422 : END DO
423 : CASE DEFAULT
424 64 : CPABORT("Unknown kpoint scheme requested")
425 : END SELECT
426 :
427 64 : END SUBROUTINE calculate_kpoints_for_bs
428 :
429 : END MODULE qs_band_structure
|
