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