Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of Overlap and Hamiltonian matrices in xTB
10 : !> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11 : !> JCTC 13, 1989-2009, (2017)
12 : !> DOI: 10.1021/acs.jctc.7b00118
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE xtb_matrices
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction
18 : USE ai_overlap, ONLY: overlap_ab
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind_set
21 : USE atprop_types, ONLY: atprop_array_init,&
22 : atprop_type
23 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
24 : gto_basis_set_type
25 : USE block_p_types, ONLY: block_p_type
26 : USE cp_blacs_env, ONLY: cp_blacs_env_type
27 : USE cp_control_types, ONLY: dft_control_type,&
28 : xtb_control_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
30 : dbcsr_create,&
31 : dbcsr_finalize,&
32 : dbcsr_get_block_p,&
33 : dbcsr_p_type
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
36 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_type
39 : USE cp_output_handling, ONLY: cp_p_file,&
40 : cp_print_key_finished_output,&
41 : cp_print_key_should_output,&
42 : cp_print_key_unit_nr
43 : USE eeq_input, ONLY: eeq_solver_type
44 : USE input_constants, ONLY: vdw_pairpot_dftd4
45 : USE input_section_types, ONLY: section_vals_val_get
46 : USE kinds, ONLY: dp
47 : USE kpoint_types, ONLY: get_kpoint_info,&
48 : kpoint_type
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE orbital_pointers, ONLY: ncoset
51 : USE particle_types, ONLY: particle_type
52 : USE qs_condnum, ONLY: overlap_condnum
53 : USE qs_dispersion_cnum, ONLY: cnumber_init,&
54 : cnumber_release,&
55 : dcnum_type
56 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
57 : USE qs_dispersion_types, ONLY: qs_dispersion_type
58 : USE qs_energy_types, ONLY: qs_energy_type
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type,&
61 : set_qs_env
62 : USE qs_force_types, ONLY: qs_force_type
63 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
64 : get_memory_usage
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : qs_kind_type
67 : USE qs_ks_types, ONLY: get_ks_env,&
68 : qs_ks_env_type,&
69 : set_ks_env
70 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
71 : neighbor_list_iterate,&
72 : neighbor_list_iterator_create,&
73 : neighbor_list_iterator_p_type,&
74 : neighbor_list_iterator_release,&
75 : neighbor_list_set_p_type
76 : USE qs_overlap, ONLY: create_sab_matrix
77 : USE qs_rho_types, ONLY: qs_rho_get,&
78 : qs_rho_type
79 : USE virial_methods, ONLY: virial_pair_force
80 : USE virial_types, ONLY: virial_type
81 : USE xtb_eeq, ONLY: xtb_eeq_calculation,&
82 : xtb_eeq_forces
83 : USE xtb_hcore, ONLY: gfn0_huckel,&
84 : gfn0_kpair,&
85 : gfn1_huckel,&
86 : gfn1_kpair
87 : USE xtb_potentials, ONLY: nonbonded_correction,&
88 : repulsive_potential,&
89 : srb_potential,&
90 : xb_interaction
91 : USE xtb_types, ONLY: get_xtb_atom_param,&
92 : xtb_atom_type
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 :
97 : PRIVATE
98 :
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_matrices'
100 :
101 : PUBLIC :: build_xtb_matrices
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param qs_env ...
108 : !> \param calculate_forces ...
109 : ! **************************************************************************************************
110 6016 : 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 6016 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
119 6016 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
120 :
121 2046 : SELECT CASE (gfn_type)
122 : CASE (0)
123 2046 : CALL build_gfn0_xtb_matrices(qs_env, calculate_forces)
124 : CASE (1)
125 3970 : 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 6016 : CPABORT("Unknown gfn_type")
130 : END SELECT
131 :
132 6016 : END SUBROUTINE build_xtb_matrices
133 :
134 : ! **************************************************************************************************
135 : !> \brief ...
136 : !> \param qs_env ...
137 : !> \param calculate_forces ...
138 : ! **************************************************************************************************
139 2046 : 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 4092 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
151 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
152 : INTEGER, DIMENSION(3) :: cell
153 2046 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
154 2046 : npgfb, nsgfa, nsgfb
155 2046 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
156 2046 : 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 4092 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges, qlagrange
163 4092 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, dqhuckel, huckel, owork
164 2046 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
165 2046 : 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 2046 : REAL(KIND=dp), DIMENSION(:), POINTER :: eeq_q, set_radius_a, set_radius_b
170 2046 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171 2046 : scon_a, scon_b, wblock, zeta, zetb
172 2046 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173 : TYPE(atprop_type), POINTER :: atprop
174 8184 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
175 : TYPE(cp_logger_type), POINTER :: logger
176 2046 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
177 2046 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
178 : TYPE(dft_control_type), POINTER :: dft_control
179 : TYPE(eeq_solver_type) :: eeq_sparam
180 2046 : 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 2046 : DIMENSION(:), POINTER :: nl_iterator
186 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187 2046 : POINTER :: sab_orb, sab_xtb_nonbond
188 2046 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
190 : TYPE(qs_energy_type), POINTER :: energy
191 2046 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
192 2046 : 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 2046 : CALL timeset(routineN, handle)
200 :
201 2046 : NULLIFY (logger, virial, atprop)
202 2046 : logger => cp_get_default_logger()
203 :
204 2046 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
205 2046 : 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 2046 : sab_orb=sab_orb)
217 :
218 2046 : nkind = SIZE(atomic_kind_set)
219 2046 : xtb_control => dft_control%qs_control%xtb_control
220 2046 : do_nonbonded = xtb_control%do_nonbonded
221 2046 : nimg = dft_control%nimages
222 2046 : nderivatives = 0
223 2046 : IF (calculate_forces) nderivatives = 1
224 2046 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
225 2046 : maxder = ncoset(nderivatives)
226 :
227 2046 : NULLIFY (particle_set)
228 2046 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
229 2046 : natom = SIZE(particle_set)
230 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
231 2046 : atom_of_kind=atom_of_kind, kind_of=kind_of)
232 :
233 2046 : IF (calculate_forces) THEN
234 68 : NULLIFY (rho, force, matrix_w)
235 : CALL get_qs_env(qs_env=qs_env, &
236 : rho=rho, matrix_w_kp=matrix_w, &
237 68 : virial=virial, force=force)
238 68 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
239 :
240 68 : IF (SIZE(matrix_p, 1) == 2) THEN
241 16 : DO img = 1, nimg
242 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
243 8 : 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 16 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
246 : END DO
247 : END IF
248 108 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
249 : END IF
250 : ! atomic energy decomposition
251 2046 : IF (atprop%energy) THEN
252 0 : CALL atprop_array_init(atprop%atecc, natom)
253 : END IF
254 :
255 2046 : NULLIFY (cell_to_index)
256 2046 : IF (nimg > 1) THEN
257 376 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
258 376 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
259 : END IF
260 :
261 : ! set up basis set lists
262 10872 : ALLOCATE (basis_set_list(nkind))
263 2046 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
264 :
265 : ! allocate overlap matrix
266 2046 : 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 2046 : sab_orb, .TRUE.)
269 2046 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
270 :
271 : ! initialize H matrix
272 2046 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
273 39138 : DO img = 1, nimg
274 37092 : ALLOCATE (matrix_h(1, img)%matrix)
275 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
276 37092 : name="HAMILTONIAN MATRIX")
277 39138 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
278 : END DO
279 2046 : 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 2046 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
285 :
286 6138 : ALLOCATE (charges(natom))
287 2046 : charges = 0.0_dp
288 2046 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
289 2046 : IF (calculate_forces) THEN
290 136 : ALLOCATE (dcharges(natom))
291 416 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
292 : END IF
293 2046 : energy%eeq = eeq_energy
294 2046 : energy%efield = ef_energy
295 :
296 2046 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
297 : ! prepare charges (needed for D4)
298 2046 : 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 2046 : energy%dispersion, calculate_forces)
311 2046 : IF (calculate_forces) THEN
312 68 : 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 2046 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
319 :
320 : ! Calculate KAB parameters and electronegativity correction
321 2046 : CALL gfn0_kpair(qs_env, kijab)
322 :
323 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
324 2046 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
325 493308 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
326 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
327 491262 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
328 491262 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
329 491262 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
330 491262 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
331 491262 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
332 491262 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
333 491262 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
334 :
335 1965048 : 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 491262 : 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 491262 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
342 :
343 491262 : IF (nimg == 1) THEN
344 : ic = 1
345 : ELSE
346 269927 : ic = cell_to_index(cell(1), cell(2), cell(3))
347 269927 : CPASSERT(ic > 0)
348 : END IF
349 :
350 491262 : icol = MAX(iatom, jatom)
351 491262 : irow = MIN(iatom, jatom)
352 491262 : NULLIFY (sblock, fblock)
353 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
354 491262 : row=irow, col=icol, BLOCK=sblock, found=found)
355 491262 : CPASSERT(found)
356 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
357 491262 : row=irow, col=icol, BLOCK=fblock, found=found)
358 491262 : CPASSERT(found)
359 :
360 491262 : IF (calculate_forces) THEN
361 28766 : NULLIFY (pblock)
362 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
363 28766 : row=irow, col=icol, block=pblock, found=found)
364 28766 : CPASSERT(ASSOCIATED(pblock))
365 28766 : NULLIFY (wblock)
366 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
367 28766 : row=irow, col=icol, block=wblock, found=found)
368 28766 : CPASSERT(ASSOCIATED(wblock))
369 115064 : DO i = 2, 4
370 86298 : NULLIFY (dsblocks(i)%block)
371 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
372 86298 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
373 115064 : CPASSERT(found)
374 : END DO
375 : END IF
376 :
377 : ! overlap
378 491262 : basis_set_a => basis_set_list(ikind)%gto_basis_set
379 491262 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
380 491262 : basis_set_b => basis_set_list(jkind)%gto_basis_set
381 491262 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
382 491262 : atom_a = atom_of_kind(iatom)
383 491262 : atom_b = atom_of_kind(jatom)
384 : ! basis ikind
385 491262 : first_sgfa => basis_set_a%first_sgf
386 491262 : la_max => basis_set_a%lmax
387 491262 : la_min => basis_set_a%lmin
388 491262 : npgfa => basis_set_a%npgf
389 491262 : nseta = basis_set_a%nset
390 491262 : nsgfa => basis_set_a%nsgf_set
391 491262 : rpgfa => basis_set_a%pgf_radius
392 491262 : set_radius_a => basis_set_a%set_radius
393 491262 : scon_a => basis_set_a%scon
394 491262 : zeta => basis_set_a%zet
395 : ! basis jkind
396 491262 : first_sgfb => basis_set_b%first_sgf
397 491262 : lb_max => basis_set_b%lmax
398 491262 : lb_min => basis_set_b%lmin
399 491262 : npgfb => basis_set_b%npgf
400 491262 : nsetb = basis_set_b%nset
401 491262 : nsgfb => basis_set_b%nsgf_set
402 491262 : rpgfb => basis_set_b%pgf_radius
403 491262 : set_radius_b => basis_set_b%set_radius
404 491262 : scon_b => basis_set_b%scon
405 491262 : zetb => basis_set_b%zet
406 :
407 491262 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
408 3930096 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
409 2456310 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
410 491262 : sint = 0.0_dp
411 :
412 1899871 : DO iset = 1, nseta
413 1408609 : ncoa = npgfa(iset)*ncoset(la_max(iset))
414 1408609 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
415 1408609 : sgfa = first_sgfa(1, iset)
416 5947231 : DO jset = 1, nsetb
417 4047360 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
418 2088629 : ncob = npgfb(jset)*ncoset(lb_max(jset))
419 2088629 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
420 2088629 : sgfb = first_sgfb(1, jset)
421 2088629 : 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 123944 : 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 1964685 : 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 2088629 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
433 2088629 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
434 3497238 : IF (calculate_forces) THEN
435 495776 : DO i = 2, 4
436 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
437 371832 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
438 495776 : 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 491262 : IF (calculate_forces) THEN
445 115064 : DO i = 1, 3
446 115064 : IF (iatom <= jatom) THEN
447 4051116 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
448 : ELSE
449 3171456 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
450 : END IF
451 : END DO
452 28766 : f1 = 2.0_dp
453 115064 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
454 115064 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
455 28766 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
456 27840 : IF (iatom == jatom) f1 = 1.0_dp
457 27840 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
458 : END IF
459 : END IF
460 : ! update S matrix
461 491262 : IF (iatom <= jatom) THEN
462 22018173 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
463 : ELSE
464 16993369 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
465 : END IF
466 491262 : IF (calculate_forces) THEN
467 115064 : DO i = 2, 4
468 115064 : IF (iatom <= jatom) THEN
469 4051116 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
470 : ELSE
471 3187608 : 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 491262 : rcovab = rcova + rcovb
478 491262 : rrab = SQRT(dr/rcovab)
479 1899871 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
480 1895231 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
481 491262 : IF (calculate_forces) THEN
482 28766 : IF (dr > 1.e-6_dp) THEN
483 28592 : drx = 0.5_dp/rrab/rcovab
484 : ELSE
485 : drx = 0.0_dp
486 : END IF
487 112664 : dpia(1:nsa) = drx*kpolya(1:nsa)
488 112572 : dpib(1:nsb) = drx*kpolyb(1:nsb)
489 : END IF
490 :
491 : ! diagonal block
492 491262 : diagblock = .FALSE.
493 491262 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
494 : !
495 : ! Eq. 10
496 : !
497 : IF (diagblock) THEN
498 30959 : DO i = 1, natorb_a
499 26103 : na = naoa(i)
500 30959 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
501 : END DO
502 : ELSE
503 4509076 : DO j = 1, natorb_b
504 4022670 : nb = naob(j)
505 38662451 : DO i = 1, natorb_a
506 34153375 : na = naoa(i)
507 34153375 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
508 38176045 : IF (iatom <= jatom) THEN
509 19223733 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
510 : ELSE
511 14929642 : 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 491262 : IF (calculate_forces) THEN
517 28766 : f0 = 1.0_dp
518 28766 : IF (irow == iatom) f0 = -1.0_dp
519 28766 : f2 = 1.0_dp
520 28766 : IF (iatom /= jatom) f2 = 2.0_dp
521 : ! Derivative wrt coordination number
522 28766 : fhua = 0.0_dp
523 28766 : fhub = 0.0_dp
524 28766 : fhud = 0.0_dp
525 28766 : fqa = 0.0_dp
526 28766 : fqb = 0.0_dp
527 28766 : IF (diagblock) THEN
528 1168 : DO i = 1, natorb_a
529 994 : la = laoa(i)
530 994 : na = naoa(i)
531 994 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
532 1168 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
533 : END DO
534 174 : dcharges(iatom) = dcharges(iatom) + fqa
535 : ELSE
536 273142 : DO j = 1, natorb_b
537 244550 : lb = laob(j)
538 244550 : nb = naob(j)
539 2399054 : DO i = 1, natorb_a
540 2125912 : la = laoa(i)
541 2125912 : na = naoa(i)
542 2125912 : hij = 0.5_dp*pia(na)*pib(nb)
543 2125912 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
544 2370462 : IF (iatom <= jatom) THEN
545 1187288 : fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
546 1187288 : fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
547 1187288 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
548 1187288 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
549 : ELSE
550 938624 : fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
551 938624 : fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
552 938624 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
553 938624 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
554 : END IF
555 : END DO
556 : END DO
557 28592 : dcharges(iatom) = dcharges(iatom) + fqa
558 28592 : dcharges(jatom) = dcharges(jatom) + fqb
559 : END IF
560 : ! iatom
561 28766 : atom_a = atom_of_kind(iatom)
562 445324 : DO i = 1, dcnum(iatom)%neighbors
563 416558 : katom = dcnum(iatom)%nlist(i)
564 416558 : kkind = kind_of(katom)
565 416558 : atom_c = atom_of_kind(katom)
566 1666232 : rik = dcnum(iatom)%rik(:, i)
567 1666232 : drk = SQRT(SUM(rik(:)**2))
568 445324 : IF (drk > 1.e-3_dp) THEN
569 1666232 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
570 1666232 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
571 1666232 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
572 1666232 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
573 1666232 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
574 1666232 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
575 416558 : IF (use_virial) THEN
576 1656556 : fdik = fdika + fdikb
577 414139 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
578 : END IF
579 : END IF
580 : END DO
581 : ! jatom
582 28766 : atom_b = atom_of_kind(jatom)
583 444486 : DO i = 1, dcnum(jatom)%neighbors
584 415720 : katom = dcnum(jatom)%nlist(i)
585 415720 : kkind = kind_of(katom)
586 415720 : atom_c = atom_of_kind(katom)
587 1662880 : rik = dcnum(jatom)%rik(:, i)
588 1662880 : drk = SQRT(SUM(rik(:)**2))
589 444486 : IF (drk > 1.e-3_dp) THEN
590 1662880 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
591 1662880 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
592 1662880 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
593 415720 : IF (use_virial) THEN
594 413391 : 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 28766 : IF (diagblock) THEN
600 174 : force_ab = 0._dp
601 : ELSE
602 28592 : n1 = SIZE(fblock, 1)
603 28592 : n2 = SIZE(fblock, 2)
604 114368 : ALLOCATE (dfblock(n1, n2))
605 28592 : dfblock = 0.0_dp
606 273142 : DO j = 1, natorb_b
607 244550 : lb = laob(j)
608 244550 : nb = naob(j)
609 2399054 : DO i = 1, natorb_a
610 2125912 : la = laoa(i)
611 2125912 : na = naoa(i)
612 2125912 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
613 2370462 : IF (iatom <= jatom) THEN
614 1187288 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
615 : ELSE
616 938624 : 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 2404438 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
621 114368 : DO ir = 1, 3
622 85776 : foab = 2.0_dp*dfp*rij(ir)/dr
623 : ! force from overlap matrix contribution to H
624 819426 : DO j = 1, natorb_b
625 733650 : lb = laob(j)
626 733650 : nb = naob(j)
627 7197162 : DO i = 1, natorb_a
628 6377736 : la = laoa(i)
629 6377736 : na = naoa(i)
630 6377736 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
631 7111386 : IF (iatom <= jatom) THEN
632 3561864 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
633 : ELSE
634 2815872 : 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 114368 : force_ab(ir) = foab
639 : END DO
640 28592 : DEALLOCATE (dfblock)
641 : END IF
642 : END IF
643 :
644 491262 : IF (calculate_forces) THEN
645 28766 : atom_a = atom_of_kind(iatom)
646 28766 : atom_b = atom_of_kind(jatom)
647 76934 : IF (irow == iatom) force_ab = -force_ab
648 115064 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
649 115064 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
650 28766 : IF (use_virial) THEN
651 27933 : f1 = 1.0_dp
652 27933 : IF (iatom == jatom) f1 = 0.5_dp
653 27933 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
654 : END IF
655 : END IF
656 :
657 2949618 : DEALLOCATE (oint, owork, sint)
658 :
659 : END DO
660 2046 : CALL neighbor_list_iterator_release(nl_iterator)
661 :
662 4092 : DO i = 1, SIZE(matrix_h, 1)
663 41184 : DO img = 1, nimg
664 37092 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
665 39138 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
666 : END DO
667 : END DO
668 :
669 : ! EEQ forces (response and direct)
670 2046 : IF (calculate_forces) THEN
671 68 : CALL para_env%sum(dcharges)
672 136 : ALLOCATE (qlagrange(natom))
673 68 : CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
674 : END IF
675 :
676 2046 : kf = xtb_control%kf
677 2046 : enscale = xtb_control%enscale
678 2046 : erep = 0.0_dp
679 2046 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
680 :
681 2046 : esrb = 0.0_dp
682 2046 : CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
683 :
684 2046 : enonbonded = 0.0_dp
685 2046 : 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 2046 : erep = erep + esrb + enonbonded
695 2046 : IF (do_nonbonded) THEN
696 0 : CALL para_env%sum(enonbonded)
697 0 : energy%xtb_nonbonded = enonbonded
698 : END IF
699 2046 : CALL para_env%sum(esrb)
700 2046 : energy%srb = esrb
701 2046 : CALL para_env%sum(erep)
702 2046 : energy%repulsive = erep
703 :
704 : ! save EEQ charges
705 2046 : NULLIFY (eeq_q)
706 2046 : CALL get_qs_env(qs_env, eeq=eeq_q)
707 2046 : IF (ASSOCIATED(eeq_q)) THEN
708 1352 : CPASSERT(SIZE(eeq_q) == natom)
709 : ELSE
710 1388 : ALLOCATE (eeq_q(natom))
711 3574 : eeq_q(1:natom) = charges(1:natom)
712 : END IF
713 2046 : CALL set_qs_env(qs_env, eeq=eeq_q)
714 :
715 : ! deallocate coordination numbers
716 2046 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
717 :
718 : ! deallocate Huckel parameters
719 2046 : DEALLOCATE (huckel)
720 2046 : IF (calculate_forces) THEN
721 68 : DEALLOCATE (dhuckel, dqhuckel)
722 : END IF
723 : ! deallocate KAB parameters
724 2046 : DEALLOCATE (kijab)
725 :
726 : ! deallocate charges
727 2046 : DEALLOCATE (charges)
728 2046 : IF (calculate_forces) THEN
729 68 : DEALLOCATE (dcharges, qlagrange)
730 : END IF
731 :
732 : ! AO matrix outputs
733 2046 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
734 :
735 2046 : DEALLOCATE (basis_set_list)
736 2046 : IF (calculate_forces) THEN
737 68 : IF (SIZE(matrix_p, 1) == 2) THEN
738 16 : 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 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
742 16 : beta_scalar=-1.0_dp)
743 : END DO
744 : END IF
745 : END IF
746 :
747 2046 : CALL timestop(handle)
748 :
749 6138 : END SUBROUTINE build_gfn0_xtb_matrices
750 :
751 : ! **************************************************************************************************
752 : !> \brief ...
753 : !> \param qs_env ...
754 : !> \param calculate_forces ...
755 : ! **************************************************************************************************
756 3970 : SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
757 :
758 : TYPE(qs_environment_type), POINTER :: qs_env
759 : LOGICAL, INTENT(IN) :: calculate_forces
760 :
761 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn1_xtb_matrices'
762 :
763 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
764 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
765 : natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
766 : nsetb, sgfa, sgfb, za, zb
767 7940 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
768 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
769 : INTEGER, DIMENSION(3) :: cell
770 3970 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
771 3970 : npgfb, nsgfa, nsgfb
772 3970 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
773 3970 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
774 : LOGICAL :: defined, diagblock, do_nonbonded, found, &
775 : use_virial, xb_inter
776 : REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, enonbonded, &
777 : enscale, erep, etaa, etab, exb, f0, &
778 : f1, fhua, fhub, fhud, foab, hij, kf, &
779 : rcova, rcovab, rcovb, rrab
780 3970 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
781 7940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
782 3970 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
783 3970 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
784 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
785 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
786 3970 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
787 3970 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
788 3970 : scon_a, scon_b, wblock, zeta, zetb
789 3970 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
790 : TYPE(atprop_type), POINTER :: atprop
791 15880 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
792 : TYPE(cp_logger_type), POINTER :: logger
793 3970 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
794 3970 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
795 : TYPE(dft_control_type), POINTER :: dft_control
796 3970 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
797 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
798 : TYPE(kpoint_type), POINTER :: kpoints
799 : TYPE(mp_para_env_type), POINTER :: para_env
800 : TYPE(neighbor_list_iterator_p_type), &
801 3970 : DIMENSION(:), POINTER :: nl_iterator
802 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
803 3970 : POINTER :: sab_orb, sab_xtb_nonbond
804 3970 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
805 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
806 : TYPE(qs_energy_type), POINTER :: energy
807 3970 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
808 3970 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
809 : TYPE(qs_ks_env_type), POINTER :: ks_env
810 : TYPE(qs_rho_type), POINTER :: rho
811 : TYPE(virial_type), POINTER :: virial
812 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
813 : TYPE(xtb_control_type), POINTER :: xtb_control
814 :
815 3970 : CALL timeset(routineN, handle)
816 :
817 3970 : NULLIFY (logger, virial, atprop)
818 3970 : logger => cp_get_default_logger()
819 :
820 3970 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
821 3970 : qs_kind_set, sab_orb, ks_env)
822 :
823 : CALL get_qs_env(qs_env=qs_env, &
824 : ks_env=ks_env, &
825 : energy=energy, &
826 : atomic_kind_set=atomic_kind_set, &
827 : qs_kind_set=qs_kind_set, &
828 : matrix_h_kp=matrix_h, &
829 : matrix_s_kp=matrix_s, &
830 : para_env=para_env, &
831 : atprop=atprop, &
832 : dft_control=dft_control, &
833 3970 : sab_orb=sab_orb)
834 :
835 3970 : nkind = SIZE(atomic_kind_set)
836 3970 : xtb_control => dft_control%qs_control%xtb_control
837 3970 : xb_inter = xtb_control%xb_interaction
838 3970 : do_nonbonded = xtb_control%do_nonbonded
839 3970 : nimg = dft_control%nimages
840 3970 : nderivatives = 0
841 3970 : IF (calculate_forces) nderivatives = 1
842 3970 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
843 3970 : maxder = ncoset(nderivatives)
844 :
845 3970 : NULLIFY (particle_set)
846 3970 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
847 3970 : natom = SIZE(particle_set)
848 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
849 3970 : atom_of_kind=atom_of_kind, kind_of=kind_of)
850 :
851 3970 : IF (calculate_forces) THEN
852 526 : NULLIFY (rho, force, matrix_w)
853 : CALL get_qs_env(qs_env=qs_env, &
854 : rho=rho, matrix_w_kp=matrix_w, &
855 526 : virial=virial, force=force)
856 526 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
857 :
858 526 : IF (SIZE(matrix_p, 1) == 2) THEN
859 440 : DO img = 1, nimg
860 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
861 418 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
862 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
863 440 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
864 : END DO
865 : END IF
866 916 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
867 : END IF
868 : ! atomic energy decomposition
869 3970 : IF (atprop%energy) THEN
870 36 : CALL atprop_array_init(atprop%atecc, natom)
871 : END IF
872 :
873 3970 : NULLIFY (cell_to_index)
874 3970 : IF (nimg > 1) THEN
875 500 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
876 500 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
877 : END IF
878 :
879 : ! set up basis set lists
880 21520 : ALLOCATE (basis_set_list(nkind))
881 3970 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
882 :
883 : ! allocate overlap matrix
884 3970 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
885 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
886 3970 : sab_orb, .TRUE.)
887 3970 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
888 :
889 : ! initialize H matrix
890 3970 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
891 56396 : DO img = 1, nimg
892 52426 : ALLOCATE (matrix_h(1, img)%matrix)
893 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
894 52426 : name="HAMILTONIAN MATRIX")
895 56396 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
896 : END DO
897 3970 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
898 :
899 : ! Calculate coordination numbers
900 : ! needed for effective atomic energy levels (Eq. 12)
901 : ! code taken from D3 dispersion energy
902 3970 : CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
903 :
904 : ! vdW Potential
905 3970 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
906 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
907 3970 : energy%dispersion, calculate_forces)
908 :
909 : ! Calculate Huckel parameters
910 3970 : CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
911 :
912 : ! Calculate KAB parameters and electronegativity correction
913 3970 : CALL gfn1_kpair(qs_env, kijab)
914 :
915 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
916 3970 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
917 1140200 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
918 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
919 1136230 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
920 1136230 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
921 1136230 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
922 1136230 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
923 1136230 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
924 1136230 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
925 1136230 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
926 :
927 4544920 : dr = SQRT(SUM(rij(:)**2))
928 :
929 : ! atomic parameters
930 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
931 1136230 : lmax=lmaxa, nshell=nsa, kpoly=kpolya)
932 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
933 1136230 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
934 :
935 1136230 : IF (nimg == 1) THEN
936 : ic = 1
937 : ELSE
938 278022 : ic = cell_to_index(cell(1), cell(2), cell(3))
939 278022 : CPASSERT(ic > 0)
940 : END IF
941 :
942 1136230 : icol = MAX(iatom, jatom)
943 1136230 : irow = MIN(iatom, jatom)
944 1136230 : NULLIFY (sblock, fblock)
945 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
946 1136230 : row=irow, col=icol, BLOCK=sblock, found=found)
947 1136230 : CPASSERT(found)
948 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
949 1136230 : row=irow, col=icol, BLOCK=fblock, found=found)
950 1136230 : CPASSERT(found)
951 :
952 1136230 : IF (calculate_forces) THEN
953 256083 : NULLIFY (pblock)
954 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
955 256083 : row=irow, col=icol, block=pblock, found=found)
956 256083 : CPASSERT(found)
957 256083 : NULLIFY (wblock)
958 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
959 256083 : row=irow, col=icol, block=wblock, found=found)
960 256083 : CPASSERT(found)
961 1024332 : DO i = 2, 4
962 768249 : NULLIFY (dsblocks(i)%block)
963 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
964 768249 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
965 1024332 : CPASSERT(found)
966 : END DO
967 : END IF
968 :
969 : ! overlap
970 1136230 : basis_set_a => basis_set_list(ikind)%gto_basis_set
971 1136230 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
972 1136230 : basis_set_b => basis_set_list(jkind)%gto_basis_set
973 1136230 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
974 1136230 : atom_a = atom_of_kind(iatom)
975 1136230 : atom_b = atom_of_kind(jatom)
976 : ! basis ikind
977 1136230 : first_sgfa => basis_set_a%first_sgf
978 1136230 : la_max => basis_set_a%lmax
979 1136230 : la_min => basis_set_a%lmin
980 1136230 : npgfa => basis_set_a%npgf
981 1136230 : nseta = basis_set_a%nset
982 1136230 : nsgfa => basis_set_a%nsgf_set
983 1136230 : rpgfa => basis_set_a%pgf_radius
984 1136230 : set_radius_a => basis_set_a%set_radius
985 1136230 : scon_a => basis_set_a%scon
986 1136230 : zeta => basis_set_a%zet
987 : ! basis jkind
988 1136230 : first_sgfb => basis_set_b%first_sgf
989 1136230 : lb_max => basis_set_b%lmax
990 1136230 : lb_min => basis_set_b%lmin
991 1136230 : npgfb => basis_set_b%npgf
992 1136230 : nsetb = basis_set_b%nset
993 1136230 : nsgfb => basis_set_b%nsgf_set
994 1136230 : rpgfb => basis_set_b%pgf_radius
995 1136230 : set_radius_b => basis_set_b%set_radius
996 1136230 : scon_b => basis_set_b%scon
997 1136230 : zetb => basis_set_b%zet
998 :
999 1136230 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1000 9089840 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1001 5681150 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
1002 1136230 : sint = 0.0_dp
1003 :
1004 3638706 : DO iset = 1, nseta
1005 2502476 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1006 2502476 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1007 2502476 : sgfa = first_sgfa(1, iset)
1008 9308932 : DO jset = 1, nsetb
1009 5670226 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1010 4436131 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1011 4436131 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1012 4436131 : sgfb = first_sgfb(1, jset)
1013 4436131 : IF (calculate_forces) THEN
1014 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1015 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1016 1015916 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1017 : ELSE
1018 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1019 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1020 3420215 : rij, sab=oint(:, :, 1))
1021 : END IF
1022 : ! Contraction
1023 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1024 4436131 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1025 4436131 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1026 6938607 : IF (calculate_forces) THEN
1027 4063664 : DO i = 2, 4
1028 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1029 3047748 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1030 4063664 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1031 : END DO
1032 : END IF
1033 : END DO
1034 : END DO
1035 : ! forces W matrix
1036 1136230 : IF (calculate_forces) THEN
1037 1024332 : DO i = 1, 3
1038 1024332 : IF (iatom <= jatom) THEN
1039 16315077 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
1040 : ELSE
1041 11622747 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
1042 : END IF
1043 : END DO
1044 256083 : f1 = 2.0_dp
1045 1024332 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1046 1024332 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1047 256083 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
1048 162215 : IF (iatom == jatom) f1 = 1.0_dp
1049 162215 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1050 : END IF
1051 : END IF
1052 : ! update S matrix
1053 1136230 : IF (iatom <= jatom) THEN
1054 17332529 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1055 : ELSE
1056 12752339 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1057 : END IF
1058 1136230 : IF (calculate_forces) THEN
1059 1024332 : DO i = 2, 4
1060 1024332 : IF (iatom <= jatom) THEN
1061 16315077 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1062 : ELSE
1063 11456607 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
1064 : END IF
1065 : END DO
1066 : END IF
1067 :
1068 : ! Calculate Pi = Pia * Pib (Eq. 11)
1069 1136230 : rcovab = rcova + rcovb
1070 1136230 : rrab = SQRT(dr/rcovab)
1071 3638706 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1072 3638711 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1073 1136230 : IF (calculate_forces) THEN
1074 256083 : IF (dr > 1.e-6_dp) THEN
1075 253738 : drx = 0.5_dp/rrab/rcovab
1076 : ELSE
1077 : drx = 0.0_dp
1078 : END IF
1079 852983 : dpia(1:nsa) = drx*kpolya(1:nsa)
1080 852985 : dpib(1:nsb) = drx*kpolyb(1:nsb)
1081 : END IF
1082 :
1083 : ! diagonal block
1084 1136230 : diagblock = .FALSE.
1085 1136230 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
1086 : !
1087 : ! Eq. 10
1088 : !
1089 : IF (diagblock) THEN
1090 68012 : DO i = 1, natorb_a
1091 51820 : na = naoa(i)
1092 68012 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1093 : END DO
1094 : ELSE
1095 5565073 : DO j = 1, natorb_b
1096 4445035 : nb = naob(j)
1097 30029565 : DO i = 1, natorb_a
1098 24464492 : na = naoa(i)
1099 24464492 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1100 28909527 : IF (iatom <= jatom) THEN
1101 14108618 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1102 : ELSE
1103 10355874 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1104 : END IF
1105 : END DO
1106 : END DO
1107 : END IF
1108 1136230 : IF (calculate_forces) THEN
1109 256083 : f0 = 1.0_dp
1110 256083 : IF (irow == iatom) f0 = -1.0_dp
1111 : ! Derivative wrt coordination number
1112 256083 : fhua = 0.0_dp
1113 256083 : fhub = 0.0_dp
1114 256083 : fhud = 0.0_dp
1115 256083 : IF (diagblock) THEN
1116 10351 : DO i = 1, natorb_a
1117 8006 : la = laoa(i)
1118 8006 : na = naoa(i)
1119 10351 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1120 : END DO
1121 : ELSE
1122 1458630 : DO j = 1, natorb_b
1123 1204892 : lb = laob(j)
1124 1204892 : nb = naob(j)
1125 9264651 : DO i = 1, natorb_a
1126 7806021 : la = laoa(i)
1127 7806021 : na = naoa(i)
1128 7806021 : hij = 0.5_dp*pia(na)*pib(nb)
1129 9010913 : IF (iatom <= jatom) THEN
1130 4600616 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1131 4600616 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1132 : ELSE
1133 3205405 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1134 3205405 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
1135 : END IF
1136 : END DO
1137 : END DO
1138 253738 : IF (iatom /= jatom) THEN
1139 215940 : fhua = 2.0_dp*fhua
1140 215940 : fhub = 2.0_dp*fhub
1141 : END IF
1142 : END IF
1143 : ! iatom
1144 256083 : atom_a = atom_of_kind(iatom)
1145 8712499 : DO i = 1, dcnum(iatom)%neighbors
1146 8456416 : katom = dcnum(iatom)%nlist(i)
1147 8456416 : kkind = kind_of(katom)
1148 8456416 : atom_c = atom_of_kind(katom)
1149 33825664 : rik = dcnum(iatom)%rik(:, i)
1150 33825664 : drk = SQRT(SUM(rik(:)**2))
1151 8712499 : IF (drk > 1.e-3_dp) THEN
1152 33825664 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1153 33825664 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1154 33825664 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1155 33825664 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1156 33825664 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1157 33825664 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1158 8456416 : IF (use_virial) THEN
1159 23653192 : fdik = fdika + fdikb
1160 5913298 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1161 : END IF
1162 : END IF
1163 : END DO
1164 : ! jatom
1165 256083 : atom_b = atom_of_kind(jatom)
1166 8709634 : DO i = 1, dcnum(jatom)%neighbors
1167 8453551 : katom = dcnum(jatom)%nlist(i)
1168 8453551 : kkind = kind_of(katom)
1169 8453551 : atom_c = atom_of_kind(katom)
1170 33814204 : rik = dcnum(jatom)%rik(:, i)
1171 33814204 : drk = SQRT(SUM(rik(:)**2))
1172 8709634 : IF (drk > 1.e-3_dp) THEN
1173 33814204 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1174 33814204 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1175 33814204 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1176 8453551 : IF (use_virial) THEN
1177 5912224 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1178 : END IF
1179 : END IF
1180 : END DO
1181 : ! force from R dendent Huckel element: Pia*Pib
1182 256083 : IF (diagblock) THEN
1183 2345 : force_ab = 0._dp
1184 : ELSE
1185 253738 : n1 = SIZE(fblock, 1)
1186 253738 : n2 = SIZE(fblock, 2)
1187 1014952 : ALLOCATE (dfblock(n1, n2))
1188 253738 : dfblock = 0.0_dp
1189 1458630 : DO j = 1, natorb_b
1190 1204892 : lb = laob(j)
1191 1204892 : nb = naob(j)
1192 9264651 : DO i = 1, natorb_a
1193 7806021 : la = laoa(i)
1194 7806021 : na = naoa(i)
1195 7806021 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1196 9010913 : IF (iatom <= jatom) THEN
1197 4600616 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1198 : ELSE
1199 3205405 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1200 : END IF
1201 : END DO
1202 : END DO
1203 9209271 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
1204 1014952 : DO ir = 1, 3
1205 761214 : foab = 2.0_dp*dfp*rij(ir)/dr
1206 : ! force from overlap matrix contribution to H
1207 4375890 : DO j = 1, natorb_b
1208 3614676 : lb = laob(j)
1209 3614676 : nb = naob(j)
1210 27793953 : DO i = 1, natorb_a
1211 23418063 : la = laoa(i)
1212 23418063 : na = naoa(i)
1213 23418063 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1214 27032739 : IF (iatom <= jatom) THEN
1215 13801848 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1216 : ELSE
1217 9616215 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1218 : END IF
1219 : END DO
1220 : END DO
1221 1014952 : force_ab(ir) = foab
1222 : END DO
1223 253738 : DEALLOCATE (dfblock)
1224 : END IF
1225 : END IF
1226 :
1227 1136230 : IF (calculate_forces) THEN
1228 256083 : atom_a = atom_of_kind(iatom)
1229 256083 : atom_b = atom_of_kind(jatom)
1230 671679 : IF (irow == iatom) force_ab = -force_ab
1231 1024332 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1232 1024332 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1233 256083 : IF (use_virial) THEN
1234 163119 : f1 = 1.0_dp
1235 163119 : IF (iatom == jatom) f1 = 0.5_dp
1236 163119 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1237 : END IF
1238 : END IF
1239 :
1240 6821350 : DEALLOCATE (oint, owork, sint)
1241 :
1242 : END DO
1243 3970 : CALL neighbor_list_iterator_release(nl_iterator)
1244 :
1245 7940 : DO i = 1, SIZE(matrix_h, 1)
1246 60366 : DO img = 1, nimg
1247 52426 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1248 56396 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1249 : END DO
1250 : END DO
1251 :
1252 3970 : kf = xtb_control%kf
1253 3970 : enscale = xtb_control%enscale
1254 3970 : erep = 0.0_dp
1255 3970 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1256 :
1257 3970 : exb = 0.0_dp
1258 3970 : IF (xb_inter) THEN
1259 3928 : CALL xb_interaction(qs_env, exb, calculate_forces)
1260 : END IF
1261 :
1262 3970 : enonbonded = 0.0_dp
1263 3970 : IF (do_nonbonded) THEN
1264 : ! nonbonded interactions
1265 34 : NULLIFY (sab_xtb_nonbond)
1266 34 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
1267 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1268 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1269 : END IF
1270 :
1271 : ! set repulsive energy
1272 3970 : erep = erep + exb + enonbonded
1273 3970 : IF (xb_inter) THEN
1274 3928 : CALL para_env%sum(exb)
1275 3928 : energy%xtb_xb_inter = exb
1276 : END IF
1277 3970 : IF (do_nonbonded) THEN
1278 34 : CALL para_env%sum(enonbonded)
1279 34 : energy%xtb_nonbonded = enonbonded
1280 : END IF
1281 3970 : CALL para_env%sum(erep)
1282 3970 : energy%repulsive = erep
1283 :
1284 : ! deallocate coordination numbers
1285 3970 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1286 :
1287 : ! deallocate Huckel parameters
1288 3970 : DEALLOCATE (huckel)
1289 3970 : IF (calculate_forces) THEN
1290 526 : DEALLOCATE (dhuckel)
1291 : END IF
1292 : ! deallocate KAB parameters
1293 3970 : DEALLOCATE (kijab)
1294 :
1295 : ! AO matrix outputs
1296 3970 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1297 :
1298 3970 : DEALLOCATE (basis_set_list)
1299 3970 : IF (calculate_forces) THEN
1300 526 : IF (SIZE(matrix_p, 1) == 2) THEN
1301 440 : DO img = 1, nimg
1302 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
1303 418 : beta_scalar=-1.0_dp)
1304 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
1305 440 : beta_scalar=-1.0_dp)
1306 : END DO
1307 : END IF
1308 : END IF
1309 :
1310 3970 : CALL timestop(handle)
1311 :
1312 11910 : END SUBROUTINE build_gfn1_xtb_matrices
1313 :
1314 : ! **************************************************************************************************
1315 : !> \brief ...
1316 : !> \param qs_env ...
1317 : !> \param matrix_h ...
1318 : !> \param matrix_s ...
1319 : !> \param calculate_forces ...
1320 : ! **************************************************************************************************
1321 12032 : SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1322 : TYPE(qs_environment_type), POINTER :: qs_env
1323 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1324 : LOGICAL, INTENT(IN) :: calculate_forces
1325 :
1326 : INTEGER :: after, i, img, iw, nimg
1327 : LOGICAL :: norml1, norml2, omit_headers, use_arnoldi
1328 : REAL(KIND=dp), DIMENSION(2) :: condnum
1329 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1330 : TYPE(cp_logger_type), POINTER :: logger
1331 : TYPE(mp_para_env_type), POINTER :: para_env
1332 :
1333 6016 : logger => cp_get_default_logger()
1334 :
1335 6016 : CALL get_qs_env(qs_env, para_env=para_env)
1336 6016 : nimg = SIZE(matrix_h, 2)
1337 6016 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1338 6016 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1339 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
1340 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
1341 0 : extension=".Log")
1342 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1343 0 : after = MIN(MAX(after, 1), 16)
1344 0 : DO img = 1, nimg
1345 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
1346 0 : output_unit=iw, omit_headers=omit_headers)
1347 : END DO
1348 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
1349 : END IF
1350 :
1351 6016 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1352 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
1353 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
1354 0 : extension=".Log")
1355 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1356 0 : after = MIN(MAX(after, 1), 16)
1357 0 : DO img = 1, nimg
1358 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
1359 0 : output_unit=iw, omit_headers=omit_headers)
1360 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1361 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
1362 0 : DO i = 2, SIZE(matrix_s, 1)
1363 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
1364 0 : output_unit=iw, omit_headers=omit_headers)
1365 : END DO
1366 : END IF
1367 : END DO
1368 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
1369 : END IF
1370 :
1371 : ! *** Overlap condition number
1372 6016 : IF (.NOT. calculate_forces) THEN
1373 5422 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1374 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1375 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1376 4 : extension=".Log")
1377 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1378 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1379 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1380 4 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1381 4 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1382 : END IF
1383 : END IF
1384 :
1385 6016 : END SUBROUTINE ao_matrix_output
1386 :
1387 : END MODULE xtb_matrices
1388 :
|
