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 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 5084 : 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 5084 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
119 5084 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
120 :
121 1646 : SELECT CASE (gfn_type)
122 : CASE (0)
123 1646 : CALL build_gfn0_xtb_matrices(qs_env, calculate_forces)
124 : CASE (1)
125 3438 : 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 5084 : CPABORT("Unknown gfn_type")
130 : END SELECT
131 :
132 5084 : END SUBROUTINE build_xtb_matrices
133 :
134 : ! **************************************************************************************************
135 : !> \brief ...
136 : !> \param qs_env ...
137 : !> \param calculate_forces ...
138 : ! **************************************************************************************************
139 1646 : 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 3292 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
151 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
152 : INTEGER, DIMENSION(3) :: cell
153 1646 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
154 1646 : npgfb, nsgfa, nsgfb
155 1646 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
156 1646 : 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 3292 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges, qlagrange
163 3292 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, dqhuckel, huckel, owork
164 1646 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
165 1646 : 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 1646 : REAL(KIND=dp), DIMENSION(:), POINTER :: eeq_q, set_radius_a, set_radius_b
170 1646 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171 1646 : scon_a, scon_b, wblock, zeta, zetb
172 1646 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173 : TYPE(atprop_type), POINTER :: atprop
174 6584 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
175 : TYPE(cp_logger_type), POINTER :: logger
176 1646 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
177 1646 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
178 : TYPE(dft_control_type), POINTER :: dft_control
179 : TYPE(eeq_solver_type) :: eeq_sparam
180 1646 : 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 1646 : DIMENSION(:), POINTER :: nl_iterator
186 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187 1646 : POINTER :: sab_orb, sab_xtb_nonbond
188 1646 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
190 : TYPE(qs_energy_type), POINTER :: energy
191 1646 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
192 1646 : 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 1646 : CALL timeset(routineN, handle)
200 :
201 1646 : NULLIFY (logger, virial, atprop)
202 1646 : logger => cp_get_default_logger()
203 :
204 1646 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
205 1646 : 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 1646 : sab_orb=sab_orb)
217 :
218 1646 : nkind = SIZE(atomic_kind_set)
219 1646 : xtb_control => dft_control%qs_control%xtb_control
220 1646 : do_nonbonded = xtb_control%do_nonbonded
221 1646 : nimg = dft_control%nimages
222 1646 : nderivatives = 0
223 1646 : IF (calculate_forces) nderivatives = 1
224 1646 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
225 1646 : maxder = ncoset(nderivatives)
226 :
227 1646 : NULLIFY (particle_set)
228 1646 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
229 1646 : natom = SIZE(particle_set)
230 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
231 1646 : atom_of_kind=atom_of_kind, kind_of=kind_of)
232 :
233 1646 : IF (calculate_forces) THEN
234 42 : NULLIFY (rho, force, matrix_w)
235 : CALL get_qs_env(qs_env=qs_env, &
236 : rho=rho, matrix_w_kp=matrix_w, &
237 42 : virial=virial, force=force)
238 42 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
239 :
240 42 : 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 74 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
249 : END IF
250 : ! atomic energy decomposition
251 1646 : IF (atprop%energy) THEN
252 0 : CALL atprop_array_init(atprop%atecc, natom)
253 : END IF
254 :
255 1646 : NULLIFY (cell_to_index)
256 1646 : IF (nimg > 1) THEN
257 142 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
258 142 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
259 : END IF
260 :
261 : ! set up basis set lists
262 9022 : ALLOCATE (basis_set_list(nkind))
263 1646 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
264 :
265 : ! allocate overlap matrix
266 1646 : 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 1646 : sab_orb, .TRUE.)
269 1646 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
270 :
271 : ! initialize H matrix
272 1646 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
273 30552 : DO img = 1, nimg
274 28906 : ALLOCATE (matrix_h(1, img)%matrix)
275 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
276 28906 : name="HAMILTONIAN MATRIX")
277 30552 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
278 : END DO
279 1646 : 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 1646 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
285 :
286 4938 : ALLOCATE (charges(natom))
287 9396 : charges = 0.0_dp
288 1646 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
289 1646 : IF (calculate_forces) THEN
290 84 : ALLOCATE (dcharges(natom))
291 242 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
292 : END IF
293 1646 : energy%eeq = eeq_energy
294 1646 : energy%efield = ef_energy
295 :
296 1646 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
297 : ! prepare charges (needed for D4)
298 1646 : 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 1646 : energy%dispersion, calculate_forces)
311 1646 : IF (calculate_forces) THEN
312 42 : 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 1646 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
319 :
320 : ! Calculate KAB parameters and electronegativity correction
321 1646 : CALL gfn0_kpair(qs_env, kijab)
322 :
323 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
324 1646 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
325 422349 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
326 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
327 420703 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
328 420703 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
329 420703 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
330 420703 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
331 420703 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
332 420703 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
333 420703 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
334 :
335 1682812 : 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 420703 : 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 420703 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
342 :
343 420703 : IF (nimg == 1) THEN
344 : ic = 1
345 : ELSE
346 199870 : ic = cell_to_index(cell(1), cell(2), cell(3))
347 199870 : CPASSERT(ic > 0)
348 : END IF
349 :
350 420703 : icol = MAX(iatom, jatom)
351 420703 : irow = MIN(iatom, jatom)
352 420703 : NULLIFY (sblock, fblock)
353 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
354 420703 : row=irow, col=icol, BLOCK=sblock, found=found)
355 420703 : CPASSERT(found)
356 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
357 420703 : row=irow, col=icol, BLOCK=fblock, found=found)
358 420703 : CPASSERT(found)
359 :
360 420703 : IF (calculate_forces) THEN
361 11584 : NULLIFY (pblock)
362 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
363 11584 : row=irow, col=icol, block=pblock, found=found)
364 11584 : CPASSERT(ASSOCIATED(pblock))
365 11584 : NULLIFY (wblock)
366 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
367 11584 : row=irow, col=icol, block=wblock, found=found)
368 11584 : CPASSERT(ASSOCIATED(wblock))
369 46336 : DO i = 2, 4
370 34752 : NULLIFY (dsblocks(i)%block)
371 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
372 34752 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
373 46336 : CPASSERT(found)
374 : END DO
375 : END IF
376 :
377 : ! overlap
378 420703 : basis_set_a => basis_set_list(ikind)%gto_basis_set
379 420703 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
380 420703 : basis_set_b => basis_set_list(jkind)%gto_basis_set
381 420703 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
382 420703 : atom_a = atom_of_kind(iatom)
383 420703 : atom_b = atom_of_kind(jatom)
384 : ! basis ikind
385 420703 : first_sgfa => basis_set_a%first_sgf
386 420703 : la_max => basis_set_a%lmax
387 420703 : la_min => basis_set_a%lmin
388 420703 : npgfa => basis_set_a%npgf
389 420703 : nseta = basis_set_a%nset
390 420703 : nsgfa => basis_set_a%nsgf_set
391 420703 : rpgfa => basis_set_a%pgf_radius
392 420703 : set_radius_a => basis_set_a%set_radius
393 420703 : scon_a => basis_set_a%scon
394 420703 : zeta => basis_set_a%zet
395 : ! basis jkind
396 420703 : first_sgfb => basis_set_b%first_sgf
397 420703 : lb_max => basis_set_b%lmax
398 420703 : lb_min => basis_set_b%lmin
399 420703 : npgfb => basis_set_b%npgf
400 420703 : nsetb = basis_set_b%nset
401 420703 : nsgfb => basis_set_b%nsgf_set
402 420703 : rpgfb => basis_set_b%pgf_radius
403 420703 : set_radius_b => basis_set_b%set_radius
404 420703 : scon_b => basis_set_b%scon
405 420703 : zetb => basis_set_b%zet
406 :
407 420703 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
408 3365624 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
409 2103515 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
410 35613059 : sint = 0.0_dp
411 :
412 1619859 : DO iset = 1, nseta
413 1199156 : ncoa = npgfa(iset)*ncoset(la_max(iset))
414 1199156 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
415 1199156 : sgfa = first_sgfa(1, iset)
416 5043308 : DO jset = 1, nsetb
417 3423449 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
418 1749537 : ncob = npgfb(jset)*ncoset(lb_max(jset))
419 1749537 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
420 1749537 : sgfb = first_sgfb(1, jset)
421 1749537 : 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 47176 : 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 1702361 : 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 1749537 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
433 1749537 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
434 2948693 : IF (calculate_forces) THEN
435 188704 : DO i = 2, 4
436 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
437 141528 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
438 188704 : 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 420703 : IF (calculate_forces) THEN
445 46336 : DO i = 1, 3
446 46336 : IF (iatom <= jatom) THEN
447 1466268 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
448 : ELSE
449 1091826 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
450 : END IF
451 : END DO
452 11584 : f1 = 2.0_dp
453 46336 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
454 46336 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
455 11584 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
456 10966 : IF (iatom == jatom) f1 = 1.0_dp
457 10966 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
458 : END IF
459 : END IF
460 : ! update S matrix
461 420703 : IF (iatom <= jatom) THEN
462 18553151 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
463 : ELSE
464 14217966 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
465 : END IF
466 420703 : IF (calculate_forces) THEN
467 46336 : DO i = 2, 4
468 46336 : IF (iatom <= jatom) THEN
469 1466268 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
470 : ELSE
471 1108062 : 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 420703 : rcovab = rcova + rcovb
478 420703 : rrab = SQRT(dr/rcovab)
479 1619859 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
480 1615219 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
481 420703 : IF (calculate_forces) THEN
482 11584 : IF (dr > 1.e-6_dp) THEN
483 11484 : drx = 0.5_dp/rrab/rcovab
484 : ELSE
485 : drx = 0.0_dp
486 : END IF
487 44044 : dpia(1:nsa) = drx*kpolya(1:nsa)
488 43952 : dpib(1:nsb) = drx*kpolyb(1:nsb)
489 : END IF
490 :
491 : ! diagonal block
492 420703 : diagblock = .FALSE.
493 420703 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
494 : !
495 : ! Eq. 10
496 : !
497 : IF (diagblock) THEN
498 23510 : DO i = 1, natorb_a
499 19635 : na = naoa(i)
500 23510 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
501 : END DO
502 : ELSE
503 3825345 : DO j = 1, natorb_b
504 3408517 : nb = naob(j)
505 32481137 : DO i = 1, natorb_a
506 28655792 : na = naoa(i)
507 28655792 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
508 32064309 : IF (iatom <= jatom) THEN
509 16195035 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
510 : ELSE
511 12460757 : 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 420703 : IF (calculate_forces) THEN
517 11584 : f0 = 1.0_dp
518 11584 : IF (irow == iatom) f0 = -1.0_dp
519 11584 : f2 = 1.0_dp
520 11584 : IF (iatom /= jatom) f2 = 2.0_dp
521 : ! Derivative wrt coordination number
522 11584 : fhua = 0.0_dp
523 11584 : fhub = 0.0_dp
524 11584 : fhud = 0.0_dp
525 11584 : fqa = 0.0_dp
526 11584 : fqb = 0.0_dp
527 11584 : IF (diagblock) THEN
528 542 : DO i = 1, natorb_a
529 442 : la = laoa(i)
530 442 : na = naoa(i)
531 442 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
532 542 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
533 : END DO
534 100 : dcharges(iatom) = dcharges(iatom) + fqa
535 : ELSE
536 102648 : DO j = 1, natorb_b
537 91164 : lb = laob(j)
538 91164 : nb = naob(j)
539 849534 : DO i = 1, natorb_a
540 746886 : la = laoa(i)
541 746886 : na = naoa(i)
542 746886 : hij = 0.5_dp*pia(na)*pib(nb)
543 746886 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
544 838050 : IF (iatom <= jatom) THEN
545 425196 : fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
546 425196 : fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
547 425196 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
548 425196 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
549 : ELSE
550 321690 : fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
551 321690 : fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
552 321690 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
553 321690 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
554 : END IF
555 : END DO
556 : END DO
557 11484 : dcharges(iatom) = dcharges(iatom) + fqa
558 11484 : dcharges(jatom) = dcharges(jatom) + fqb
559 : END IF
560 : ! iatom
561 11584 : atom_a = atom_of_kind(iatom)
562 158260 : DO i = 1, dcnum(iatom)%neighbors
563 146676 : katom = dcnum(iatom)%nlist(i)
564 146676 : kkind = kind_of(katom)
565 146676 : atom_c = atom_of_kind(katom)
566 586704 : rik = dcnum(iatom)%rik(:, i)
567 586704 : drk = SQRT(SUM(rik(:)**2))
568 158260 : IF (drk > 1.e-3_dp) THEN
569 586704 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
570 586704 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
571 586704 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
572 586704 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
573 586704 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
574 586704 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
575 146676 : IF (use_virial) THEN
576 583032 : fdik = fdika + fdikb
577 145758 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
578 : END IF
579 : END IF
580 : END DO
581 : ! jatom
582 11584 : atom_b = atom_of_kind(jatom)
583 157422 : DO i = 1, dcnum(jatom)%neighbors
584 145838 : katom = dcnum(jatom)%nlist(i)
585 145838 : kkind = kind_of(katom)
586 145838 : atom_c = atom_of_kind(katom)
587 583352 : rik = dcnum(jatom)%rik(:, i)
588 583352 : drk = SQRT(SUM(rik(:)**2))
589 157422 : IF (drk > 1.e-3_dp) THEN
590 583352 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
591 583352 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
592 583352 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
593 145838 : IF (use_virial) THEN
594 145010 : 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 11584 : IF (diagblock) THEN
600 100 : force_ab = 0._dp
601 : ELSE
602 11484 : n1 = SIZE(fblock, 1)
603 11484 : n2 = SIZE(fblock, 2)
604 45936 : ALLOCATE (dfblock(n1, n2))
605 854946 : dfblock = 0.0_dp
606 102648 : DO j = 1, natorb_b
607 91164 : lb = laob(j)
608 91164 : nb = naob(j)
609 849534 : DO i = 1, natorb_a
610 746886 : la = laoa(i)
611 746886 : na = naoa(i)
612 746886 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
613 838050 : IF (iatom <= jatom) THEN
614 425196 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
615 : ELSE
616 321690 : 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 854946 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
621 45936 : DO ir = 1, 3
622 34452 : foab = 2.0_dp*dfp*rij(ir)/dr
623 : ! force from overlap matrix contribution to H
624 307944 : DO j = 1, natorb_b
625 273492 : lb = laob(j)
626 273492 : nb = naob(j)
627 2548602 : DO i = 1, natorb_a
628 2240658 : la = laoa(i)
629 2240658 : na = naoa(i)
630 2240658 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
631 2514150 : IF (iatom <= jatom) THEN
632 1275588 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
633 : ELSE
634 965070 : 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 45936 : force_ab(ir) = foab
639 : END DO
640 11484 : DEALLOCATE (dfblock)
641 : END IF
642 : END IF
643 :
644 420703 : IF (calculate_forces) THEN
645 11584 : atom_a = atom_of_kind(iatom)
646 11584 : atom_b = atom_of_kind(jatom)
647 31156 : IF (irow == iatom) force_ab = -force_ab
648 46336 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
649 46336 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
650 11584 : IF (use_virial) THEN
651 11002 : f1 = 1.0_dp
652 11002 : IF (iatom == jatom) f1 = 0.5_dp
653 11002 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
654 : END IF
655 : END IF
656 :
657 2525864 : DEALLOCATE (oint, owork, sint)
658 :
659 : END DO
660 1646 : CALL neighbor_list_iterator_release(nl_iterator)
661 :
662 3292 : DO i = 1, SIZE(matrix_h, 1)
663 32198 : DO img = 1, nimg
664 28906 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
665 30552 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
666 : END DO
667 : END DO
668 :
669 : ! EEQ forces (response and direct)
670 1646 : IF (calculate_forces) THEN
671 42 : CALL para_env%sum(dcharges)
672 84 : ALLOCATE (qlagrange(natom))
673 42 : CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
674 : END IF
675 :
676 1646 : kf = xtb_control%kf
677 1646 : enscale = xtb_control%enscale
678 1646 : erep = 0.0_dp
679 1646 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
680 :
681 1646 : esrb = 0.0_dp
682 1646 : CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
683 :
684 1646 : enonbonded = 0.0_dp
685 1646 : 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 1646 : erep = erep + esrb + enonbonded
695 1646 : IF (do_nonbonded) THEN
696 0 : CALL para_env%sum(enonbonded)
697 0 : energy%xtb_nonbonded = enonbonded
698 : END IF
699 1646 : CALL para_env%sum(esrb)
700 1646 : energy%srb = esrb
701 1646 : CALL para_env%sum(erep)
702 1646 : energy%repulsive = erep
703 :
704 : ! save EEQ charges
705 1646 : NULLIFY (eeq_q)
706 1646 : CALL get_qs_env(qs_env, eeq=eeq_q)
707 1646 : IF (ASSOCIATED(eeq_q)) THEN
708 972 : CPASSERT(SIZE(eeq_q) == natom)
709 : ELSE
710 1348 : ALLOCATE (eeq_q(natom))
711 3432 : eeq_q(1:natom) = charges(1:natom)
712 : END IF
713 1646 : CALL set_qs_env(qs_env, eeq=eeq_q)
714 :
715 : ! deallocate coordination numbers
716 1646 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
717 :
718 : ! deallocate Huckel parameters
719 1646 : DEALLOCATE (huckel)
720 1646 : IF (calculate_forces) THEN
721 42 : DEALLOCATE (dhuckel, dqhuckel)
722 : END IF
723 : ! deallocate KAB parameters
724 1646 : DEALLOCATE (kijab)
725 :
726 : ! deallocate charges
727 1646 : DEALLOCATE (charges)
728 1646 : IF (calculate_forces) THEN
729 42 : DEALLOCATE (dcharges, qlagrange)
730 : END IF
731 :
732 : ! AO matrix outputs
733 1646 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
734 :
735 1646 : DEALLOCATE (basis_set_list)
736 1646 : IF (calculate_forces) THEN
737 42 : 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 1646 : CALL timestop(handle)
746 :
747 4938 : END SUBROUTINE build_gfn0_xtb_matrices
748 :
749 : ! **************************************************************************************************
750 : !> \brief ...
751 : !> \param qs_env ...
752 : !> \param calculate_forces ...
753 : ! **************************************************************************************************
754 3438 : 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 6876 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
766 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
767 : INTEGER, DIMENSION(3) :: cell
768 3438 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
769 3438 : npgfb, nsgfa, nsgfb
770 3438 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
771 3438 : 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 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
779 6876 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
780 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
781 3438 : 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 3438 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
785 3438 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
786 3438 : scon_a, scon_b, wblock, zeta, zetb
787 3438 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
788 : TYPE(atprop_type), POINTER :: atprop
789 13752 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
790 : TYPE(cp_logger_type), POINTER :: logger
791 3438 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
792 3438 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
793 : TYPE(dft_control_type), POINTER :: dft_control
794 3438 : 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 3438 : DIMENSION(:), POINTER :: nl_iterator
800 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
801 3438 : POINTER :: sab_orb, sab_xtb_nonbond
802 3438 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
803 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
804 : TYPE(qs_energy_type), POINTER :: energy
805 3438 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
806 3438 : 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 3438 : CALL timeset(routineN, handle)
814 :
815 3438 : NULLIFY (logger, virial, atprop)
816 3438 : logger => cp_get_default_logger()
817 :
818 3438 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
819 3438 : 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 3438 : sab_orb=sab_orb)
832 :
833 3438 : nkind = SIZE(atomic_kind_set)
834 3438 : xtb_control => dft_control%qs_control%xtb_control
835 3438 : xb_inter = xtb_control%xb_interaction
836 3438 : do_nonbonded = xtb_control%do_nonbonded
837 3438 : nimg = dft_control%nimages
838 3438 : nderivatives = 0
839 3438 : IF (calculate_forces) nderivatives = 1
840 3438 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
841 3438 : maxder = ncoset(nderivatives)
842 :
843 3438 : NULLIFY (particle_set)
844 3438 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
845 3438 : natom = SIZE(particle_set)
846 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
847 3438 : atom_of_kind=atom_of_kind, kind_of=kind_of)
848 :
849 3438 : IF (calculate_forces) THEN
850 444 : NULLIFY (rho, force, matrix_w)
851 : CALL get_qs_env(qs_env=qs_env, &
852 : rho=rho, matrix_w_kp=matrix_w, &
853 444 : virial=virial, force=force)
854 444 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
855 :
856 444 : 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 834 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
865 : END IF
866 : ! atomic energy decomposition
867 3438 : IF (atprop%energy) THEN
868 36 : CALL atprop_array_init(atprop%atecc, natom)
869 : END IF
870 :
871 3438 : NULLIFY (cell_to_index)
872 3438 : IF (nimg > 1) THEN
873 336 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
874 336 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
875 : END IF
876 :
877 : ! set up basis set lists
878 18904 : ALLOCATE (basis_set_list(nkind))
879 3438 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
880 :
881 : ! allocate overlap matrix
882 3438 : 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 3438 : sab_orb, .TRUE.)
885 3438 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
886 :
887 : ! initialize H matrix
888 3438 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
889 50512 : DO img = 1, nimg
890 47074 : ALLOCATE (matrix_h(1, img)%matrix)
891 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
892 47074 : name="HAMILTONIAN MATRIX")
893 50512 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
894 : END DO
895 3438 : 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 3438 : CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
901 :
902 : ! vdW Potential
903 3438 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
904 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
905 3438 : energy%dispersion, calculate_forces)
906 :
907 : ! Calculate Huckel parameters
908 3438 : CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
909 :
910 : ! Calculate KAB parameters and electronegativity correction
911 3438 : CALL gfn1_kpair(qs_env, kijab)
912 :
913 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
914 3438 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
915 952422 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
916 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
917 948984 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
918 948984 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
919 948984 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
920 948984 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
921 948984 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
922 948984 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
923 948984 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
924 :
925 3795936 : 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 948984 : 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 948984 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
932 :
933 948984 : IF (nimg == 1) THEN
934 : ic = 1
935 : ELSE
936 241506 : ic = cell_to_index(cell(1), cell(2), cell(3))
937 241506 : CPASSERT(ic > 0)
938 : END IF
939 :
940 948984 : icol = MAX(iatom, jatom)
941 948984 : irow = MIN(iatom, jatom)
942 948984 : NULLIFY (sblock, fblock)
943 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
944 948984 : row=irow, col=icol, BLOCK=sblock, found=found)
945 948984 : CPASSERT(found)
946 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
947 948984 : row=irow, col=icol, BLOCK=fblock, found=found)
948 948984 : CPASSERT(found)
949 :
950 948984 : IF (calculate_forces) THEN
951 193867 : NULLIFY (pblock)
952 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
953 193867 : row=irow, col=icol, block=pblock, found=found)
954 193867 : CPASSERT(found)
955 193867 : NULLIFY (wblock)
956 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
957 193867 : row=irow, col=icol, block=wblock, found=found)
958 193867 : CPASSERT(found)
959 775468 : DO i = 2, 4
960 581601 : NULLIFY (dsblocks(i)%block)
961 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
962 581601 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
963 775468 : CPASSERT(found)
964 : END DO
965 : END IF
966 :
967 : ! overlap
968 948984 : basis_set_a => basis_set_list(ikind)%gto_basis_set
969 948984 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
970 948984 : basis_set_b => basis_set_list(jkind)%gto_basis_set
971 948984 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
972 948984 : atom_a = atom_of_kind(iatom)
973 948984 : atom_b = atom_of_kind(jatom)
974 : ! basis ikind
975 948984 : first_sgfa => basis_set_a%first_sgf
976 948984 : la_max => basis_set_a%lmax
977 948984 : la_min => basis_set_a%lmin
978 948984 : npgfa => basis_set_a%npgf
979 948984 : nseta = basis_set_a%nset
980 948984 : nsgfa => basis_set_a%nsgf_set
981 948984 : rpgfa => basis_set_a%pgf_radius
982 948984 : set_radius_a => basis_set_a%set_radius
983 948984 : scon_a => basis_set_a%scon
984 948984 : zeta => basis_set_a%zet
985 : ! basis jkind
986 948984 : first_sgfb => basis_set_b%first_sgf
987 948984 : lb_max => basis_set_b%lmax
988 948984 : lb_min => basis_set_b%lmin
989 948984 : npgfb => basis_set_b%npgf
990 948984 : nsetb = basis_set_b%nset
991 948984 : nsgfb => basis_set_b%nsgf_set
992 948984 : rpgfb => basis_set_b%pgf_radius
993 948984 : set_radius_b => basis_set_b%set_radius
994 948984 : scon_b => basis_set_b%scon
995 948984 : zetb => basis_set_b%zet
996 :
997 948984 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
998 7591872 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
999 4744920 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
1000 34052672 : sint = 0.0_dp
1001 :
1002 2931985 : DO iset = 1, nseta
1003 1983001 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1004 1983001 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1005 1983001 : sgfa = first_sgfa(1, iset)
1006 7151082 : DO jset = 1, nsetb
1007 4219097 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1008 3581445 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1009 3581445 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1010 3581445 : sgfb = first_sgfb(1, jset)
1011 3581445 : 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 738061 : 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 2843384 : 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 3581445 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1023 3581445 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1024 5564446 : IF (calculate_forces) THEN
1025 2952244 : DO i = 2, 4
1026 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1027 2214183 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1028 2952244 : 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 948984 : IF (calculate_forces) THEN
1035 775468 : DO i = 1, 3
1036 775468 : IF (iatom <= jatom) THEN
1037 8167620 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
1038 : ELSE
1039 6154986 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
1040 : END IF
1041 : END DO
1042 193867 : f1 = 2.0_dp
1043 775468 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1044 775468 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1045 193867 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
1046 100045 : IF (iatom == jatom) f1 = 1.0_dp
1047 100045 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1048 : END IF
1049 : END IF
1050 : ! update S matrix
1051 948984 : IF (iatom <= jatom) THEN
1052 9383324 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1053 : ELSE
1054 7370208 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1055 : END IF
1056 948984 : IF (calculate_forces) THEN
1057 775468 : DO i = 2, 4
1058 775468 : IF (iatom <= jatom) THEN
1059 8167620 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1060 : ELSE
1061 6114060 : 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 948984 : rcovab = rcova + rcovb
1068 948984 : rrab = SQRT(dr/rcovab)
1069 2931985 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1070 2931990 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1071 948984 : IF (calculate_forces) THEN
1072 193867 : IF (dr > 1.e-6_dp) THEN
1073 191787 : drx = 0.5_dp/rrab/rcovab
1074 : ELSE
1075 : drx = 0.0_dp
1076 : END IF
1077 617188 : dpia(1:nsa) = drx*kpolya(1:nsa)
1078 617190 : dpib(1:nsb) = drx*kpolyb(1:nsb)
1079 : END IF
1080 :
1081 : ! diagonal block
1082 948984 : diagblock = .FALSE.
1083 948984 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
1084 : !
1085 : ! Eq. 10
1086 : !
1087 : IF (diagblock) THEN
1088 60553 : DO i = 1, natorb_a
1089 45603 : na = naoa(i)
1090 60553 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1091 : END DO
1092 : ELSE
1093 3919853 : DO j = 1, natorb_b
1094 2985819 : nb = naob(j)
1095 16630245 : DO i = 1, natorb_a
1096 12710392 : na = naoa(i)
1097 12710392 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1098 15696211 : IF (iatom <= jatom) THEN
1099 7112778 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1100 : ELSE
1101 5597614 : 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 948984 : IF (calculate_forces) THEN
1107 193867 : f0 = 1.0_dp
1108 193867 : IF (irow == iatom) f0 = -1.0_dp
1109 : ! Derivative wrt coordination number
1110 193867 : fhua = 0.0_dp
1111 193867 : fhub = 0.0_dp
1112 193867 : fhud = 0.0_dp
1113 193867 : IF (diagblock) THEN
1114 8408 : DO i = 1, natorb_a
1115 6328 : la = laoa(i)
1116 6328 : na = naoa(i)
1117 8408 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1118 : END DO
1119 : ELSE
1120 903752 : DO j = 1, natorb_b
1121 711965 : lb = laob(j)
1122 711965 : nb = naob(j)
1123 4740616 : DO i = 1, natorb_a
1124 3836864 : la = laoa(i)
1125 3836864 : na = naoa(i)
1126 3836864 : hij = 0.5_dp*pia(na)*pib(nb)
1127 4548829 : IF (iatom <= jatom) THEN
1128 2206609 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1129 2206609 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1130 : ELSE
1131 1630255 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1132 1630255 : 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 191787 : IF (iatom /= jatom) THEN
1137 172335 : fhua = 2.0_dp*fhua
1138 172335 : fhub = 2.0_dp*fhub
1139 : END IF
1140 : END IF
1141 : ! iatom
1142 193867 : atom_a = atom_of_kind(iatom)
1143 5044297 : DO i = 1, dcnum(iatom)%neighbors
1144 4850430 : katom = dcnum(iatom)%nlist(i)
1145 4850430 : kkind = kind_of(katom)
1146 4850430 : atom_c = atom_of_kind(katom)
1147 19401720 : rik = dcnum(iatom)%rik(:, i)
1148 19401720 : drk = SQRT(SUM(rik(:)**2))
1149 5044297 : IF (drk > 1.e-3_dp) THEN
1150 19401720 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1151 19401720 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1152 19401720 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1153 19401720 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1154 19401720 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1155 19401720 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1156 4850430 : IF (use_virial) THEN
1157 9227644 : fdik = fdika + fdikb
1158 2306911 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1159 : END IF
1160 : END IF
1161 : END DO
1162 : ! jatom
1163 193867 : atom_b = atom_of_kind(jatom)
1164 5041492 : DO i = 1, dcnum(jatom)%neighbors
1165 4847625 : katom = dcnum(jatom)%nlist(i)
1166 4847625 : kkind = kind_of(katom)
1167 4847625 : atom_c = atom_of_kind(katom)
1168 19390500 : rik = dcnum(jatom)%rik(:, i)
1169 19390500 : drk = SQRT(SUM(rik(:)**2))
1170 5041492 : IF (drk > 1.e-3_dp) THEN
1171 19390500 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1172 19390500 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1173 19390500 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1174 4847625 : IF (use_virial) THEN
1175 2305897 : 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 193867 : IF (diagblock) THEN
1181 2080 : force_ab = 0._dp
1182 : ELSE
1183 191787 : n1 = SIZE(fblock, 1)
1184 191787 : n2 = SIZE(fblock, 2)
1185 767148 : ALLOCATE (dfblock(n1, n2))
1186 4726974 : dfblock = 0.0_dp
1187 903752 : DO j = 1, natorb_b
1188 711965 : lb = laob(j)
1189 711965 : nb = naob(j)
1190 4740616 : DO i = 1, natorb_a
1191 3836864 : la = laoa(i)
1192 3836864 : na = naoa(i)
1193 3836864 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1194 4548829 : IF (iatom <= jatom) THEN
1195 2206609 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1196 : ELSE
1197 1630255 : 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 4726974 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
1202 767148 : DO ir = 1, 3
1203 575361 : foab = 2.0_dp*dfp*rij(ir)/dr
1204 : ! force from overlap matrix contribution to H
1205 2711256 : DO j = 1, natorb_b
1206 2135895 : lb = laob(j)
1207 2135895 : nb = naob(j)
1208 14221848 : DO i = 1, natorb_a
1209 11510592 : la = laoa(i)
1210 11510592 : na = naoa(i)
1211 11510592 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1212 13646487 : IF (iatom <= jatom) THEN
1213 6619827 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1214 : ELSE
1215 4890765 : 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 767148 : force_ab(ir) = foab
1220 : END DO
1221 191787 : DEALLOCATE (dfblock)
1222 : END IF
1223 : END IF
1224 :
1225 948984 : IF (calculate_forces) THEN
1226 193867 : atom_a = atom_of_kind(iatom)
1227 193867 : atom_b = atom_of_kind(jatom)
1228 500332 : IF (irow == iatom) force_ab = -force_ab
1229 775468 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1230 775468 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1231 193867 : IF (use_virial) THEN
1232 100685 : f1 = 1.0_dp
1233 100685 : IF (iatom == jatom) f1 = 0.5_dp
1234 100685 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1235 : END IF
1236 : END IF
1237 :
1238 5697342 : DEALLOCATE (oint, owork, sint)
1239 :
1240 : END DO
1241 3438 : CALL neighbor_list_iterator_release(nl_iterator)
1242 :
1243 6876 : DO i = 1, SIZE(matrix_h, 1)
1244 53950 : DO img = 1, nimg
1245 47074 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1246 50512 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1247 : END DO
1248 : END DO
1249 :
1250 3438 : kf = xtb_control%kf
1251 3438 : enscale = xtb_control%enscale
1252 3438 : erep = 0.0_dp
1253 3438 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1254 :
1255 3438 : exb = 0.0_dp
1256 3438 : IF (xb_inter) THEN
1257 3438 : CALL xb_interaction(qs_env, exb, calculate_forces)
1258 : END IF
1259 :
1260 3438 : enonbonded = 0.0_dp
1261 3438 : 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 3438 : erep = erep + exb + enonbonded
1271 3438 : IF (xb_inter) THEN
1272 3438 : CALL para_env%sum(exb)
1273 3438 : energy%xtb_xb_inter = exb
1274 : END IF
1275 3438 : IF (do_nonbonded) THEN
1276 34 : CALL para_env%sum(enonbonded)
1277 34 : energy%xtb_nonbonded = enonbonded
1278 : END IF
1279 3438 : CALL para_env%sum(erep)
1280 3438 : energy%repulsive = erep
1281 :
1282 : ! deallocate coordination numbers
1283 3438 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1284 :
1285 : ! deallocate Huckel parameters
1286 3438 : DEALLOCATE (huckel)
1287 3438 : IF (calculate_forces) THEN
1288 444 : DEALLOCATE (dhuckel)
1289 : END IF
1290 : ! deallocate KAB parameters
1291 3438 : DEALLOCATE (kijab)
1292 :
1293 : ! AO matrix outputs
1294 3438 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1295 :
1296 3438 : DEALLOCATE (basis_set_list)
1297 3438 : IF (calculate_forces) THEN
1298 444 : 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 3438 : CALL timestop(handle)
1307 :
1308 10314 : 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 10168 : 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 5084 : logger => cp_get_default_logger()
1330 :
1331 5084 : CALL get_qs_env(qs_env, para_env=para_env)
1332 5084 : nimg = SIZE(matrix_h, 2)
1333 5084 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1334 5084 : 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 5084 : 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 5084 : IF (.NOT. calculate_forces) THEN
1369 4598 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1370 : "DFT%PRINT%OVERLAP_CONDITION") .NE. 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 5084 : END SUBROUTINE ao_matrix_output
1382 :
1383 : END MODULE xtb_matrices
1384 :
|
