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 28419 : 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 28419 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
82 28419 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
83 :
84 2576 : SELECT CASE (gfn_type)
85 : CASE (0)
86 2576 : CPASSERT(.NOT. PRESENT(ext_ks_matrix))
87 2576 : CALL build_gfn0_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
88 : CASE (1)
89 25843 : 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 28419 : CPABORT("Unknown gfn_type")
94 : END SELECT
95 :
96 28418 : END SUBROUTINE build_xtb_ks_matrix
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param qs_env ...
101 : !> \param calculate_forces ...
102 : !> \param just_energy ...
103 : ! **************************************************************************************************
104 2576 : 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 2576 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(cp_logger_type), POINTER :: logger
115 2576 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
116 2576 : 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 2576 : 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 2576 : CALL timeset(routineN, handle)
127 :
128 : MARK_USED(calculate_forces)
129 :
130 2576 : NULLIFY (dft_control, logger, scf_section, ks_env, ks_matrix, rho, &
131 2576 : energy)
132 2576 : CPASSERT(ASSOCIATED(qs_env))
133 :
134 2576 : logger => cp_get_default_logger()
135 2576 : 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 2576 : energy=energy)
146 :
147 2576 : energy%hartree = 0.0_dp
148 2576 : energy%qmmm_el = 0.0_dp
149 :
150 2576 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
151 2576 : nspins = dft_control%nspins
152 2576 : nimg = dft_control%nimages
153 2576 : CPASSERT(ASSOCIATED(matrix_h))
154 7728 : CPASSERT(SIZE(ks_matrix) > 0)
155 :
156 5480 : DO ispin = 1, nspins
157 70376 : DO img = 1, nimg
158 : ! copy the core matrix into the fock matrix
159 67800 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
160 : END DO
161 : END DO
162 :
163 2576 : 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 2576 : energy%repulsive + energy%dispersion + energy%kTS
183 :
184 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
185 2576 : extension=".scfLog")
186 2576 : 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 2576 : "PRINT%DETAILED_ENERGY")
207 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
208 2576 : 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 2576 : CALL timestop(handle)
225 :
226 2576 : 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 25843 : 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 25843 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
251 25843 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, charges
252 25843 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
253 : TYPE(cp_logger_type), POINTER :: logger
254 25843 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs, p_matrix
255 25843 : 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 25843 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
260 : TYPE(qs_energy_type), POINTER :: energy
261 25843 : 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 25843 : CALL timeset(routineN, handle)
269 :
270 25843 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
271 25843 : ks_matrix, rho, energy)
272 25843 : CPASSERT(ASSOCIATED(qs_env))
273 :
274 25843 : logger => cp_get_default_logger()
275 25843 : 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 25843 : energy=energy)
287 :
288 25843 : 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 25843 : energy%hartree = 0.0_dp
296 25843 : energy%qmmm_el = 0.0_dp
297 25843 : energy%efield = 0.0_dp
298 :
299 25843 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
300 25843 : nspins = dft_control%nspins
301 25843 : nimg = dft_control%nimages
302 25843 : CPASSERT(ASSOCIATED(matrix_h))
303 25843 : CPASSERT(ASSOCIATED(rho))
304 77529 : CPASSERT(SIZE(ks_matrix) > 0)
305 :
306 53882 : DO ispin = 1, nspins
307 418981 : DO img = 1, nimg
308 : ! copy the core matrix into the fock matrix
309 393138 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
310 : END DO
311 : END DO
312 :
313 25843 : 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 18799 : do_efield = .FALSE.
318 : END IF
319 :
320 25843 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
321 : ! Mulliken charges
322 23737 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
323 23737 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
324 23737 : natom = SIZE(particle_set)
325 118685 : ALLOCATE (mcharge(natom), charges(natom, 5))
326 1239412 : charges = 0.0_dp
327 23737 : nkind = SIZE(atomic_kind_set)
328 23737 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
329 94948 : ALLOCATE (aocg(nsgf, natom))
330 1213687 : aocg = 0.0_dp
331 23737 : IF (nimg > 1) THEN
332 2344 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
333 : ELSE
334 21393 : p_matrix => matrix_p(:, 1)
335 21393 : s_matrix => matrix_s(1, 1)%matrix
336 21393 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
337 : END IF
338 86710 : DO ikind = 1, nkind
339 62973 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
340 62973 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
341 62973 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
342 369081 : DO iatom = 1, na
343 219398 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
344 1316388 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
345 999475 : DO is = 1, natorb
346 717104 : ns = lao(is) + 1
347 936502 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
348 : END DO
349 : END DO
350 : END DO
351 23737 : DEALLOCATE (aocg)
352 : ! charge mixing
353 23737 : IF (dft_control%qs_control%do_ls_scf) THEN
354 : !
355 : ELSE
356 23525 : 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_damping=dft_control%qs_control%xtb_control%tblite_mixer_damping, &
361 23525 : tblite_mixer_memory=qs_env%scf_control%max_scf)
362 : END IF
363 :
364 266872 : DO iatom = 1, natom
365 1342231 : mcharge(iatom) = SUM(charges(iatom, :))
366 : END DO
367 : END IF
368 :
369 25843 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
370 : CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
371 23737 : calculate_forces, just_energy)
372 : END IF
373 :
374 25843 : IF (do_efield) THEN
375 7044 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
376 : END IF
377 :
378 25843 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
379 23737 : IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
380 21811 : pass_check = .TRUE.
381 79888 : DO ikind = 1, nkind
382 58077 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
383 58077 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
384 58077 : CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
385 343737 : DO iatom = 1, na
386 205772 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
387 205772 : achrg = mcharge(atom_a)
388 263849 : IF (ABS(achrg) > chmax) THEN
389 3 : IF (iounit > 0) THEN
390 3 : WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
391 6 : " Kind Atom=", ikind, atom_a, " Limit=", chmax, " Charge=", achrg
392 : END IF
393 : pass_check = .FALSE.
394 : END IF
395 : END DO
396 : END DO
397 21811 : IF (.NOT. pass_check) THEN
398 : CALL cp_warn(__LOCATION__, "Atomic charges outside chemical range were detected."// &
399 : " Switch-off CHECK_ATOMIC_CHARGES keyword in the &xTB section"// &
400 1 : " if you want to force to continue the calculation.")
401 1 : CPABORT("xTB Charges")
402 : END IF
403 : END IF
404 : END IF
405 :
406 25842 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
407 23736 : DEALLOCATE (mcharge, charges)
408 : END IF
409 :
410 25842 : IF (qs_env%qmmm) THEN
411 5146 : CPASSERT(SIZE(ks_matrix, 2) == 1)
412 10292 : DO ispin = 1, nspins
413 : ! If QM/MM sumup the 1el Hamiltonian
414 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
415 5146 : 1.0_dp, 1.0_dp)
416 5146 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
417 : ! Compute QM/MM Energy
418 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
419 5146 : matrix_p1(ispin)%matrix, qmmm_el)
420 10292 : energy%qmmm_el = energy%qmmm_el + qmmm_el
421 : END DO
422 5146 : pc_ener = qs_env%ks_qmmm_env%pc_ener
423 5146 : energy%qmmm_el = energy%qmmm_el + pc_ener
424 : END IF
425 :
426 : energy%total = energy%core + energy%hartree + energy%efield + energy%qmmm_el + &
427 25842 : energy%repulsive + energy%dispersion + energy%dftb3 + energy%kTS
428 :
429 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
430 25842 : extension=".scfLog")
431 25842 : IF (iounit > 0) THEN
432 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
433 0 : "Repulsive pair potential energy: ", energy%repulsive, &
434 0 : "Zeroth order Hamiltonian energy: ", energy%core, &
435 0 : "Charge fluctuation energy: ", energy%hartree, &
436 0 : "London dispersion energy: ", energy%dispersion
437 0 : IF (dft_control%qs_control%xtb_control%xb_interaction) &
438 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
439 0 : "Correction for halogen bonding: ", energy%xtb_xb_inter
440 0 : IF (dft_control%qs_control%xtb_control%do_nonbonded) &
441 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
442 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
443 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
444 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
445 0 : "Electric field interaction energy: ", energy%efield
446 : END IF
447 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
448 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
449 0 : IF (qs_env%qmmm) THEN
450 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
451 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
452 : END IF
453 : END IF
454 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
455 25842 : "PRINT%DETAILED_ENERGY")
456 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
457 25842 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
458 7762 : CPASSERT(SIZE(ks_matrix, 2) == 1)
459 : BLOCK
460 7762 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
461 7762 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
462 15572 : DO ispin = 1, SIZE(mo_derivs)
463 7810 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
464 7810 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
465 0 : CPABORT("")
466 : END IF
467 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
468 15572 : 0.0_dp, mo_derivs(ispin)%matrix)
469 : END DO
470 : END BLOCK
471 : END IF
472 :
473 25842 : CALL timestop(handle)
474 :
475 51685 : END SUBROUTINE build_gfn1_xtb_ks_matrix
476 :
477 : END MODULE xtb_ks_matrix
|
