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 6288 : 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 6288 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
119 6288 : 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 4242 : 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 6288 : CPABORT("Unknown gfn_type")
130 : END SELECT
131 :
132 6288 : 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 : eeq_sparam = xtb_control%eeq_sparam
221 2046 : do_nonbonded = xtb_control%do_nonbonded
222 2046 : nimg = dft_control%nimages
223 2046 : nderivatives = 0
224 2046 : IF (calculate_forces) nderivatives = 1
225 2046 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
226 2046 : maxder = ncoset(nderivatives)
227 :
228 2046 : NULLIFY (particle_set)
229 2046 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
230 2046 : natom = SIZE(particle_set)
231 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
232 2046 : atom_of_kind=atom_of_kind, kind_of=kind_of)
233 :
234 2046 : IF (calculate_forces) THEN
235 68 : NULLIFY (rho, force, matrix_w)
236 : CALL get_qs_env(qs_env=qs_env, &
237 : rho=rho, matrix_w_kp=matrix_w, &
238 68 : virial=virial, force=force)
239 68 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
240 :
241 68 : IF (SIZE(matrix_p, 1) == 2) THEN
242 16 : DO img = 1, nimg
243 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
244 8 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
245 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
246 16 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
247 : END DO
248 : END IF
249 108 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
250 : END IF
251 : ! atomic energy decomposition
252 2046 : IF (atprop%energy) THEN
253 0 : CALL atprop_array_init(atprop%atecc, natom)
254 : END IF
255 :
256 2046 : NULLIFY (cell_to_index)
257 2046 : IF (nimg > 1) THEN
258 376 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
259 376 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
260 : END IF
261 :
262 : ! set up basis set lists
263 10872 : ALLOCATE (basis_set_list(nkind))
264 2046 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
265 :
266 : ! allocate overlap matrix
267 2046 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
268 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
269 2046 : sab_orb, .TRUE.)
270 2046 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
271 :
272 : ! initialize H matrix
273 2046 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
274 39142 : DO img = 1, nimg
275 37096 : ALLOCATE (matrix_h(1, img)%matrix)
276 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
277 37096 : name="HAMILTONIAN MATRIX")
278 39142 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
279 : END DO
280 2046 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
281 :
282 : ! Calculate coordination numbers
283 : ! needed for effective atomic energy levels
284 : ! code taken from D3 dispersion energy
285 2046 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
286 :
287 6138 : ALLOCATE (charges(natom))
288 2046 : charges = 0.0_dp
289 2046 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
290 2046 : IF (calculate_forces) THEN
291 136 : ALLOCATE (dcharges(natom))
292 416 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
293 : END IF
294 2046 : energy%eeq = eeq_energy
295 2046 : energy%efield = ef_energy
296 :
297 2046 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
298 : ! prepare charges (needed for D4)
299 2046 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
300 858 : dispersion_env%ext_charges = .TRUE.
301 858 : IF (ASSOCIATED(dispersion_env%charges)) DEALLOCATE (dispersion_env%charges)
302 1716 : ALLOCATE (dispersion_env%charges(natom))
303 4320 : dispersion_env%charges = charges
304 858 : IF (calculate_forces) THEN
305 12 : IF (ASSOCIATED(dispersion_env%dcharges)) DEALLOCATE (dispersion_env%dcharges)
306 24 : ALLOCATE (dispersion_env%dcharges(natom))
307 60 : dispersion_env%dcharges = 0.0_dp
308 : END IF
309 : END IF
310 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
311 2046 : energy%dispersion, calculate_forces)
312 2046 : IF (calculate_forces) THEN
313 68 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
314 60 : dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
315 : END IF
316 : END IF
317 :
318 : ! Calculate Huckel parameters
319 2046 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
320 :
321 : ! Calculate KAB parameters and electronegativity correction
322 2046 : CALL gfn0_kpair(qs_env, kijab)
323 :
324 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
325 2046 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
326 493308 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
327 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
328 491262 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
329 491262 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
330 491262 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
331 491262 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
332 491262 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
333 491262 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
334 491262 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
335 :
336 1965048 : dr = SQRT(SUM(rij(:)**2))
337 :
338 : ! atomic parameters
339 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
340 491262 : lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
341 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
342 491262 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
343 :
344 491262 : IF (nimg == 1) THEN
345 : ic = 1
346 : ELSE
347 269927 : ic = cell_to_index(cell(1), cell(2), cell(3))
348 269927 : CPASSERT(ic > 0)
349 : END IF
350 :
351 491262 : icol = MAX(iatom, jatom)
352 491262 : irow = MIN(iatom, jatom)
353 491262 : NULLIFY (sblock, fblock)
354 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
355 491262 : row=irow, col=icol, BLOCK=sblock, found=found)
356 491262 : CPASSERT(found)
357 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
358 491262 : row=irow, col=icol, BLOCK=fblock, found=found)
359 491262 : CPASSERT(found)
360 :
361 491262 : IF (calculate_forces) THEN
362 28766 : NULLIFY (pblock)
363 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
364 28766 : row=irow, col=icol, block=pblock, found=found)
365 28766 : CPASSERT(ASSOCIATED(pblock))
366 28766 : NULLIFY (wblock)
367 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
368 28766 : row=irow, col=icol, block=wblock, found=found)
369 28766 : CPASSERT(ASSOCIATED(wblock))
370 115064 : DO i = 2, 4
371 86298 : NULLIFY (dsblocks(i)%block)
372 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
373 86298 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
374 115064 : CPASSERT(found)
375 : END DO
376 : END IF
377 :
378 : ! overlap
379 491262 : basis_set_a => basis_set_list(ikind)%gto_basis_set
380 491262 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
381 491262 : basis_set_b => basis_set_list(jkind)%gto_basis_set
382 491262 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
383 491262 : atom_a = atom_of_kind(iatom)
384 491262 : atom_b = atom_of_kind(jatom)
385 : ! basis ikind
386 491262 : first_sgfa => basis_set_a%first_sgf
387 491262 : la_max => basis_set_a%lmax
388 491262 : la_min => basis_set_a%lmin
389 491262 : npgfa => basis_set_a%npgf
390 491262 : nseta = basis_set_a%nset
391 491262 : nsgfa => basis_set_a%nsgf_set
392 491262 : rpgfa => basis_set_a%pgf_radius
393 491262 : set_radius_a => basis_set_a%set_radius
394 491262 : scon_a => basis_set_a%scon
395 491262 : zeta => basis_set_a%zet
396 : ! basis jkind
397 491262 : first_sgfb => basis_set_b%first_sgf
398 491262 : lb_max => basis_set_b%lmax
399 491262 : lb_min => basis_set_b%lmin
400 491262 : npgfb => basis_set_b%npgf
401 491262 : nsetb = basis_set_b%nset
402 491262 : nsgfb => basis_set_b%nsgf_set
403 491262 : rpgfb => basis_set_b%pgf_radius
404 491262 : set_radius_b => basis_set_b%set_radius
405 491262 : scon_b => basis_set_b%scon
406 491262 : zetb => basis_set_b%zet
407 :
408 491262 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
409 3930096 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
410 2456310 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
411 491262 : sint = 0.0_dp
412 :
413 1899871 : DO iset = 1, nseta
414 1408609 : ncoa = npgfa(iset)*ncoset(la_max(iset))
415 1408609 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
416 1408609 : sgfa = first_sgfa(1, iset)
417 5947231 : DO jset = 1, nsetb
418 4047360 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
419 2088629 : ncob = npgfb(jset)*ncoset(lb_max(jset))
420 2088629 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
421 2088629 : sgfb = first_sgfb(1, jset)
422 2088629 : IF (calculate_forces) THEN
423 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
424 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
425 123944 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
426 : ELSE
427 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
428 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
429 1964685 : rij, sab=oint(:, :, 1))
430 : END IF
431 : ! Contraction
432 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
433 2088629 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
434 2088629 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
435 3497238 : IF (calculate_forces) THEN
436 495776 : DO i = 2, 4
437 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
438 371832 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
439 495776 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
440 : END DO
441 : END IF
442 : END DO
443 : END DO
444 : ! forces W matrix
445 491262 : IF (calculate_forces) THEN
446 115064 : DO i = 1, 3
447 115064 : IF (iatom <= jatom) THEN
448 4051116 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
449 : ELSE
450 3171456 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
451 : END IF
452 : END DO
453 28766 : f1 = 2.0_dp
454 115064 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
455 115064 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
456 28766 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
457 27840 : IF (iatom == jatom) f1 = 1.0_dp
458 27840 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
459 : END IF
460 : END IF
461 : ! update S matrix
462 491262 : IF (iatom <= jatom) THEN
463 22018173 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
464 : ELSE
465 16993369 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
466 : END IF
467 491262 : IF (calculate_forces) THEN
468 115064 : DO i = 2, 4
469 115064 : IF (iatom <= jatom) THEN
470 4051116 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
471 : ELSE
472 3187608 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
473 : END IF
474 : END DO
475 : END IF
476 :
477 : ! Calculate Pi = Pia * Pib (Eq. 11)
478 491262 : rcovab = rcova + rcovb
479 491262 : rrab = SQRT(dr/rcovab)
480 1899871 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
481 1895231 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
482 491262 : IF (calculate_forces) THEN
483 28766 : IF (dr > 1.e-6_dp) THEN
484 28592 : drx = 0.5_dp/rrab/rcovab
485 : ELSE
486 : drx = 0.0_dp
487 : END IF
488 112664 : dpia(1:nsa) = drx*kpolya(1:nsa)
489 112572 : dpib(1:nsb) = drx*kpolyb(1:nsb)
490 : END IF
491 :
492 : ! diagonal block
493 491262 : diagblock = .FALSE.
494 491262 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
495 : !
496 : ! Eq. 10
497 : !
498 : IF (diagblock) THEN
499 30959 : DO i = 1, natorb_a
500 26103 : na = naoa(i)
501 30959 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
502 : END DO
503 : ELSE
504 4509076 : DO j = 1, natorb_b
505 4022670 : nb = naob(j)
506 38662451 : DO i = 1, natorb_a
507 34153375 : na = naoa(i)
508 34153375 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
509 38176045 : IF (iatom <= jatom) THEN
510 19223733 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
511 : ELSE
512 14929642 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
513 : END IF
514 : END DO
515 : END DO
516 : END IF
517 491262 : IF (calculate_forces) THEN
518 28766 : f0 = 1.0_dp
519 28766 : IF (irow == iatom) f0 = -1.0_dp
520 28766 : f2 = 1.0_dp
521 28766 : IF (iatom /= jatom) f2 = 2.0_dp
522 : ! Derivative wrt coordination number
523 28766 : fhua = 0.0_dp
524 28766 : fhub = 0.0_dp
525 28766 : fhud = 0.0_dp
526 28766 : fqa = 0.0_dp
527 28766 : fqb = 0.0_dp
528 28766 : IF (diagblock) THEN
529 1168 : DO i = 1, natorb_a
530 994 : la = laoa(i)
531 994 : na = naoa(i)
532 994 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
533 1168 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
534 : END DO
535 174 : dcharges(iatom) = dcharges(iatom) + fqa
536 : ELSE
537 273142 : DO j = 1, natorb_b
538 244550 : lb = laob(j)
539 244550 : nb = naob(j)
540 2399054 : DO i = 1, natorb_a
541 2125912 : la = laoa(i)
542 2125912 : na = naoa(i)
543 2125912 : hij = 0.5_dp*pia(na)*pib(nb)
544 2125912 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
545 2370462 : IF (iatom <= jatom) THEN
546 1187288 : fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
547 1187288 : fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
548 1187288 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
549 1187288 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
550 : ELSE
551 938624 : fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
552 938624 : fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
553 938624 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
554 938624 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
555 : END IF
556 : END DO
557 : END DO
558 28592 : dcharges(iatom) = dcharges(iatom) + fqa
559 28592 : dcharges(jatom) = dcharges(jatom) + fqb
560 : END IF
561 : ! iatom
562 28766 : atom_a = atom_of_kind(iatom)
563 445324 : DO i = 1, dcnum(iatom)%neighbors
564 416558 : katom = dcnum(iatom)%nlist(i)
565 416558 : kkind = kind_of(katom)
566 416558 : atom_c = atom_of_kind(katom)
567 1666232 : rik = dcnum(iatom)%rik(:, i)
568 1666232 : drk = SQRT(SUM(rik(:)**2))
569 445324 : IF (drk > 1.e-3_dp) THEN
570 1666232 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
571 1666232 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
572 1666232 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
573 1666232 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
574 1666232 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
575 1666232 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
576 416558 : IF (use_virial) THEN
577 1656556 : fdik = fdika + fdikb
578 414139 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
579 : END IF
580 : END IF
581 : END DO
582 : ! jatom
583 28766 : atom_b = atom_of_kind(jatom)
584 444486 : DO i = 1, dcnum(jatom)%neighbors
585 415720 : katom = dcnum(jatom)%nlist(i)
586 415720 : kkind = kind_of(katom)
587 415720 : atom_c = atom_of_kind(katom)
588 1662880 : rik = dcnum(jatom)%rik(:, i)
589 1662880 : drk = SQRT(SUM(rik(:)**2))
590 444486 : IF (drk > 1.e-3_dp) THEN
591 1662880 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
592 1662880 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
593 1662880 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
594 415720 : IF (use_virial) THEN
595 413391 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
596 : END IF
597 : END IF
598 : END DO
599 : ! force from R dendent Huckel element: Pia*Pib
600 28766 : IF (diagblock) THEN
601 174 : force_ab = 0._dp
602 : ELSE
603 28592 : n1 = SIZE(fblock, 1)
604 28592 : n2 = SIZE(fblock, 2)
605 114368 : ALLOCATE (dfblock(n1, n2))
606 28592 : dfblock = 0.0_dp
607 273142 : DO j = 1, natorb_b
608 244550 : lb = laob(j)
609 244550 : nb = naob(j)
610 2399054 : DO i = 1, natorb_a
611 2125912 : la = laoa(i)
612 2125912 : na = naoa(i)
613 2125912 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
614 2370462 : IF (iatom <= jatom) THEN
615 1187288 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
616 : ELSE
617 938624 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
618 : END IF
619 : END DO
620 : END DO
621 2404438 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
622 114368 : DO ir = 1, 3
623 85776 : foab = 2.0_dp*dfp*rij(ir)/dr
624 : ! force from overlap matrix contribution to H
625 819426 : DO j = 1, natorb_b
626 733650 : lb = laob(j)
627 733650 : nb = naob(j)
628 7197162 : DO i = 1, natorb_a
629 6377736 : la = laoa(i)
630 6377736 : na = naoa(i)
631 6377736 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
632 7111386 : IF (iatom <= jatom) THEN
633 3561864 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
634 : ELSE
635 2815872 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
636 : END IF
637 : END DO
638 : END DO
639 114368 : force_ab(ir) = foab
640 : END DO
641 28592 : DEALLOCATE (dfblock)
642 : END IF
643 : END IF
644 :
645 491262 : IF (calculate_forces) THEN
646 28766 : atom_a = atom_of_kind(iatom)
647 28766 : atom_b = atom_of_kind(jatom)
648 76934 : IF (irow == iatom) force_ab = -force_ab
649 115064 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
650 115064 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
651 28766 : IF (use_virial) THEN
652 27933 : f1 = 1.0_dp
653 27933 : IF (iatom == jatom) f1 = 0.5_dp
654 27933 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
655 : END IF
656 : END IF
657 :
658 2949618 : DEALLOCATE (oint, owork, sint)
659 :
660 : END DO
661 2046 : CALL neighbor_list_iterator_release(nl_iterator)
662 :
663 4092 : DO i = 1, SIZE(matrix_h, 1)
664 41188 : DO img = 1, nimg
665 37096 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
666 39142 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
667 : END DO
668 : END DO
669 :
670 : ! EEQ forces (response and direct)
671 2046 : IF (calculate_forces) THEN
672 68 : CALL para_env%sum(dcharges)
673 136 : ALLOCATE (qlagrange(natom))
674 68 : CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
675 : END IF
676 :
677 2046 : kf = xtb_control%kf
678 2046 : enscale = xtb_control%enscale
679 2046 : erep = 0.0_dp
680 2046 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
681 :
682 2046 : esrb = 0.0_dp
683 2046 : CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
684 :
685 2046 : enonbonded = 0.0_dp
686 2046 : IF (do_nonbonded) THEN
687 : ! nonbonded interactions
688 0 : NULLIFY (sab_xtb_nonbond)
689 0 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
690 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
691 0 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
692 : END IF
693 :
694 : ! set repulsive energy
695 2046 : erep = erep + esrb + enonbonded
696 2046 : IF (do_nonbonded) THEN
697 0 : CALL para_env%sum(enonbonded)
698 0 : energy%xtb_nonbonded = enonbonded
699 : END IF
700 2046 : CALL para_env%sum(esrb)
701 2046 : energy%srb = esrb
702 2046 : CALL para_env%sum(erep)
703 2046 : energy%repulsive = erep
704 :
705 : ! save EEQ charges
706 2046 : NULLIFY (eeq_q)
707 2046 : CALL get_qs_env(qs_env, eeq=eeq_q)
708 2046 : IF (ASSOCIATED(eeq_q)) THEN
709 1352 : CPASSERT(SIZE(eeq_q) == natom)
710 : ELSE
711 1388 : ALLOCATE (eeq_q(natom))
712 3574 : eeq_q(1:natom) = charges(1:natom)
713 : END IF
714 2046 : CALL set_qs_env(qs_env, eeq=eeq_q)
715 :
716 : ! deallocate coordination numbers
717 2046 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
718 :
719 : ! deallocate Huckel parameters
720 2046 : DEALLOCATE (huckel)
721 2046 : IF (calculate_forces) THEN
722 68 : DEALLOCATE (dhuckel, dqhuckel)
723 : END IF
724 : ! deallocate KAB parameters
725 2046 : DEALLOCATE (kijab)
726 :
727 : ! deallocate charges
728 2046 : DEALLOCATE (charges)
729 2046 : IF (calculate_forces) THEN
730 68 : DEALLOCATE (dcharges, qlagrange)
731 : END IF
732 :
733 : ! AO matrix outputs
734 2046 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
735 :
736 2046 : DEALLOCATE (basis_set_list)
737 2046 : IF (calculate_forces) THEN
738 68 : IF (SIZE(matrix_p, 1) == 2) THEN
739 16 : DO img = 1, nimg
740 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
741 8 : beta_scalar=-1.0_dp)
742 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
743 16 : beta_scalar=-1.0_dp)
744 : END DO
745 : END IF
746 : END IF
747 :
748 2046 : CALL timestop(handle)
749 :
750 6138 : END SUBROUTINE build_gfn0_xtb_matrices
751 :
752 : ! **************************************************************************************************
753 : !> \brief ...
754 : !> \param qs_env ...
755 : !> \param calculate_forces ...
756 : ! **************************************************************************************************
757 4242 : SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
758 :
759 : TYPE(qs_environment_type), POINTER :: qs_env
760 : LOGICAL, INTENT(IN) :: calculate_forces
761 :
762 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn1_xtb_matrices'
763 :
764 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
765 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
766 : natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
767 : nsetb, sgfa, sgfb, za, zb
768 8484 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
769 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
770 : INTEGER, DIMENSION(3) :: cell
771 4242 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
772 4242 : npgfb, nsgfa, nsgfb
773 4242 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
774 4242 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
775 : LOGICAL :: defined, diagblock, do_nonbonded, found, &
776 : use_virial, xb_inter
777 : REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, enonbonded, &
778 : enscale, erep, etaa, etab, exb, f0, &
779 : f1, fhua, fhub, fhud, foab, hij, kf, &
780 : rcova, rcovab, rcovb, rrab
781 4242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
782 8484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
783 4242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
784 4242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
785 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
786 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
787 4242 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
788 4242 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
789 4242 : scon_a, scon_b, wblock, zeta, zetb
790 4242 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
791 : TYPE(atprop_type), POINTER :: atprop
792 16968 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
793 : TYPE(cp_logger_type), POINTER :: logger
794 4242 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
795 4242 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
796 : TYPE(dft_control_type), POINTER :: dft_control
797 4242 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
798 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
799 : TYPE(kpoint_type), POINTER :: kpoints
800 : TYPE(mp_para_env_type), POINTER :: para_env
801 : TYPE(neighbor_list_iterator_p_type), &
802 4242 : DIMENSION(:), POINTER :: nl_iterator
803 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
804 4242 : POINTER :: sab_orb, sab_xtb_nonbond
805 4242 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
806 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
807 : TYPE(qs_energy_type), POINTER :: energy
808 4242 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
809 4242 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
810 : TYPE(qs_ks_env_type), POINTER :: ks_env
811 : TYPE(qs_rho_type), POINTER :: rho
812 : TYPE(virial_type), POINTER :: virial
813 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
814 : TYPE(xtb_control_type), POINTER :: xtb_control
815 :
816 4242 : CALL timeset(routineN, handle)
817 :
818 4242 : NULLIFY (logger, virial, atprop)
819 4242 : logger => cp_get_default_logger()
820 :
821 4242 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
822 4242 : qs_kind_set, sab_orb, ks_env)
823 :
824 : CALL get_qs_env(qs_env=qs_env, &
825 : ks_env=ks_env, &
826 : energy=energy, &
827 : atomic_kind_set=atomic_kind_set, &
828 : qs_kind_set=qs_kind_set, &
829 : matrix_h_kp=matrix_h, &
830 : matrix_s_kp=matrix_s, &
831 : para_env=para_env, &
832 : atprop=atprop, &
833 : dft_control=dft_control, &
834 4242 : sab_orb=sab_orb)
835 :
836 4242 : nkind = SIZE(atomic_kind_set)
837 4242 : xtb_control => dft_control%qs_control%xtb_control
838 4242 : xb_inter = xtb_control%xb_interaction
839 4242 : do_nonbonded = xtb_control%do_nonbonded
840 4242 : nimg = dft_control%nimages
841 4242 : nderivatives = 0
842 4242 : IF (calculate_forces) nderivatives = 1
843 4242 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
844 4242 : maxder = ncoset(nderivatives)
845 :
846 4242 : NULLIFY (particle_set)
847 4242 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
848 4242 : natom = SIZE(particle_set)
849 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
850 4242 : atom_of_kind=atom_of_kind, kind_of=kind_of)
851 :
852 4242 : IF (calculate_forces) THEN
853 564 : NULLIFY (rho, force, matrix_w)
854 : CALL get_qs_env(qs_env=qs_env, &
855 : rho=rho, matrix_w_kp=matrix_w, &
856 564 : virial=virial, force=force)
857 564 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
858 :
859 564 : IF (SIZE(matrix_p, 1) == 2) THEN
860 752 : DO img = 1, nimg
861 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
862 710 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
863 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
864 752 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
865 : END DO
866 : END IF
867 974 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
868 : END IF
869 : ! atomic energy decomposition
870 4242 : IF (atprop%energy) THEN
871 36 : CALL atprop_array_init(atprop%atecc, natom)
872 : END IF
873 :
874 4242 : NULLIFY (cell_to_index)
875 4242 : IF (nimg > 1) THEN
876 550 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
877 550 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
878 : END IF
879 :
880 : ! set up basis set lists
881 23024 : ALLOCATE (basis_set_list(nkind))
882 4242 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
883 :
884 : ! allocate overlap matrix
885 4242 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
886 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
887 4242 : sab_orb, .TRUE.)
888 4242 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
889 :
890 : ! initialize H matrix
891 4242 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
892 60348 : DO img = 1, nimg
893 56106 : ALLOCATE (matrix_h(1, img)%matrix)
894 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
895 56106 : name="HAMILTONIAN MATRIX")
896 60348 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
897 : END DO
898 4242 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
899 :
900 : ! Calculate coordination numbers
901 : ! needed for effective atomic energy levels (Eq. 12)
902 : ! code taken from D3 dispersion energy
903 4242 : CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
904 :
905 : ! vdW Potential
906 4242 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
907 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
908 4242 : energy%dispersion, calculate_forces)
909 :
910 : ! Calculate Huckel parameters
911 4242 : CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
912 :
913 : ! Calculate KAB parameters and electronegativity correction
914 4242 : CALL gfn1_kpair(qs_env, kijab)
915 :
916 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
917 4242 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
918 1162992 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
919 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
920 1158750 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
921 1158750 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
922 1158750 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
923 1158750 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
924 1158750 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
925 1158750 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
926 1158750 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
927 :
928 4635000 : dr = SQRT(SUM(rij(:)**2))
929 :
930 : ! atomic parameters
931 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
932 1158750 : lmax=lmaxa, nshell=nsa, kpoly=kpolya)
933 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
934 1158750 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
935 :
936 1158750 : IF (nimg == 1) THEN
937 : ic = 1
938 : ELSE
939 287872 : ic = cell_to_index(cell(1), cell(2), cell(3))
940 287872 : CPASSERT(ic > 0)
941 : END IF
942 :
943 1158750 : icol = MAX(iatom, jatom)
944 1158750 : irow = MIN(iatom, jatom)
945 1158750 : NULLIFY (sblock, fblock)
946 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
947 1158750 : row=irow, col=icol, BLOCK=sblock, found=found)
948 1158750 : CPASSERT(found)
949 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
950 1158750 : row=irow, col=icol, BLOCK=fblock, found=found)
951 1158750 : CPASSERT(found)
952 :
953 1158750 : IF (calculate_forces) THEN
954 258799 : NULLIFY (pblock)
955 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
956 258799 : row=irow, col=icol, block=pblock, found=found)
957 258799 : CPASSERT(found)
958 258799 : NULLIFY (wblock)
959 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
960 258799 : row=irow, col=icol, block=wblock, found=found)
961 258799 : CPASSERT(found)
962 1035196 : DO i = 2, 4
963 776397 : NULLIFY (dsblocks(i)%block)
964 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
965 776397 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
966 1035196 : CPASSERT(found)
967 : END DO
968 : END IF
969 :
970 : ! overlap
971 1158750 : basis_set_a => basis_set_list(ikind)%gto_basis_set
972 1158750 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
973 1158750 : basis_set_b => basis_set_list(jkind)%gto_basis_set
974 1158750 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
975 1158750 : atom_a = atom_of_kind(iatom)
976 1158750 : atom_b = atom_of_kind(jatom)
977 : ! basis ikind
978 1158750 : first_sgfa => basis_set_a%first_sgf
979 1158750 : la_max => basis_set_a%lmax
980 1158750 : la_min => basis_set_a%lmin
981 1158750 : npgfa => basis_set_a%npgf
982 1158750 : nseta = basis_set_a%nset
983 1158750 : nsgfa => basis_set_a%nsgf_set
984 1158750 : rpgfa => basis_set_a%pgf_radius
985 1158750 : set_radius_a => basis_set_a%set_radius
986 1158750 : scon_a => basis_set_a%scon
987 1158750 : zeta => basis_set_a%zet
988 : ! basis jkind
989 1158750 : first_sgfb => basis_set_b%first_sgf
990 1158750 : lb_max => basis_set_b%lmax
991 1158750 : lb_min => basis_set_b%lmin
992 1158750 : npgfb => basis_set_b%npgf
993 1158750 : nsetb = basis_set_b%nset
994 1158750 : nsgfb => basis_set_b%nsgf_set
995 1158750 : rpgfb => basis_set_b%pgf_radius
996 1158750 : set_radius_b => basis_set_b%set_radius
997 1158750 : scon_b => basis_set_b%scon
998 1158750 : zetb => basis_set_b%zet
999 :
1000 1158750 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1001 9270000 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1002 5793750 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
1003 1158750 : sint = 0.0_dp
1004 :
1005 3708577 : DO iset = 1, nseta
1006 2549827 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1007 2549827 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1008 2549827 : sgfa = first_sgfa(1, iset)
1009 9480340 : DO jset = 1, nsetb
1010 5771763 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1011 4523918 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1012 4523918 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1013 4523918 : sgfb = first_sgfb(1, jset)
1014 4523918 : IF (calculate_forces) THEN
1015 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1016 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1017 1028585 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1018 : ELSE
1019 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1020 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1021 3495333 : rij, sab=oint(:, :, 1))
1022 : END IF
1023 : ! Contraction
1024 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1025 4523918 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1026 4523918 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1027 7073745 : IF (calculate_forces) THEN
1028 4114340 : DO i = 2, 4
1029 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1030 3085755 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1031 4114340 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1032 : END DO
1033 : END IF
1034 : END DO
1035 : END DO
1036 : ! forces W matrix
1037 1158750 : IF (calculate_forces) THEN
1038 1035196 : DO i = 1, 3
1039 1035196 : IF (iatom <= jatom) THEN
1040 16537518 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
1041 : ELSE
1042 11758932 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
1043 : END IF
1044 : END DO
1045 258799 : f1 = 2.0_dp
1046 1035196 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1047 1035196 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1048 258799 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
1049 164110 : IF (iatom == jatom) f1 = 1.0_dp
1050 164110 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1051 : END IF
1052 : END IF
1053 : ! update S matrix
1054 1158750 : IF (iatom <= jatom) THEN
1055 17592067 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1056 : ELSE
1057 12917773 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1058 : END IF
1059 1158750 : IF (calculate_forces) THEN
1060 1035196 : DO i = 2, 4
1061 1035196 : IF (iatom <= jatom) THEN
1062 16537518 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1063 : ELSE
1064 11591007 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
1065 : END IF
1066 : END DO
1067 : END IF
1068 :
1069 : ! Calculate Pi = Pia * Pib (Eq. 11)
1070 1158750 : rcovab = rcova + rcovb
1071 1158750 : rrab = SQRT(dr/rcovab)
1072 3708577 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1073 3708566 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1074 1158750 : IF (calculate_forces) THEN
1075 258799 : IF (dr > 1.e-6_dp) THEN
1076 256395 : drx = 0.5_dp/rrab/rcovab
1077 : ELSE
1078 : drx = 0.0_dp
1079 : END IF
1080 862245 : dpia(1:nsa) = drx*kpolya(1:nsa)
1081 862246 : dpib(1:nsb) = drx*kpolyb(1:nsb)
1082 : END IF
1083 :
1084 : ! diagonal block
1085 1158750 : diagblock = .FALSE.
1086 1158750 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
1087 : !
1088 : ! Eq. 10
1089 : !
1090 : IF (diagblock) THEN
1091 70582 : DO i = 1, natorb_a
1092 53785 : na = naoa(i)
1093 70582 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1094 : END DO
1095 : ELSE
1096 5659996 : DO j = 1, natorb_b
1097 4518043 : nb = naob(j)
1098 30450680 : DO i = 1, natorb_a
1099 24790684 : na = naoa(i)
1100 24790684 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1101 29308727 : IF (iatom <= jatom) THEN
1102 14304788 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1103 : ELSE
1104 10485896 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1105 : END IF
1106 : END DO
1107 : END DO
1108 : END IF
1109 1158750 : IF (calculate_forces) THEN
1110 258799 : f0 = 1.0_dp
1111 258799 : IF (irow == iatom) f0 = -1.0_dp
1112 : ! Derivative wrt coordination number
1113 258799 : fhua = 0.0_dp
1114 258799 : fhub = 0.0_dp
1115 258799 : fhud = 0.0_dp
1116 258799 : IF (diagblock) THEN
1117 10631 : DO i = 1, natorb_a
1118 8227 : la = laoa(i)
1119 8227 : na = naoa(i)
1120 10631 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1121 : END DO
1122 : ELSE
1123 1475463 : DO j = 1, natorb_b
1124 1219068 : lb = laob(j)
1125 1219068 : nb = naob(j)
1126 9382744 : DO i = 1, natorb_a
1127 7907281 : la = laoa(i)
1128 7907281 : na = naoa(i)
1129 7907281 : hij = 0.5_dp*pia(na)*pib(nb)
1130 9126349 : IF (iatom <= jatom) THEN
1131 4663312 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1132 4663312 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1133 : ELSE
1134 3243969 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1135 3243969 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
1136 : END IF
1137 : END DO
1138 : END DO
1139 256395 : IF (iatom /= jatom) THEN
1140 217285 : fhua = 2.0_dp*fhua
1141 217285 : fhub = 2.0_dp*fhub
1142 : END IF
1143 : END IF
1144 : ! iatom
1145 258799 : atom_a = atom_of_kind(iatom)
1146 8758876 : DO i = 1, dcnum(iatom)%neighbors
1147 8500077 : katom = dcnum(iatom)%nlist(i)
1148 8500077 : kkind = kind_of(katom)
1149 8500077 : atom_c = atom_of_kind(katom)
1150 34000308 : rik = dcnum(iatom)%rik(:, i)
1151 34000308 : drk = SQRT(SUM(rik(:)**2))
1152 8758876 : IF (drk > 1.e-3_dp) THEN
1153 34000308 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1154 34000308 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1155 34000308 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1156 34000308 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1157 34000308 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1158 34000308 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1159 8500077 : IF (use_virial) THEN
1160 23804424 : fdik = fdika + fdikb
1161 5951106 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1162 : END IF
1163 : END IF
1164 : END DO
1165 : ! jatom
1166 258799 : atom_b = atom_of_kind(jatom)
1167 8755735 : DO i = 1, dcnum(jatom)%neighbors
1168 8496936 : katom = dcnum(jatom)%nlist(i)
1169 8496936 : kkind = kind_of(katom)
1170 8496936 : atom_c = atom_of_kind(katom)
1171 33987744 : rik = dcnum(jatom)%rik(:, i)
1172 33987744 : drk = SQRT(SUM(rik(:)**2))
1173 8755735 : IF (drk > 1.e-3_dp) THEN
1174 33987744 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1175 33987744 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1176 33987744 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1177 8496936 : IF (use_virial) THEN
1178 5949864 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1179 : END IF
1180 : END IF
1181 : END DO
1182 : ! force from R dendent Huckel element: Pia*Pib
1183 258799 : IF (diagblock) THEN
1184 2404 : force_ab = 0._dp
1185 : ELSE
1186 256395 : n1 = SIZE(fblock, 1)
1187 256395 : n2 = SIZE(fblock, 2)
1188 1025580 : ALLOCATE (dfblock(n1, n2))
1189 256395 : dfblock = 0.0_dp
1190 1475463 : DO j = 1, natorb_b
1191 1219068 : lb = laob(j)
1192 1219068 : nb = naob(j)
1193 9382744 : DO i = 1, natorb_a
1194 7907281 : la = laoa(i)
1195 7907281 : na = naoa(i)
1196 7907281 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1197 9126349 : IF (iatom <= jatom) THEN
1198 4663312 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1199 : ELSE
1200 3243969 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1201 : END IF
1202 : END DO
1203 : END DO
1204 9326769 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
1205 1025580 : DO ir = 1, 3
1206 769185 : foab = 2.0_dp*dfp*rij(ir)/dr
1207 : ! force from overlap matrix contribution to H
1208 4426389 : DO j = 1, natorb_b
1209 3657204 : lb = laob(j)
1210 3657204 : nb = naob(j)
1211 28148232 : DO i = 1, natorb_a
1212 23721843 : la = laoa(i)
1213 23721843 : na = naoa(i)
1214 23721843 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1215 27379047 : IF (iatom <= jatom) THEN
1216 13989936 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1217 : ELSE
1218 9731907 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1219 : END IF
1220 : END DO
1221 : END DO
1222 1025580 : force_ab(ir) = foab
1223 : END DO
1224 256395 : DEALLOCATE (dfblock)
1225 : END IF
1226 : END IF
1227 :
1228 1158750 : IF (calculate_forces) THEN
1229 258799 : atom_a = atom_of_kind(iatom)
1230 258799 : atom_b = atom_of_kind(jatom)
1231 679414 : IF (irow == iatom) force_ab = -force_ab
1232 1035196 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1233 1035196 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1234 258799 : IF (use_virial) THEN
1235 165040 : f1 = 1.0_dp
1236 165040 : IF (iatom == jatom) f1 = 0.5_dp
1237 165040 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1238 : END IF
1239 : END IF
1240 :
1241 6956742 : DEALLOCATE (oint, owork, sint)
1242 :
1243 : END DO
1244 4242 : CALL neighbor_list_iterator_release(nl_iterator)
1245 :
1246 8484 : DO i = 1, SIZE(matrix_h, 1)
1247 64590 : DO img = 1, nimg
1248 56106 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1249 60348 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1250 : END DO
1251 : END DO
1252 :
1253 4242 : kf = xtb_control%kf
1254 4242 : enscale = xtb_control%enscale
1255 4242 : erep = 0.0_dp
1256 4242 : CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1257 :
1258 4242 : exb = 0.0_dp
1259 4242 : IF (xb_inter) THEN
1260 4168 : CALL xb_interaction(qs_env, exb, calculate_forces)
1261 : END IF
1262 :
1263 4242 : enonbonded = 0.0_dp
1264 4242 : IF (do_nonbonded) THEN
1265 : ! nonbonded interactions
1266 34 : NULLIFY (sab_xtb_nonbond)
1267 34 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
1268 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1269 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1270 : END IF
1271 :
1272 : ! set repulsive energy
1273 4242 : erep = erep + exb + enonbonded
1274 4242 : IF (xb_inter) THEN
1275 4168 : CALL para_env%sum(exb)
1276 4168 : energy%xtb_xb_inter = exb
1277 : END IF
1278 4242 : IF (do_nonbonded) THEN
1279 34 : CALL para_env%sum(enonbonded)
1280 34 : energy%xtb_nonbonded = enonbonded
1281 : END IF
1282 4242 : CALL para_env%sum(erep)
1283 4242 : energy%repulsive = erep
1284 :
1285 : ! deallocate coordination numbers
1286 4242 : CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1287 :
1288 : ! deallocate Huckel parameters
1289 4242 : DEALLOCATE (huckel)
1290 4242 : IF (calculate_forces) THEN
1291 564 : DEALLOCATE (dhuckel)
1292 : END IF
1293 : ! deallocate KAB parameters
1294 4242 : DEALLOCATE (kijab)
1295 :
1296 : ! AO matrix outputs
1297 4242 : CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1298 :
1299 4242 : DEALLOCATE (basis_set_list)
1300 4242 : IF (calculate_forces) THEN
1301 564 : IF (SIZE(matrix_p, 1) == 2) THEN
1302 752 : DO img = 1, nimg
1303 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
1304 710 : beta_scalar=-1.0_dp)
1305 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
1306 752 : beta_scalar=-1.0_dp)
1307 : END DO
1308 : END IF
1309 : END IF
1310 :
1311 4242 : CALL timestop(handle)
1312 :
1313 12726 : END SUBROUTINE build_gfn1_xtb_matrices
1314 :
1315 : ! **************************************************************************************************
1316 : !> \brief ...
1317 : !> \param qs_env ...
1318 : !> \param matrix_h ...
1319 : !> \param matrix_s ...
1320 : !> \param calculate_forces ...
1321 : ! **************************************************************************************************
1322 12576 : SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1323 : TYPE(qs_environment_type), POINTER :: qs_env
1324 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1325 : LOGICAL, INTENT(IN) :: calculate_forces
1326 :
1327 : INTEGER :: after, i, img, iw, nimg
1328 : LOGICAL :: norml1, norml2, omit_headers, use_arnoldi
1329 : REAL(KIND=dp), DIMENSION(2) :: condnum
1330 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1331 : TYPE(cp_logger_type), POINTER :: logger
1332 : TYPE(mp_para_env_type), POINTER :: para_env
1333 :
1334 6288 : logger => cp_get_default_logger()
1335 :
1336 6288 : CALL get_qs_env(qs_env, para_env=para_env)
1337 6288 : nimg = SIZE(matrix_h, 2)
1338 6288 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1339 6288 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1340 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
1341 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
1342 0 : extension=".Log")
1343 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1344 0 : after = MIN(MAX(after, 1), 16)
1345 0 : DO img = 1, nimg
1346 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
1347 0 : output_unit=iw, omit_headers=omit_headers)
1348 : END DO
1349 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
1350 : END IF
1351 :
1352 6288 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1353 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
1354 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
1355 0 : extension=".Log")
1356 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1357 0 : after = MIN(MAX(after, 1), 16)
1358 0 : DO img = 1, nimg
1359 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
1360 0 : output_unit=iw, omit_headers=omit_headers)
1361 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1362 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
1363 0 : DO i = 2, SIZE(matrix_s, 1)
1364 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
1365 0 : output_unit=iw, omit_headers=omit_headers)
1366 : END DO
1367 : END IF
1368 : END DO
1369 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
1370 : END IF
1371 :
1372 : ! *** Overlap condition number
1373 6288 : IF (.NOT. calculate_forces) THEN
1374 5656 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1375 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1376 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1377 4 : extension=".Log")
1378 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1379 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1380 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1381 4 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1382 4 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1383 : END IF
1384 : END IF
1385 :
1386 6288 : END SUBROUTINE ao_matrix_output
1387 :
1388 : END MODULE xtb_matrices
|
