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 Overlap and Hamiltonian matrices 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_matrices
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction
18 : USE ai_overlap, ONLY: overlap_ab
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind_set
21 : USE atprop_types, ONLY: atprop_array_init,&
22 : atprop_type
23 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
24 : gto_basis_set_type
25 : USE block_p_types, ONLY: block_p_type
26 : USE cp_blacs_env, ONLY: cp_blacs_env_type
27 : USE cp_control_types, ONLY: dft_control_type,&
28 : xtb_control_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
30 : dbcsr_create,&
31 : dbcsr_finalize,&
32 : dbcsr_get_block_p,&
33 : dbcsr_p_type
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
36 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_type
39 : USE cp_output_handling, ONLY: cp_p_file,&
40 : cp_print_key_finished_output,&
41 : cp_print_key_should_output,&
42 : cp_print_key_unit_nr
43 : USE eeq_input, ONLY: eeq_solver_type
44 : USE input_constants, ONLY: vdw_pairpot_dftd4
45 : USE input_section_types, ONLY: section_vals_val_get
46 : USE kinds, ONLY: dp
47 : USE kpoint_types, ONLY: get_kpoint_info,&
48 : kpoint_type
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE orbital_pointers, ONLY: ncoset
51 : USE particle_types, ONLY: particle_type
52 : USE qs_condnum, ONLY: overlap_condnum
53 : USE qs_dispersion_cnum, ONLY: cnumber_init,&
54 : cnumber_release,&
55 : dcnum_type
56 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
57 : USE qs_dispersion_types, ONLY: qs_dispersion_type
58 : USE qs_energy_types, ONLY: qs_energy_type
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type,&
61 : set_qs_env
62 : USE qs_force_types, ONLY: qs_force_type
63 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
64 : get_memory_usage
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : qs_kind_type
67 : USE qs_ks_types, ONLY: get_ks_env,&
68 : qs_ks_env_type,&
69 : set_ks_env
70 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
71 : neighbor_list_iterate,&
72 : neighbor_list_iterator_create,&
73 : neighbor_list_iterator_p_type,&
74 : neighbor_list_iterator_release,&
75 : neighbor_list_set_p_type
76 : USE qs_overlap, ONLY: create_sab_matrix
77 : USE qs_rho_types, ONLY: qs_rho_get,&
78 : qs_rho_type
79 : USE virial_methods, ONLY: virial_pair_force
80 : USE virial_types, ONLY: virial_type
81 : USE xtb_eeq, ONLY: xtb_eeq_calculation,&
82 : xtb_eeq_forces
83 : USE xtb_hcore, ONLY: gfn0_huckel,&
84 : gfn0_kpair,&
85 : gfn1_huckel,&
86 : gfn1_kpair
87 : USE xtb_potentials, ONLY: nonbonded_correction,&
88 : repulsive_potential,&
89 : srb_potential,&
90 : xb_interaction
91 : USE xtb_types, ONLY: get_xtb_atom_param,&
92 : xtb_atom_type
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 :
97 : PRIVATE
98 :
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_matrices'
100 :
101 : PUBLIC :: build_xtb_matrices
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param qs_env ...
108 : !> \param calculate_forces ...
109 : ! **************************************************************************************************
110 5299 : SUBROUTINE build_xtb_matrices(qs_env, calculate_forces)
111 :
112 : TYPE(qs_environment_type), POINTER :: qs_env
113 : LOGICAL, INTENT(IN) :: calculate_forces
114 :
115 : INTEGER :: gfn_type
116 : TYPE(dft_control_type), POINTER :: dft_control
117 :
118 5299 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
119 5299 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
120 :
121 1680 : SELECT CASE (gfn_type)
122 : CASE (0)
123 1680 : CALL build_gfn0_xtb_matrices(qs_env, calculate_forces)
124 : CASE (1)
125 3619 : CALL build_gfn1_xtb_matrices(qs_env, calculate_forces)
126 : CASE (2)
127 0 : CPABORT("gfn_type = 2 not yet available")
128 : CASE DEFAULT
129 5299 : CPABORT("Unknown gfn_type")
130 : END SELECT
131 :
132 5299 : END SUBROUTINE build_xtb_matrices
133 :
134 : ! **************************************************************************************************
135 : !> \brief ...
136 : !> \param qs_env ...
137 : !> \param calculate_forces ...
138 : ! **************************************************************************************************
139 1680 : SUBROUTINE build_gfn0_xtb_matrices(qs_env, calculate_forces)
140 :
141 : TYPE(qs_environment_type), POINTER :: qs_env
142 : LOGICAL, INTENT(IN) :: calculate_forces
143 :
144 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn0_xtb_matrices'
145 :
146 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
147 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
148 : natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
149 : nsetb, sgfa, sgfb, za, zb
150 3360 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
151 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
152 : INTEGER, DIMENSION(3) :: cell
153 1680 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
154 1680 : npgfb, nsgfa, nsgfb
155 1680 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
156 1680 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
157 : LOGICAL :: defined, diagblock, do_nonbonded, found, &
158 : use_virial
159 : REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, eeq_energy, ef_energy, enonbonded, enscale, erep, &
160 : esrb, etaa, etab, f0, f1, f2, fhua, fhub, fhud, foab, fqa, fqb, hij, kf, qlambda, rcova, &
161 : rcovab, rcovb, rrab
162 3360 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges, qlagrange
163 3360 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, dqhuckel, huckel, owork
164 1680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
165 1680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
166 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
167 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kpolya, kpolyb, &
168 : pia, pib
169 1680 : REAL(KIND=dp), DIMENSION(:), POINTER :: eeq_q, set_radius_a, set_radius_b
170 1680 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171 1680 : scon_a, scon_b, wblock, zeta, zetb
172 1680 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173 : TYPE(atprop_type), POINTER :: atprop
174 6720 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
175 : TYPE(cp_logger_type), POINTER :: logger
176 1680 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
177 1680 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
178 : TYPE(dft_control_type), POINTER :: dft_control
179 : TYPE(eeq_solver_type) :: eeq_sparam
180 1680 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
181 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
182 : TYPE(kpoint_type), POINTER :: kpoints
183 : TYPE(mp_para_env_type), POINTER :: para_env
184 : TYPE(neighbor_list_iterator_p_type), &
185 1680 : DIMENSION(:), POINTER :: nl_iterator
186 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187 1680 : POINTER :: sab_orb, sab_xtb_nonbond
188 1680 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
190 : TYPE(qs_energy_type), POINTER :: energy
191 1680 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
192 1680 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
193 : TYPE(qs_ks_env_type), POINTER :: ks_env
194 : TYPE(qs_rho_type), POINTER :: rho
195 : TYPE(virial_type), POINTER :: virial
196 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
197 : TYPE(xtb_control_type), POINTER :: xtb_control
198 :
199 1680 : CALL timeset(routineN, handle)
200 :
201 1680 : NULLIFY (logger, virial, atprop)
202 1680 : logger => cp_get_default_logger()
203 :
204 1680 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
205 1680 : qs_kind_set, sab_orb, ks_env)
206 : CALL get_qs_env(qs_env=qs_env, &
207 : ks_env=ks_env, &
208 : energy=energy, &
209 : atomic_kind_set=atomic_kind_set, &
210 : qs_kind_set=qs_kind_set, &
211 : matrix_h_kp=matrix_h, &
212 : matrix_s_kp=matrix_s, &
213 : para_env=para_env, &
214 : atprop=atprop, &
215 : dft_control=dft_control, &
216 1680 : sab_orb=sab_orb)
217 :
218 1680 : nkind = SIZE(atomic_kind_set)
219 1680 : xtb_control => dft_control%qs_control%xtb_control
220 1680 : do_nonbonded = xtb_control%do_nonbonded
221 1680 : nimg = dft_control%nimages
222 1680 : nderivatives = 0
223 1680 : IF (calculate_forces) nderivatives = 1
224 1680 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
225 1680 : maxder = ncoset(nderivatives)
226 :
227 1680 : NULLIFY (particle_set)
228 1680 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
229 1680 : natom = SIZE(particle_set)
230 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
231 1680 : atom_of_kind=atom_of_kind, kind_of=kind_of)
232 :
233 1680 : IF (calculate_forces) THEN
234 52 : NULLIFY (rho, force, matrix_w)
235 : CALL get_qs_env(qs_env=qs_env, &
236 : rho=rho, matrix_w_kp=matrix_w, &
237 52 : virial=virial, force=force)
238 52 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
239 :
240 52 : IF (SIZE(matrix_p, 1) == 2) THEN
241 8 : DO img = 1, nimg
242 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
243 4 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
244 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
245 8 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
246 : END DO
247 : END IF
248 84 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
249 : END IF
250 : ! atomic energy decomposition
251 1680 : IF (atprop%energy) THEN
252 0 : CALL atprop_array_init(atprop%atecc, natom)
253 : END IF
254 :
255 1680 : NULLIFY (cell_to_index)
256 1680 : IF (nimg > 1) THEN
257 176 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
258 176 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
259 : END IF
260 :
261 : ! set up basis set lists
262 9158 : ALLOCATE (basis_set_list(nkind))
263 1680 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
264 :
265 : ! allocate overlap matrix
266 1680 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
267 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
268 1680 : sab_orb, .TRUE.)
269 1680 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
270 :
271 : ! initialize H matrix
272 1680 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
273 37042 : DO img = 1, nimg
274 35362 : ALLOCATE (matrix_h(1, img)%matrix)
275 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
276 35362 : name="HAMILTONIAN MATRIX")
277 37042 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
278 : END DO
279 1680 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
280 :
281 : ! Calculate coordination numbers
282 : ! needed for effective atomic energy levels
283 : ! code taken from D3 dispersion energy
284 1680 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
285 :
286 5040 : ALLOCATE (charges(natom))
287 9702 : charges = 0.0_dp
288 1680 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
289 1680 : IF (calculate_forces) THEN
290 104 : ALLOCATE (dcharges(natom))
291 332 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
292 : END IF
293 1680 : energy%eeq = eeq_energy
294 1680 : energy%efield = ef_energy
295 :
296 1680 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
297 : ! prepare charges (needed for D4)
298 1680 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
299 858 : dispersion_env%ext_charges = .TRUE.
300 858 : IF (ASSOCIATED(dispersion_env%charges)) DEALLOCATE (dispersion_env%charges)
301 1716 : ALLOCATE (dispersion_env%charges(natom))
302 4320 : dispersion_env%charges = charges
303 858 : IF (calculate_forces) THEN
304 12 : IF (ASSOCIATED(dispersion_env%dcharges)) DEALLOCATE (dispersion_env%dcharges)
305 24 : ALLOCATE (dispersion_env%dcharges(natom))
306 60 : dispersion_env%dcharges = 0.0_dp
307 : END IF
308 : END IF
309 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
310 1680 : energy%dispersion, calculate_forces)
311 1680 : IF (calculate_forces) THEN
312 52 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
313 60 : dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
314 : END IF
315 : END IF
316 :
317 : ! Calculate Huckel parameters
318 1680 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
319 :
320 : ! Calculate KAB parameters and electronegativity correction
321 1680 : CALL gfn0_kpair(qs_env, kijab)
322 :
323 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
324 1680 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
325 479095 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
326 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
327 477415 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
328 477415 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
329 477415 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
330 477415 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
331 477415 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
332 477415 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
333 477415 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
334 :
335 1909660 : dr = SQRT(SUM(rij(:)**2))
336 :
337 : ! atomic parameters
338 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
339 477415 : lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
340 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
341 477415 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
342 :
343 477415 : IF (nimg == 1) THEN
344 : ic = 1
345 : ELSE
346 256582 : ic = cell_to_index(cell(1), cell(2), cell(3))
347 256582 : CPASSERT(ic > 0)
348 : END IF
349 :
350 477415 : icol = MAX(iatom, jatom)
351 477415 : irow = MIN(iatom, jatom)
352 477415 : NULLIFY (sblock, fblock)
353 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
354 477415 : row=irow, col=icol, BLOCK=sblock, found=found)
355 477415 : CPASSERT(found)
356 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
357 477415 : row=irow, col=icol, BLOCK=fblock, found=found)
358 477415 : CPASSERT(found)
359 :
360 477415 : IF (calculate_forces) THEN
361 28264 : NULLIFY (pblock)
362 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
363 28264 : row=irow, col=icol, block=pblock, found=found)
364 28264 : CPASSERT(ASSOCIATED(pblock))
365 28264 : NULLIFY (wblock)
366 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
367 28264 : row=irow, col=icol, block=wblock, found=found)
368 28264 : CPASSERT(ASSOCIATED(wblock))
369 113056 : DO i = 2, 4
370 84792 : NULLIFY (dsblocks(i)%block)
371 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
372 84792 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
373 113056 : CPASSERT(found)
374 : END DO
375 : END IF
376 :
377 : ! overlap
378 477415 : basis_set_a => basis_set_list(ikind)%gto_basis_set
379 477415 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
380 477415 : basis_set_b => basis_set_list(jkind)%gto_basis_set
381 477415 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
382 477415 : atom_a = atom_of_kind(iatom)
383 477415 : atom_b = atom_of_kind(jatom)
384 : ! basis ikind
385 477415 : first_sgfa => basis_set_a%first_sgf
386 477415 : la_max => basis_set_a%lmax
387 477415 : la_min => basis_set_a%lmin
388 477415 : npgfa => basis_set_a%npgf
389 477415 : nseta = basis_set_a%nset
390 477415 : nsgfa => basis_set_a%nsgf_set
391 477415 : rpgfa => basis_set_a%pgf_radius
392 477415 : set_radius_a => basis_set_a%set_radius
393 477415 : scon_a => basis_set_a%scon
394 477415 : zeta => basis_set_a%zet
395 : ! basis jkind
396 477415 : first_sgfb => basis_set_b%first_sgf
397 477415 : lb_max => basis_set_b%lmax
398 477415 : lb_min => basis_set_b%lmin
399 477415 : npgfb => basis_set_b%npgf
400 477415 : nsetb = basis_set_b%nset
401 477415 : nsgfb => basis_set_b%nsgf_set
402 477415 : rpgfb => basis_set_b%pgf_radius
403 477415 : set_radius_b => basis_set_b%set_radius
404 477415 : scon_b => basis_set_b%scon
405 477415 : zetb => basis_set_b%zet
406 :
407 477415 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
408 3819320 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
409 2387075 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
410 45384203 : sint = 0.0_dp
411 :
412 1846707 : DO iset = 1, nseta
413 1369292 : ncoa = npgfa(iset)*ncoset(la_max(iset))
414 1369292 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
415 1369292 : sgfa = first_sgfa(1, iset)
416 5780564 : DO jset = 1, nsetb
417 3933857 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
418 1999913 : ncob = npgfb(jset)*ncoset(lb_max(jset))
419 1999913 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
420 1999913 : sgfb = first_sgfb(1, jset)
421 1999913 : IF (calculate_forces) THEN
422 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
423 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
424 120816 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
425 : ELSE
426 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
427 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
428 1879097 : rij, sab=oint(:, :, 1))
429 : END IF
430 : ! Contraction
431 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
432 1999913 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
433 1999913 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
434 3369205 : IF (calculate_forces) THEN
435 483264 : DO i = 2, 4
436 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
437 362448 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
438 483264 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
439 : END DO
440 : END IF
441 : END DO
442 : END DO
443 : ! forces W matrix
444 477415 : IF (calculate_forces) THEN
445 113056 : DO i = 1, 3
446 113056 : IF (iatom <= jatom) THEN
447 3988788 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
448 : ELSE
449 3122946 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
450 : END IF
451 : END DO
452 28264 : f1 = 2.0_dp
453 113056 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
454 113056 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
455 28264 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
456 27606 : IF (iatom == jatom) f1 = 1.0_dp
457 27606 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
458 : END IF
459 : END IF
460 : ! update S matrix
461 477415 : IF (iatom <= jatom) THEN
462 21412007 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
463 : ELSE
464 16519902 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
465 : END IF
466 477415 : IF (calculate_forces) THEN
467 113056 : DO i = 2, 4
468 113056 : IF (iatom <= jatom) THEN
469 3988788 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
470 : ELSE
471 3139182 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
472 : END IF
473 : END DO
474 : END IF
475 :
476 : ! Calculate Pi = Pia * Pib (Eq. 11)
477 477415 : rcovab = rcova + rcovb
478 477415 : rrab = SQRT(dr/rcovab)
479 1846707 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
480 1842067 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
481 477415 : IF (calculate_forces) THEN
482 28264 : IF (dr > 1.e-6_dp) THEN
483 28124 : drx = 0.5_dp/rrab/rcovab
484 : ELSE
485 : drx = 0.0_dp
486 : END IF
487 110764 : dpia(1:nsa) = drx*kpolya(1:nsa)
488 110672 : dpib(1:nsb) = drx*kpolyb(1:nsb)
489 : END IF
490 :
491 : ! diagonal block
492 477415 : diagblock = .FALSE.
493 477415 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
494 : !
495 : ! Eq. 10
496 : !
497 : IF (diagblock) THEN
498 24870 : DO i = 1, natorb_a
499 20859 : na = naoa(i)
500 24870 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
501 : END DO
502 : ELSE
503 4391105 : DO j = 1, natorb_b
504 3917701 : nb = naob(j)
505 37629553 : DO i = 1, natorb_a
506 33238448 : na = naoa(i)
507 33238448 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
508 37156149 : IF (iatom <= jatom) THEN
509 18728715 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
510 : ELSE
511 14509733 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
512 : END IF
513 : END DO
514 : END DO
515 : END IF
516 477415 : IF (calculate_forces) THEN
517 28264 : f0 = 1.0_dp
518 28264 : IF (irow == iatom) f0 = -1.0_dp
519 28264 : f2 = 1.0_dp
520 28264 : IF (iatom /= jatom) f2 = 2.0_dp
521 : ! Derivative wrt coordination number
522 28264 : fhua = 0.0_dp
523 28264 : fhub = 0.0_dp
524 28264 : fhud = 0.0_dp
525 28264 : fqa = 0.0_dp
526 28264 : fqb = 0.0_dp
527 28264 : IF (diagblock) THEN
528 942 : DO i = 1, natorb_a
529 802 : la = laoa(i)
530 802 : na = naoa(i)
531 802 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
532 942 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
533 : END DO
534 140 : dcharges(iatom) = dcharges(iatom) + fqa
535 : ELSE
536 269048 : DO j = 1, natorb_b
537 240924 : lb = laob(j)
538 240924 : nb = naob(j)
539 2363774 : DO i = 1, natorb_a
540 2094726 : la = laoa(i)
541 2094726 : na = naoa(i)
542 2094726 : hij = 0.5_dp*pia(na)*pib(nb)
543 2094726 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
544 2335650 : IF (iatom <= jatom) THEN
545 1170396 : fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
546 1170396 : fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
547 1170396 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
548 1170396 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
549 : ELSE
550 924330 : fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
551 924330 : fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
552 924330 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
553 924330 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
554 : END IF
555 : END DO
556 : END DO
557 28124 : dcharges(iatom) = dcharges(iatom) + fqa
558 28124 : dcharges(jatom) = dcharges(jatom) + fqb
559 : END IF
560 : ! iatom
561 28264 : atom_a = atom_of_kind(iatom)
562 441820 : DO i = 1, dcnum(iatom)%neighbors
563 413556 : katom = dcnum(iatom)%nlist(i)
564 413556 : kkind = kind_of(katom)
565 413556 : atom_c = atom_of_kind(katom)
566 1654224 : rik = dcnum(iatom)%rik(:, i)
567 1654224 : drk = SQRT(SUM(rik(:)**2))
568 441820 : IF (drk > 1.e-3_dp) THEN
569 1654224 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
570 1654224 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
571 1654224 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
572 1654224 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
573 1654224 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
574 1654224 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
575 413556 : IF (use_virial) THEN
576 1650552 : fdik = fdika + fdikb
577 412638 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
578 : END IF
579 : END IF
580 : END DO
581 : ! jatom
582 28264 : atom_b = atom_of_kind(jatom)
583 440982 : DO i = 1, dcnum(jatom)%neighbors
584 412718 : katom = dcnum(jatom)%nlist(i)
585 412718 : kkind = kind_of(katom)
586 412718 : atom_c = atom_of_kind(katom)
587 1650872 : rik = dcnum(jatom)%rik(:, i)
588 1650872 : drk = SQRT(SUM(rik(:)**2))
589 440982 : IF (drk > 1.e-3_dp) THEN
590 1650872 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
591 1650872 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
592 1650872 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
593 412718 : IF (use_virial) THEN
594 411890 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
595 : END IF
596 : END IF
597 : END DO
598 : ! force from R dendent Huckel element: Pia*Pib
599 28264 : IF (diagblock) THEN
600 140 : force_ab = 0._dp
601 : ELSE
602 28124 : n1 = SIZE(fblock, 1)
603 28124 : n2 = SIZE(fblock, 2)
604 112496 : ALLOCATE (dfblock(n1, n2))
605 2369186 : dfblock = 0.0_dp
606 269048 : DO j = 1, natorb_b
607 240924 : lb = laob(j)
608 240924 : nb = naob(j)
609 2363774 : DO i = 1, natorb_a
610 2094726 : la = laoa(i)
611 2094726 : na = naoa(i)
612 2094726 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
613 2335650 : IF (iatom <= jatom) THEN
614 1170396 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
615 : ELSE
616 924330 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
617 : END IF
618 : END DO
619 : END DO
620 2369186 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
621 112496 : DO ir = 1, 3
622 84372 : foab = 2.0_dp*dfp*rij(ir)/dr
623 : ! force from overlap matrix contribution to H
624 807144 : DO j = 1, natorb_b
625 722772 : lb = laob(j)
626 722772 : nb = naob(j)
627 7091322 : DO i = 1, natorb_a
628 6284178 : la = laoa(i)
629 6284178 : na = naoa(i)
630 6284178 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
631 7006950 : IF (iatom <= jatom) THEN
632 3511188 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
633 : ELSE
634 2772990 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
635 : END IF
636 : END DO
637 : END DO
638 112496 : force_ab(ir) = foab
639 : END DO
640 28124 : DEALLOCATE (dfblock)
641 : END IF
642 : END IF
643 :
644 477415 : IF (calculate_forces) THEN
645 28264 : atom_a = atom_of_kind(iatom)
646 28264 : atom_b = atom_of_kind(jatom)
647 75556 : IF (irow == iatom) force_ab = -force_ab
648 113056 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
649 113056 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
650 28264 : IF (use_virial) THEN
651 27682 : f1 = 1.0_dp
652 27682 : IF (iatom == jatom) f1 = 0.5_dp
653 27682 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
654 : END IF
655 : END IF
656 :
657 2866170 : DEALLOCATE (oint, owork, sint)
658 :
659 : END DO
660 1680 : CALL neighbor_list_iterator_release(nl_iterator)
661 :
662 3360 : DO i = 1, SIZE(matrix_h, 1)
663 38722 : DO img = 1, nimg
664 35362 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
665 37042 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
666 : END DO
667 : END DO
668 :
669 : ! EEQ forces (response and direct)
670 1680 : IF (calculate_forces) THEN
671 52 : CALL para_env%sum(dcharges)
672 104 : ALLOCATE (qlagrange(natom))
673 52 : CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
674 : END IF
675 :
676 1680 : kf = xtb_control%kf
677 1680 : enscale = xtb_control%enscale
678 1680 : erep = 0.0_dp
679 1680 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
680 :
681 1680 : esrb = 0.0_dp
682 1680 : CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
683 :
684 1680 : enonbonded = 0.0_dp
685 1680 : IF (do_nonbonded) THEN
686 : ! nonbonded interactions
687 0 : NULLIFY (sab_xtb_nonbond)
688 0 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
689 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
690 0 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
691 : END IF
692 :
693 : ! set repulsive energy
694 1680 : erep = erep + esrb + enonbonded
695 1680 : IF (do_nonbonded) THEN
696 0 : CALL para_env%sum(enonbonded)
697 0 : energy%xtb_nonbonded = enonbonded
698 : END IF
699 1680 : CALL para_env%sum(esrb)
700 1680 : energy%srb = esrb
701 1680 : CALL para_env%sum(erep)
702 1680 : energy%repulsive = erep
703 :
704 : ! save EEQ charges
705 1680 : NULLIFY (eeq_q)
706 1680 : CALL get_qs_env(qs_env, eeq=eeq_q)
707 1680 : IF (ASSOCIATED(eeq_q)) THEN
708 996 : CPASSERT(SIZE(eeq_q) == natom)
709 : ELSE
710 1368 : ALLOCATE (eeq_q(natom))
711 3522 : eeq_q(1:natom) = charges(1:natom)
712 : END IF
713 1680 : CALL set_qs_env(qs_env, eeq=eeq_q)
714 :
715 : ! deallocate coordination numbers
716 1680 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
717 :
718 : ! deallocate Huckel parameters
719 1680 : DEALLOCATE (huckel)
720 1680 : IF (calculate_forces) THEN
721 52 : DEALLOCATE (dhuckel, dqhuckel)
722 : END IF
723 : ! deallocate KAB parameters
724 1680 : DEALLOCATE (kijab)
725 :
726 : ! deallocate charges
727 1680 : DEALLOCATE (charges)
728 1680 : IF (calculate_forces) THEN
729 52 : DEALLOCATE (dcharges, qlagrange)
730 : END IF
731 :
732 : ! AO matrix outputs
733 1680 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
734 :
735 1680 : DEALLOCATE (basis_set_list)
736 1680 : IF (calculate_forces) THEN
737 52 : IF (SIZE(matrix_p, 1) == 2) THEN
738 8 : DO img = 1, nimg
739 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
740 8 : beta_scalar=-1.0_dp)
741 : END DO
742 : END IF
743 : END IF
744 :
745 1680 : CALL timestop(handle)
746 :
747 5040 : END SUBROUTINE build_gfn0_xtb_matrices
748 :
749 : ! **************************************************************************************************
750 : !> \brief ...
751 : !> \param qs_env ...
752 : !> \param calculate_forces ...
753 : ! **************************************************************************************************
754 3619 : SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
755 :
756 : TYPE(qs_environment_type), POINTER :: qs_env
757 : LOGICAL, INTENT(IN) :: calculate_forces
758 :
759 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn1_xtb_matrices'
760 :
761 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
762 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
763 : natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
764 : nsetb, sgfa, sgfb, za, zb
765 7238 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
766 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
767 : INTEGER, DIMENSION(3) :: cell
768 3619 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
769 3619 : npgfb, nsgfa, nsgfb
770 3619 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
771 3619 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
772 : LOGICAL :: defined, diagblock, do_nonbonded, found, &
773 : use_virial, xb_inter
774 : REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, enonbonded, &
775 : enscale, erep, etaa, etab, exb, f0, &
776 : f1, fhua, fhub, fhud, foab, hij, kf, &
777 : rcova, rcovab, rcovb, rrab
778 3619 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
779 7238 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
780 3619 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
781 3619 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
782 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
783 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
784 3619 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
785 3619 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
786 3619 : scon_a, scon_b, wblock, zeta, zetb
787 3619 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
788 : TYPE(atprop_type), POINTER :: atprop
789 14476 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
790 : TYPE(cp_logger_type), POINTER :: logger
791 3619 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
792 3619 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
793 : TYPE(dft_control_type), POINTER :: dft_control
794 3619 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
795 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
796 : TYPE(kpoint_type), POINTER :: kpoints
797 : TYPE(mp_para_env_type), POINTER :: para_env
798 : TYPE(neighbor_list_iterator_p_type), &
799 3619 : DIMENSION(:), POINTER :: nl_iterator
800 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
801 3619 : POINTER :: sab_orb, sab_xtb_nonbond
802 3619 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
803 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
804 : TYPE(qs_energy_type), POINTER :: energy
805 3619 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
806 3619 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
807 : TYPE(qs_ks_env_type), POINTER :: ks_env
808 : TYPE(qs_rho_type), POINTER :: rho
809 : TYPE(virial_type), POINTER :: virial
810 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
811 : TYPE(xtb_control_type), POINTER :: xtb_control
812 :
813 3619 : CALL timeset(routineN, handle)
814 :
815 3619 : NULLIFY (logger, virial, atprop)
816 3619 : logger => cp_get_default_logger()
817 :
818 3619 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
819 3619 : qs_kind_set, sab_orb, ks_env)
820 :
821 : CALL get_qs_env(qs_env=qs_env, &
822 : ks_env=ks_env, &
823 : energy=energy, &
824 : atomic_kind_set=atomic_kind_set, &
825 : qs_kind_set=qs_kind_set, &
826 : matrix_h_kp=matrix_h, &
827 : matrix_s_kp=matrix_s, &
828 : para_env=para_env, &
829 : atprop=atprop, &
830 : dft_control=dft_control, &
831 3619 : sab_orb=sab_orb)
832 :
833 3619 : nkind = SIZE(atomic_kind_set)
834 3619 : xtb_control => dft_control%qs_control%xtb_control
835 3619 : xb_inter = xtb_control%xb_interaction
836 3619 : do_nonbonded = xtb_control%do_nonbonded
837 3619 : nimg = dft_control%nimages
838 3619 : nderivatives = 0
839 3619 : IF (calculate_forces) nderivatives = 1
840 3619 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
841 3619 : maxder = ncoset(nderivatives)
842 :
843 3619 : NULLIFY (particle_set)
844 3619 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
845 3619 : natom = SIZE(particle_set)
846 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
847 3619 : atom_of_kind=atom_of_kind, kind_of=kind_of)
848 :
849 3619 : IF (calculate_forces) THEN
850 494 : NULLIFY (rho, force, matrix_w)
851 : CALL get_qs_env(qs_env=qs_env, &
852 : rho=rho, matrix_w_kp=matrix_w, &
853 494 : virial=virial, force=force)
854 494 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
855 :
856 494 : IF (SIZE(matrix_p, 1) == 2) THEN
857 432 : DO img = 1, nimg
858 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
859 414 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
860 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
861 432 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
862 : END DO
863 : END IF
864 876 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
865 : END IF
866 : ! atomic energy decomposition
867 3619 : IF (atprop%energy) THEN
868 36 : CALL atprop_array_init(atprop%atecc, natom)
869 : END IF
870 :
871 3619 : NULLIFY (cell_to_index)
872 3619 : IF (nimg > 1) THEN
873 340 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
874 340 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
875 : END IF
876 :
877 : ! set up basis set lists
878 19814 : ALLOCATE (basis_set_list(nkind))
879 3619 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
880 :
881 : ! allocate overlap matrix
882 3619 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
883 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
884 3619 : sab_orb, .TRUE.)
885 3619 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
886 :
887 : ! initialize H matrix
888 3619 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
889 51550 : DO img = 1, nimg
890 47931 : ALLOCATE (matrix_h(1, img)%matrix)
891 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
892 47931 : name="HAMILTONIAN MATRIX")
893 51550 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
894 : END DO
895 3619 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
896 :
897 : ! Calculate coordination numbers
898 : ! needed for effective atomic energy levels (Eq. 12)
899 : ! code taken from D3 dispersion energy
900 3619 : CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
901 :
902 : ! vdW Potential
903 3619 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
904 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
905 3619 : energy%dispersion, calculate_forces)
906 :
907 : ! Calculate Huckel parameters
908 3619 : CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
909 :
910 : ! Calculate KAB parameters and electronegativity correction
911 3619 : CALL gfn1_kpair(qs_env, kijab)
912 :
913 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
914 3619 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
915 1087558 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
916 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
917 1083939 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
918 1083939 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
919 1083939 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
920 1083939 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
921 1083939 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
922 1083939 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
923 1083939 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
924 :
925 4335756 : dr = SQRT(SUM(rij(:)**2))
926 :
927 : ! atomic parameters
928 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
929 1083939 : lmax=lmaxa, nshell=nsa, kpoly=kpolya)
930 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
931 1083939 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
932 :
933 1083939 : IF (nimg == 1) THEN
934 : ic = 1
935 : ELSE
936 247218 : ic = cell_to_index(cell(1), cell(2), cell(3))
937 247218 : CPASSERT(ic > 0)
938 : END IF
939 :
940 1083939 : icol = MAX(iatom, jatom)
941 1083939 : irow = MIN(iatom, jatom)
942 1083939 : NULLIFY (sblock, fblock)
943 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
944 1083939 : row=irow, col=icol, BLOCK=sblock, found=found)
945 1083939 : CPASSERT(found)
946 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
947 1083939 : row=irow, col=icol, BLOCK=fblock, found=found)
948 1083939 : CPASSERT(found)
949 :
950 1083939 : IF (calculate_forces) THEN
951 242583 : NULLIFY (pblock)
952 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
953 242583 : row=irow, col=icol, block=pblock, found=found)
954 242583 : CPASSERT(found)
955 242583 : NULLIFY (wblock)
956 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
957 242583 : row=irow, col=icol, block=wblock, found=found)
958 242583 : CPASSERT(found)
959 970332 : DO i = 2, 4
960 727749 : NULLIFY (dsblocks(i)%block)
961 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
962 727749 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
963 970332 : CPASSERT(found)
964 : END DO
965 : END IF
966 :
967 : ! overlap
968 1083939 : basis_set_a => basis_set_list(ikind)%gto_basis_set
969 1083939 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
970 1083939 : basis_set_b => basis_set_list(jkind)%gto_basis_set
971 1083939 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
972 1083939 : atom_a = atom_of_kind(iatom)
973 1083939 : atom_b = atom_of_kind(jatom)
974 : ! basis ikind
975 1083939 : first_sgfa => basis_set_a%first_sgf
976 1083939 : la_max => basis_set_a%lmax
977 1083939 : la_min => basis_set_a%lmin
978 1083939 : npgfa => basis_set_a%npgf
979 1083939 : nseta = basis_set_a%nset
980 1083939 : nsgfa => basis_set_a%nsgf_set
981 1083939 : rpgfa => basis_set_a%pgf_radius
982 1083939 : set_radius_a => basis_set_a%set_radius
983 1083939 : scon_a => basis_set_a%scon
984 1083939 : zeta => basis_set_a%zet
985 : ! basis jkind
986 1083939 : first_sgfb => basis_set_b%first_sgf
987 1083939 : lb_max => basis_set_b%lmax
988 1083939 : lb_min => basis_set_b%lmin
989 1083939 : npgfb => basis_set_b%npgf
990 1083939 : nsetb = basis_set_b%nset
991 1083939 : nsgfb => basis_set_b%nsgf_set
992 1083939 : rpgfb => basis_set_b%pgf_radius
993 1083939 : set_radius_b => basis_set_b%set_radius
994 1083939 : scon_b => basis_set_b%scon
995 1083939 : zetb => basis_set_b%zet
996 :
997 1083939 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
998 8671512 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
999 5419695 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
1000 53061798 : sint = 0.0_dp
1001 :
1002 3432273 : DO iset = 1, nseta
1003 2348334 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1004 2348334 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1005 2348334 : sgfa = first_sgfa(1, iset)
1006 8645535 : DO jset = 1, nsetb
1007 5213262 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1008 4178179 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1009 4178179 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1010 4178179 : sgfb = first_sgfb(1, jset)
1011 4178179 : IF (calculate_forces) THEN
1012 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1013 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1014 952840 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1015 : ELSE
1016 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1017 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1018 3225339 : rij, sab=oint(:, :, 1))
1019 : END IF
1020 : ! Contraction
1021 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1022 4178179 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1023 4178179 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1024 6526513 : IF (calculate_forces) THEN
1025 3811360 : DO i = 2, 4
1026 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1027 2858520 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1028 3811360 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1029 : END DO
1030 : END IF
1031 : END DO
1032 : END DO
1033 : ! forces W matrix
1034 1083939 : IF (calculate_forces) THEN
1035 970332 : DO i = 1, 3
1036 970332 : IF (iatom <= jatom) THEN
1037 14172753 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
1038 : ELSE
1039 10114419 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
1040 : END IF
1041 : END DO
1042 242583 : f1 = 2.0_dp
1043 970332 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1044 970332 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1045 242583 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
1046 148835 : IF (iatom == jatom) f1 = 1.0_dp
1047 148835 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1048 : END IF
1049 : END IF
1050 : ! update S matrix
1051 1083939 : IF (iatom <= jatom) THEN
1052 14727940 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1053 : ELSE
1054 10819659 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1055 : END IF
1056 1083939 : IF (calculate_forces) THEN
1057 970332 : DO i = 2, 4
1058 970332 : IF (iatom <= jatom) THEN
1059 14172753 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1060 : ELSE
1061 9948399 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
1062 : END IF
1063 : END DO
1064 : END IF
1065 :
1066 : ! Calculate Pi = Pia * Pib (Eq. 11)
1067 1083939 : rcovab = rcova + rcovb
1068 1083939 : rrab = SQRT(dr/rcovab)
1069 3432273 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1070 3432278 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1071 1083939 : IF (calculate_forces) THEN
1072 242583 : IF (dr > 1.e-6_dp) THEN
1073 240298 : drx = 0.5_dp/rrab/rcovab
1074 : ELSE
1075 : drx = 0.0_dp
1076 : END IF
1077 799127 : dpia(1:nsa) = drx*kpolya(1:nsa)
1078 799129 : dpib(1:nsb) = drx*kpolyb(1:nsb)
1079 : END IF
1080 :
1081 : ! diagonal block
1082 1083939 : diagblock = .FALSE.
1083 1083939 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
1084 : !
1085 : ! Eq. 10
1086 : !
1087 : IF (diagblock) THEN
1088 64988 : DO i = 1, natorb_a
1089 49386 : na = naoa(i)
1090 64988 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1091 : END DO
1092 : ELSE
1093 5062814 : DO j = 1, natorb_b
1094 3994477 : nb = naob(j)
1095 25509218 : DO i = 1, natorb_a
1096 20446404 : na = naoa(i)
1097 20446404 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1098 24440881 : IF (iatom <= jatom) THEN
1099 11808962 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1100 : ELSE
1101 8637442 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1102 : END IF
1103 : END DO
1104 : END DO
1105 : END IF
1106 1083939 : IF (calculate_forces) THEN
1107 242583 : f0 = 1.0_dp
1108 242583 : IF (irow == iatom) f0 = -1.0_dp
1109 : ! Derivative wrt coordination number
1110 242583 : fhua = 0.0_dp
1111 242583 : fhub = 0.0_dp
1112 242583 : fhud = 0.0_dp
1113 242583 : IF (diagblock) THEN
1114 9903 : DO i = 1, natorb_a
1115 7618 : la = laoa(i)
1116 7618 : na = naoa(i)
1117 9903 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1118 : END DO
1119 : ELSE
1120 1325006 : DO j = 1, natorb_b
1121 1084708 : lb = laob(j)
1122 1084708 : nb = naob(j)
1123 8051323 : DO i = 1, natorb_a
1124 6726317 : la = laoa(i)
1125 6726317 : na = naoa(i)
1126 6726317 : hij = 0.5_dp*pia(na)*pib(nb)
1127 7811025 : IF (iatom <= jatom) THEN
1128 3968304 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1129 3968304 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1130 : ELSE
1131 2758013 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1132 2758013 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
1133 : END IF
1134 : END DO
1135 : END DO
1136 240298 : IF (iatom /= jatom) THEN
1137 207280 : fhua = 2.0_dp*fhua
1138 207280 : fhub = 2.0_dp*fhub
1139 : END IF
1140 : END IF
1141 : ! iatom
1142 242583 : atom_a = atom_of_kind(iatom)
1143 7931791 : DO i = 1, dcnum(iatom)%neighbors
1144 7689208 : katom = dcnum(iatom)%nlist(i)
1145 7689208 : kkind = kind_of(katom)
1146 7689208 : atom_c = atom_of_kind(katom)
1147 30756832 : rik = dcnum(iatom)%rik(:, i)
1148 30756832 : drk = SQRT(SUM(rik(:)**2))
1149 7931791 : IF (drk > 1.e-3_dp) THEN
1150 30756832 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1151 30756832 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1152 30756832 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1153 30756832 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1154 30756832 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1155 30756832 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1156 7689208 : IF (use_virial) THEN
1157 20584936 : fdik = fdika + fdikb
1158 5146234 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1159 : END IF
1160 : END IF
1161 : END DO
1162 : ! jatom
1163 242583 : atom_b = atom_of_kind(jatom)
1164 7928926 : DO i = 1, dcnum(jatom)%neighbors
1165 7686343 : katom = dcnum(jatom)%nlist(i)
1166 7686343 : kkind = kind_of(katom)
1167 7686343 : atom_c = atom_of_kind(katom)
1168 30745372 : rik = dcnum(jatom)%rik(:, i)
1169 30745372 : drk = SQRT(SUM(rik(:)**2))
1170 7928926 : IF (drk > 1.e-3_dp) THEN
1171 30745372 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1172 30745372 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1173 30745372 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1174 7686343 : IF (use_virial) THEN
1175 5145160 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1176 : END IF
1177 : END IF
1178 : END DO
1179 : ! force from R dendent Huckel element: Pia*Pib
1180 242583 : IF (diagblock) THEN
1181 2285 : force_ab = 0._dp
1182 : ELSE
1183 240298 : n1 = SIZE(fblock, 1)
1184 240298 : n2 = SIZE(fblock, 2)
1185 961192 : ALLOCATE (dfblock(n1, n2))
1186 7995983 : dfblock = 0.0_dp
1187 1325006 : DO j = 1, natorb_b
1188 1084708 : lb = laob(j)
1189 1084708 : nb = naob(j)
1190 8051323 : DO i = 1, natorb_a
1191 6726317 : la = laoa(i)
1192 6726317 : na = naoa(i)
1193 6726317 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1194 7811025 : IF (iatom <= jatom) THEN
1195 3968304 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1196 : ELSE
1197 2758013 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1198 : END IF
1199 : END DO
1200 : END DO
1201 7995983 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
1202 961192 : DO ir = 1, 3
1203 720894 : foab = 2.0_dp*dfp*rij(ir)/dr
1204 : ! force from overlap matrix contribution to H
1205 3975018 : DO j = 1, natorb_b
1206 3254124 : lb = laob(j)
1207 3254124 : nb = naob(j)
1208 24153969 : DO i = 1, natorb_a
1209 20178951 : la = laoa(i)
1210 20178951 : na = naoa(i)
1211 20178951 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1212 23433075 : IF (iatom <= jatom) THEN
1213 11904912 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1214 : ELSE
1215 8274039 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1216 : END IF
1217 : END DO
1218 : END DO
1219 961192 : force_ab(ir) = foab
1220 : END DO
1221 240298 : DEALLOCATE (dfblock)
1222 : END IF
1223 : END IF
1224 :
1225 1083939 : IF (calculate_forces) THEN
1226 242583 : atom_a = atom_of_kind(iatom)
1227 242583 : atom_b = atom_of_kind(jatom)
1228 634383 : IF (irow == iatom) force_ab = -force_ab
1229 970332 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1230 970332 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1231 242583 : IF (use_virial) THEN
1232 149691 : f1 = 1.0_dp
1233 149691 : IF (iatom == jatom) f1 = 0.5_dp
1234 149691 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1235 : END IF
1236 : END IF
1237 :
1238 6507253 : DEALLOCATE (oint, owork, sint)
1239 :
1240 : END DO
1241 3619 : CALL neighbor_list_iterator_release(nl_iterator)
1242 :
1243 7238 : DO i = 1, SIZE(matrix_h, 1)
1244 55169 : DO img = 1, nimg
1245 47931 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1246 51550 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1247 : END DO
1248 : END DO
1249 :
1250 3619 : kf = xtb_control%kf
1251 3619 : enscale = xtb_control%enscale
1252 3619 : erep = 0.0_dp
1253 3619 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1254 :
1255 3619 : exb = 0.0_dp
1256 3619 : IF (xb_inter) THEN
1257 3577 : CALL xb_interaction(qs_env, exb, calculate_forces)
1258 : END IF
1259 :
1260 3619 : enonbonded = 0.0_dp
1261 3619 : IF (do_nonbonded) THEN
1262 : ! nonbonded interactions
1263 34 : NULLIFY (sab_xtb_nonbond)
1264 34 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
1265 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1266 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1267 : END IF
1268 :
1269 : ! set repulsive energy
1270 3619 : erep = erep + exb + enonbonded
1271 3619 : IF (xb_inter) THEN
1272 3577 : CALL para_env%sum(exb)
1273 3577 : energy%xtb_xb_inter = exb
1274 : END IF
1275 3619 : IF (do_nonbonded) THEN
1276 34 : CALL para_env%sum(enonbonded)
1277 34 : energy%xtb_nonbonded = enonbonded
1278 : END IF
1279 3619 : CALL para_env%sum(erep)
1280 3619 : energy%repulsive = erep
1281 :
1282 : ! deallocate coordination numbers
1283 3619 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1284 :
1285 : ! deallocate Huckel parameters
1286 3619 : DEALLOCATE (huckel)
1287 3619 : IF (calculate_forces) THEN
1288 494 : DEALLOCATE (dhuckel)
1289 : END IF
1290 : ! deallocate KAB parameters
1291 3619 : DEALLOCATE (kijab)
1292 :
1293 : ! AO matrix outputs
1294 3619 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1295 :
1296 3619 : DEALLOCATE (basis_set_list)
1297 3619 : IF (calculate_forces) THEN
1298 494 : IF (SIZE(matrix_p, 1) == 2) THEN
1299 432 : DO img = 1, nimg
1300 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
1301 432 : beta_scalar=-1.0_dp)
1302 : END DO
1303 : END IF
1304 : END IF
1305 :
1306 3619 : CALL timestop(handle)
1307 :
1308 10857 : END SUBROUTINE build_gfn1_xtb_matrices
1309 :
1310 : ! **************************************************************************************************
1311 : !> \brief ...
1312 : !> \param qs_env ...
1313 : !> \param matrix_h ...
1314 : !> \param matrix_s ...
1315 : !> \param calculate_forces ...
1316 : ! **************************************************************************************************
1317 10598 : SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1318 : TYPE(qs_environment_type), POINTER :: qs_env
1319 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1320 : LOGICAL, INTENT(IN) :: calculate_forces
1321 :
1322 : INTEGER :: after, i, img, iw, nimg
1323 : LOGICAL :: norml1, norml2, omit_headers, use_arnoldi
1324 : REAL(KIND=dp), DIMENSION(2) :: condnum
1325 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1326 : TYPE(cp_logger_type), POINTER :: logger
1327 : TYPE(mp_para_env_type), POINTER :: para_env
1328 :
1329 5299 : logger => cp_get_default_logger()
1330 :
1331 5299 : CALL get_qs_env(qs_env, para_env=para_env)
1332 5299 : nimg = SIZE(matrix_h, 2)
1333 5299 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1334 5299 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1335 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
1336 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
1337 0 : extension=".Log")
1338 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1339 0 : after = MIN(MAX(after, 1), 16)
1340 0 : DO img = 1, nimg
1341 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
1342 0 : output_unit=iw, omit_headers=omit_headers)
1343 : END DO
1344 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
1345 : END IF
1346 :
1347 5299 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1348 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
1349 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
1350 0 : extension=".Log")
1351 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1352 0 : after = MIN(MAX(after, 1), 16)
1353 0 : DO img = 1, nimg
1354 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
1355 0 : output_unit=iw, omit_headers=omit_headers)
1356 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1357 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
1358 0 : DO i = 2, SIZE(matrix_s, 1)
1359 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
1360 0 : output_unit=iw, omit_headers=omit_headers)
1361 : END DO
1362 : END IF
1363 : END DO
1364 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
1365 : END IF
1366 :
1367 : ! *** Overlap condition number
1368 5299 : IF (.NOT. calculate_forces) THEN
1369 4753 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1370 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1371 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1372 4 : extension=".Log")
1373 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1374 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1375 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1376 4 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1377 4 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1378 : END IF
1379 : END IF
1380 :
1381 5299 : END SUBROUTINE ao_matrix_output
1382 :
1383 : END MODULE xtb_matrices
1384 :
|
