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