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 Interface to Wannier90 code
10 : !> \par History
11 : !> 06.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_wannier90
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell
18 : USE cp_blacs_env, ONLY: cp_blacs_env_type
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
21 : dbcsr_deallocate_matrix,&
22 : dbcsr_p_type,&
23 : dbcsr_set,&
24 : dbcsr_type,&
25 : dbcsr_type_antisymmetric,&
26 : dbcsr_type_symmetric
27 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
28 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
29 : dbcsr_deallocate_matrix_set
30 : USE cp_files, ONLY: close_file,&
31 : open_file
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
36 : cp_fm_create,&
37 : cp_fm_get_element,&
38 : cp_fm_release,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
41 : cp_logger_type
42 : USE input_section_types, ONLY: section_vals_get,&
43 : section_vals_get_subs_vals,&
44 : section_vals_type,&
45 : section_vals_val_get
46 : USE kinds, ONLY: default_string_length,&
47 : dp
48 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
49 : kpoint_init_cell_index,&
50 : kpoint_initialize_mo_set,&
51 : kpoint_initialize_mos,&
52 : rskp_transform
53 : USE kpoint_types, ONLY: get_kpoint_info,&
54 : kpoint_create,&
55 : kpoint_env_type,&
56 : kpoint_release,&
57 : kpoint_type
58 : USE machine, ONLY: m_timestamp,&
59 : timestamp_length
60 : USE mathconstants, ONLY: twopi
61 : USE message_passing, ONLY: mp_para_env_type
62 : USE parallel_gemm_api, ONLY: parallel_gemm
63 : USE particle_types, ONLY: particle_type
64 : USE physcon, ONLY: angstrom,&
65 : evolt
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_env_release,&
68 : qs_environment_type
69 : USE qs_gamma2kp, ONLY: create_kp_from_gamma
70 : USE qs_mo_types, ONLY: get_mo_set,&
71 : mo_set_type
72 : USE qs_moments, ONLY: build_berry_kpoint_matrix
73 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
74 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
75 : USE qs_scf_types, ONLY: qs_scf_env_type
76 : USE scf_control_types, ONLY: scf_control_type
77 : USE wannier90, ONLY: wannier_setup
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 : PRIVATE
82 :
83 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wannier90'
84 :
85 : TYPE berry_matrix_type
86 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: sinmat => NULL(), cosmat => NULL()
87 : END TYPE berry_matrix_type
88 :
89 : PUBLIC :: wannier90_interface
90 :
91 : ! **************************************************************************************************
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief ...
97 : !> \param input ...
98 : !> \param logger ...
99 : !> \param qs_env ...
100 : ! **************************************************************************************************
101 20102 : SUBROUTINE wannier90_interface(input, logger, qs_env)
102 : TYPE(section_vals_type), POINTER :: input
103 : TYPE(cp_logger_type), POINTER :: logger
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 :
106 : CHARACTER(len=*), PARAMETER :: routineN = 'wannier90_interface'
107 :
108 : INTEGER :: handle, iw
109 : LOGICAL :: explicit
110 : TYPE(section_vals_type), POINTER :: w_input
111 :
112 : !--------------------------------------------------------------------------------------------!
113 :
114 10051 : CALL timeset(routineN, handle)
115 : w_input => section_vals_get_subs_vals(section_vals=input, &
116 10051 : subsection_name="DFT%PRINT%WANNIER90")
117 10051 : CALL section_vals_get(w_input, explicit=explicit)
118 10051 : IF (explicit) THEN
119 :
120 0 : iw = cp_logger_get_default_io_unit(logger)
121 :
122 0 : IF (iw > 0) THEN
123 : WRITE (iw, '(/,T2,A)') &
124 0 : '!-----------------------------------------------------------------------------!'
125 0 : WRITE (iw, '(T32,A)') "Interface to Wannier90"
126 : WRITE (iw, '(T2,A)') &
127 0 : '!-----------------------------------------------------------------------------!'
128 : END IF
129 :
130 0 : CALL wannier90_files(qs_env, w_input, iw)
131 :
132 0 : IF (iw > 0) THEN
133 : WRITE (iw, '(/,T2,A)') &
134 0 : '!--------------------------------End of Wannier90-----------------------------!'
135 : END IF
136 : END IF
137 10051 : CALL timestop(handle)
138 :
139 10051 : END SUBROUTINE wannier90_interface
140 :
141 : ! **************************************************************************************************
142 : !> \brief ...
143 : !> \param qs_env ...
144 : !> \param input ...
145 : !> \param iw ...
146 : ! **************************************************************************************************
147 0 : SUBROUTINE wannier90_files(qs_env, input, iw)
148 : TYPE(qs_environment_type), POINTER :: qs_env
149 : TYPE(section_vals_type), POINTER :: input
150 : INTEGER, INTENT(IN) :: iw
151 :
152 : INTEGER, PARAMETER :: num_nnmax = 12
153 :
154 : CHARACTER(len=2) :: asym
155 0 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
156 : CHARACTER(len=default_string_length) :: filename, seed_name
157 : CHARACTER(LEN=timestamp_length) :: timestamp
158 : INTEGER :: i, i_rep, ib, ib1, ib2, ibs, ik, ik2, ikk, ikpgr, ispin, iunit, ix, iy, iz, k, &
159 : n_rep, nadd, nao, nbs, nexcl, nkp, nmo, nntot, nspins, num_atoms, num_bands, &
160 : num_bands_tot, num_kpts, num_wann
161 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: exclude_bands
162 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nblist, nnlist
163 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: nncell
164 : INTEGER, DIMENSION(2) :: kp_range
165 : INTEGER, DIMENSION(3) :: mp_grid
166 0 : INTEGER, DIMENSION(:), POINTER :: invals
167 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
168 : LOGICAL :: diis_step, do_kpoints, gamma_only, &
169 : my_kpgrp, mygrp, spinors
170 : REAL(KIND=dp) :: cmmn, ksign, rmmn
171 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
172 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, b_latt, kpt_latt
173 : REAL(KIND=dp), DIMENSION(3) :: bvec
174 : REAL(KIND=dp), DIMENSION(3, 3) :: real_lattice, recip_lattice
175 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
176 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
177 0 : TYPE(berry_matrix_type), DIMENSION(:), POINTER :: berry_matrix
178 : TYPE(cell_type), POINTER :: cell
179 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
180 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_mmn, matrix_struct_work
181 : TYPE(cp_fm_type) :: fm_tmp, mmn_imag, mmn_real
182 0 : TYPE(cp_fm_type), DIMENSION(2) :: fmk1, fmk2
183 : TYPE(cp_fm_type), POINTER :: fmdummy, fmi, fmr
184 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
185 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix
186 : TYPE(dft_control_type), POINTER :: dft_control
187 : TYPE(kpoint_env_type), POINTER :: kp
188 : TYPE(kpoint_type), POINTER :: kpoint
189 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
190 : TYPE(mp_para_env_type), POINTER :: para_env
191 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
192 0 : POINTER :: sab_nl
193 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
194 : TYPE(qs_environment_type), POINTER :: qs_env_kp
195 : TYPE(qs_scf_env_type), POINTER :: scf_env
196 : TYPE(scf_control_type), POINTER :: scf_control
197 :
198 : !--------------------------------------------------------------------------------------------!
199 :
200 : ! add code for exclude_bands and projectors
201 :
202 : ! generate all arrays needed for the setup call
203 0 : CALL section_vals_val_get(input, "SEED_NAME", c_val=seed_name)
204 0 : CALL section_vals_val_get(input, "MP_GRID", i_vals=invals)
205 0 : CALL section_vals_val_get(input, "WANNIER_FUNCTIONS", i_val=num_wann)
206 0 : CALL section_vals_val_get(input, "ADDED_MOS", i_val=nadd)
207 0 : mp_grid(1:3) = invals(1:3)
208 0 : num_kpts = mp_grid(1)*mp_grid(2)*mp_grid(3)
209 : ! excluded bands
210 0 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", n_rep_val=n_rep)
211 0 : nexcl = 0
212 0 : DO i_rep = 1, n_rep
213 0 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
214 0 : nexcl = nexcl + SIZE(invals)
215 : END DO
216 0 : IF (nexcl > 0) THEN
217 0 : ALLOCATE (exclude_bands(nexcl))
218 0 : nexcl = 0
219 0 : DO i_rep = 1, n_rep
220 0 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
221 0 : exclude_bands(nexcl + 1:nexcl + SIZE(invals)) = invals(:)
222 0 : nexcl = nexcl + SIZE(invals)
223 : END DO
224 : END IF
225 : !
226 : ! lattice -> Angstrom
227 0 : CALL get_qs_env(qs_env, cell=cell)
228 0 : CALL get_cell(cell, h=real_lattice, h_inv=recip_lattice)
229 0 : real_lattice(1:3, 1:3) = angstrom*real_lattice(1:3, 1:3)
230 0 : recip_lattice(1:3, 1:3) = (twopi/angstrom)*TRANSPOSE(recip_lattice(1:3, 1:3))
231 : ! k-points
232 0 : ALLOCATE (kpt_latt(3, num_kpts))
233 0 : CALL get_qs_env(qs_env, particle_set=particle_set)
234 0 : NULLIFY (kpoint)
235 0 : CALL kpoint_create(kpoint)
236 0 : kpoint%kp_scheme = "MONKHORST-PACK"
237 0 : kpoint%symmetry = .FALSE.
238 0 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
239 0 : kpoint%verbose = .FALSE.
240 0 : kpoint%full_grid = .TRUE.
241 0 : kpoint%eps_geo = 1.0e-6_dp
242 0 : kpoint%use_real_wfn = .FALSE.
243 0 : kpoint%parallel_group_size = 0
244 0 : i = 0
245 0 : DO ix = 0, mp_grid(1) - 1
246 0 : DO iy = 0, mp_grid(2) - 1
247 0 : DO iz = 0, mp_grid(3) - 1
248 0 : i = i + 1
249 0 : kpt_latt(1, i) = REAL(ix, KIND=dp)/REAL(mp_grid(1), KIND=dp)
250 0 : kpt_latt(2, i) = REAL(iy, KIND=dp)/REAL(mp_grid(2), KIND=dp)
251 0 : kpt_latt(3, i) = REAL(iz, KIND=dp)/REAL(mp_grid(3), KIND=dp)
252 : END DO
253 : END DO
254 : END DO
255 0 : kpoint%nkp = num_kpts
256 0 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
257 0 : kpoint%wkp(:) = 1._dp/REAL(num_kpts, KIND=dp)
258 0 : DO i = 1, num_kpts
259 0 : kpoint%xkp(1:3, i) = (angstrom/twopi)*MATMUL(recip_lattice, kpt_latt(:, i))
260 : END DO
261 : ! number of bands in calculation
262 0 : CALL get_qs_env(qs_env, mos=mos)
263 0 : CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=num_bands_tot)
264 0 : num_bands_tot = MIN(nao, num_bands_tot + nadd)
265 0 : num_bands = num_wann
266 0 : num_atoms = SIZE(particle_set)
267 0 : ALLOCATE (atoms_cart(3, num_atoms))
268 0 : ALLOCATE (atom_symbols(num_atoms))
269 0 : DO i = 1, num_atoms
270 0 : atoms_cart(1:3, i) = particle_set(i)%r(1:3)
271 0 : CALL get_atomic_kind(particle_set(i)%atomic_kind, element_symbol=asym)
272 0 : atom_symbols(i) = asym
273 : END DO
274 0 : gamma_only = .FALSE.
275 0 : spinors = .FALSE.
276 : ! output
277 0 : ALLOCATE (nnlist(num_kpts, num_nnmax))
278 0 : ALLOCATE (nncell(3, num_kpts, num_nnmax))
279 0 : nnlist = 0
280 0 : nncell = 0
281 0 : nntot = 0
282 :
283 0 : IF (iw > 0) THEN
284 : ! setup
285 : CALL wannier_setup(mp_grid, num_kpts, real_lattice, recip_lattice, &
286 0 : kpt_latt, nntot, nnlist, nncell, iw)
287 : END IF
288 :
289 0 : CALL get_qs_env(qs_env, para_env=para_env)
290 0 : CALL para_env%sum(nntot)
291 0 : CALL para_env%sum(nnlist)
292 0 : CALL para_env%sum(nncell)
293 :
294 0 : IF (para_env%is_source()) THEN
295 : ! Write the Wannier90 input file "seed_name.win"
296 0 : WRITE (filename, '(A,A)') TRIM(seed_name), ".win"
297 0 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
298 : !
299 0 : CALL m_timestamp(timestamp)
300 0 : WRITE (iunit, "(A)") "! Wannier90 input file generated by CP2K "
301 0 : WRITE (iunit, "(A,/)") "! Creation date "//timestamp
302 : !
303 0 : WRITE (iunit, "(A,I5)") "num_wann = ", num_wann
304 0 : IF (num_bands /= num_wann) THEN
305 0 : WRITE (iunit, "(A,I5)") "num_bands = ", num_bands
306 : END IF
307 0 : WRITE (iunit, "(/,A,/)") "length_unit = bohr "
308 0 : WRITE (iunit, "(/,A,/)") "! System"
309 0 : WRITE (iunit, "(/,A)") "begin unit_cell_cart"
310 0 : WRITE (iunit, "(A)") "bohr"
311 0 : DO i = 1, 3
312 0 : WRITE (iunit, "(3F12.6)") cell%hmat(i, 1:3)
313 : END DO
314 0 : WRITE (iunit, "(A,/)") "end unit_cell_cart"
315 0 : WRITE (iunit, "(/,A)") "begin atoms_cart"
316 0 : DO i = 1, num_atoms
317 0 : WRITE (iunit, "(A,3F15.10)") atom_symbols(i), atoms_cart(1:3, i)
318 : END DO
319 0 : WRITE (iunit, "(A,/)") "end atoms_cart"
320 0 : WRITE (iunit, "(/,A,/)") "! Kpoints"
321 0 : WRITE (iunit, "(/,A,3I6/)") "mp_grid = ", mp_grid(1:3)
322 0 : WRITE (iunit, "(A)") "begin kpoints"
323 0 : DO i = 1, num_kpts
324 0 : WRITE (iunit, "(3F12.6)") kpt_latt(1:3, i)
325 : END DO
326 0 : WRITE (iunit, "(A)") "end kpoints"
327 0 : CALL close_file(iunit)
328 : ELSE
329 0 : iunit = -1
330 : END IF
331 :
332 : ! calculate bands
333 0 : NULLIFY (qs_env_kp)
334 0 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
335 0 : IF (do_kpoints) THEN
336 : ! we already do kpoints
337 0 : qs_env_kp => qs_env
338 : ELSE
339 : ! we start from gamma point only
340 0 : ALLOCATE (qs_env_kp)
341 0 : CALL create_kp_from_gamma(qs_env, qs_env_kp)
342 : END IF
343 0 : IF (iw > 0) THEN
344 0 : WRITE (unit=iw, FMT="(/,T2,A)") "Start K-Point Calculation ..."
345 : END IF
346 0 : CALL get_qs_env(qs_env=qs_env_kp, para_env=para_env, blacs_env=blacs_env)
347 0 : CALL kpoint_env_initialize(kpoint, para_env, blacs_env)
348 0 : CALL kpoint_initialize_mos(kpoint, mos, nadd)
349 0 : CALL kpoint_initialize_mo_set(kpoint)
350 : !
351 0 : CALL get_qs_env(qs_env=qs_env_kp, sab_orb=sab_nl, dft_control=dft_control)
352 0 : CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
353 : !
354 : CALL get_qs_env(qs_env=qs_env_kp, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
355 0 : scf_env=scf_env, scf_control=scf_control)
356 0 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
357 : !
358 0 : IF (iw > 0) THEN
359 0 : WRITE (iw, '(T69,A)') "... Finished"
360 : END IF
361 : !
362 : ! Calculate and print Overlaps
363 : !
364 0 : IF (para_env%is_source()) THEN
365 0 : WRITE (filename, '(A,A)') TRIM(seed_name), ".mmn"
366 0 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
367 0 : CALL m_timestamp(timestamp)
368 0 : WRITE (iunit, "(A)") "! Wannier90 file generated by CP2K "//timestamp
369 0 : WRITE (iunit, "(3I8)") num_bands, num_kpts, nntot
370 : ELSE
371 0 : iunit = -1
372 : END IF
373 : ! create a list of unique b vectors and a table of pointers
374 : ! nblist(ik,i) -> +/- b_latt(1:3,x)
375 0 : ALLOCATE (nblist(num_kpts, nntot))
376 0 : ALLOCATE (b_latt(3, num_kpts*nntot))
377 0 : nblist = 0
378 0 : nbs = 0
379 0 : DO ik = 1, num_kpts
380 0 : DO i = 1, nntot
381 0 : bvec(1:3) = kpt_latt(1:3, nnlist(ik, i)) - kpt_latt(1:3, ik) + nncell(1:3, ik, i)
382 0 : ibs = 0
383 0 : DO k = 1, nbs
384 0 : IF (SUM(ABS(bvec(1:3) - b_latt(1:3, k))) < 1.e-6_dp) THEN
385 : ibs = k
386 : EXIT
387 : END IF
388 0 : IF (SUM(ABS(bvec(1:3) + b_latt(1:3, k))) < 1.e-6_dp) THEN
389 0 : ibs = -k
390 0 : EXIT
391 : END IF
392 : END DO
393 0 : IF (ibs /= 0) THEN
394 : ! old lattice vector
395 0 : nblist(ik, i) = ibs
396 : ELSE
397 : ! new lattice vector
398 0 : nbs = nbs + 1
399 0 : b_latt(1:3, nbs) = bvec(1:3)
400 0 : nblist(ik, i) = nbs
401 : END IF
402 : END DO
403 : END DO
404 : ! calculate all the operator matrices (a|bvec|b)
405 0 : ALLOCATE (berry_matrix(nbs))
406 0 : DO i = 1, nbs
407 0 : NULLIFY (berry_matrix(i)%cosmat)
408 0 : NULLIFY (berry_matrix(i)%sinmat)
409 0 : bvec(1:3) = twopi*MATMUL(TRANSPOSE(cell%h_inv(1:3, 1:3)), b_latt(1:3, i))
410 : CALL build_berry_kpoint_matrix(qs_env_kp, berry_matrix(i)%cosmat, &
411 0 : berry_matrix(i)%sinmat, bvec)
412 : END DO
413 : ! work matrices for MOs (all group)
414 0 : kp => kpoint%kp_env(1)%kpoint_env
415 0 : CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
416 0 : NULLIFY (matrix_struct_work)
417 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, &
418 : ncol_global=nmo, &
419 : para_env=para_env, &
420 0 : context=blacs_env)
421 0 : CALL cp_fm_create(fm_tmp, matrix_struct_work)
422 0 : DO i = 1, 2
423 0 : CALL cp_fm_create(fmk1(i), matrix_struct_work)
424 0 : CALL cp_fm_create(fmk2(i), matrix_struct_work)
425 : END DO
426 : ! work matrices for Mmn(k,b) integrals
427 0 : NULLIFY (matrix_struct_mmn)
428 : CALL cp_fm_struct_create(matrix_struct_mmn, nrow_global=nmo, &
429 : ncol_global=nmo, &
430 : para_env=para_env, &
431 0 : context=blacs_env)
432 0 : CALL cp_fm_create(mmn_real, matrix_struct_mmn)
433 0 : CALL cp_fm_create(mmn_imag, matrix_struct_mmn)
434 : ! allocate some work matrices
435 0 : ALLOCATE (rmatrix, cmatrix)
436 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
437 0 : matrix_type=dbcsr_type_symmetric)
438 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
439 0 : matrix_type=dbcsr_type_antisymmetric)
440 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
441 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
442 : !
443 0 : CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
444 0 : NULLIFY (fmdummy)
445 0 : nspins = dft_control%nspins
446 0 : DO ispin = 1, nspins
447 : ! loop over all k-points
448 0 : DO ik = 1, num_kpts
449 : ! get the MO coefficients for this k-point
450 0 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
451 : IF (my_kpgrp) THEN
452 0 : ikk = ik - kpoint%kp_range(1) + 1
453 0 : kp => kpoint%kp_env(ikk)%kpoint_env
454 0 : CPASSERT(SIZE(kp%mos, 1) == 2)
455 0 : fmr => kp%mos(1, ispin)%mo_coeff
456 0 : fmi => kp%mos(2, ispin)%mo_coeff
457 0 : CALL cp_fm_copy_general(fmr, fmk1(1), para_env)
458 0 : CALL cp_fm_copy_general(fmi, fmk1(2), para_env)
459 : ELSE
460 0 : NULLIFY (fmr, fmi, kp)
461 0 : CALL cp_fm_copy_general(fmdummy, fmk1(1), para_env)
462 0 : CALL cp_fm_copy_general(fmdummy, fmk1(2), para_env)
463 : END IF
464 : ! loop over all connected neighbors
465 0 : DO i = 1, nntot
466 : ! get the MO coefficients for the connected k-point
467 0 : ik2 = nnlist(ik, i)
468 0 : mygrp = (ik2 >= kpoint%kp_range(1) .AND. ik2 <= kpoint%kp_range(2))
469 : IF (mygrp) THEN
470 0 : ikk = ik2 - kpoint%kp_range(1) + 1
471 0 : kp => kpoint%kp_env(ikk)%kpoint_env
472 0 : CPASSERT(SIZE(kp%mos, 1) == 2)
473 0 : fmr => kp%mos(1, ispin)%mo_coeff
474 0 : fmi => kp%mos(2, ispin)%mo_coeff
475 0 : CALL cp_fm_copy_general(fmr, fmk2(1), para_env)
476 0 : CALL cp_fm_copy_general(fmi, fmk2(2), para_env)
477 : ELSE
478 0 : NULLIFY (fmr, fmi, kp)
479 0 : CALL cp_fm_copy_general(fmdummy, fmk2(1), para_env)
480 0 : CALL cp_fm_copy_general(fmdummy, fmk2(2), para_env)
481 : END IF
482 : !
483 : ! transfer realspace overlaps to connected k-point
484 0 : ibs = nblist(ik, i)
485 0 : ksign = SIGN(1.0_dp, REAL(ibs, KIND=dp))
486 0 : ibs = ABS(ibs)
487 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
488 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
489 : CALL rskp_transform(rmatrix, cmatrix, rsmat=berry_matrix(ibs)%cosmat, ispin=1, &
490 : xkp=kpoint%xkp(1:3, ik2), cell_to_index=cell_to_index, sab_nl=sab_nl, &
491 0 : is_complex=.FALSE., rs_sign=ksign)
492 : CALL rskp_transform(cmatrix, rmatrix, rsmat=berry_matrix(ibs)%sinmat, ispin=1, &
493 : xkp=kpoint%xkp(1:3, ik2), cell_to_index=cell_to_index, sab_nl=sab_nl, &
494 0 : is_complex=.TRUE., rs_sign=ksign)
495 : !
496 : ! calculate M_(mn)^(k,b)
497 0 : CALL cp_dbcsr_sm_fm_multiply(rmatrix, fmk2(1), fm_tmp, nmo)
498 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(1), fm_tmp, 0.0_dp, mmn_real)
499 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(2), fm_tmp, 0.0_dp, mmn_imag)
500 0 : CALL cp_dbcsr_sm_fm_multiply(rmatrix, fmk2(2), fm_tmp, nmo)
501 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(1), fm_tmp, 1.0_dp, mmn_imag)
502 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(2), fm_tmp, -1.0_dp, mmn_real)
503 0 : CALL cp_dbcsr_sm_fm_multiply(cmatrix, fmk2(1), fm_tmp, nmo)
504 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(1), fm_tmp, 1.0_dp, mmn_imag)
505 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(2), fm_tmp, -1.0_dp, mmn_real)
506 0 : CALL cp_dbcsr_sm_fm_multiply(cmatrix, fmk2(2), fm_tmp, nmo)
507 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(1), fm_tmp, -1.0_dp, mmn_real)
508 0 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, fmk1(2), fm_tmp, -1.0_dp, mmn_imag)
509 : !
510 : ! write to output file
511 0 : IF (para_env%is_source()) THEN
512 0 : WRITE (iunit, "(2I8,3I5)") ik, ik2, nncell(1:3, ik, i)
513 : END IF
514 0 : DO ib2 = 1, nmo
515 0 : DO ib1 = 1, nmo
516 0 : CALL cp_fm_get_element(mmn_real, ib1, ib2, rmmn)
517 0 : CALL cp_fm_get_element(mmn_imag, ib1, ib2, cmmn)
518 0 : IF (para_env%is_source()) THEN
519 0 : WRITE (iunit, "(2E30.14)") rmmn, cmmn
520 : END IF
521 : END DO
522 : END DO
523 : !
524 : END DO
525 : END DO
526 : END DO
527 0 : DO i = 1, nbs
528 0 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%cosmat)
529 0 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%sinmat)
530 : END DO
531 0 : DEALLOCATE (berry_matrix)
532 0 : CALL cp_fm_struct_release(matrix_struct_work)
533 0 : DO i = 1, 2
534 0 : CALL cp_fm_release(fmk1(i))
535 0 : CALL cp_fm_release(fmk2(i))
536 : END DO
537 0 : CALL cp_fm_release(fm_tmp)
538 0 : CALL cp_fm_struct_release(matrix_struct_mmn)
539 0 : CALL cp_fm_release(mmn_real)
540 0 : CALL cp_fm_release(mmn_imag)
541 0 : CALL dbcsr_deallocate_matrix(rmatrix)
542 0 : CALL dbcsr_deallocate_matrix(cmatrix)
543 : !
544 0 : IF (para_env%is_source()) THEN
545 0 : CALL close_file(iunit)
546 : END IF
547 : !
548 : ! Calculate and print Projections
549 : !
550 : ! Print eigenvalues
551 0 : nspins = dft_control%nspins
552 0 : kp => kpoint%kp_env(1)%kpoint_env
553 0 : CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
554 0 : ALLOCATE (eigval(nmo))
555 0 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp)
556 0 : IF (para_env%is_source()) THEN
557 0 : WRITE (filename, '(A,A)') TRIM(seed_name), ".eig"
558 0 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
559 : ELSE
560 0 : iunit = -1
561 : END IF
562 : !
563 0 : DO ik = 1, nkp
564 0 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
565 0 : DO ispin = 1, nspins
566 0 : IF (my_kpgrp) THEN
567 0 : ikpgr = ik - kp_range(1) + 1
568 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
569 0 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
570 0 : eigval(1:nmo) = eigenvalues(1:nmo)
571 : ELSE
572 0 : eigval(1:nmo) = 0.0_dp
573 : END IF
574 0 : CALL kpoint%para_env_inter_kp%sum(eigval)
575 0 : eigval(1:nmo) = eigval(1:nmo)*evolt
576 : ! output
577 0 : IF (iunit > 0) THEN
578 0 : DO ib = 1, nmo
579 0 : WRITE (iunit, "(2I8,F24.14)") ib, ik, eigval(ib)
580 : END DO
581 : END IF
582 : END DO
583 : END DO
584 0 : IF (para_env%is_source()) THEN
585 0 : CALL close_file(iunit)
586 : END IF
587 : !
588 : ! clean up
589 0 : DEALLOCATE (kpt_latt, atoms_cart, atom_symbols, eigval)
590 0 : DEALLOCATE (nnlist, nncell)
591 0 : DEALLOCATE (nblist, b_latt)
592 0 : IF (nexcl > 0) THEN
593 0 : DEALLOCATE (exclude_bands)
594 : END IF
595 0 : IF (do_kpoints) THEN
596 0 : NULLIFY (qs_env_kp)
597 : ELSE
598 0 : CALL qs_env_release(qs_env_kp)
599 0 : DEALLOCATE (qs_env_kp)
600 : NULLIFY (qs_env_kp)
601 : END IF
602 :
603 0 : CALL kpoint_release(kpoint)
604 :
605 0 : END SUBROUTINE wannier90_files
606 :
607 0 : END MODULE qs_wannier90
|
