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