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 Utility subroutine for qs energy calculation
10 : !> \par History
11 : !> 11.2016 split out from qs_energy_utils
12 : !> \author MK (29.10.2002)
13 : ! **************************************************************************************************
14 : MODULE qs_energy_init
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 : dbcsr_p_type,&
18 : dbcsr_set,&
19 : dbcsr_type
20 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
21 : USE efield_utils, ONLY: calculate_ecore_efield
22 : USE input_constants, ONLY: kg_tnadd_atomic,&
23 : kg_tnadd_embed,&
24 : kg_tnadd_embed_ri,&
25 : kg_tnadd_none
26 : USE input_section_types, ONLY: section_vals_type
27 : USE kg_environment, ONLY: kg_build_neighborlist,&
28 : kg_build_subsets
29 : USE kg_environment_types, ONLY: kg_environment_type
30 : USE kinds, ONLY: dp
31 : USE kpoint_methods, ONLY: kpoint_init_cell_index
32 : USE kpoint_types, ONLY: kpoint_type,&
33 : set_kpoint_info
34 : USE lri_environment_methods, ONLY: build_lri_matrices,&
35 : calculate_lri_integrals
36 : USE lri_environment_types, ONLY: lri_environment_type
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE molecule_types, ONLY: molecule_of_atom,&
39 : molecule_type
40 : USE optimize_embedding_potential, ONLY: given_embed_pot
41 : USE qs_core_energies, ONLY: calculate_ecore_overlap,&
42 : calculate_ecore_self
43 : USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix
44 : USE qs_dftb_dispersion, ONLY: calculate_dftb_dispersion
45 : USE qs_dftb_matrices, ONLY: build_dftb_matrices
46 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
47 : USE qs_dispersion_types, ONLY: qs_dispersion_type
48 : USE qs_energy_types, ONLY: qs_energy_type
49 : USE qs_environment_types, ONLY: get_qs_env,&
50 : qs_environment_type
51 : USE qs_external_density, ONLY: external_read_density
52 : USE qs_external_potential, ONLY: external_c_potential,&
53 : external_e_potential
54 : USE qs_gcp_method, ONLY: calculate_gcp_pairpot
55 : USE qs_gcp_types, ONLY: qs_gcp_type
56 : USE qs_ks_methods, ONLY: qs_ks_allocate_basics
57 : USE qs_ks_types, ONLY: get_ks_env,&
58 : qs_ks_env_type,&
59 : set_ks_env
60 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
61 : USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
62 : USE qs_update_s_mstruct, ONLY: qs_env_update_s_mstruct
63 : USE ri_environment_methods, ONLY: build_ri_matrices
64 : USE se_core_core, ONLY: se_core_core_interaction
65 : USE se_core_matrix, ONLY: build_se_core_matrix
66 : USE tblite_interface, ONLY: build_tblite_matrices
67 : USE xtb_matrices, ONLY: build_xtb_matrices
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : ! *** Global parameters ***
75 :
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_init'
77 :
78 : PUBLIC :: qs_energies_init
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief Refactoring of qs_energies_scf. Driver routine for the initial
84 : !> setup and calculations for a qs energy calculation
85 : !> \param qs_env ...
86 : !> \param calc_forces ...
87 : !> \par History
88 : !> 05.2013 created [Florian Schiffmann]
89 : ! **************************************************************************************************
90 :
91 22374 : SUBROUTINE qs_energies_init(qs_env, calc_forces)
92 : TYPE(qs_environment_type), POINTER :: qs_env
93 : LOGICAL, INTENT(IN) :: calc_forces
94 :
95 : INTEGER :: img, ispin, nimg, nspin
96 : LOGICAL :: has_unit_metric, ks_is_complex, &
97 : molecule_only
98 22374 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_w
99 : TYPE(dbcsr_type), POINTER :: matrix
100 : TYPE(dft_control_type), POINTER :: dft_control
101 : TYPE(qs_ks_env_type), POINTER :: ks_env
102 :
103 22374 : NULLIFY (ks_env, matrix_w, matrix_s, dft_control)
104 :
105 22374 : CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env)
106 22374 : IF (dft_control%qs_control%do_kg) THEN
107 112 : molecule_only = .TRUE.
108 112 : CALL qs_energies_init_kg(qs_env)
109 : ELSE
110 22262 : molecule_only = .FALSE.
111 : END IF
112 22374 : CALL qs_energies_init_hamiltonians(qs_env, calc_forces, molecule_only)
113 22374 : CALL get_ks_env(ks_env, complex_ks=ks_is_complex)
114 22374 : CALL qs_ks_allocate_basics(qs_env, is_complex=ks_is_complex)
115 :
116 : ! if need forces allocate energy weighted density matrices
117 22374 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
118 22374 : IF (calc_forces .AND. .NOT. has_unit_metric) THEN
119 : CALL get_qs_env(qs_env, &
120 : ks_env=ks_env, &
121 6367 : matrix_s_kp=matrix_s)
122 6367 : nspin = dft_control%nspins
123 6367 : nimg = dft_control%nimages
124 6367 : matrix => matrix_s(1, 1)%matrix
125 6367 : CALL dbcsr_allocate_matrix_set(matrix_w, nspin, nimg)
126 13496 : DO ispin = 1, nspin
127 45359 : DO img = 1, nimg
128 31863 : ALLOCATE (matrix_w(ispin, img)%matrix)
129 31863 : CALL dbcsr_copy(matrix_w(ispin, img)%matrix, matrix, name="W MATRIX")
130 38992 : CALL dbcsr_set(matrix_w(ispin, img)%matrix, 0.0_dp)
131 : END DO
132 : END DO
133 6367 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w)
134 : END IF
135 :
136 22374 : END SUBROUTINE qs_energies_init
137 :
138 : ! **************************************************************************************************
139 : !> \brief Refactoring of qs_energies_scf. Puts initialization of the Kim-Gordon
140 : !> settings into separate subroutine
141 : !> \param qs_env ...
142 : !> \par History
143 : !> 05.2013 created [Florian Schiffmann]
144 : ! **************************************************************************************************
145 :
146 112 : SUBROUTINE qs_energies_init_kg(qs_env)
147 : TYPE(qs_environment_type), POINTER :: qs_env
148 :
149 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_init_kg'
150 :
151 : INTEGER :: handle, isubset, natom
152 : TYPE(dft_control_type), POINTER :: dft_control
153 : TYPE(kg_environment_type), POINTER :: kg_env
154 112 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
155 : TYPE(mp_para_env_type), POINTER :: para_env
156 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
157 112 : POINTER :: soo_list, soo_list1
158 :
159 112 : CALL timeset(routineN, handle)
160 :
161 112 : CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
162 112 : CPASSERT(dft_control%qs_control%do_kg)
163 :
164 112 : kg_env => qs_env%kg_env
165 :
166 : ! get the set of molecules
167 112 : CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set, natom=natom)
168 112 : kg_env%natom = natom
169 : ! store set of molecules in kg_env
170 112 : kg_env%molecule_set => molecule_set
171 : ! build the (new) full neighborlist
172 112 : CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%sab_orb_full)
173 :
174 112 : IF (.NOT. ALLOCATED(kg_env%atom_to_molecule)) THEN
175 192 : ALLOCATE (kg_env%atom_to_molecule(natom))
176 : ! get the mapping from atoms to molecules
177 64 : CALL molecule_of_atom(molecule_set, atom_to_mol=kg_env%atom_to_molecule)
178 : END IF
179 :
180 192 : SELECT CASE (kg_env%tnadd_method)
181 : CASE (kg_tnadd_embed)
182 : ! allocate the subset list
183 80 : IF (.NOT. ASSOCIATED(kg_env%subset_of_mol)) THEN
184 144 : ALLOCATE (kg_env%subset_of_mol(SIZE(molecule_set)))
185 : END IF
186 : !
187 80 : CALL kg_build_subsets(kg_env, para_env)
188 : !
189 246 : DO isubset = 1, kg_env%nsubsets
190 : ! build the (new) molecular neighborlist of the current subset
191 : CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%subset(isubset)%sab_orb, molecular=.TRUE., &
192 246 : subset_of_mol=kg_env%subset_of_mol, current_subset=isubset)
193 : END DO
194 : CASE (kg_tnadd_embed_ri)
195 : ! should be deleted as soon as atomic grids work
196 : ! allocate the subset list
197 2 : IF (.NOT. ASSOCIATED(kg_env%subset_of_mol)) THEN
198 6 : ALLOCATE (kg_env%subset_of_mol(SIZE(molecule_set)))
199 : END IF
200 : !
201 2 : CALL kg_build_subsets(kg_env, para_env)
202 : !
203 6 : DO isubset = 1, kg_env%nsubsets
204 : ! build the (new) molecular neighborlist of the current subset
205 : CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%subset(isubset)%sab_orb, molecular=.TRUE., &
206 6 : subset_of_mol=kg_env%subset_of_mol, current_subset=isubset)
207 : END DO
208 : !
209 : ! LRI neighborlist
210 2 : NULLIFY (soo_list)
211 2 : CALL kg_build_neighborlist(qs_env, sab_orb=soo_list, molecular=.TRUE.)
212 2 : kg_env%lri_env%soo_list => soo_list
213 2 : CALL calculate_lri_integrals(kg_env%lri_env, qs_env)
214 2 : IF (qs_env%energy_correction) THEN
215 0 : NULLIFY (soo_list1)
216 0 : CALL kg_build_neighborlist(qs_env, sab_orb=soo_list1, molecular=.TRUE.)
217 0 : kg_env%lri_env1%soo_list => soo_list1
218 0 : CALL calculate_lri_integrals(kg_env%lri_env1, qs_env)
219 : END IF
220 :
221 : ! Atomic grids
222 : CASE (kg_tnadd_atomic)
223 : ! build the A-C list for the nonadditive kinetic energy potential
224 22 : CALL kg_build_neighborlist(qs_env, sac_kin=kg_env%sac_kin)
225 : CASE (kg_tnadd_none)
226 : ! nothing to do
227 : CASE DEFAULT
228 112 : CPABORT("KG:TNADD METHOD")
229 : END SELECT
230 :
231 112 : CALL timestop(handle)
232 :
233 112 : END SUBROUTINE qs_energies_init_kg
234 :
235 : ! **************************************************************************************************
236 : !> \brief Refactoring of qs_energies_scf. Moves computation of the different
237 : !> core hamiltonians into separate subroutine
238 : !> \param qs_env QS environment
239 : !> \param calc_forces Calculate forces
240 : !> \param molecule_only restrict neighbor list to molecules
241 : !> \par History
242 : !> 05.2013 created [Florian Schiffmann]
243 : !> 08.2014 Kpoints [JGH]
244 : ! **************************************************************************************************
245 :
246 22374 : SUBROUTINE qs_energies_init_hamiltonians(qs_env, calc_forces, molecule_only)
247 : TYPE(qs_environment_type), POINTER :: qs_env
248 : LOGICAL, INTENT(IN) :: calc_forces
249 : LOGICAL :: molecule_only
250 :
251 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_init_hamiltonians'
252 :
253 : INTEGER :: handle
254 : LOGICAL :: do_kpoints
255 : TYPE(dft_control_type), POINTER :: dft_control
256 : TYPE(kpoint_type), POINTER :: kpoints
257 : TYPE(lri_environment_type), POINTER :: lri_env
258 : TYPE(mp_para_env_type), POINTER :: para_env
259 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
260 22374 : POINTER :: sab_nl, sab_nl_nosym
261 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
262 : TYPE(qs_energy_type), POINTER :: energy
263 : TYPE(qs_gcp_type), POINTER :: gcp_env
264 : TYPE(section_vals_type), POINTER :: input
265 :
266 22374 : CALL timeset(routineN, handle)
267 :
268 : CALL get_qs_env(qs_env, &
269 : input=input, &
270 : dft_control=dft_control, &
271 : para_env=para_env, &
272 : kpoints=kpoints, &
273 22374 : do_kpoints=do_kpoints)
274 :
275 : ! create neighbor lists for standard use in QS
276 : CALL build_qs_neighbor_lists(qs_env, para_env, molecular=molecule_only, &
277 22374 : force_env_section=input)
278 :
279 : ! calculate cell index for k-point calculations
280 22374 : IF (do_kpoints) THEN
281 926 : CALL get_qs_env(qs_env, sab_kp=sab_nl, sab_kp_nosym=sab_nl_nosym)
282 926 : CALL kpoint_init_cell_index(kpoints, sab_nl, para_env, dft_control)
283 926 : CALL set_kpoint_info(kpoints, sab_nl_nosym=sab_nl_nosym)
284 : END IF
285 22374 : IF (dft_control%qs_control%cdft) THEN
286 326 : IF (.NOT. (dft_control%qs_control%cdft_control%external_control)) &
287 292 : dft_control%qs_control%cdft_control%need_pot = .TRUE.
288 326 : IF (ASSOCIATED(dft_control%qs_control%cdft_control%group)) THEN
289 : ! In case CDFT weight function was built beforehand (in mixed force_eval)
290 326 : IF (ASSOCIATED(dft_control%qs_control%cdft_control%group(1)%weight)) &
291 114 : dft_control%qs_control%cdft_control%need_pot = .FALSE.
292 : END IF
293 : END IF
294 :
295 : ! Calculate the overlap and the core Hamiltonian integral matrix
296 22374 : IF (dft_control%qs_control%semi_empirical) THEN
297 : CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
298 4336 : calculate_forces=.FALSE.)
299 4336 : CALL qs_env_update_s_mstruct(qs_env)
300 4336 : CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.FALSE.)
301 4336 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env, energy=energy)
302 4336 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, energy%dispersion, calc_forces)
303 18038 : ELSEIF (dft_control%qs_control%dftb) THEN
304 : CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
305 2098 : calculate_forces=.FALSE.)
306 : CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
307 2098 : calculate_forces=.FALSE.)
308 2098 : CALL qs_env_update_s_mstruct(qs_env)
309 15940 : ELSEIF (dft_control%qs_control%xtb) THEN
310 5340 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
311 656 : CALL build_tblite_matrices(qs_env=qs_env, calculate_forces=.FALSE.)
312 : ELSE
313 4684 : CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.FALSE.)
314 : END IF
315 5340 : CALL qs_env_update_s_mstruct(qs_env)
316 : ELSE
317 10600 : CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.FALSE.)
318 10600 : CALL qs_env_update_s_mstruct(qs_env)
319 10600 : CALL calculate_ecore_self(qs_env)
320 10600 : CALL calculate_ecore_efield(qs_env, calculate_forces=.FALSE.)
321 10600 : CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.FALSE.)
322 : !swap external_e_potential before external_c_potential, to ensure
323 : !that external potential on grid is loaded before calculating energy of cores
324 10600 : CALL external_e_potential(qs_env)
325 10600 : IF (.NOT. dft_control%qs_control%gapw) THEN
326 9218 : CALL external_c_potential(qs_env, calculate_forces=.FALSE.)
327 : END IF
328 : ! LRIGPW/RIGPW matrices
329 10600 : IF (dft_control%qs_control%lrigpw) THEN
330 58 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
331 58 : CALL build_lri_matrices(lri_env, qs_env)
332 10542 : ELSE IF (dft_control%qs_control%rigpw) THEN
333 0 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
334 0 : CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.FALSE.)
335 : END IF
336 :
337 : ! ZMP addition to read external density
338 10600 : CALL external_read_density(qs_env)
339 :
340 : ! Add possible pair potential dispersion energy - Evaluate first so we can print
341 : ! energy info at the end of the SCF
342 10600 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env, energy=energy)
343 10600 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, energy%dispersion, calc_forces)
344 : ! Add possible pair potential gCP energy - Evaluate first so we can print
345 : ! energy info at the end of the SCF
346 10600 : CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env, energy=energy)
347 10600 : IF (ASSOCIATED(gcp_env)) THEN
348 10600 : CALL calculate_gcp_pairpot(qs_env, gcp_env, energy%gcp, calc_forces)
349 : END IF
350 :
351 : END IF
352 : ! Embedding potential
353 22374 : IF (dft_control%qs_control%dfet_embedded) THEN
354 2 : dft_control%apply_embed_pot = .TRUE.
355 2 : CALL given_embed_pot(qs_env)
356 : END IF
357 : ! Matrix embedding potential
358 22374 : IF (dft_control%qs_control%dmfet_embedded) THEN
359 0 : dft_control%apply_dmfet_pot = .TRUE.
360 : END IF
361 :
362 22374 : CALL timestop(handle)
363 :
364 22374 : END SUBROUTINE qs_energies_init_hamiltonians
365 :
366 : END MODULE qs_energy_init
|
