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