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