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 MAO's and analyze wavefunctions
10 : !> \par History
11 : !> 03.2016 created [JGH]
12 : !> 12.2016 split into four modules [JGH]
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE mao_basis
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
18 : gto_basis_set_type
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
21 : dbcsr_distribution_type,&
22 : dbcsr_p_type,&
23 : dbcsr_type_no_symmetry
24 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks
25 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
26 : dbcsr_deallocate_matrix_set
27 : USE kinds, ONLY: dp
28 : USE mao_methods, ONLY: mao_build_q
29 : USE mao_optimizer, ONLY: mao_optimize
30 : USE particle_methods, ONLY: get_particle_set
31 : USE particle_types, ONLY: particle_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_kind_types, ONLY: get_qs_kind,&
35 : qs_kind_type
36 : USE qs_ks_types, ONLY: get_ks_env,&
37 : qs_ks_env_type
38 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
39 : release_neighbor_list_sets
40 : USE qs_neighbor_lists, ONLY: setup_neighbor_list
41 : USE qs_overlap, ONLY: build_overlap_matrix_simple
42 : USE qs_rho_types, ONLY: qs_rho_get,&
43 : qs_rho_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 : PRIVATE
48 :
49 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_basis'
50 :
51 : PUBLIC :: mao_generate_basis
52 :
53 : ! **************************************************************************************************
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param qs_env ...
60 : !> \param mao_coef ...
61 : !> \param ref_basis_set ...
62 : !> \param pmat_external ...
63 : !> \param smat_external ...
64 : !> \param molecular ...
65 : !> \param max_iter ...
66 : !> \param eps_grad ...
67 : !> \param nmao_external ...
68 : !> \param eps1_mao ...
69 : !> \param iolevel ...
70 : !> \param unit_nr ...
71 : ! **************************************************************************************************
72 4 : SUBROUTINE mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, &
73 0 : molecular, max_iter, eps_grad, nmao_external, &
74 : eps1_mao, iolevel, unit_nr)
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
77 : CHARACTER(len=*), OPTIONAL :: ref_basis_set
78 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
79 : POINTER :: pmat_external, smat_external
80 : LOGICAL, INTENT(IN), OPTIONAL :: molecular
81 : INTEGER, INTENT(IN), OPTIONAL :: max_iter
82 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_grad
83 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: nmao_external
84 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps1_mao
85 : INTEGER, INTENT(IN), OPTIONAL :: iolevel, unit_nr
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'mao_generate_basis'
88 :
89 : CHARACTER(len=10) :: mao_basis_set
90 : INTEGER :: handle, iab, iatom, ikind, iolev, ispin, &
91 : iw, mao_max_iter, natom, nbas, &
92 : nimages, nkind, nmao, nspin
93 4 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
94 : LOGICAL :: do_nmao_external, molecule
95 : REAL(KIND=dp) :: electra(2), eps1, eps_filter, eps_fun, &
96 : mao_eps_grad
97 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
98 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q, matrix_smm, matrix_smo
99 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
100 : TYPE(dft_control_type), POINTER :: dft_control
101 4 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
102 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
103 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
104 4 : POINTER :: smm_list, smo_list
105 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107 : TYPE(qs_kind_type), POINTER :: qs_kind
108 : TYPE(qs_ks_env_type), POINTER :: ks_env
109 : TYPE(qs_rho_type), POINTER :: rho
110 :
111 4 : CALL timeset(routineN, handle)
112 :
113 : ! k-points?
114 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
115 4 : nimages = dft_control%nimages
116 4 : CPASSERT(nimages == 1)
117 :
118 : ! output
119 4 : IF (PRESENT(unit_nr)) THEN
120 4 : iw = unit_nr
121 : ELSE
122 0 : iw = -1
123 : END IF
124 4 : IF (PRESENT(iolevel)) THEN
125 4 : iolev = iolevel
126 : ELSE
127 0 : iolev = 1
128 : END IF
129 4 : IF (iolevel == 0) THEN
130 0 : iw = -1
131 : END IF
132 :
133 : ! molecules
134 4 : IF (PRESENT(molecular)) THEN
135 4 : molecule = molecular
136 : ELSE
137 0 : molecule = .FALSE.
138 : END IF
139 :
140 : ! iterations
141 4 : IF (PRESENT(max_iter)) THEN
142 4 : mao_max_iter = max_iter
143 : ELSE
144 0 : mao_max_iter = 0
145 : END IF
146 :
147 : ! threshold
148 4 : IF (PRESENT(eps_grad)) THEN
149 4 : mao_eps_grad = eps_grad
150 : ELSE
151 0 : mao_eps_grad = 0.00001_dp
152 : END IF
153 4 : eps_fun = 10._dp*mao_eps_grad
154 :
155 4 : do_nmao_external = .FALSE.
156 : ! external number of MAOs per atom
157 4 : IF (PRESENT(nmao_external)) THEN
158 0 : do_nmao_external = .TRUE.
159 : END IF
160 :
161 : ! mao_threshold
162 4 : IF (PRESENT(eps1_mao)) THEN
163 4 : eps1 = eps1_mao
164 : ELSE
165 0 : eps1 = 1000._dp
166 : END IF
167 :
168 4 : IF (iw > 0) THEN
169 2 : WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
170 2 : WRITE (UNIT=iw, FMT="(T37,A)") "MAO BASIS"
171 2 : WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
172 : END IF
173 :
174 : ! Reference basis set
175 4 : IF (PRESENT(ref_basis_set)) THEN
176 4 : mao_basis_set = ref_basis_set
177 : ELSE
178 0 : mao_basis_set = "ORB"
179 : END IF
180 :
181 4 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
182 4 : nkind = SIZE(qs_kind_set)
183 36 : ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
184 12 : DO ikind = 1, nkind
185 8 : qs_kind => qs_kind_set(ikind)
186 8 : NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
187 8 : NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
188 8 : NULLIFY (basis_set_a, basis_set_b)
189 8 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="ORB")
190 8 : IF (ASSOCIATED(basis_set_a)) orb_basis_set_list(ikind)%gto_basis_set => basis_set_a
191 8 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_b, basis_type=mao_basis_set)
192 12 : IF (ASSOCIATED(basis_set_b)) mao_basis_set_list(ikind)%gto_basis_set => basis_set_b
193 : END DO
194 4 : IF (iw > 0) THEN
195 6 : DO ikind = 1, nkind
196 6 : IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
197 : WRITE (UNIT=iw, FMT="(T2,A,I4)") &
198 0 : "WARNING: No MAO basis set associated with Kind ", ikind
199 : ELSE
200 4 : IF (do_nmao_external) THEN
201 0 : nmao = nmao_external(ikind)
202 : ELSE
203 4 : CALL get_qs_kind(qs_kind_set(ikind), mao=nmao)
204 : END IF
205 4 : nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
206 : WRITE (UNIT=iw, FMT="(T2,A,I4,T30,A,I2,T56,A,I10)") &
207 4 : "MAO basis set Kind ", ikind, " MinNum of MAO:", nmao, " Number of BSF:", nbas
208 : END IF
209 : END DO
210 : END IF
211 :
212 : ! neighbor lists
213 4 : NULLIFY (smm_list, smo_list)
214 4 : CALL setup_neighbor_list(smm_list, mao_basis_set_list, molecular=molecule, qs_env=qs_env)
215 : CALL setup_neighbor_list(smo_list, mao_basis_set_list, orb_basis_set_list, &
216 4 : molecular=molecule, qs_env=qs_env)
217 :
218 : ! overlap matrices
219 4 : NULLIFY (matrix_smm, matrix_smo)
220 4 : CALL get_qs_env(qs_env, ks_env=ks_env)
221 : CALL build_overlap_matrix_simple(ks_env, matrix_smm, &
222 4 : mao_basis_set_list, mao_basis_set_list, smm_list)
223 : CALL build_overlap_matrix_simple(ks_env, matrix_smo, &
224 4 : mao_basis_set_list, orb_basis_set_list, smo_list)
225 :
226 : ! get reference density matrix and overlap matrix
227 4 : IF (PRESENT(pmat_external)) THEN
228 0 : matrix_p => pmat_external
229 : ELSE
230 4 : CALL get_qs_env(qs_env, rho=rho)
231 4 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
232 : END IF
233 4 : IF (PRESENT(smat_external)) THEN
234 0 : matrix_s => smat_external
235 : ELSE
236 4 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
237 : END IF
238 :
239 4 : nspin = SIZE(matrix_p, 1)
240 4 : eps_filter = 0.0_dp
241 : ! Q matrix
242 4 : CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter)
243 :
244 : ! MAO matrices
245 4 : CALL get_qs_env(qs_env=qs_env, natom=natom)
246 4 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
247 4 : NULLIFY (mao_coef)
248 4 : CALL dbcsr_allocate_matrix_set(mao_coef, nspin)
249 20 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
250 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
251 4 : basis=mao_basis_set_list)
252 4 : IF (do_nmao_external) THEN
253 0 : DO iatom = 1, natom
254 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
255 0 : kind_number=ikind)
256 0 : col_blk_sizes(iatom) = nmao_external(ikind)
257 : END DO
258 : ELSE
259 4 : CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
260 : END IF
261 :
262 : ! check if MAOs have been specified
263 28 : DO iab = 1, natom
264 24 : IF (col_blk_sizes(iab) < 0) &
265 4 : CPABORT("Minimum number of MAOs has to be specified in KIND section for all elements")
266 : END DO
267 8 : DO ispin = 1, nspin
268 : ! coefficients
269 4 : ALLOCATE (mao_coef(ispin)%matrix)
270 : CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, &
271 : name="MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
272 4 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
273 8 : CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
274 : END DO
275 4 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
276 :
277 : ! optimize MAOs
278 : CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, mao_max_iter, mao_eps_grad, eps1, &
279 4 : iolev, iw)
280 :
281 : ! Deallocate the neighbor list structure
282 4 : CALL release_neighbor_list_sets(smm_list)
283 4 : CALL release_neighbor_list_sets(smo_list)
284 :
285 4 : DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
286 :
287 4 : IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm)
288 4 : IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo)
289 4 : IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q)
290 :
291 4 : CALL timestop(handle)
292 :
293 4 : END SUBROUTINE mao_generate_basis
294 :
295 : END MODULE mao_basis
|
