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 Calculate localized minimal basis
10 : !> \par History
11 : !> 12.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE minbas_methods
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_iterator_blocks_left, &
19 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
20 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks
22 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 : copy_fm_to_dbcsr,&
24 : dbcsr_allocate_matrix_set,&
25 : dbcsr_deallocate_matrix_set
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
27 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
28 : cp_fm_power
29 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30 : cp_fm_struct_release,&
31 : cp_fm_struct_type
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_get_diag,&
34 : cp_fm_release,&
35 : cp_fm_to_fm_submat,&
36 : cp_fm_type
37 : USE kinds, ONLY: dp
38 : USE mao_basis, ONLY: mao_generate_basis
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE parallel_gemm_api, ONLY: parallel_gemm
41 : USE particle_methods, ONLY: get_particle_set
42 : USE particle_types, ONLY: particle_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kind_types, ONLY: qs_kind_type
46 : USE qs_ks_types, ONLY: get_ks_env,&
47 : qs_ks_env_type
48 : USE qs_mo_types, ONLY: get_mo_set,&
49 : mo_set_type
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_methods'
56 :
57 : PUBLIC :: minbas_calculation
58 :
59 : ! **************************************************************************************************
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief ...
65 : !> \param qs_env ...
66 : !> \param mos ...
67 : !> \param quambo ...
68 : !> \param mao ...
69 : !> \param iounit ...
70 : !> \param full_ortho ...
71 : !> \param eps_filter ...
72 : ! **************************************************************************************************
73 0 : SUBROUTINE minbas_calculation(qs_env, mos, quambo, mao, iounit, full_ortho, eps_filter)
74 : TYPE(qs_environment_type), POINTER :: qs_env
75 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quambo
77 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
78 : POINTER :: mao
79 : INTEGER, INTENT(IN), OPTIONAL :: iounit
80 : LOGICAL, INTENT(IN), OPTIONAL :: full_ortho
81 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
82 :
83 : CHARACTER(len=*), PARAMETER :: routineN = 'minbas_calculation'
84 :
85 : INTEGER :: handle, homo, i, iab, ispin, nao, natom, &
86 : ndep, nmao, nmo, nmx, np, np1, nspin, &
87 : nvirt, unit_nr
88 0 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
89 : LOGICAL :: do_minbas, my_full_ortho
90 : REAL(KIND=dp) :: my_eps_filter
91 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dval, dvalo, dvalv, eigval
92 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
93 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c, &
94 : fm_struct_d, fm_struct_e
95 : TYPE(cp_fm_type) :: fm1, fm2, fm3, fm4, fm5, fm6, fma, fmb, &
96 : fmwork
97 : TYPE(cp_fm_type), POINTER :: fm_mos
98 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
99 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
100 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
101 : TYPE(dbcsr_type) :: smao, sortho
102 : TYPE(dbcsr_type), POINTER :: smat
103 : TYPE(dft_control_type), POINTER :: dft_control
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107 : TYPE(qs_ks_env_type), POINTER :: ks_env
108 :
109 0 : CALL timeset(routineN, handle)
110 :
111 0 : IF (PRESENT(iounit)) THEN
112 0 : unit_nr = iounit
113 : ELSE
114 : unit_nr = -1
115 : END IF
116 :
117 0 : IF (PRESENT(full_ortho)) THEN
118 0 : my_full_ortho = full_ortho
119 : ELSE
120 : my_full_ortho = .FALSE.
121 : END IF
122 :
123 0 : IF (PRESENT(eps_filter)) THEN
124 0 : my_eps_filter = eps_filter
125 : ELSE
126 0 : my_eps_filter = 1.0e-10_dp
127 : END IF
128 :
129 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
130 0 : nspin = dft_control%nspins
131 :
132 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
133 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
134 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
135 0 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
136 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
137 0 : CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
138 0 : nmao = SUM(col_blk_sizes)
139 : ! check if MAOs have been specified
140 0 : DO iab = 1, natom
141 0 : IF (col_blk_sizes(iab) < 0) &
142 0 : CPABORT("Number of MAOs has to be specified in KIND section for all elements")
143 : END DO
144 0 : CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=nmo)
145 :
146 0 : IF (unit_nr > 0) THEN
147 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Atomic Basis Set Functions :', nao
148 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Minimal Basis Set Functions :', nmao
149 0 : IF (nspin == 1) THEN
150 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Molecular Orbitals available :', nmo
151 : ELSE
152 0 : DO ispin = 1, nspin
153 0 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
154 : WRITE (unit_nr, '(T2,A,i2,T71,I10)') &
155 0 : 'Total Number of Molecular Orbitals available for Spin ', ispin, nmx
156 : END DO
157 : END IF
158 : END IF
159 0 : CPASSERT(nmao <= nao)
160 0 : DO ispin = 1, nspin
161 0 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
162 0 : IF (nmx /= nmo) EXIT
163 : END DO
164 0 : do_minbas = .TRUE.
165 0 : IF (nmao > nmo) THEN
166 0 : IF (unit_nr > 0) THEN
167 0 : WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
168 : END IF
169 : do_minbas = .FALSE.
170 0 : ELSEIF (nmo /= nmx) THEN
171 0 : IF (unit_nr > 0) THEN
172 0 : WRITE (unit_nr, '(T2,A)') 'Different Number of Alpha and Beta MOs'
173 0 : WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
174 : END IF
175 : do_minbas = .FALSE.
176 : ELSE
177 0 : IF (nao > nmo) THEN
178 0 : IF (unit_nr > 0) THEN
179 0 : WRITE (unit_nr, '(T2,A)') 'WARNING: Only a subset of MOs is available: Analysis depends on MOs'
180 : END IF
181 : END IF
182 : END IF
183 :
184 0 : IF (do_minbas) THEN
185 : ! initialize QUAMBOs
186 0 : NULLIFY (quambo)
187 0 : CALL dbcsr_allocate_matrix_set(quambo, nspin)
188 0 : DO ispin = 1, nspin
189 : ! coeficients
190 0 : ALLOCATE (quambo(ispin)%matrix)
191 : CALL dbcsr_create(matrix=quambo(ispin)%matrix, &
192 : name="QUAMBO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
193 0 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
194 : END DO
195 :
196 : ! initialize MAOs
197 : ! optimize MAOs (mao_coef is allocated in the routine)
198 0 : CALL mao_generate_basis(qs_env, mao_coef)
199 :
200 : ! sortho (nmao x nmao)
201 : CALL dbcsr_create(sortho, name="SORTHO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
202 0 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
203 0 : CALL dbcsr_reserve_diag_blocks(matrix=sortho)
204 :
205 0 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
206 :
207 : ! temporary FM matrices
208 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
209 0 : NULLIFY (fm_struct_a, fm_struct_b)
210 : CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, &
211 0 : para_env=para_env, context=blacs_env)
212 : CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmo, ncol_global=nmao, &
213 0 : para_env=para_env, context=blacs_env)
214 0 : CALL cp_fm_create(fm1, fm_struct_a)
215 0 : CALL cp_fm_create(fm2, fm_struct_b)
216 0 : CALL cp_fm_create(fma, fm_struct_b)
217 0 : CALL cp_fm_create(fmb, fm_struct_b)
218 :
219 0 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
220 0 : smat => matrix_s(1, 1)%matrix
221 0 : DO ispin = 1, nspin
222 :
223 : ! SMAO = Overlap*MAOs
224 0 : CALL dbcsr_create(smao, name="S*MAO", template=mao_coef(1)%matrix)
225 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef(ispin)%matrix, 0.0_dp, smao)
226 : ! a(nj)* = C(vn)(T) * SMAO(vj)
227 0 : CALL copy_dbcsr_to_fm(smao, fm1)
228 0 : CALL get_mo_set(mos(ispin), mo_coeff=fm_mos)
229 0 : CALL parallel_gemm("T", "N", nmo, nmao, nao, 1.0_dp, fm_mos, fm1, 0.0_dp, fm2)
230 0 : CALL dbcsr_release(smao)
231 0 : CALL get_mo_set(mo_set=mos(ispin), homo=homo)
232 0 : IF (unit_nr > 0) THEN
233 0 : WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Occupied Valence Set', 'Spin ', ispin, homo
234 : END IF
235 0 : nvirt = nmo - homo
236 0 : NULLIFY (fm_struct_c)
237 : CALL cp_fm_struct_create(fm_struct_c, nrow_global=nvirt, ncol_global=nvirt, &
238 0 : para_env=para_env, context=blacs_env)
239 0 : CALL cp_fm_create(fm3, fm_struct_c)
240 0 : CALL cp_fm_create(fm4, fm_struct_c)
241 : ! B(vw) = a(vj)* * a(wj)*
242 : CALL parallel_gemm("N", "T", nvirt, nvirt, nmao, 1.0_dp, fm2, fm2, 0.0_dp, fm3, &
243 0 : a_first_row=homo + 1, b_first_row=homo + 1)
244 0 : ALLOCATE (eigval(nvirt))
245 0 : CALL choose_eigv_solver(fm3, fm4, eigval)
246 : ! SVD(B) -> select p largest eigenvalues and vectors
247 0 : np = nmao - homo
248 0 : np1 = nvirt - np + 1
249 0 : IF (unit_nr > 0) THEN
250 0 : WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Virtual Valence Set', 'Spin ', ispin, np
251 : END IF
252 : ! R(vw) = SUM_p T(vp)*T(wp)
253 : CALL parallel_gemm("N", "T", nvirt, nvirt, np, 1.0_dp, fm4, fm4, 0.0_dp, fm3, &
254 0 : a_first_col=np1, b_first_col=np1)
255 : !
256 0 : ALLOCATE (dval(nmao), dvalo(nmao), dvalv(nmao))
257 0 : NULLIFY (fm_struct_d)
258 : CALL cp_fm_struct_create(fm_struct_d, nrow_global=nvirt, ncol_global=nmao, &
259 0 : para_env=para_env, context=blacs_env)
260 0 : CALL cp_fm_create(fm5, fm_struct_d)
261 0 : NULLIFY (fm_struct_e)
262 : CALL cp_fm_struct_create(fm_struct_e, nrow_global=nmao, ncol_global=nmao, &
263 0 : para_env=para_env, context=blacs_env)
264 0 : CALL cp_fm_create(fm6, fm_struct_e)
265 : ! D(j) = SUM_n (a(nj)*)^2 + SUM_vw R(vw) * a(vj)* * a(wj)*
266 : CALL parallel_gemm("N", "N", nvirt, nmao, nvirt, 1.0_dp, fm3, fm2, 0.0_dp, fm5, &
267 0 : b_first_row=homo + 1)
268 : CALL parallel_gemm("T", "N", nmao, nmao, nvirt, 1.0_dp, fm2, fm5, 0.0_dp, fm6, &
269 0 : a_first_row=homo + 1)
270 0 : CALL cp_fm_get_diag(fm6, dvalv(1:nmao))
271 0 : CALL parallel_gemm("T", "N", nmao, nmao, homo, 1.0_dp, fm2, fm2, 0.0_dp, fm6)
272 0 : CALL cp_fm_get_diag(fm6, dvalo(1:nmao))
273 0 : DO i = 1, nmao
274 0 : dval(i) = 1.0_dp/SQRT(dvalo(i) + dvalv(i))
275 : END DO
276 : ! scale intermediate expansion
277 0 : CALL cp_fm_to_fm_submat(fm2, fma, homo, nmao, 1, 1, 1, 1)
278 0 : CALL cp_fm_to_fm_submat(fm5, fma, nvirt, nmao, 1, 1, homo + 1, 1)
279 0 : CALL cp_fm_column_scale(fma, dval)
280 : ! Orthogonalization
281 0 : CALL cp_fm_create(fmwork, fm_struct_e)
282 0 : CALL parallel_gemm("T", "N", nmao, nmao, nmo, 1.0_dp, fma, fma, 0.0_dp, fm6)
283 0 : IF (my_full_ortho) THEN
284 : ! full orthogonalization
285 0 : CALL cp_fm_power(fm6, fmwork, -0.5_dp, 1.0e-12_dp, ndep)
286 0 : IF (ndep > 0 .AND. unit_nr > 0) THEN
287 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Warning: linear dependent basis ', ndep
288 : END IF
289 0 : CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
290 : ELSE
291 : ! orthogonalize on-atom blocks
292 0 : CALL copy_fm_to_dbcsr(fm6, sortho, keep_sparsity=.TRUE.)
293 0 : CALL diag_sqrt_invert(sortho)
294 0 : CALL copy_dbcsr_to_fm(sortho, fm6)
295 0 : CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
296 : END IF
297 : ! store as QUAMBO
298 0 : CALL parallel_gemm("N", "N", nao, nmao, nmo, 1.0_dp, fm_mos, fmb, 0.0_dp, fm1)
299 0 : CALL copy_fm_to_dbcsr(fm1, quambo(ispin)%matrix)
300 0 : CALL dbcsr_filter(quambo(ispin)%matrix, my_eps_filter)
301 : !
302 0 : DEALLOCATE (eigval, dval, dvalo, dvalv)
303 0 : CALL cp_fm_release(fm3)
304 0 : CALL cp_fm_release(fm4)
305 0 : CALL cp_fm_release(fm5)
306 0 : CALL cp_fm_release(fm6)
307 0 : CALL cp_fm_release(fmwork)
308 0 : CALL cp_fm_struct_release(fm_struct_c)
309 0 : CALL cp_fm_struct_release(fm_struct_d)
310 0 : CALL cp_fm_struct_release(fm_struct_e)
311 :
312 : END DO
313 :
314 : ! clean up
315 0 : CALL cp_fm_release(fm1)
316 0 : CALL cp_fm_release(fm2)
317 0 : CALL cp_fm_release(fma)
318 0 : CALL cp_fm_release(fmb)
319 0 : CALL cp_fm_struct_release(fm_struct_a)
320 0 : CALL cp_fm_struct_release(fm_struct_b)
321 0 : CALL dbcsr_release(sortho)
322 :
323 : ! return MAOs if requested
324 0 : IF (PRESENT(mao)) THEN
325 0 : mao => mao_coef
326 : ELSE
327 0 : CALL dbcsr_deallocate_matrix_set(mao_coef)
328 : END IF
329 :
330 : ELSE
331 0 : NULLIFY (quambo)
332 : END IF
333 :
334 0 : CALL timestop(handle)
335 :
336 0 : END SUBROUTINE minbas_calculation
337 :
338 : ! **************************************************************************************************
339 : !> \brief ...
340 : !> \param sortho ...
341 : ! **************************************************************************************************
342 0 : SUBROUTINE diag_sqrt_invert(sortho)
343 : TYPE(dbcsr_type) :: sortho
344 :
345 : INTEGER :: i, iatom, info, jatom, lwork, n
346 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
347 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
348 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock
349 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
350 :
351 0 : CALL dbcsr_iterator_start(dbcsr_iter, sortho)
352 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
353 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock)
354 0 : CPASSERT(iatom == jatom)
355 0 : n = SIZE(sblock, 1)
356 0 : lwork = MAX(n*n, 100)
357 0 : ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
358 0 : amat(1:n, 1:n) = sblock(1:n, 1:n)
359 0 : info = 0
360 0 : CALL dsyev("V", "U", n, amat, n, w, work, lwork, info)
361 0 : CPASSERT(info == 0)
362 0 : w(1:n) = 1._dp/SQRT(w(1:n))
363 0 : DO i = 1, n
364 0 : bmat(1:n, i) = amat(1:n, i)*w(i)
365 : END DO
366 0 : sblock(1:n, 1:n) = MATMUL(amat, TRANSPOSE(bmat))
367 0 : DEALLOCATE (amat, bmat, w, work)
368 : END DO
369 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
370 :
371 0 : END SUBROUTINE diag_sqrt_invert
372 :
373 : END MODULE minbas_methods
|
