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 Calculation of KS matrix in xTB
10 : !> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11 : !> JCTC 13, 1989-2009, (2017)
12 : !> DOI: 10.1021/acs.jctc.7b00118
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE xtb_ks_matrix
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
20 : dbcsr_copy,&
21 : dbcsr_multiply,&
22 : dbcsr_p_type,&
23 : dbcsr_type
24 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_get_default_io_unit,&
27 : cp_logger_type
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_unit_nr
30 : USE efield_tb_methods, ONLY: efield_tb_matrix
31 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
32 : section_vals_type
33 : USE kinds, ONLY: dp
34 : USE message_passing, ONLY: mp_para_env_type
35 : USE mulliken, ONLY: ao_charges
36 : USE particle_types, ONLY: particle_type
37 : USE qs_charge_mixing, ONLY: charge_mixing
38 : USE qs_energy_types, ONLY: qs_energy_type
39 : USE qs_environment_types, ONLY: get_qs_env,&
40 : qs_environment_type
41 : USE qs_kind_types, ONLY: get_qs_kind,&
42 : get_qs_kind_set,&
43 : qs_kind_type
44 : USE qs_ks_types, ONLY: qs_ks_env_type
45 : USE qs_mo_types, ONLY: get_mo_set,&
46 : mo_set_type
47 : USE qs_rho_types, ONLY: qs_rho_get,&
48 : qs_rho_type
49 : USE qs_scf_types, ONLY: qs_scf_env_type
50 : USE xtb_coulomb, ONLY: build_xtb_coulomb
51 : USE xtb_types, ONLY: get_xtb_atom_param,&
52 : xtb_atom_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ks_matrix'
60 :
61 : PUBLIC :: build_xtb_ks_matrix
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief ...
67 : !> \param qs_env ...
68 : !> \param calculate_forces ...
69 : !> \param just_energy ...
70 : !> \param ext_ks_matrix ...
71 : ! **************************************************************************************************
72 37082 : SUBROUTINE build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
73 : TYPE(qs_environment_type), POINTER :: qs_env
74 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
75 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
76 : POINTER :: ext_ks_matrix
77 :
78 : INTEGER :: gfn_type
79 : TYPE(dft_control_type), POINTER :: dft_control
80 :
81 37082 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
82 37082 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
83 :
84 3274 : SELECT CASE (gfn_type)
85 : CASE (0)
86 3274 : CPASSERT(.NOT. PRESENT(ext_ks_matrix))
87 3274 : CALL build_gfn0_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
88 : CASE (1)
89 33808 : CALL build_gfn1_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
90 : CASE (2)
91 0 : CPABORT("gfn_type = 2 not yet available")
92 : CASE DEFAULT
93 37082 : CPABORT("Unknown gfn_type")
94 : END SELECT
95 :
96 37082 : END SUBROUTINE build_xtb_ks_matrix
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param qs_env ...
101 : !> \param calculate_forces ...
102 : !> \param just_energy ...
103 : ! **************************************************************************************************
104 3274 : SUBROUTINE build_gfn0_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
107 :
108 : CHARACTER(len=*), PARAMETER :: routineN = 'build_gfn0_xtb_ks_matrix'
109 :
110 : INTEGER :: handle, img, iounit, ispin, natom, nimg, &
111 : nspins
112 : REAL(KIND=dp) :: pc_ener, qmmm_el
113 3274 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(cp_logger_type), POINTER :: logger
115 3274 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
116 3274 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h
117 : TYPE(dbcsr_type), POINTER :: mo_coeff
118 : TYPE(dft_control_type), POINTER :: dft_control
119 : TYPE(mp_para_env_type), POINTER :: para_env
120 : TYPE(qs_energy_type), POINTER :: energy
121 3274 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
122 : TYPE(qs_ks_env_type), POINTER :: ks_env
123 : TYPE(qs_rho_type), POINTER :: rho
124 : TYPE(section_vals_type), POINTER :: scf_section
125 :
126 3274 : CALL timeset(routineN, handle)
127 :
128 : MARK_USED(calculate_forces)
129 :
130 3274 : NULLIFY (dft_control, logger, scf_section, ks_env, ks_matrix, rho, &
131 3274 : energy)
132 3274 : CPASSERT(ASSOCIATED(qs_env))
133 :
134 3274 : logger => cp_get_default_logger()
135 3274 : iounit = cp_logger_get_default_io_unit(logger)
136 :
137 : CALL get_qs_env(qs_env, &
138 : dft_control=dft_control, &
139 : atomic_kind_set=atomic_kind_set, &
140 : qs_kind_set=qs_kind_set, &
141 : matrix_h_kp=matrix_h, &
142 : para_env=para_env, &
143 : ks_env=ks_env, &
144 : matrix_ks_kp=ks_matrix, &
145 3274 : energy=energy)
146 :
147 3274 : energy%hartree = 0.0_dp
148 3274 : energy%qmmm_el = 0.0_dp
149 :
150 3274 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
151 3274 : nspins = dft_control%nspins
152 3274 : nimg = dft_control%nimages
153 3274 : CPASSERT(ASSOCIATED(matrix_h))
154 9822 : CPASSERT(SIZE(ks_matrix) > 0)
155 :
156 7034 : DO ispin = 1, nspins
157 75402 : DO img = 1, nimg
158 : ! copy the core matrix into the fock matrix
159 72128 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
160 : END DO
161 : END DO
162 :
163 3274 : IF (qs_env%qmmm) THEN
164 0 : CPABORT("gfn0 QMMM NYA")
165 0 : CALL get_qs_env(qs_env=qs_env, rho=rho, natom=natom)
166 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
167 0 : DO ispin = 1, nspins
168 : ! If QM/MM sumup the 1el Hamiltonian
169 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
170 0 : 1.0_dp, 1.0_dp)
171 0 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
172 : ! Compute QM/MM Energy
173 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
174 0 : matrix_p1(ispin)%matrix, qmmm_el)
175 0 : energy%qmmm_el = energy%qmmm_el + qmmm_el
176 : END DO
177 0 : pc_ener = qs_env%ks_qmmm_env%pc_ener
178 0 : energy%qmmm_el = energy%qmmm_el + pc_ener
179 : END IF
180 :
181 : energy%total = energy%core + energy%eeq + energy%efield + energy%qmmm_el + &
182 3274 : energy%repulsive + energy%dispersion + energy%kTS
183 :
184 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
185 3274 : extension=".scfLog")
186 3274 : IF (iounit > 0) THEN
187 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
188 0 : "Repulsive pair potential energy: ", energy%repulsive, &
189 0 : "SRB Correction energy: ", energy%srb, &
190 0 : "Zeroth order Hamiltonian energy: ", energy%core, &
191 0 : "Charge equilibration energy: ", energy%eeq, &
192 0 : "London dispersion energy: ", energy%dispersion
193 0 : IF (dft_control%qs_control%xtb_control%do_nonbonded) THEN
194 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
195 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
196 : END IF
197 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
198 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
199 0 : "Electric field interaction energy: ", energy%efield
200 : END IF
201 0 : IF (qs_env%qmmm) THEN
202 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
203 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
204 : END IF
205 : END IF
206 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
207 3274 : "PRINT%DETAILED_ENERGY")
208 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
209 3274 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
210 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
211 : BLOCK
212 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
213 0 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
214 0 : DO ispin = 1, SIZE(mo_derivs)
215 0 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
216 0 : CPASSERT(mo_array(ispin)%use_mo_coeff_b)
217 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
218 0 : 0.0_dp, mo_derivs(ispin)%matrix)
219 : END DO
220 : END BLOCK
221 : END IF
222 :
223 3274 : CALL timestop(handle)
224 :
225 3274 : END SUBROUTINE build_gfn0_xtb_ks_matrix
226 :
227 : ! **************************************************************************************************
228 : !> \brief ...
229 : !> \param qs_env ...
230 : !> \param calculate_forces ...
231 : !> \param just_energy ...
232 : !> \param ext_ks_matrix ...
233 : ! **************************************************************************************************
234 33808 : SUBROUTINE build_gfn1_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
235 : TYPE(qs_environment_type), POINTER :: qs_env
236 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
237 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
238 : POINTER :: ext_ks_matrix
239 :
240 : CHARACTER(len=*), PARAMETER :: routineN = 'build_gfn1_xtb_ks_matrix'
241 :
242 : INTEGER :: atom_a, handle, iatom, ikind, img, &
243 : iounit, is, ispin, na, natom, natorb, &
244 : nimg, nkind, ns, nsgf, nspins
245 : INTEGER, DIMENSION(25) :: lao
246 : INTEGER, DIMENSION(5) :: occ
247 : LOGICAL :: do_efield, pass_check
248 : REAL(KIND=dp) :: achrg, chmax, pc_ener, qmmm_el
249 33808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
250 33808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, charges
251 33808 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
252 : TYPE(cp_logger_type), POINTER :: logger
253 33808 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs, p_matrix
254 33808 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
255 : TYPE(dbcsr_type), POINTER :: mo_coeff, s_matrix
256 : TYPE(dft_control_type), POINTER :: dft_control
257 : TYPE(mp_para_env_type), POINTER :: para_env
258 33808 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
259 : TYPE(qs_energy_type), POINTER :: energy
260 33808 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
261 : TYPE(qs_ks_env_type), POINTER :: ks_env
262 : TYPE(qs_rho_type), POINTER :: rho
263 : TYPE(qs_scf_env_type), POINTER :: scf_env
264 : TYPE(section_vals_type), POINTER :: scf_section
265 : TYPE(xtb_atom_type), POINTER :: xtb_kind
266 :
267 33808 : CALL timeset(routineN, handle)
268 :
269 33808 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
270 33808 : ks_matrix, rho, energy)
271 33808 : CPASSERT(ASSOCIATED(qs_env))
272 :
273 33808 : logger => cp_get_default_logger()
274 33808 : iounit = cp_logger_get_default_io_unit(logger)
275 :
276 : CALL get_qs_env(qs_env, &
277 : dft_control=dft_control, &
278 : atomic_kind_set=atomic_kind_set, &
279 : qs_kind_set=qs_kind_set, &
280 : matrix_h_kp=matrix_h, &
281 : para_env=para_env, &
282 : ks_env=ks_env, &
283 : matrix_ks_kp=ks_matrix, &
284 : rho=rho, &
285 33808 : energy=energy)
286 :
287 33808 : IF (PRESENT(ext_ks_matrix)) THEN
288 : ! remap pointer to allow for non-kpoint external ks matrix
289 : ! ext_ks_matrix is used in linear response code
290 16 : ns = SIZE(ext_ks_matrix)
291 16 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
292 : END IF
293 :
294 33808 : energy%hartree = 0.0_dp
295 33808 : energy%qmmm_el = 0.0_dp
296 33808 : energy%efield = 0.0_dp
297 :
298 33808 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
299 33808 : nspins = dft_control%nspins
300 33808 : nimg = dft_control%nimages
301 33808 : CPASSERT(ASSOCIATED(matrix_h))
302 33808 : CPASSERT(ASSOCIATED(rho))
303 101424 : CPASSERT(SIZE(ks_matrix) > 0)
304 :
305 73194 : DO ispin = 1, nspins
306 613332 : DO img = 1, nimg
307 : ! copy the core matrix into the fock matrix
308 579524 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
309 : END DO
310 : END DO
311 :
312 33808 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
313 : dft_control%apply_efield_field) THEN
314 : do_efield = .TRUE.
315 : ELSE
316 26390 : do_efield = .FALSE.
317 : END IF
318 :
319 33808 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
320 : ! Mulliken charges
321 31702 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
322 31702 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
323 31702 : natom = SIZE(particle_set)
324 158510 : ALLOCATE (mcharge(natom), charges(natom, 5))
325 31702 : charges = 0.0_dp
326 31702 : nkind = SIZE(atomic_kind_set)
327 31702 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
328 126808 : ALLOCATE (aocg(nsgf, natom))
329 31702 : aocg = 0.0_dp
330 31702 : IF (nimg > 1) THEN
331 6532 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
332 : ELSE
333 25170 : p_matrix => matrix_p(:, 1)
334 25170 : s_matrix => matrix_s(1, 1)%matrix
335 25170 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
336 : END IF
337 112112 : DO ikind = 1, nkind
338 80410 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
339 80410 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
340 80410 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
341 445408 : DO iatom = 1, na
342 252886 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
343 1517316 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
344 1154592 : DO is = 1, natorb
345 821296 : ns = lao(is) + 1
346 1074182 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
347 : END DO
348 : END DO
349 : END DO
350 31702 : DEALLOCATE (aocg)
351 : ! charge mixing
352 31702 : IF (dft_control%qs_control%do_ls_scf) THEN
353 : !
354 : ELSE
355 31490 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
356 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
357 : charges, para_env, scf_env%iter_count, &
358 : scc_mixer=dft_control%qs_control%xtb_control%tblite_scc_mixer, &
359 : tblite_mixer_iterations= &
360 : dft_control%qs_control%xtb_control%tblite_mixer_iterations, &
361 : tblite_mixer_damping=dft_control%qs_control%xtb_control%tblite_mixer_damping, &
362 : tblite_mixer_memory=dft_control%qs_control%xtb_control%tblite_mixer_memory, &
363 : tblite_mixer_omega0=dft_control%qs_control%xtb_control%tblite_mixer_omega0, &
364 : tblite_mixer_min_weight= &
365 : dft_control%qs_control%xtb_control%tblite_mixer_min_weight, &
366 : tblite_mixer_max_weight= &
367 : dft_control%qs_control%xtb_control%tblite_mixer_max_weight, &
368 : tblite_mixer_weight_factor= &
369 31490 : dft_control%qs_control%xtb_control%tblite_mixer_weight_factor)
370 : END IF
371 :
372 316290 : DO iatom = 1, natom
373 1551124 : mcharge(iatom) = SUM(charges(iatom, :))
374 : END DO
375 : END IF
376 33808 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .AND. &
377 : (.NOT. dft_control%qs_control%xtb_control%do_tblite)) THEN
378 : CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
379 31702 : calculate_forces, just_energy)
380 : END IF
381 :
382 33808 : IF (do_efield) THEN
383 7418 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
384 : END IF
385 :
386 33808 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
387 31702 : IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
388 26684 : pass_check = .TRUE.
389 94192 : DO ikind = 1, nkind
390 67508 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
391 67508 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
392 67508 : CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
393 382514 : DO iatom = 1, na
394 220814 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
395 220814 : achrg = mcharge(atom_a)
396 288322 : IF (ABS(achrg) > chmax) THEN
397 0 : IF (iounit > 0) THEN
398 0 : WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
399 0 : " Kind Atom=", ikind, atom_a, " Limit=", chmax, " Charge=", achrg
400 : END IF
401 : pass_check = .FALSE.
402 : END IF
403 : END DO
404 : END DO
405 26684 : IF (.NOT. pass_check) THEN
406 : CALL cp_warn(__LOCATION__, "Atomic charges outside chemical range were detected."// &
407 : " Switch-off CHECK_ATOMIC_CHARGES keyword in the &xTB section"// &
408 0 : " if you want to force to continue the calculation.")
409 0 : CPABORT("xTB Charges")
410 : END IF
411 : END IF
412 : END IF
413 :
414 33808 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
415 31702 : DEALLOCATE (mcharge, charges)
416 : END IF
417 :
418 33808 : IF (qs_env%qmmm) THEN
419 5146 : CPASSERT(SIZE(ks_matrix, 2) == 1)
420 10292 : DO ispin = 1, nspins
421 : ! If QM/MM sumup the 1el Hamiltonian
422 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
423 5146 : 1.0_dp, 1.0_dp)
424 5146 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
425 : ! Compute QM/MM Energy
426 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
427 5146 : matrix_p1(ispin)%matrix, qmmm_el)
428 10292 : energy%qmmm_el = energy%qmmm_el + qmmm_el
429 : END DO
430 5146 : pc_ener = qs_env%ks_qmmm_env%pc_ener
431 5146 : energy%qmmm_el = energy%qmmm_el + pc_ener
432 : END IF
433 :
434 : energy%total = energy%core + energy%repulsive + &
435 : energy%hartree + energy%xtb_spinpol + energy%efield + &
436 33808 : energy%qmmm_el + energy%dispersion + energy%dftb3 + energy%kTS
437 :
438 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
439 33808 : extension=".scfLog")
440 33808 : IF (iounit > 0) THEN
441 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
442 0 : "Repulsive pair potential energy: ", energy%repulsive, &
443 0 : "Zeroth order Hamiltonian energy: ", energy%core, &
444 0 : "Charge fluctuation energy: ", energy%hartree, &
445 0 : "London dispersion energy: ", energy%dispersion
446 0 : IF (dft_control%qs_control%xtb_control%do_spinpol) THEN
447 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
448 0 : "Spin polarisation correction: ", energy%xtb_spinpol
449 : END IF
450 0 : IF (dft_control%qs_control%xtb_control%xb_interaction) THEN
451 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
452 0 : "Correction for halogen bonding: ", energy%xtb_xb_inter
453 : END IF
454 0 : IF (dft_control%qs_control%xtb_control%do_nonbonded) THEN
455 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
456 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
457 : END IF
458 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
459 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
460 0 : "Electric field interaction energy: ", energy%efield
461 : END IF
462 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
463 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
464 0 : IF (qs_env%qmmm) THEN
465 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
466 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
467 : END IF
468 : END IF
469 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
470 33808 : "PRINT%DETAILED_ENERGY")
471 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
472 33808 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
473 7976 : CPASSERT(SIZE(ks_matrix, 2) == 1)
474 : BLOCK
475 7976 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
476 7976 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
477 16216 : DO ispin = 1, SIZE(mo_derivs)
478 8240 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
479 8240 : CPASSERT(mo_array(ispin)%use_mo_coeff_b)
480 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
481 16216 : 0.0_dp, mo_derivs(ispin)%matrix)
482 : END DO
483 : END BLOCK
484 : END IF
485 :
486 33808 : CALL timestop(handle)
487 :
488 67616 : END SUBROUTINE build_gfn1_xtb_ks_matrix
489 :
490 : END MODULE xtb_ks_matrix
|
