Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate localized minimal basis and analyze wavefunctions
10 : !> \par History
11 : !> 12.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE minbas_wfn_analysis
15 : USE atomic_charges, ONLY: print_atomic_charges,&
16 : print_bond_orders
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE bibliography, ONLY: Lu2004,&
19 : cite_reference
20 : USE cell_types, ONLY: cell_type
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_block_p, &
25 : dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
26 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
27 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
28 : dbcsr_type_symmetric
29 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
30 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
31 : cp_dbcsr_plus_fm_fm_t,&
32 : cp_dbcsr_sm_fm_multiply,&
33 : dbcsr_allocate_matrix_set,&
34 : dbcsr_deallocate_matrix_set
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_get_diag,&
41 : cp_fm_release,&
42 : cp_fm_to_fm,&
43 : cp_fm_type
44 : USE cp_log_handling, ONLY: cp_get_default_logger,&
45 : cp_logger_type
46 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
47 : cp_print_key_unit_nr
48 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
49 : USE input_section_types, ONLY: section_get_ivals,&
50 : section_vals_get,&
51 : section_vals_get_subs_vals,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE iterate_matrix, ONLY: invert_Hotelling
55 : USE kinds, ONLY: default_path_length,&
56 : dp
57 : USE message_passing, ONLY: mp_para_env_type
58 : USE minbas_methods, ONLY: minbas_calculation
59 : USE molden_utils, ONLY: write_mos_molden
60 : USE mulliken, ONLY: compute_bond_order,&
61 : mulliken_charges
62 : USE parallel_gemm_api, ONLY: parallel_gemm
63 : USE particle_list_types, ONLY: particle_list_type
64 : USE particle_methods, ONLY: get_particle_set
65 : USE particle_types, ONLY: particle_type
66 : USE pw_env_types, ONLY: pw_env_get,&
67 : pw_env_type
68 : USE pw_pool_types, ONLY: pw_pool_type
69 : USE pw_types, ONLY: pw_c1d_gs_type,&
70 : pw_r3d_rs_type
71 : USE qs_collocate_density, ONLY: calculate_wavefunction
72 : USE qs_environment_types, ONLY: get_qs_env,&
73 : qs_environment_type
74 : USE qs_kind_types, ONLY: qs_kind_type
75 : USE qs_ks_types, ONLY: get_ks_env,&
76 : qs_ks_env_type
77 : USE qs_mo_methods, ONLY: make_basis_lowdin
78 : USE qs_mo_types, ONLY: allocate_mo_set,&
79 : deallocate_mo_set,&
80 : get_mo_set,&
81 : mo_set_type,&
82 : set_mo_set
83 : USE qs_subsys_types, ONLY: qs_subsys_get,&
84 : qs_subsys_type
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_wfn_analysis'
91 :
92 : PUBLIC :: minbas_analysis
93 :
94 : ! **************************************************************************************************
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param qs_env ...
101 : !> \param input_section ...
102 : !> \param unit_nr ...
103 : ! **************************************************************************************************
104 28 : SUBROUTINE minbas_analysis(qs_env, input_section, unit_nr)
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : TYPE(section_vals_type), POINTER :: input_section
107 : INTEGER, INTENT(IN) :: unit_nr
108 :
109 : CHARACTER(len=*), PARAMETER :: routineN = 'minbas_analysis'
110 :
111 : INTEGER :: handle, homo, i, ispin, nao, natom, &
112 : nimages, nmao, nmo, nspin
113 28 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: ecount
114 28 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
115 : LOGICAL :: do_bondorder, explicit, full_ortho, occeq
116 : REAL(KIND=dp) :: alpha, amax, eps_filter, filter_eps, &
117 : trace
118 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: border, fnorm, mcharge, prmao
119 28 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
120 : TYPE(cell_type), POINTER :: cell
121 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
122 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c
123 : TYPE(cp_fm_type) :: fm1, fm2, fm3, fm4
124 : TYPE(cp_fm_type), POINTER :: fm_mos
125 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
126 28 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, pqmat, quambo, sqmat
127 28 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
128 : TYPE(dbcsr_type) :: psmat, sinv, smao, smaox, spmat
129 : TYPE(dbcsr_type), POINTER :: smat
130 : TYPE(dft_control_type), POINTER :: dft_control
131 28 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mbas
132 28 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
133 : TYPE(mp_para_env_type), POINTER :: para_env
134 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
135 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
136 : TYPE(qs_ks_env_type), POINTER :: ks_env
137 : TYPE(section_vals_type), POINTER :: molden_section
138 :
139 : ! only do MINBAS analysis if explicitly requested
140 28 : CALL section_vals_get(input_section, explicit=explicit)
141 28 : IF (.NOT. explicit) RETURN
142 :
143 : ! k-points?
144 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
145 0 : nspin = dft_control%nspins
146 0 : nimages = dft_control%nimages
147 0 : IF (nimages > 1) THEN
148 0 : IF (unit_nr > 0) THEN
149 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
150 0 : "K-Points: Localized Minimal Basis Analysis not available."
151 : END IF
152 : END IF
153 0 : IF (nimages > 1) RETURN
154 :
155 0 : CALL timeset(routineN, handle)
156 :
157 0 : IF (unit_nr > 0) THEN
158 0 : WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
159 0 : WRITE (UNIT=unit_nr, FMT="(T26,A)") "LOCALIZED MINIMAL BASIS ANALYSIS"
160 0 : WRITE (UNIT=unit_nr, FMT="(T18,A)") "W.C. Lu et al, J. Chem. Phys. 120, 2629 (2004)"
161 0 : WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
162 : END IF
163 0 : CALL cite_reference(Lu2004)
164 :
165 : ! input options
166 0 : CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter)
167 0 : CALL section_vals_val_get(input_section, "FULL_ORTHOGONALIZATION", l_val=full_ortho)
168 0 : CALL section_vals_val_get(input_section, "BOND_ORDER", l_val=do_bondorder)
169 :
170 : ! generate MAOs and QUAMBOs
171 0 : CALL get_qs_env(qs_env, mos=mos)
172 0 : NULLIFY (quambo, mao_coef)
173 : CALL minbas_calculation(qs_env, mos, quambo, mao=mao_coef, iounit=unit_nr, &
174 0 : full_ortho=full_ortho, eps_filter=eps_filter)
175 0 : IF (ASSOCIATED(quambo)) THEN
176 0 : CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=nmo)
177 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
178 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
179 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
180 0 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
181 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
182 0 : CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
183 0 : nmao = SUM(col_blk_sizes)
184 :
185 0 : NULLIFY (pqmat, sqmat)
186 0 : CALL dbcsr_allocate_matrix_set(sqmat, nspin)
187 0 : CALL dbcsr_allocate_matrix_set(pqmat, nspin)
188 0 : DO ispin = 1, nspin
189 0 : ALLOCATE (sqmat(ispin)%matrix)
190 : CALL dbcsr_create(matrix=sqmat(ispin)%matrix, &
191 : name="SQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
192 0 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
193 0 : ALLOCATE (pqmat(ispin)%matrix)
194 : CALL dbcsr_create(matrix=pqmat(ispin)%matrix, &
195 : name="PQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
196 0 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
197 : END DO
198 0 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
199 :
200 : ! Start wfn analysis
201 0 : IF (unit_nr > 0) THEN
202 0 : WRITE (unit_nr, '(/,T2,A)') 'Localized Minimal Basis Wavefunction Analysis'
203 : END IF
204 :
205 : ! localization of basis
206 0 : DO ispin = 1, nspin
207 0 : amax = dbcsr_get_occupation(quambo(ispin)%matrix)
208 0 : IF (unit_nr > 0) THEN
209 : WRITE (unit_nr, '(/,T2,A,I2,T69,F10.3,A2,/)') &
210 0 : 'Occupation of Basis Function Representation (Spin) ', ispin, amax*100._dp, ' %'
211 : END IF
212 : END DO
213 :
214 0 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
215 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
216 : CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, &
217 0 : para_env=para_env, context=blacs_env)
218 0 : CALL cp_fm_create(fm1, fm_struct_a)
219 : CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmao, ncol_global=nmo, &
220 0 : para_env=para_env, context=blacs_env)
221 0 : CALL cp_fm_create(fm2, fm_struct_b)
222 0 : CALL cp_fm_create(fm3, fm_struct_b)
223 : CALL cp_fm_struct_create(fm_struct_c, nrow_global=nmo, ncol_global=nmo, &
224 0 : para_env=para_env, context=blacs_env)
225 0 : CALL cp_fm_create(fm4, fm_struct_c)
226 0 : ALLOCATE (fnorm(nmo, nspin), ecount(natom, 3, nspin), prmao(natom, nspin))
227 0 : ecount = 0
228 0 : prmao = 0.0_dp
229 0 : DO ispin = 1, nspin
230 0 : CALL dbcsr_create(smao, name="S*QM", template=mao_coef(1)%matrix)
231 0 : smat => matrix_s(1, 1)%matrix
232 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, quambo(ispin)%matrix, 0.0_dp, smao)
233 : ! calculate atomic extend of basis
234 0 : CALL pm_extend(quambo(ispin)%matrix, smao, ecount(:, :, ispin))
235 0 : CALL dbcsr_create(sinv, name="QM*S*QM", template=sqmat(ispin)%matrix)
236 0 : CALL dbcsr_multiply("T", "N", 1.0_dp, quambo(ispin)%matrix, smao, 0.0_dp, sqmat(ispin)%matrix)
237 : ! atomic MAO projection
238 0 : CALL project_mao(mao_coef(ispin)%matrix, smao, sqmat(ispin)%matrix, prmao(:, ispin))
239 : ! invert overlap
240 0 : CALL invert_Hotelling(sinv, sqmat(ispin)%matrix, 1.e-6_dp, silent=.TRUE.)
241 0 : CALL dbcsr_create(smaox, name="S*QM*SINV", template=smao)
242 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, smao, sinv, 0.0_dp, smaox)
243 0 : CALL copy_dbcsr_to_fm(smaox, fm1)
244 0 : CALL get_mo_set(mos(ispin), mo_coeff=fm_mos, homo=homo)
245 0 : CALL parallel_gemm("T", "N", nmao, nmo, nao, 1.0_dp, fm1, fm_mos, 0.0_dp, fm2)
246 0 : CALL cp_dbcsr_sm_fm_multiply(sqmat(ispin)%matrix, fm2, fm3, nmo)
247 0 : CALL parallel_gemm("T", "N", nmo, nmo, nmao, 1.0_dp, fm2, fm3, 0.0_dp, fm4)
248 0 : CALL cp_fm_get_diag(fm4, fnorm(1:nmo, ispin))
249 : ! fm2 are the projected MOs (in MAO basis); orthogonalize the occupied subspace
250 0 : CALL make_basis_lowdin(vmatrix=fm2, ncol=homo, matrix_s=sqmat(ispin)%matrix)
251 : ! pmat
252 0 : CALL get_mo_set(mos(ispin), occupation_numbers=occupation_numbers, maxocc=alpha)
253 0 : occeq = ALL(occupation_numbers(1:homo) == alpha)
254 0 : CALL dbcsr_copy(pqmat(ispin)%matrix, sqmat(ispin)%matrix)
255 0 : CALL dbcsr_set(pqmat(ispin)%matrix, 0.0_dp)
256 0 : IF (occeq) THEN
257 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, &
258 0 : ncol=homo, alpha=alpha, keep_sparsity=.FALSE.)
259 : ELSE
260 0 : CALL cp_fm_to_fm(fm2, fm3)
261 0 : CALL cp_fm_column_scale(fm3, occupation_numbers(1:homo))
262 0 : alpha = 1.0_dp
263 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, &
264 0 : matrix_g=fm3, ncol=homo, alpha=alpha, keep_sparsity=.TRUE.)
265 : END IF
266 :
267 0 : CALL dbcsr_release(smao)
268 0 : CALL dbcsr_release(smaox)
269 0 : CALL dbcsr_release(sinv)
270 : END DO
271 : ! Basis extension
272 0 : CALL para_env%sum(ecount)
273 0 : IF (unit_nr > 0) THEN
274 0 : IF (nspin == 1) THEN
275 0 : WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 '
276 0 : DO i = 1, natom
277 0 : WRITE (unit_nr, '(T2,I8,T20,I10,T40,I10,T60,I10)') i, ecount(i, 1:3, 1)
278 : END DO
279 : ELSE
280 0 : WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 '
281 0 : DO i = 1, natom
282 : WRITE (unit_nr, '(T2,I8,T20,2I6,T40,2I6,T60,2I6)') &
283 0 : i, ecount(i, 1, 1:2), ecount(i, 2, 1:2), ecount(i, 3, 1:2)
284 : END DO
285 : END IF
286 : END IF
287 : ! MAO projection
288 0 : CALL para_env%sum(prmao)
289 0 : IF (unit_nr > 0) THEN
290 0 : DO ispin = 1, nspin
291 0 : WRITE (unit_nr, '(/,T2,A,I2)') 'Projection on same atom MAO orbitals: Spin ', ispin
292 0 : DO i = 1, natom, 2
293 0 : IF (i < natom) THEN
294 : WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6,T42,A,I8,T60,A,F10.6)') &
295 0 : " Atom:", i, "Projection:", prmao(i, ispin), " Atom:", i + 1, "Projection:", prmao(i + 1, ispin)
296 : ELSE
297 0 : WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6)') " Atom:", i, "Projection:", prmao(i, ispin)
298 : END IF
299 : END DO
300 : END DO
301 : END IF
302 : ! MO expansion completness
303 0 : DO ispin = 1, nspin
304 0 : CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo)
305 0 : IF (unit_nr > 0) THEN
306 0 : WRITE (unit_nr, '(/,T2,A,I2)') 'MO expansion in Localized Minimal Basis: Spin ', ispin
307 0 : WRITE (unit_nr, '(T64,A)') 'Occupied Orbitals'
308 0 : WRITE (unit_nr, '(8F10.6)') fnorm(1:homo, ispin)
309 0 : WRITE (unit_nr, '(T65,A)') 'Virtual Orbitals'
310 0 : WRITE (unit_nr, '(8F10.6)') fnorm(homo + 1:nmo, ispin)
311 : END IF
312 : END DO
313 : ! Mulliken population
314 0 : IF (unit_nr > 0) THEN
315 0 : WRITE (unit_nr, '(/,T2,A)') 'Mulliken Population Analysis '
316 : END IF
317 0 : ALLOCATE (mcharge(natom, nspin))
318 0 : DO ispin = 1, nspin
319 0 : CALL dbcsr_dot(pqmat(ispin)%matrix, sqmat(ispin)%matrix, trace)
320 0 : IF (unit_nr > 0) THEN
321 0 : WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(PS) Spin ', ispin, trace
322 : END IF
323 0 : CALL mulliken_charges(pqmat(ispin)%matrix, sqmat(ispin)%matrix, para_env, mcharge(:, ispin))
324 : END DO
325 : CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Minimal Basis Mulliken Charges", &
326 0 : electronic_charges=mcharge)
327 : ! Mayer bond orders
328 0 : IF (do_bondorder) THEN
329 0 : ALLOCATE (border(natom, natom))
330 0 : border = 0.0_dp
331 0 : CALL dbcsr_create(psmat, name="PS", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
332 0 : CALL dbcsr_create(spmat, name="SP", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
333 0 : filter_eps = 1.e-6_dp
334 0 : DO ispin = 1, nspin
335 : CALL dbcsr_multiply("N", "N", 1.0_dp, pqmat(ispin)%matrix, sqmat(ispin)%matrix, 0.0_dp, psmat, &
336 0 : filter_eps=filter_eps)
337 : CALL dbcsr_multiply("N", "N", 1.0_dp, sqmat(ispin)%matrix, pqmat(ispin)%matrix, 0.0_dp, spmat, &
338 0 : filter_eps=filter_eps)
339 0 : CALL compute_bond_order(psmat, spmat, border)
340 : END DO
341 0 : CALL para_env%sum(border)
342 0 : border = border*REAL(nspin, KIND=dp)
343 0 : CALL dbcsr_release(psmat)
344 0 : CALL dbcsr_release(spmat)
345 0 : CALL print_bond_orders(particle_set, unit_nr, border)
346 0 : DEALLOCATE (border)
347 : END IF
348 :
349 : ! for printing purposes we now copy the QUAMBOs into MO format
350 0 : ALLOCATE (mbas(nspin))
351 0 : DO ispin = 1, nspin
352 0 : CALL allocate_mo_set(mbas(ispin), nao, nmao, nmao, 0.0_dp, 1.0_dp, 0.0_dp)
353 0 : CALL set_mo_set(mbas(ispin), homo=nmao)
354 0 : ALLOCATE (mbas(ispin)%eigenvalues(nmao))
355 0 : mbas(ispin)%eigenvalues = 0.0_dp
356 0 : ALLOCATE (mbas(ispin)%occupation_numbers(nmao))
357 0 : mbas(ispin)%occupation_numbers = 1.0_dp
358 0 : CALL cp_fm_create(mbas(ispin)%mo_coeff, fm_struct_a)
359 0 : CALL copy_dbcsr_to_fm(quambo(ispin)%matrix, mbas(ispin)%mo_coeff)
360 : END DO
361 :
362 : ! Print basis functions: cube files
363 0 : DO ispin = 1, nspin
364 0 : CALL get_mo_set(mbas(ispin), mo_coeff=fm_mos)
365 0 : CALL post_minbas_cubes(qs_env, input_section, fm_mos, ispin)
366 : END DO
367 : ! Print basis functions: molden format
368 0 : CALL get_qs_env(qs_env, cell=cell)
369 0 : molden_section => section_vals_get_subs_vals(input_section, "MINBAS_MOLDEN")
370 0 : CALL write_mos_molden(mbas, qs_kind_set, particle_set, molden_section, cell=cell)
371 0 : DO ispin = 1, nspin
372 0 : CALL deallocate_mo_set(mbas(ispin))
373 : END DO
374 0 : DEALLOCATE (mbas)
375 :
376 0 : DEALLOCATE (fnorm, ecount, prmao, mcharge)
377 0 : CALL cp_fm_release(fm1)
378 0 : CALL cp_fm_release(fm2)
379 0 : CALL cp_fm_release(fm3)
380 0 : CALL cp_fm_release(fm4)
381 0 : CALL cp_fm_struct_release(fm_struct_a)
382 0 : CALL cp_fm_struct_release(fm_struct_b)
383 0 : CALL cp_fm_struct_release(fm_struct_c)
384 :
385 : ! clean up
386 0 : CALL dbcsr_deallocate_matrix_set(sqmat)
387 0 : CALL dbcsr_deallocate_matrix_set(pqmat)
388 0 : CALL dbcsr_deallocate_matrix_set(mao_coef)
389 0 : CALL dbcsr_deallocate_matrix_set(quambo)
390 :
391 : END IF
392 :
393 0 : IF (unit_nr > 0) THEN
394 : WRITE (unit_nr, '(/,T2,A)') &
395 0 : '!--------------------------END OF MINBAS ANALYSIS-----------------------------!'
396 : END IF
397 :
398 0 : CALL timestop(handle)
399 :
400 56 : END SUBROUTINE minbas_analysis
401 :
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param quambo ...
405 : !> \param smao ...
406 : !> \param ecount ...
407 : ! **************************************************************************************************
408 0 : SUBROUTINE pm_extend(quambo, smao, ecount)
409 : TYPE(dbcsr_type) :: quambo, smao
410 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: ecount
411 :
412 : INTEGER :: iatom, jatom, n
413 : LOGICAL :: found
414 : REAL(KIND=dp) :: wij
415 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: qblock, sblock
416 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
417 :
418 0 : CALL dbcsr_iterator_start(dbcsr_iter, quambo)
419 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
420 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock)
421 0 : CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, BLOCK=sblock, found=found)
422 0 : IF (found) THEN
423 0 : n = SIZE(qblock, 2)
424 0 : wij = ABS(SUM(qblock*sblock))/REAL(n, KIND=dp)
425 0 : IF (wij > 0.1_dp) THEN
426 0 : ecount(jatom, 1) = ecount(jatom, 1) + 1
427 0 : ELSEIF (wij > 0.01_dp) THEN
428 0 : ecount(jatom, 2) = ecount(jatom, 2) + 1
429 0 : ELSEIF (wij > 0.001_dp) THEN
430 0 : ecount(jatom, 3) = ecount(jatom, 3) + 1
431 : END IF
432 : END IF
433 : END DO
434 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
435 :
436 0 : END SUBROUTINE pm_extend
437 :
438 : ! **************************************************************************************************
439 : !> \brief ...
440 : !> \param mao ...
441 : !> \param smao ...
442 : !> \param sovl ...
443 : !> \param prmao ...
444 : ! **************************************************************************************************
445 0 : SUBROUTINE project_mao(mao, smao, sovl, prmao)
446 : TYPE(dbcsr_type) :: mao, smao, sovl
447 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: prmao
448 :
449 : INTEGER :: i, iatom, jatom, n
450 : LOGICAL :: found
451 : REAL(KIND=dp) :: wi
452 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: qblock, sblock, so
453 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
454 :
455 0 : CALL dbcsr_iterator_start(dbcsr_iter, mao)
456 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
457 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock)
458 0 : CPASSERT(iatom == jatom)
459 0 : CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, BLOCK=sblock, found=found)
460 0 : IF (found) THEN
461 0 : CALL dbcsr_get_block_p(matrix=sovl, row=iatom, col=jatom, BLOCK=so, found=found)
462 0 : n = SIZE(qblock, 2)
463 0 : DO i = 1, n
464 0 : wi = SUM(qblock(:, i)*sblock(:, i))
465 0 : prmao(iatom) = prmao(iatom) + wi/so(i, i)
466 : END DO
467 : END IF
468 : END DO
469 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
470 :
471 0 : END SUBROUTINE project_mao
472 :
473 : ! **************************************************************************************************
474 : !> \brief Computes and prints the Cube Files for the minimal basis set
475 : !> \param qs_env ...
476 : !> \param print_section ...
477 : !> \param minbas_coeff ...
478 : !> \param ispin ...
479 : ! **************************************************************************************************
480 0 : SUBROUTINE post_minbas_cubes(qs_env, print_section, minbas_coeff, ispin)
481 : TYPE(qs_environment_type), POINTER :: qs_env
482 : TYPE(section_vals_type), POINTER :: print_section
483 : TYPE(cp_fm_type), INTENT(IN) :: minbas_coeff
484 : INTEGER, INTENT(IN) :: ispin
485 :
486 : CHARACTER(LEN=default_path_length) :: filename, title
487 : INTEGER :: i, i_rep, ivec, iw, j, n_rep, natom
488 0 : INTEGER, DIMENSION(:), POINTER :: blk_sizes, first_bas, ilist, stride
489 : LOGICAL :: explicit, mpi_io
490 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
491 : TYPE(cell_type), POINTER :: cell
492 : TYPE(cp_logger_type), POINTER :: logger
493 : TYPE(dft_control_type), POINTER :: dft_control
494 : TYPE(particle_list_type), POINTER :: particles
495 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
496 : TYPE(pw_c1d_gs_type) :: wf_g
497 : TYPE(pw_env_type), POINTER :: pw_env
498 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
499 : TYPE(pw_r3d_rs_type) :: wf_r
500 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
501 : TYPE(qs_subsys_type), POINTER :: subsys
502 : TYPE(section_vals_type), POINTER :: minbas_section
503 :
504 0 : minbas_section => section_vals_get_subs_vals(print_section, "MINBAS_CUBE")
505 0 : CALL section_vals_get(minbas_section, explicit=explicit)
506 0 : IF (.NOT. explicit) RETURN
507 :
508 0 : logger => cp_get_default_logger()
509 0 : stride => section_get_ivals(print_section, "MINBAS_CUBE%STRIDE")
510 :
511 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
512 0 : subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
513 0 : CALL qs_subsys_get(subsys, particles=particles)
514 :
515 0 : CALL get_qs_env(qs_env=qs_env, natom=natom)
516 0 : ALLOCATE (blk_sizes(natom), first_bas(0:natom))
517 0 : CALL get_particle_set(particle_set, qs_kind_set, nmao=blk_sizes)
518 0 : first_bas(0) = 0
519 0 : DO i = 1, natom
520 0 : first_bas(i) = first_bas(i - 1) + blk_sizes(i)
521 : END DO
522 :
523 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
524 0 : CALL auxbas_pw_pool%create_pw(wf_r)
525 0 : CALL auxbas_pw_pool%create_pw(wf_g)
526 :
527 : ! loop over list of atoms
528 0 : CALL section_vals_val_get(minbas_section, "ATOM_LIST", n_rep_val=n_rep)
529 0 : IF (n_rep == 0) THEN
530 0 : DO i = 1, natom
531 0 : DO ivec = first_bas(i - 1) + 1, first_bas(i)
532 0 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin
533 0 : WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", i, " spin ", ispin
534 0 : mpi_io = .TRUE.
535 : iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", &
536 : middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., &
537 0 : mpi_io=mpi_io)
538 : CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
539 0 : cell, dft_control, particle_set, pw_env)
540 0 : CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
541 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io)
542 : END DO
543 : END DO
544 : ELSE
545 0 : DO i_rep = 1, n_rep
546 0 : CALL section_vals_val_get(minbas_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=ilist)
547 0 : DO i = 1, SIZE(ilist, 1)
548 0 : j = ilist(i)
549 0 : DO ivec = first_bas(j - 1) + 1, first_bas(j)
550 0 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin
551 0 : WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", j, " spin ", ispin
552 0 : mpi_io = .TRUE.
553 : iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", &
554 : middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., &
555 0 : mpi_io=mpi_io)
556 : CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
557 0 : cell, dft_control, particle_set, pw_env)
558 0 : CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
559 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io)
560 : END DO
561 : END DO
562 : END DO
563 : END IF
564 0 : DEALLOCATE (blk_sizes, first_bas)
565 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
566 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
567 :
568 0 : END SUBROUTINE post_minbas_cubes
569 :
570 : END MODULE minbas_wfn_analysis
|
