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 34880 : 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 34880 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
82 34880 : 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 31606 : 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 34880 : CPABORT("Unknown gfn_type")
94 : END SELECT
95 :
96 34880 : 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) &
194 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
195 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
196 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
197 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
198 0 : "Electric field interaction energy: ", energy%efield
199 : END IF
200 0 : IF (qs_env%qmmm) THEN
201 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
202 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
203 : END IF
204 : END IF
205 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
206 3274 : "PRINT%DETAILED_ENERGY")
207 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
208 3274 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
209 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
210 : BLOCK
211 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
212 0 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
213 0 : DO ispin = 1, SIZE(mo_derivs)
214 0 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
215 0 : CPASSERT(mo_array(ispin)%use_mo_coeff_b)
216 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
217 0 : 0.0_dp, mo_derivs(ispin)%matrix)
218 : END DO
219 : END BLOCK
220 : END IF
221 :
222 3274 : CALL timestop(handle)
223 :
224 3274 : END SUBROUTINE build_gfn0_xtb_ks_matrix
225 :
226 : ! **************************************************************************************************
227 : !> \brief ...
228 : !> \param qs_env ...
229 : !> \param calculate_forces ...
230 : !> \param just_energy ...
231 : !> \param ext_ks_matrix ...
232 : ! **************************************************************************************************
233 31606 : SUBROUTINE build_gfn1_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
234 : TYPE(qs_environment_type), POINTER :: qs_env
235 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
236 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
237 : POINTER :: ext_ks_matrix
238 :
239 : CHARACTER(len=*), PARAMETER :: routineN = 'build_gfn1_xtb_ks_matrix'
240 :
241 : INTEGER :: atom_a, handle, iatom, ikind, img, &
242 : iounit, is, ispin, na, natom, natorb, &
243 : nimg, nkind, ns, nsgf, nspins
244 : INTEGER, DIMENSION(25) :: lao
245 : INTEGER, DIMENSION(5) :: occ
246 : LOGICAL :: do_efield, pass_check
247 : REAL(KIND=dp) :: achrg, chmax, pc_ener, qmmm_el
248 31606 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
249 31606 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, charges
250 31606 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
251 : TYPE(cp_logger_type), POINTER :: logger
252 31606 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs, p_matrix
253 31606 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
254 : TYPE(dbcsr_type), POINTER :: mo_coeff, s_matrix
255 : TYPE(dft_control_type), POINTER :: dft_control
256 : TYPE(mp_para_env_type), POINTER :: para_env
257 31606 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
258 : TYPE(qs_energy_type), POINTER :: energy
259 31606 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
260 : TYPE(qs_ks_env_type), POINTER :: ks_env
261 : TYPE(qs_rho_type), POINTER :: rho
262 : TYPE(qs_scf_env_type), POINTER :: scf_env
263 : TYPE(section_vals_type), POINTER :: scf_section
264 : TYPE(xtb_atom_type), POINTER :: xtb_kind
265 :
266 31606 : CALL timeset(routineN, handle)
267 :
268 31606 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
269 31606 : ks_matrix, rho, energy)
270 31606 : CPASSERT(ASSOCIATED(qs_env))
271 :
272 31606 : logger => cp_get_default_logger()
273 31606 : iounit = cp_logger_get_default_io_unit(logger)
274 :
275 : CALL get_qs_env(qs_env, &
276 : dft_control=dft_control, &
277 : atomic_kind_set=atomic_kind_set, &
278 : qs_kind_set=qs_kind_set, &
279 : matrix_h_kp=matrix_h, &
280 : para_env=para_env, &
281 : ks_env=ks_env, &
282 : matrix_ks_kp=ks_matrix, &
283 : rho=rho, &
284 31606 : energy=energy)
285 :
286 31606 : IF (PRESENT(ext_ks_matrix)) THEN
287 : ! remap pointer to allow for non-kpoint external ks matrix
288 : ! ext_ks_matrix is used in linear response code
289 16 : ns = SIZE(ext_ks_matrix)
290 16 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
291 : END IF
292 :
293 31606 : energy%hartree = 0.0_dp
294 31606 : energy%qmmm_el = 0.0_dp
295 31606 : energy%efield = 0.0_dp
296 :
297 31606 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
298 31606 : nspins = dft_control%nspins
299 31606 : nimg = dft_control%nimages
300 31606 : CPASSERT(ASSOCIATED(matrix_h))
301 31606 : CPASSERT(ASSOCIATED(rho))
302 94818 : CPASSERT(SIZE(ks_matrix) > 0)
303 :
304 66598 : DO ispin = 1, nspins
305 501322 : DO img = 1, nimg
306 : ! copy the core matrix into the fock matrix
307 469716 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
308 : END DO
309 : END DO
310 :
311 31606 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
312 : dft_control%apply_efield_field) THEN
313 : do_efield = .TRUE.
314 : ELSE
315 24564 : do_efield = .FALSE.
316 : END IF
317 :
318 31606 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
319 : ! Mulliken charges
320 29500 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
321 29500 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
322 29500 : natom = SIZE(particle_set)
323 147500 : ALLOCATE (mcharge(natom), charges(natom, 5))
324 29500 : charges = 0.0_dp
325 29500 : nkind = SIZE(atomic_kind_set)
326 29500 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
327 118000 : ALLOCATE (aocg(nsgf, natom))
328 29500 : aocg = 0.0_dp
329 29500 : IF (nimg > 1) THEN
330 5810 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
331 : ELSE
332 23690 : p_matrix => matrix_p(:, 1)
333 23690 : s_matrix => matrix_s(1, 1)%matrix
334 23690 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
335 : END IF
336 103528 : DO ikind = 1, nkind
337 74028 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
338 74028 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
339 74028 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
340 414590 : DO iatom = 1, na
341 237034 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
342 1422204 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
343 1082982 : DO is = 1, natorb
344 771920 : ns = lao(is) + 1
345 1008954 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
346 : END DO
347 : END DO
348 : END DO
349 29500 : DEALLOCATE (aocg)
350 : ! charge mixing
351 29500 : IF (dft_control%qs_control%do_ls_scf) THEN
352 : !
353 : ELSE
354 29288 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
355 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
356 : charges, para_env, scf_env%iter_count, &
357 : scc_mixer=dft_control%qs_control%xtb_control%tblite_scc_mixer, &
358 : tblite_mixer_iterations= &
359 : dft_control%qs_control%xtb_control%tblite_mixer_iterations, &
360 : tblite_mixer_damping=dft_control%qs_control%xtb_control%tblite_mixer_damping, &
361 : tblite_mixer_memory=dft_control%qs_control%xtb_control%tblite_mixer_memory, &
362 : tblite_mixer_omega0=dft_control%qs_control%xtb_control%tblite_mixer_omega0, &
363 : tblite_mixer_min_weight= &
364 : dft_control%qs_control%xtb_control%tblite_mixer_min_weight, &
365 : tblite_mixer_max_weight= &
366 : dft_control%qs_control%xtb_control%tblite_mixer_max_weight, &
367 : tblite_mixer_weight_factor= &
368 29288 : dft_control%qs_control%xtb_control%tblite_mixer_weight_factor)
369 : END IF
370 :
371 296034 : DO iatom = 1, natom
372 1453810 : mcharge(iatom) = SUM(charges(iatom, :))
373 : END DO
374 : END IF
375 31606 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .AND. &
376 : (.NOT. dft_control%qs_control%xtb_control%do_tblite)) THEN
377 : CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
378 29500 : calculate_forces, just_energy)
379 : END IF
380 :
381 31606 : IF (do_efield) THEN
382 7042 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
383 : END IF
384 :
385 31606 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
386 29500 : IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
387 26368 : pass_check = .TRUE.
388 93100 : DO ikind = 1, nkind
389 66732 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
390 66732 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
391 66732 : CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
392 379562 : DO iatom = 1, na
393 219730 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
394 219730 : achrg = mcharge(atom_a)
395 286462 : IF (ABS(achrg) > chmax) THEN
396 0 : IF (iounit > 0) THEN
397 0 : WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
398 0 : " Kind Atom=", ikind, atom_a, " Limit=", chmax, " Charge=", achrg
399 : END IF
400 : pass_check = .FALSE.
401 : END IF
402 : END DO
403 : END DO
404 26368 : IF (.NOT. pass_check) THEN
405 : CALL cp_warn(__LOCATION__, "Atomic charges outside chemical range were detected."// &
406 : " Switch-off CHECK_ATOMIC_CHARGES keyword in the &xTB section"// &
407 0 : " if you want to force to continue the calculation.")
408 0 : CPABORT("xTB Charges")
409 : END IF
410 : END IF
411 : END IF
412 :
413 31606 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
414 29500 : DEALLOCATE (mcharge, charges)
415 : END IF
416 :
417 31606 : IF (qs_env%qmmm) THEN
418 5146 : CPASSERT(SIZE(ks_matrix, 2) == 1)
419 10292 : DO ispin = 1, nspins
420 : ! If QM/MM sumup the 1el Hamiltonian
421 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
422 5146 : 1.0_dp, 1.0_dp)
423 5146 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
424 : ! Compute QM/MM Energy
425 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
426 5146 : matrix_p1(ispin)%matrix, qmmm_el)
427 10292 : energy%qmmm_el = energy%qmmm_el + qmmm_el
428 : END DO
429 5146 : pc_ener = qs_env%ks_qmmm_env%pc_ener
430 5146 : energy%qmmm_el = energy%qmmm_el + pc_ener
431 : END IF
432 :
433 : energy%total = energy%core + energy%hartree + energy%efield + energy%qmmm_el + &
434 31606 : energy%repulsive + energy%dispersion + energy%dftb3 + energy%kTS
435 :
436 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
437 31606 : extension=".scfLog")
438 31606 : IF (iounit > 0) THEN
439 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
440 0 : "Repulsive pair potential energy: ", energy%repulsive, &
441 0 : "Zeroth order Hamiltonian energy: ", energy%core, &
442 0 : "Charge fluctuation energy: ", energy%hartree, &
443 0 : "London dispersion energy: ", energy%dispersion
444 0 : IF (dft_control%qs_control%xtb_control%xb_interaction) &
445 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
446 0 : "Correction for halogen bonding: ", energy%xtb_xb_inter
447 0 : IF (dft_control%qs_control%xtb_control%do_nonbonded) &
448 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
449 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
450 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
451 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
452 0 : "Electric field interaction energy: ", energy%efield
453 : END IF
454 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
455 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
456 0 : IF (qs_env%qmmm) THEN
457 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
458 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
459 : END IF
460 : END IF
461 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
462 31606 : "PRINT%DETAILED_ENERGY")
463 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
464 31606 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
465 7760 : CPASSERT(SIZE(ks_matrix, 2) == 1)
466 : BLOCK
467 7760 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
468 7760 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
469 15568 : DO ispin = 1, SIZE(mo_derivs)
470 7808 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
471 7808 : CPASSERT(mo_array(ispin)%use_mo_coeff_b)
472 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
473 15568 : 0.0_dp, mo_derivs(ispin)%matrix)
474 : END DO
475 : END BLOCK
476 : END IF
477 :
478 31606 : CALL timestop(handle)
479 :
480 63212 : END SUBROUTINE build_gfn1_xtb_ks_matrix
481 :
482 : END MODULE xtb_ks_matrix
|
