Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of Overlap and Hamiltonian matrices in xTB
10 : !> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11 : !> JCTC 13, 1989-2009, (2017)
12 : !> DOI: 10.1021/acs.jctc.7b00118
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE xtb_matrices
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction
18 : USE ai_overlap, ONLY: overlap_ab
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind,&
21 : get_atomic_kind_set
22 : USE atprop_types, ONLY: atprop_array_init,&
23 : atprop_type
24 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
25 : gto_basis_set_type
26 : USE block_p_types, ONLY: block_p_type
27 : USE cp_blacs_env, ONLY: cp_blacs_env_type
28 : USE cp_control_types, ONLY: dft_control_type,&
29 : xtb_control_type
30 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_get_default_io_unit,&
36 : cp_logger_type,&
37 : cp_to_string
38 : USE cp_output_handling, ONLY: cp_p_file,&
39 : cp_print_key_finished_output,&
40 : cp_print_key_should_output,&
41 : cp_print_key_unit_nr
42 : USE dbcsr_api, ONLY: &
43 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_finalize, dbcsr_get_block_p, &
44 : dbcsr_multiply, dbcsr_p_type, dbcsr_type
45 : USE efield_tb_methods, ONLY: efield_tb_matrix
46 : USE ewald_environment_types, ONLY: ewald_env_get,&
47 : ewald_environment_type
48 : USE fparser, ONLY: evalfd,&
49 : finalizef
50 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
51 : section_vals_type,&
52 : section_vals_val_get
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE kpoint_types, ONLY: get_kpoint_info,&
56 : kpoint_type
57 : USE message_passing, ONLY: mp_para_env_type
58 : USE mulliken, ONLY: ao_charges
59 : USE orbital_pointers, ONLY: ncoset
60 : USE pair_potential, ONLY: init_genpot
61 : USE pair_potential_types, ONLY: not_initialized,&
62 : pair_potential_p_type,&
63 : pair_potential_pp_create,&
64 : pair_potential_pp_release,&
65 : pair_potential_pp_type,&
66 : pair_potential_single_clean,&
67 : pair_potential_single_copy,&
68 : pair_potential_single_type
69 : USE pair_potential_util, ONLY: ener_pot
70 : USE particle_types, ONLY: particle_type
71 : USE qs_charge_mixing, ONLY: charge_mixing
72 : USE qs_condnum, ONLY: overlap_condnum
73 : USE qs_dispersion_pairpot, ONLY: d3_cnumber,&
74 : dcnum_distribute,&
75 : dcnum_type
76 : USE qs_dispersion_types, ONLY: qs_dispersion_type
77 : USE qs_energy_types, ONLY: qs_energy_type
78 : USE qs_environment_types, ONLY: get_qs_env,&
79 : qs_environment_type
80 : USE qs_force_types, ONLY: qs_force_type
81 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
82 : get_memory_usage
83 : USE qs_kind_types, ONLY: get_qs_kind,&
84 : get_qs_kind_set,&
85 : qs_kind_type
86 : USE qs_ks_types, ONLY: get_ks_env,&
87 : qs_ks_env_type,&
88 : set_ks_env
89 : USE qs_mo_types, ONLY: get_mo_set,&
90 : mo_set_type
91 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
92 : neighbor_list_iterate,&
93 : neighbor_list_iterator_create,&
94 : neighbor_list_iterator_p_type,&
95 : neighbor_list_iterator_release,&
96 : neighbor_list_set_p_type
97 : USE qs_overlap, ONLY: create_sab_matrix
98 : USE qs_rho_types, ONLY: qs_rho_get,&
99 : qs_rho_type
100 : USE qs_scf_types, ONLY: qs_scf_env_type
101 : USE string_utilities, ONLY: compress,&
102 : uppercase
103 : USE virial_methods, ONLY: virial_pair_force
104 : USE virial_types, ONLY: virial_type
105 : USE xtb_coulomb, ONLY: build_xtb_coulomb
106 : USE xtb_parameters, ONLY: xtb_set_kab
107 : USE xtb_types, ONLY: get_xtb_atom_param,&
108 : xtb_atom_type
109 : #include "./base/base_uses.f90"
110 :
111 : IMPLICIT NONE
112 :
113 : TYPE neighbor_atoms_type
114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: coord
115 : REAL(KIND=dp), DIMENSION(:), POINTER :: rab
116 : INTEGER, DIMENSION(:), POINTER :: katom
117 : END TYPE neighbor_atoms_type
118 :
119 : PRIVATE
120 :
121 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_matrices'
122 :
123 : PUBLIC :: build_xtb_matrices, build_xtb_ks_matrix, xtb_hab_force
124 :
125 : CONTAINS
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param qs_env ...
130 : !> \param para_env ...
131 : !> \param calculate_forces ...
132 : ! **************************************************************************************************
133 3184 : SUBROUTINE build_xtb_matrices(qs_env, para_env, calculate_forces)
134 :
135 : TYPE(qs_environment_type), POINTER :: qs_env
136 : TYPE(mp_para_env_type), POINTER :: para_env
137 : LOGICAL, INTENT(IN) :: calculate_forces
138 :
139 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_xtb_matrices'
140 :
141 : INTEGER :: after, atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, &
142 : iset, iw, j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, maxs, &
143 : n1, n2, na, natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, &
144 : nsb, nseta, nsetb, nshell, sgfa, sgfb, za, zat, zb
145 6368 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
146 : INTEGER, DIMENSION(25) :: laoa, laob, lval, naoa, naob
147 : INTEGER, DIMENSION(3) :: cell
148 3184 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
149 3184 : npgfb, nsgfa, nsgfb
150 3184 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
151 3184 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
152 : LOGICAL :: defined, diagblock, do_nonbonded, floating_a, found, ghost_a, norml1, norml2, &
153 : omit_headers, use_arnoldi, use_virial, xb_inter
154 3184 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: floating, ghost
155 : REAL(KIND=dp) :: alphaa, alphab, derepij, dfp, dhij, dr, drk, drx, ena, enb, enonbonded, &
156 : erep, erepij, etaa, etab, exb, f0, f1, fen, fhua, fhub, fhud, foab, hij, k2sh, kab, kcnd, &
157 : kcnp, kcns, kd, ken, kf, kg, kia, kjb, kp, ks, ksp, kx2, kxr, rcova, rcovab, rcovb, rrab, &
158 : zneffa, zneffb
159 6368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers, kx
160 6368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, huckel, kcnlk, owork, rcab
161 3184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
162 3184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
163 : REAL(KIND=dp), DIMENSION(0:3) :: kl
164 : REAL(KIND=dp), DIMENSION(2) :: condnum
165 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, force_rr, &
166 : rij, rik
167 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kappaa, kappab, &
168 : kpolya, kpolyb, pia, pib
169 3184 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
170 3184 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171 3184 : scon_a, scon_b, wblock, zeta, zetb
172 3184 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173 : TYPE(atprop_type), POINTER :: atprop
174 12736 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
175 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
176 : TYPE(cp_logger_type), POINTER :: logger
177 3184 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
178 3184 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
179 : TYPE(dft_control_type), POINTER :: dft_control
180 3184 : 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(neighbor_atoms_type), ALLOCATABLE, &
184 3184 : DIMENSION(:) :: neighbor_atoms
185 : TYPE(neighbor_list_iterator_p_type), &
186 3184 : DIMENSION(:), POINTER :: nl_iterator
187 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
188 3184 : POINTER :: sab_orb, sab_xb, sab_xtb_nonbond
189 3184 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
190 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
191 : TYPE(qs_energy_type), POINTER :: energy
192 3184 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
193 3184 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
194 : TYPE(qs_ks_env_type), POINTER :: ks_env
195 : TYPE(qs_rho_type), POINTER :: rho
196 : TYPE(virial_type), POINTER :: virial
197 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
198 : TYPE(xtb_control_type), POINTER :: xtb_control
199 :
200 : !, na1
201 :
202 3184 : CALL timeset(routineN, handle)
203 :
204 3184 : NULLIFY (logger, virial, atprop)
205 3184 : logger => cp_get_default_logger()
206 :
207 3184 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
208 3184 : qs_kind_set, sab_orb, ks_env)
209 :
210 : CALL get_qs_env(qs_env=qs_env, &
211 : ks_env=ks_env, &
212 : energy=energy, &
213 : atomic_kind_set=atomic_kind_set, &
214 : qs_kind_set=qs_kind_set, &
215 : matrix_h_kp=matrix_h, &
216 : matrix_s_kp=matrix_s, &
217 : atprop=atprop, &
218 : dft_control=dft_control, &
219 3184 : sab_orb=sab_orb)
220 :
221 3184 : nkind = SIZE(atomic_kind_set)
222 3184 : xtb_control => dft_control%qs_control%xtb_control
223 3184 : xb_inter = xtb_control%xb_interaction
224 3184 : do_nonbonded = xtb_control%do_nonbonded
225 3184 : nimg = dft_control%nimages
226 3184 : nderivatives = 0
227 3184 : IF (calculate_forces) nderivatives = 1
228 3184 : IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
229 3184 : maxder = ncoset(nderivatives)
230 :
231 3184 : IF (xb_inter) energy%xtb_xb_inter = 0.0_dp
232 3184 : IF (do_nonbonded) energy%xtb_nonbonded = 0.0_dp
233 :
234 : ! global parameters (Table 2 from Ref.)
235 3184 : ks = xtb_control%ks
236 3184 : kp = xtb_control%kp
237 3184 : kd = xtb_control%kd
238 3184 : ksp = xtb_control%ksp
239 3184 : k2sh = xtb_control%k2sh
240 3184 : kg = xtb_control%kg
241 3184 : kf = xtb_control%kf
242 3184 : kcns = xtb_control%kcns
243 3184 : kcnp = xtb_control%kcnp
244 3184 : kcnd = xtb_control%kcnd
245 3184 : ken = xtb_control%ken
246 3184 : kxr = xtb_control%kxr
247 3184 : kx2 = xtb_control%kx2
248 :
249 3184 : NULLIFY (particle_set)
250 3184 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
251 3184 : natom = SIZE(particle_set)
252 3184 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
253 3184 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
254 :
255 3184 : IF (calculate_forces) THEN
256 428 : NULLIFY (rho, force, matrix_w)
257 : CALL get_qs_env(qs_env=qs_env, &
258 : rho=rho, matrix_w_kp=matrix_w, &
259 428 : virial=virial, force=force)
260 428 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
261 :
262 428 : IF (SIZE(matrix_p, 1) == 2) THEN
263 432 : DO img = 1, nimg
264 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
265 414 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
266 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
267 432 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
268 : END DO
269 : END IF
270 806 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
271 : END IF
272 : ! atomic energy decomposition
273 3184 : IF (atprop%energy) THEN
274 36 : CALL atprop_array_init(atprop%atecc, natom)
275 : END IF
276 :
277 3184 : NULLIFY (cell_to_index)
278 3184 : IF (nimg > 1) THEN
279 336 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
280 336 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
281 : END IF
282 :
283 : ! set up basis set lists
284 17380 : ALLOCATE (basis_set_list(nkind))
285 3184 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
286 :
287 : ! allocate overlap matrix
288 3184 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
289 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
290 3184 : sab_orb, .TRUE.)
291 3184 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
292 :
293 : ! initialize H matrix
294 3184 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
295 50004 : DO img = 1, nimg
296 46820 : ALLOCATE (matrix_h(1, img)%matrix)
297 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
298 46820 : name="HAMILTONIAN MATRIX")
299 50004 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
300 : END DO
301 3184 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
302 :
303 : ! Calculate coordination numbers
304 : ! needed for effective atomic energy levels (Eq. 12)
305 : ! code taken from D3 dispersion energy
306 9552 : ALLOCATE (cnumbers(natom))
307 31458 : cnumbers = 0._dp
308 3184 : IF (calculate_forces) THEN
309 1284 : ALLOCATE (dcnum(natom))
310 4218 : dcnum(:)%neighbors = 0
311 4218 : DO iatom = 1, natom
312 4218 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
313 : END DO
314 : ELSE
315 2756 : ALLOCATE (dcnum(1))
316 : END IF
317 :
318 22288 : ALLOCATE (ghost(nkind), floating(nkind), atomnumber(nkind))
319 11012 : DO ikind = 1, nkind
320 7828 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
321 7828 : CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_a, floating=floating_a)
322 7828 : ghost(ikind) = ghost_a
323 7828 : floating(ikind) = floating_a
324 18840 : atomnumber(ikind) = za
325 : END DO
326 3184 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
327 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
328 3184 : calculate_forces, .FALSE.)
329 3184 : DEALLOCATE (ghost, floating, atomnumber)
330 3184 : CALL para_env%sum(cnumbers)
331 3184 : IF (calculate_forces) THEN
332 428 : CALL dcnum_distribute(dcnum, para_env)
333 : END IF
334 :
335 : ! Calculate Huckel parameters
336 : ! Eq 12
337 : ! huckel(nshell,natom)
338 9552 : ALLOCATE (kcnlk(0:3, nkind))
339 11012 : DO ikind = 1, nkind
340 7828 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
341 11012 : IF (metal(za)) THEN
342 250 : kcnlk(0:3, ikind) = 0.0_dp
343 7778 : ELSEIF (early3d(za)) THEN
344 0 : kcnlk(0, ikind) = kcns
345 0 : kcnlk(1, ikind) = kcnp
346 0 : kcnlk(2, ikind) = 0.005_dp
347 0 : kcnlk(3, ikind) = 0.0_dp
348 : ELSE
349 7778 : kcnlk(0, ikind) = kcns
350 7778 : kcnlk(1, ikind) = kcnp
351 7778 : kcnlk(2, ikind) = kcnd
352 7778 : kcnlk(3, ikind) = 0.0_dp
353 : END IF
354 : END DO
355 9552 : ALLOCATE (huckel(5, natom))
356 3184 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
357 31458 : DO iatom = 1, natom
358 28274 : ikind = kind_of(iatom)
359 28274 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
360 28274 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
361 169644 : huckel(:, iatom) = 0.0_dp
362 116834 : DO i = 1, nshell
363 85376 : huckel(i, iatom) = hena(i)*(1._dp + kcnlk(lval(i), ikind)*cnumbers(iatom))
364 : END DO
365 : END DO
366 :
367 : ! Calculate KAB parameters and Huckel parameters and electronegativity correction
368 3184 : kl(0) = ks
369 3184 : kl(1) = kp
370 3184 : kl(2) = kd
371 3184 : kl(3) = 0.0_dp
372 19104 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
373 595064 : kijab = 0.0_dp
374 11012 : DO ikind = 1, nkind
375 7828 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
376 7828 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
377 7828 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
378 7826 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
379 40562 : DO jkind = 1, nkind
380 21724 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
381 21724 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
382 21724 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
383 21722 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
384 : ! get Fen = (1+ken*deltaEN^2)
385 21722 : fen = 1.0_dp + ken*(ena - enb)**2
386 : ! Kab
387 21722 : kab = xtb_set_kab(za, zb, xtb_control)
388 128010 : DO j = 1, natorb_b
389 76736 : lb = laob(j)
390 76736 : nb = naob(j)
391 390334 : DO i = 1, natorb_a
392 291874 : la = laoa(i)
393 291874 : na = naoa(i)
394 291874 : kia = kl(la)
395 291874 : kjb = kl(lb)
396 291874 : IF (zb == 1 .AND. nb == 2) kjb = k2sh
397 291874 : IF (za == 1 .AND. na == 2) kia = k2sh
398 368610 : IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
399 48328 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
400 : ELSE
401 243546 : IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
402 83676 : kijab(i, j, ikind, jkind) = ksp*kab*fen
403 : ELSE
404 159870 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)*kab*fen
405 : END IF
406 : END IF
407 : END DO
408 : END DO
409 : END DO
410 : END DO
411 :
412 3184 : IF (xb_inter) THEN
413 : ! list of neighbor atoms for XB term
414 9552 : ALLOCATE (neighbor_atoms(nkind))
415 11012 : DO ikind = 1, nkind
416 7828 : NULLIFY (neighbor_atoms(ikind)%coord)
417 7828 : NULLIFY (neighbor_atoms(ikind)%rab)
418 7828 : NULLIFY (neighbor_atoms(ikind)%katom)
419 7828 : CALL get_atomic_kind(atomic_kind_set(ikind), z=zat, natom=na)
420 11012 : IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
421 174 : ALLOCATE (neighbor_atoms(ikind)%coord(3, na))
422 290 : neighbor_atoms(ikind)%coord(1:3, 1:na) = 0.0_dp
423 174 : ALLOCATE (neighbor_atoms(ikind)%rab(na))
424 116 : neighbor_atoms(ikind)%rab(1:na) = HUGE(0.0_dp)
425 174 : ALLOCATE (neighbor_atoms(ikind)%katom(na))
426 116 : neighbor_atoms(ikind)%katom(1:na) = 0
427 : END IF
428 : END DO
429 : ! kx parameters
430 9552 : ALLOCATE (kx(nkind))
431 11012 : DO ikind = 1, nkind
432 7828 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
433 11012 : CALL get_xtb_atom_param(xtb_atom_a, kx=kx(ikind))
434 : END DO
435 : !
436 12736 : ALLOCATE (rcab(nkind, nkind))
437 11012 : DO ikind = 1, nkind
438 7828 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
439 7828 : CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova)
440 32740 : DO jkind = 1, nkind
441 21728 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
442 21728 : CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb)
443 29556 : rcab(ikind, jkind) = kxr*(rcova + rcovb)
444 : END DO
445 : END DO
446 : !
447 3184 : CALL get_qs_env(qs_env=qs_env, sab_xb=sab_xb)
448 : END IF
449 :
450 : ! initialize repulsion energy
451 3184 : erep = 0._dp
452 :
453 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
454 3184 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
455 943739 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
456 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
457 940555 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
458 940555 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
459 940555 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
460 940555 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
461 940555 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
462 940555 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
463 940555 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
464 :
465 3762220 : dr = SQRT(SUM(rij(:)**2))
466 :
467 : ! neighbor atom for XB term
468 940555 : IF (xb_inter .AND. (dr > 1.e-3_dp)) THEN
469 926113 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
470 250 : atom_a = atom_of_kind(iatom)
471 250 : IF (dr < neighbor_atoms(ikind)%rab(atom_a)) THEN
472 91 : neighbor_atoms(ikind)%rab(atom_a) = dr
473 364 : neighbor_atoms(ikind)%coord(1:3, atom_a) = rij(1:3)
474 91 : neighbor_atoms(ikind)%katom(atom_a) = jatom
475 : END IF
476 : END IF
477 926113 : IF (ASSOCIATED(neighbor_atoms(jkind)%rab)) THEN
478 287 : atom_b = atom_of_kind(jatom)
479 287 : IF (dr < neighbor_atoms(jkind)%rab(atom_b)) THEN
480 94 : neighbor_atoms(jkind)%rab(atom_b) = dr
481 376 : neighbor_atoms(jkind)%coord(1:3, atom_b) = -rij(1:3)
482 94 : neighbor_atoms(jkind)%katom(atom_b) = iatom
483 : END IF
484 : END IF
485 : END IF
486 :
487 : ! atomic parameters
488 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
489 : lmax=lmaxa, nshell=nsa, alpha=alphaa, zneff=zneffa, kpoly=kpolya, &
490 940555 : kappa=kappaa, hen=hena)
491 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
492 : lmax=lmaxb, nshell=nsb, alpha=alphab, zneff=zneffb, kpoly=kpolyb, &
493 940555 : kappa=kappab, hen=henb)
494 :
495 940555 : IF (nimg == 1) THEN
496 : ic = 1
497 : ELSE
498 241506 : ic = cell_to_index(cell(1), cell(2), cell(3))
499 241506 : CPASSERT(ic > 0)
500 : END IF
501 :
502 940555 : icol = MAX(iatom, jatom)
503 940555 : irow = MIN(iatom, jatom)
504 940555 : NULLIFY (sblock, fblock)
505 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
506 940555 : row=irow, col=icol, BLOCK=sblock, found=found)
507 940555 : CPASSERT(found)
508 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
509 940555 : row=irow, col=icol, BLOCK=fblock, found=found)
510 940555 : CPASSERT(found)
511 :
512 940555 : IF (calculate_forces) THEN
513 193498 : NULLIFY (pblock)
514 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
515 193498 : row=irow, col=icol, block=pblock, found=found)
516 193498 : CPASSERT(ASSOCIATED(pblock))
517 193498 : NULLIFY (wblock)
518 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
519 193498 : row=irow, col=icol, block=wblock, found=found)
520 193498 : CPASSERT(ASSOCIATED(wblock))
521 773992 : DO i = 2, 4
522 580494 : NULLIFY (dsblocks(i)%block)
523 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
524 580494 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
525 773992 : CPASSERT(found)
526 : END DO
527 : END IF
528 :
529 : ! overlap
530 940555 : basis_set_a => basis_set_list(ikind)%gto_basis_set
531 940555 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
532 940555 : basis_set_b => basis_set_list(jkind)%gto_basis_set
533 940555 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
534 940555 : atom_a = atom_of_kind(iatom)
535 940555 : atom_b = atom_of_kind(jatom)
536 : ! basis ikind
537 940555 : first_sgfa => basis_set_a%first_sgf
538 940555 : la_max => basis_set_a%lmax
539 940555 : la_min => basis_set_a%lmin
540 940555 : npgfa => basis_set_a%npgf
541 940555 : nseta = basis_set_a%nset
542 940555 : nsgfa => basis_set_a%nsgf_set
543 940555 : rpgfa => basis_set_a%pgf_radius
544 940555 : set_radius_a => basis_set_a%set_radius
545 940555 : scon_a => basis_set_a%scon
546 940555 : zeta => basis_set_a%zet
547 : ! basis jkind
548 940555 : first_sgfb => basis_set_b%first_sgf
549 940555 : lb_max => basis_set_b%lmax
550 940555 : lb_min => basis_set_b%lmin
551 940555 : npgfb => basis_set_b%npgf
552 940555 : nsetb = basis_set_b%nset
553 940555 : nsgfb => basis_set_b%nsgf_set
554 940555 : rpgfb => basis_set_b%pgf_radius
555 940555 : set_radius_b => basis_set_b%set_radius
556 940555 : scon_b => basis_set_b%scon
557 940555 : zetb => basis_set_b%zet
558 :
559 940555 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
560 7524440 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
561 4702775 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
562 33934923 : sint = 0.0_dp
563 :
564 2906698 : DO iset = 1, nseta
565 1966143 : ncoa = npgfa(iset)*ncoset(la_max(iset))
566 1966143 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
567 1966143 : sgfa = first_sgfa(1, iset)
568 7092079 : DO jset = 1, nsetb
569 4185381 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
570 3550035 : ncob = npgfb(jset)*ncoset(lb_max(jset))
571 3550035 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
572 3550035 : sgfb = first_sgfb(1, jset)
573 3550035 : IF (calculate_forces) THEN
574 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
575 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
576 736677 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
577 : ELSE
578 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
579 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
580 2813358 : rij, sab=oint(:, :, 1))
581 : END IF
582 : ! Contraction
583 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
584 3550035 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
585 3550035 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
586 5516178 : IF (calculate_forces) THEN
587 2946708 : DO i = 2, 4
588 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
589 2210031 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
590 2946708 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
591 : END DO
592 : END IF
593 : END DO
594 : END DO
595 : ! forces W matrix
596 940555 : IF (calculate_forces) THEN
597 773992 : DO i = 1, 3
598 773992 : IF (iatom <= jatom) THEN
599 8160408 : force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
600 : ELSE
601 6149439 : force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
602 : END IF
603 : END DO
604 193498 : f1 = 2.0_dp
605 773992 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
606 773992 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
607 193498 : IF (use_virial .AND. dr > 1.e-3_dp) THEN
608 99969 : IF (iatom == jatom) f1 = 1.0_dp
609 99969 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
610 99969 : IF (atprop%stress) THEN
611 64044 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_ab, rij)
612 64044 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*f1, force_ab, rij)
613 : END IF
614 : END IF
615 : END IF
616 : ! update S matrix
617 940555 : IF (iatom <= jatom) THEN
618 9328748 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
619 : ELSE
620 7331711 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
621 : END IF
622 940555 : IF (calculate_forces) THEN
623 773992 : DO i = 2, 4
624 773992 : IF (iatom <= jatom) THEN
625 8160408 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
626 : ELSE
627 6108969 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
628 : END IF
629 : END DO
630 : END IF
631 :
632 : ! Calculate Pi = Pia * Pib (Eq. 11)
633 940555 : rcovab = rcova + rcovb
634 940555 : rrab = SQRT(dr/rcovab)
635 2906698 : DO i = 1, nsa
636 2906698 : pia(i) = 1._dp + kpolya(i)*rrab
637 : END DO
638 2906703 : DO i = 1, nsb
639 2906703 : pib(i) = 1._dp + kpolyb(i)*rrab
640 : END DO
641 940555 : IF (calculate_forces) THEN
642 193498 : IF (dr > 1.e-6_dp) THEN
643 191450 : drx = 0.5_dp/rrab/rcovab
644 : ELSE
645 : drx = 0.0_dp
646 : END IF
647 616081 : dpia(1:nsa) = drx*kpolya(1:nsa)
648 616083 : dpib(1:nsb) = drx*kpolyb(1:nsb)
649 : END IF
650 :
651 : ! diagonal block
652 940555 : diagblock = .FALSE.
653 940555 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
654 : !
655 : ! Eq. 10
656 : !
657 : IF (diagblock) THEN
658 58521 : DO i = 1, natorb_a
659 44079 : na = naoa(i)
660 58521 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
661 : END DO
662 : ELSE
663 3890236 : DO j = 1, natorb_b
664 2964123 : nb = naob(j)
665 16540796 : DO i = 1, natorb_a
666 12650560 : na = naoa(i)
667 12650560 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
668 15614683 : IF (iatom <= jatom) THEN
669 7080170 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
670 : ELSE
671 5570390 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
672 : END IF
673 : END DO
674 : END DO
675 : END IF
676 940555 : IF (calculate_forces) THEN
677 193498 : f0 = 1.0_dp
678 193498 : IF (irow == iatom) f0 = -1.0_dp
679 : ! Derivative wrt coordination number
680 193498 : fhua = 0.0_dp
681 193498 : fhub = 0.0_dp
682 193498 : fhud = 0.0_dp
683 193498 : IF (diagblock) THEN
684 8280 : DO i = 1, natorb_a
685 6232 : la = laoa(i)
686 6232 : na = naoa(i)
687 8280 : fhud = fhud + pblock(i, i)*kcnlk(la, ikind)*hena(na)
688 : END DO
689 : ELSE
690 902491 : DO j = 1, natorb_b
691 711041 : lb = laob(j)
692 711041 : nb = naob(j)
693 4736811 : DO i = 1, natorb_a
694 3834320 : la = laoa(i)
695 3834320 : na = naoa(i)
696 3834320 : hij = 0.5_dp*pia(na)*pib(nb)
697 4545361 : IF (iatom <= jatom) THEN
698 2205265 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(la, ikind)*hena(na)
699 2205265 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(lb, jkind)*henb(nb)
700 : ELSE
701 1629055 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(la, ikind)*hena(na)
702 1629055 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(lb, jkind)*henb(nb)
703 : END IF
704 : END DO
705 : END DO
706 191450 : IF (iatom /= jatom) THEN
707 172110 : fhua = 2.0_dp*fhua
708 172110 : fhub = 2.0_dp*fhub
709 : END IF
710 : END IF
711 : ! iatom
712 193498 : atom_a = atom_of_kind(iatom)
713 10886952 : DO i = 1, dcnum(iatom)%neighbors
714 10693454 : katom = dcnum(iatom)%nlist(i)
715 10693454 : kkind = kind_of(katom)
716 10693454 : atom_c = atom_of_kind(katom)
717 42773816 : rik = dcnum(iatom)%rik(:, i)
718 42773816 : drk = SQRT(SUM(rik(:)**2))
719 10886952 : IF (drk > 1.e-3_dp) THEN
720 42773816 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
721 42773816 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
722 42773816 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
723 42773816 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
724 42773816 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
725 42773816 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
726 10693454 : IF (use_virial) THEN
727 20544948 : fdik = fdika + fdikb
728 5136237 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
729 5136237 : IF (atprop%stress) THEN
730 1477368 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
731 1477368 : CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fdik, rik)
732 : END IF
733 : END IF
734 : END IF
735 : END DO
736 : ! jatom
737 193498 : atom_b = atom_of_kind(jatom)
738 10881412 : DO i = 1, dcnum(jatom)%neighbors
739 10687914 : katom = dcnum(jatom)%nlist(i)
740 10687914 : kkind = kind_of(katom)
741 10687914 : atom_c = atom_of_kind(katom)
742 42751656 : rik = dcnum(jatom)%rik(:, i)
743 42751656 : drk = SQRT(SUM(rik(:)**2))
744 10881412 : IF (drk > 1.e-3_dp) THEN
745 42751656 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
746 42751656 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
747 42751656 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
748 10687914 : IF (use_virial) THEN
749 5134568 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
750 5134568 : IF (atprop%stress) THEN
751 1476450 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
752 1476450 : CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fdik, rik)
753 : END IF
754 : END IF
755 : END IF
756 : END DO
757 193498 : IF (diagblock) THEN
758 2048 : force_ab = 0._dp
759 : ELSE
760 : ! force from R dendent Huckel element
761 191450 : n1 = SIZE(fblock, 1)
762 191450 : n2 = SIZE(fblock, 2)
763 765800 : ALLOCATE (dfblock(n1, n2))
764 4723321 : dfblock = 0.0_dp
765 902491 : DO j = 1, natorb_b
766 711041 : lb = laob(j)
767 711041 : nb = naob(j)
768 4736811 : DO i = 1, natorb_a
769 3834320 : la = laoa(i)
770 3834320 : na = naoa(i)
771 3834320 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
772 4545361 : IF (iatom <= jatom) THEN
773 2205265 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
774 : ELSE
775 1629055 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
776 : END IF
777 : END DO
778 : END DO
779 4723321 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
780 765800 : DO ir = 1, 3
781 574350 : foab = 2.0_dp*dfp*rij(ir)/dr
782 : ! force from overlap matrix contribution to H
783 2707473 : DO j = 1, natorb_b
784 2133123 : lb = laob(j)
785 2133123 : nb = naob(j)
786 14210433 : DO i = 1, natorb_a
787 11502960 : la = laoa(i)
788 11502960 : na = naoa(i)
789 11502960 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
790 13636083 : IF (iatom <= jatom) THEN
791 6615795 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
792 : ELSE
793 4887165 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
794 : END IF
795 : END DO
796 : END DO
797 765800 : force_ab(ir) = foab
798 : END DO
799 191450 : DEALLOCATE (dfblock)
800 : END IF
801 : END IF
802 :
803 940555 : IF (calculate_forces) THEN
804 193498 : atom_a = atom_of_kind(iatom)
805 193498 : atom_b = atom_of_kind(jatom)
806 499303 : IF (irow == iatom) force_ab = -force_ab
807 773992 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
808 773992 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
809 193498 : IF (use_virial) THEN
810 100601 : f1 = 1.0_dp
811 100601 : IF (iatom == jatom) f1 = 0.5_dp
812 100601 : CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
813 100601 : IF (atprop%stress) THEN
814 64548 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_ab, rij)
815 64548 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*f1, force_ab, rij)
816 : END IF
817 : END IF
818 : END IF
819 :
820 940555 : DEALLOCATE (oint, owork, sint)
821 :
822 : ! repulsive potential
823 5646514 : IF (dr > 0.001_dp) THEN
824 926113 : erepij = zneffa*zneffb/dr*EXP(-SQRT(alphaa*alphab)*dr**kf)
825 926113 : erep = erep + erepij
826 926113 : IF (atprop%energy) THEN
827 128088 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
828 128088 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
829 : END IF
830 926113 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
831 191450 : derepij = -(1.0_dp/dr + SQRT(alphaa*alphab)*kf*dr**(kf - 1.0_dp))*erepij
832 191450 : force_rr(1) = derepij*rij(1)/dr
833 191450 : force_rr(2) = derepij*rij(2)/dr
834 191450 : force_rr(3) = derepij*rij(3)/dr
835 191450 : atom_a = atom_of_kind(iatom)
836 191450 : atom_b = atom_of_kind(jatom)
837 765800 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
838 765800 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
839 191450 : IF (use_virial) THEN
840 99969 : f1 = 1.0_dp
841 99969 : IF (iatom == jatom) f1 = 0.5_dp
842 99969 : CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
843 99969 : IF (atprop%stress) THEN
844 64044 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_rr, rij)
845 64044 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*f1, force_rr, rij)
846 : END IF
847 : END IF
848 : END IF
849 : END IF
850 :
851 : END DO
852 3184 : CALL neighbor_list_iterator_release(nl_iterator)
853 :
854 6368 : DO i = 1, SIZE(matrix_h, 1)
855 53188 : DO img = 1, nimg
856 46820 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
857 50004 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
858 : END DO
859 : END DO
860 :
861 3184 : exb = 0.0_dp
862 3184 : IF (xb_inter) THEN
863 3184 : CALL xb_neighbors(neighbor_atoms, para_env)
864 : CALL xb_interaction(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
865 3184 : calculate_forces, use_virial, force, virial, atprop)
866 : END IF
867 :
868 3184 : enonbonded = 0.0_dp
869 3184 : IF (do_nonbonded) THEN
870 : ! nonbonded interactions
871 34 : NULLIFY (sab_xtb_nonbond)
872 34 : CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
873 : CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
874 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
875 : END IF
876 : ! set repulsive energy
877 3184 : erep = erep + exb + enonbonded
878 3184 : IF (xb_inter) THEN
879 3184 : CALL para_env%sum(exb)
880 3184 : energy%xtb_xb_inter = exb
881 : END IF
882 3184 : IF (do_nonbonded) THEN
883 34 : CALL para_env%sum(enonbonded)
884 34 : energy%xtb_nonbonded = enonbonded
885 : END IF
886 3184 : CALL para_env%sum(erep)
887 3184 : energy%repulsive = erep
888 :
889 : ! deallocate coordination numbers
890 3184 : DEALLOCATE (cnumbers)
891 3184 : IF (calculate_forces) THEN
892 4218 : DO iatom = 1, natom
893 4218 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
894 : END DO
895 : END IF
896 3184 : DEALLOCATE (dcnum)
897 :
898 : ! deallocate Huckel parameters
899 3184 : DEALLOCATE (huckel)
900 : ! deallocate KAB parameters
901 3184 : DEALLOCATE (kijab)
902 :
903 3184 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
904 3184 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
905 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
906 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
907 0 : extension=".Log")
908 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
909 0 : after = MIN(MAX(after, 1), 16)
910 0 : DO img = 1, nimg
911 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
912 0 : output_unit=iw, omit_headers=omit_headers)
913 : END DO
914 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
915 : END IF
916 :
917 3184 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
918 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
919 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
920 0 : extension=".Log")
921 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
922 0 : after = MIN(MAX(after, 1), 16)
923 0 : DO img = 1, nimg
924 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
925 0 : output_unit=iw, omit_headers=omit_headers)
926 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
927 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
928 0 : DO i = 2, SIZE(matrix_s, 1)
929 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
930 0 : output_unit=iw, omit_headers=omit_headers)
931 : END DO
932 : END IF
933 : END DO
934 0 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
935 : END IF
936 :
937 : ! *** Overlap condition number
938 3184 : IF (.NOT. calculate_forces) THEN
939 2756 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
940 : "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
941 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
942 4 : extension=".Log")
943 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
944 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
945 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
946 4 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
947 4 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
948 : END IF
949 : END IF
950 :
951 3184 : DEALLOCATE (basis_set_list)
952 3184 : IF (calculate_forces) THEN
953 428 : IF (SIZE(matrix_p, 1) == 2) THEN
954 432 : DO img = 1, nimg
955 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
956 432 : beta_scalar=-1.0_dp)
957 : END DO
958 : END IF
959 : END IF
960 :
961 3184 : IF (xb_inter) THEN
962 11012 : DO ikind = 1, nkind
963 7828 : IF (ASSOCIATED(neighbor_atoms(ikind)%coord)) THEN
964 58 : DEALLOCATE (neighbor_atoms(ikind)%coord)
965 : END IF
966 7828 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
967 58 : DEALLOCATE (neighbor_atoms(ikind)%rab)
968 : END IF
969 11012 : IF (ASSOCIATED(neighbor_atoms(ikind)%katom)) THEN
970 58 : DEALLOCATE (neighbor_atoms(ikind)%katom)
971 : END IF
972 : END DO
973 3184 : DEALLOCATE (neighbor_atoms)
974 3184 : DEALLOCATE (kx, rcab)
975 : END IF
976 :
977 3184 : CALL timestop(handle)
978 :
979 15920 : END SUBROUTINE build_xtb_matrices
980 :
981 : ! **************************************************************************************************
982 : !> \brief ...
983 : !> \param qs_env ...
984 : !> \param p_matrix ...
985 : ! **************************************************************************************************
986 16 : SUBROUTINE xtb_hab_force(qs_env, p_matrix)
987 :
988 : TYPE(qs_environment_type), POINTER :: qs_env
989 : TYPE(dbcsr_type), POINTER :: p_matrix
990 :
991 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xtb_hab_force'
992 :
993 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
994 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, maxs, n1, n2, &
995 : na, natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, &
996 : nseta, nsetb, nshell, sgfa, sgfb, za, zb
997 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
998 : INTEGER, DIMENSION(25) :: laoa, laob, lval, naoa, naob
999 : INTEGER, DIMENSION(3) :: cell
1000 16 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1001 16 : npgfb, nsgfa, nsgfb
1002 16 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1003 : LOGICAL :: defined, diagblock, floating_a, found, &
1004 : ghost_a, use_virial
1005 16 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: floating, ghost
1006 : REAL(KIND=dp) :: alphaa, alphab, dfp, dhij, dr, drk, drx, ena, enb, etaa, etab, f0, fen, &
1007 : fhua, fhub, fhud, foab, hij, k2sh, kab, kcnd, kcnp, kcns, kd, ken, kf, kg, kia, kjb, kp, &
1008 : ks, ksp, kx2, kxr, rcova, rcovab, rcovb, rrab, zneffa, zneffb
1009 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
1010 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, huckel, kcnlk, owork
1011 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
1012 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
1013 : REAL(KIND=dp), DIMENSION(0:3) :: kl
1014 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
1015 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kappaa, kappab, &
1016 : kpolya, kpolyb, pia, pib
1017 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1018 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
1019 16 : scon_a, scon_b, zeta, zetb
1020 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1021 64 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
1022 : TYPE(cp_logger_type), POINTER :: logger
1023 16 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1024 16 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
1025 : TYPE(dft_control_type), POINTER :: dft_control
1026 16 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1027 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1028 : TYPE(mp_para_env_type), POINTER :: para_env
1029 : TYPE(neighbor_list_iterator_p_type), &
1030 16 : DIMENSION(:), POINTER :: nl_iterator
1031 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1032 16 : POINTER :: sab_orb
1033 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1034 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1035 16 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1036 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1037 : TYPE(qs_ks_env_type), POINTER :: ks_env
1038 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
1039 : TYPE(xtb_control_type), POINTER :: xtb_control
1040 :
1041 16 : CALL timeset(routineN, handle)
1042 :
1043 16 : NULLIFY (logger)
1044 16 : logger => cp_get_default_logger()
1045 :
1046 16 : NULLIFY (matrix_h, matrix_s, atomic_kind_set, qs_kind_set, sab_orb)
1047 :
1048 : CALL get_qs_env(qs_env=qs_env, &
1049 : atomic_kind_set=atomic_kind_set, &
1050 : qs_kind_set=qs_kind_set, &
1051 : dft_control=dft_control, &
1052 : para_env=para_env, &
1053 16 : sab_orb=sab_orb)
1054 :
1055 16 : nkind = SIZE(atomic_kind_set)
1056 16 : xtb_control => dft_control%qs_control%xtb_control
1057 16 : nimg = dft_control%nimages
1058 16 : nderivatives = 1
1059 16 : maxder = ncoset(nderivatives)
1060 :
1061 : ! global parameters (Table 2 from Ref.)
1062 16 : ks = xtb_control%ks
1063 16 : kp = xtb_control%kp
1064 16 : kd = xtb_control%kd
1065 16 : ksp = xtb_control%ksp
1066 16 : k2sh = xtb_control%k2sh
1067 16 : kg = xtb_control%kg
1068 16 : kf = xtb_control%kf
1069 16 : kcns = xtb_control%kcns
1070 16 : kcnp = xtb_control%kcnp
1071 16 : kcnd = xtb_control%kcnd
1072 16 : ken = xtb_control%ken
1073 16 : kxr = xtb_control%kxr
1074 16 : kx2 = xtb_control%kx2
1075 :
1076 16 : NULLIFY (particle_set)
1077 16 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
1078 16 : natom = SIZE(particle_set)
1079 16 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1080 16 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
1081 :
1082 16 : NULLIFY (force)
1083 16 : CALL get_qs_env(qs_env=qs_env, force=force)
1084 16 : use_virial = .FALSE.
1085 16 : CPASSERT(nimg == 1)
1086 :
1087 : ! set up basis set lists
1088 86 : ALLOCATE (basis_set_list(nkind))
1089 16 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1090 :
1091 : ! allocate overlap matrix
1092 16 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
1093 16 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
1094 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
1095 16 : sab_orb, .TRUE.)
1096 : ! initialize H matrix
1097 16 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
1098 32 : DO img = 1, nimg
1099 16 : ALLOCATE (matrix_h(1, img)%matrix)
1100 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
1101 16 : name="HAMILTONIAN MATRIX")
1102 32 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
1103 : END DO
1104 :
1105 : ! Calculate coordination numbers
1106 : ! needed for effective atomic energy levels (Eq. 12)
1107 : ! code taken from D3 dispersion energy
1108 48 : ALLOCATE (cnumbers(natom))
1109 256 : cnumbers = 0._dp
1110 48 : ALLOCATE (dcnum(natom))
1111 256 : dcnum(:)%neighbors = 0
1112 256 : DO iatom = 1, natom
1113 256 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
1114 : END DO
1115 :
1116 112 : ALLOCATE (ghost(nkind), floating(nkind), atomnumber(nkind))
1117 54 : DO ikind = 1, nkind
1118 38 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1119 38 : CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_a, floating=floating_a)
1120 38 : ghost(ikind) = ghost_a
1121 38 : floating(ikind) = floating_a
1122 92 : atomnumber(ikind) = za
1123 : END DO
1124 16 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
1125 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
1126 16 : .TRUE., .FALSE.)
1127 16 : DEALLOCATE (ghost, floating, atomnumber)
1128 16 : CALL para_env%sum(cnumbers)
1129 16 : CALL dcnum_distribute(dcnum, para_env)
1130 :
1131 : ! Calculate Huckel parameters
1132 : ! Eq 12
1133 : ! huckel(nshell,natom)
1134 48 : ALLOCATE (kcnlk(0:3, nkind))
1135 54 : DO ikind = 1, nkind
1136 38 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1137 54 : IF (metal(za)) THEN
1138 0 : kcnlk(0:3, ikind) = 0.0_dp
1139 38 : ELSEIF (early3d(za)) THEN
1140 0 : kcnlk(0, ikind) = kcns
1141 0 : kcnlk(1, ikind) = kcnp
1142 0 : kcnlk(2, ikind) = 0.005_dp
1143 0 : kcnlk(3, ikind) = 0.0_dp
1144 : ELSE
1145 38 : kcnlk(0, ikind) = kcns
1146 38 : kcnlk(1, ikind) = kcnp
1147 38 : kcnlk(2, ikind) = kcnd
1148 38 : kcnlk(3, ikind) = 0.0_dp
1149 : END IF
1150 : END DO
1151 48 : ALLOCATE (huckel(5, natom))
1152 16 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1153 256 : DO iatom = 1, natom
1154 240 : ikind = kind_of(iatom)
1155 240 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1156 240 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
1157 1440 : huckel(:, iatom) = 0.0_dp
1158 976 : DO i = 1, nshell
1159 720 : huckel(i, iatom) = hena(i)*(1._dp + kcnlk(lval(i), ikind)*cnumbers(iatom))
1160 : END DO
1161 : END DO
1162 :
1163 : ! Calculate KAB parameters and Huckel parameters and electronegativity correction
1164 16 : kl(0) = ks
1165 16 : kl(1) = kp
1166 16 : kl(2) = kd
1167 16 : kl(3) = 0.0_dp
1168 96 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
1169 2028 : kijab = 0.0_dp
1170 54 : DO ikind = 1, nkind
1171 38 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1172 38 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
1173 38 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
1174 38 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
1175 186 : DO jkind = 1, nkind
1176 94 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1177 94 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
1178 94 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
1179 94 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
1180 : ! get Fen = (1+ken*deltaEN^2)
1181 94 : fen = 1.0_dp + ken*(ena - enb)**2
1182 : ! Kab
1183 94 : kab = xtb_set_kab(za, zb, xtb_control)
1184 526 : DO j = 1, natorb_b
1185 300 : lb = laob(j)
1186 300 : nb = naob(j)
1187 1354 : DO i = 1, natorb_a
1188 960 : la = laoa(i)
1189 960 : na = naoa(i)
1190 960 : kia = kl(la)
1191 960 : kjb = kl(lb)
1192 960 : IF (zb == 1 .AND. nb == 2) kjb = k2sh
1193 960 : IF (za == 1 .AND. na == 2) kia = k2sh
1194 1260 : IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
1195 224 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
1196 : ELSE
1197 736 : IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
1198 336 : kijab(i, j, ikind, jkind) = ksp*kab*fen
1199 : ELSE
1200 400 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)*kab*fen
1201 : END IF
1202 : END IF
1203 : END DO
1204 : END DO
1205 : END DO
1206 : END DO
1207 :
1208 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
1209 16 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1210 13003 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1211 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1212 12987 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
1213 12987 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1214 12987 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
1215 12987 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
1216 12987 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1217 12987 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
1218 12987 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
1219 :
1220 51948 : dr = SQRT(SUM(rij(:)**2))
1221 :
1222 : ! atomic parameters
1223 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
1224 : lmax=lmaxa, nshell=nsa, alpha=alphaa, zneff=zneffa, kpoly=kpolya, &
1225 12987 : kappa=kappaa, hen=hena)
1226 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
1227 : lmax=lmaxb, nshell=nsb, alpha=alphab, zneff=zneffb, kpoly=kpolyb, &
1228 12987 : kappa=kappab, hen=henb)
1229 :
1230 12987 : ic = 1
1231 :
1232 12987 : icol = MAX(iatom, jatom)
1233 12987 : irow = MIN(iatom, jatom)
1234 12987 : NULLIFY (sblock, fblock)
1235 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1236 12987 : row=irow, col=icol, BLOCK=sblock, found=found)
1237 12987 : CPASSERT(found)
1238 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
1239 12987 : row=irow, col=icol, BLOCK=fblock, found=found)
1240 12987 : CPASSERT(found)
1241 :
1242 12987 : NULLIFY (pblock)
1243 : CALL dbcsr_get_block_p(matrix=p_matrix, &
1244 12987 : row=irow, col=icol, block=pblock, found=found)
1245 12987 : CPASSERT(ASSOCIATED(pblock))
1246 51948 : DO i = 2, 4
1247 38961 : NULLIFY (dsblocks(i)%block)
1248 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
1249 38961 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
1250 51948 : CPASSERT(found)
1251 : END DO
1252 :
1253 : ! overlap
1254 12987 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1255 12987 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1256 12987 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1257 12987 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1258 12987 : atom_a = atom_of_kind(iatom)
1259 12987 : atom_b = atom_of_kind(jatom)
1260 : ! basis ikind
1261 12987 : first_sgfa => basis_set_a%first_sgf
1262 12987 : la_max => basis_set_a%lmax
1263 12987 : la_min => basis_set_a%lmin
1264 12987 : npgfa => basis_set_a%npgf
1265 12987 : nseta = basis_set_a%nset
1266 12987 : nsgfa => basis_set_a%nsgf_set
1267 12987 : rpgfa => basis_set_a%pgf_radius
1268 12987 : set_radius_a => basis_set_a%set_radius
1269 12987 : scon_a => basis_set_a%scon
1270 12987 : zeta => basis_set_a%zet
1271 : ! basis jkind
1272 12987 : first_sgfb => basis_set_b%first_sgf
1273 12987 : lb_max => basis_set_b%lmax
1274 12987 : lb_min => basis_set_b%lmin
1275 12987 : npgfb => basis_set_b%npgf
1276 12987 : nsetb = basis_set_b%nset
1277 12987 : nsgfb => basis_set_b%nsgf_set
1278 12987 : rpgfb => basis_set_b%pgf_radius
1279 12987 : set_radius_b => basis_set_b%set_radius
1280 12987 : scon_b => basis_set_b%scon
1281 12987 : zetb => basis_set_b%zet
1282 :
1283 12987 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1284 103896 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1285 64935 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
1286 515591 : sint = 0.0_dp
1287 :
1288 38961 : DO iset = 1, nseta
1289 25974 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1290 25974 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1291 25974 : sgfa = first_sgfa(1, iset)
1292 90909 : DO jset = 1, nsetb
1293 51948 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1294 45815 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1295 45815 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1296 45815 : sgfb = first_sgfb(1, jset)
1297 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1298 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1299 45815 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1300 : ! Contraction
1301 255049 : DO i = 1, 4
1302 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1303 183260 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1304 235208 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1305 : END DO
1306 : END DO
1307 : END DO
1308 : ! update S matrix
1309 12987 : IF (iatom <= jatom) THEN
1310 61008 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1311 : ELSE
1312 59865 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1313 : END IF
1314 51948 : DO i = 2, 4
1315 51948 : IF (iatom <= jatom) THEN
1316 183024 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1317 : ELSE
1318 179595 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
1319 : END IF
1320 : END DO
1321 :
1322 : ! Calculate Pi = Pia * Pib (Eq. 11)
1323 12987 : rcovab = rcova + rcovb
1324 12987 : rrab = SQRT(dr/rcovab)
1325 38961 : DO i = 1, nsa
1326 38961 : pia(i) = 1._dp + kpolya(i)*rrab
1327 : END DO
1328 38961 : DO i = 1, nsb
1329 38961 : pib(i) = 1._dp + kpolyb(i)*rrab
1330 : END DO
1331 12987 : IF (dr > 1.e-6_dp) THEN
1332 12867 : drx = 0.5_dp/rrab/rcovab
1333 : ELSE
1334 : drx = 0.0_dp
1335 : END IF
1336 38961 : dpia(1:nsa) = drx*kpolya(1:nsa)
1337 38961 : dpib(1:nsb) = drx*kpolyb(1:nsb)
1338 :
1339 : ! diagonal block
1340 12987 : diagblock = .FALSE.
1341 12987 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
1342 : !
1343 : ! Eq. 10
1344 : !
1345 : IF (diagblock) THEN
1346 444 : DO i = 1, natorb_a
1347 324 : na = naoa(i)
1348 444 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1349 : END DO
1350 : ELSE
1351 44819 : DO j = 1, natorb_b
1352 31952 : nb = naob(j)
1353 124223 : DO i = 1, natorb_a
1354 79404 : na = naoa(i)
1355 79404 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1356 111356 : IF (iatom <= jatom) THEN
1357 39648 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1358 : ELSE
1359 39756 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1360 : END IF
1361 : END DO
1362 : END DO
1363 : END IF
1364 :
1365 12987 : f0 = 1.0_dp
1366 12987 : IF (irow == iatom) f0 = -1.0_dp
1367 : ! Derivative wrt coordination number
1368 12987 : fhua = 0.0_dp
1369 12987 : fhub = 0.0_dp
1370 12987 : fhud = 0.0_dp
1371 12987 : IF (diagblock) THEN
1372 444 : DO i = 1, natorb_a
1373 324 : la = laoa(i)
1374 324 : na = naoa(i)
1375 444 : fhud = fhud + pblock(i, i)*kcnlk(la, ikind)*hena(na)
1376 : END DO
1377 : ELSE
1378 44819 : DO j = 1, natorb_b
1379 31952 : lb = laob(j)
1380 31952 : nb = naob(j)
1381 124223 : DO i = 1, natorb_a
1382 79404 : la = laoa(i)
1383 79404 : na = naoa(i)
1384 79404 : hij = 0.5_dp*pia(na)*pib(nb)
1385 111356 : IF (iatom <= jatom) THEN
1386 39648 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(la, ikind)*hena(na)
1387 39648 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(lb, jkind)*henb(nb)
1388 : ELSE
1389 39756 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(la, ikind)*hena(na)
1390 39756 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(lb, jkind)*henb(nb)
1391 : END IF
1392 : END DO
1393 : END DO
1394 12867 : IF (iatom /= jatom) THEN
1395 12825 : fhua = 2.0_dp*fhua
1396 12825 : fhub = 2.0_dp*fhub
1397 : END IF
1398 : END IF
1399 : ! iatom
1400 12987 : atom_a = atom_of_kind(iatom)
1401 305732 : DO i = 1, dcnum(iatom)%neighbors
1402 292745 : katom = dcnum(iatom)%nlist(i)
1403 292745 : kkind = kind_of(katom)
1404 292745 : atom_c = atom_of_kind(katom)
1405 1170980 : rik = dcnum(iatom)%rik(:, i)
1406 1170980 : drk = SQRT(SUM(rik(:)**2))
1407 305732 : IF (drk > 1.e-3_dp) THEN
1408 1170980 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1409 1170980 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1410 1170980 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1411 1170980 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1412 1170980 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1413 1170980 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1414 : END IF
1415 : END DO
1416 : ! jatom
1417 12987 : atom_b = atom_of_kind(jatom)
1418 304187 : DO i = 1, dcnum(jatom)%neighbors
1419 291200 : katom = dcnum(jatom)%nlist(i)
1420 291200 : kkind = kind_of(katom)
1421 291200 : atom_c = atom_of_kind(katom)
1422 1164800 : rik = dcnum(jatom)%rik(:, i)
1423 1164800 : drk = SQRT(SUM(rik(:)**2))
1424 304187 : IF (drk > 1.e-3_dp) THEN
1425 1164800 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1426 1164800 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1427 1164800 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1428 : END IF
1429 : END DO
1430 12987 : IF (diagblock) THEN
1431 120 : force_ab = 0._dp
1432 : ELSE
1433 : ! force from R dendent Huckel element
1434 12867 : n1 = SIZE(fblock, 1)
1435 12867 : n2 = SIZE(fblock, 2)
1436 51468 : ALLOCATE (dfblock(n1, n2))
1437 119445 : dfblock = 0.0_dp
1438 44819 : DO j = 1, natorb_b
1439 31952 : lb = laob(j)
1440 31952 : nb = naob(j)
1441 124223 : DO i = 1, natorb_a
1442 79404 : la = laoa(i)
1443 79404 : na = naoa(i)
1444 79404 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1445 111356 : IF (iatom <= jatom) THEN
1446 39648 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1447 : ELSE
1448 39756 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1449 : END IF
1450 : END DO
1451 : END DO
1452 119445 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
1453 51468 : DO ir = 1, 3
1454 38601 : foab = 2.0_dp*dfp*rij(ir)/dr
1455 : ! force from overlap matrix contribution to H
1456 134457 : DO j = 1, natorb_b
1457 95856 : lb = laob(j)
1458 95856 : nb = naob(j)
1459 372669 : DO i = 1, natorb_a
1460 238212 : la = laoa(i)
1461 238212 : na = naoa(i)
1462 238212 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1463 334068 : IF (iatom <= jatom) THEN
1464 118944 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1465 : ELSE
1466 119268 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1467 : END IF
1468 : END DO
1469 : END DO
1470 51468 : force_ab(ir) = foab
1471 : END DO
1472 12867 : DEALLOCATE (dfblock)
1473 : END IF
1474 :
1475 12987 : atom_a = atom_of_kind(iatom)
1476 12987 : atom_b = atom_of_kind(jatom)
1477 32565 : IF (irow == iatom) force_ab = -force_ab
1478 51948 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1479 51948 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1480 :
1481 77938 : DEALLOCATE (oint, owork, sint)
1482 :
1483 : END DO
1484 16 : CALL neighbor_list_iterator_release(nl_iterator)
1485 :
1486 32 : DO i = 1, SIZE(matrix_h, 1)
1487 48 : DO img = 1, nimg
1488 16 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1489 32 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1490 : END DO
1491 : END DO
1492 16 : CALL dbcsr_deallocate_matrix_set(matrix_s)
1493 16 : CALL dbcsr_deallocate_matrix_set(matrix_h)
1494 :
1495 : ! deallocate coordination numbers
1496 16 : DEALLOCATE (cnumbers)
1497 256 : DO iatom = 1, natom
1498 256 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
1499 : END DO
1500 16 : DEALLOCATE (dcnum)
1501 :
1502 : ! deallocate Huckel parameters
1503 16 : DEALLOCATE (huckel)
1504 : ! deallocate KAB parameters
1505 16 : DEALLOCATE (kijab)
1506 :
1507 16 : DEALLOCATE (basis_set_list)
1508 :
1509 16 : CALL timestop(handle)
1510 :
1511 64 : END SUBROUTINE xtb_hab_force
1512 :
1513 : ! **************************************************************************************************
1514 : !> \brief ...
1515 : !> \param qs_env ...
1516 : !> \param calculate_forces ...
1517 : !> \param just_energy ...
1518 : !> \param ext_ks_matrix ...
1519 : ! **************************************************************************************************
1520 23020 : SUBROUTINE build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
1521 : TYPE(qs_environment_type), POINTER :: qs_env
1522 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
1523 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
1524 : POINTER :: ext_ks_matrix
1525 :
1526 : CHARACTER(len=*), PARAMETER :: routineN = 'build_xtb_ks_matrix'
1527 :
1528 : INTEGER :: atom_a, handle, iatom, ikind, img, &
1529 : iounit, is, ispin, na, natom, natorb, &
1530 : nimg, nkind, ns, nsgf, nspins
1531 : INTEGER, DIMENSION(25) :: lao
1532 : INTEGER, DIMENSION(5) :: occ
1533 : LOGICAL :: do_efield, pass_check
1534 : REAL(KIND=dp) :: achrg, chmax, pc_ener, qmmm_el
1535 23020 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
1536 23020 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, charges
1537 23020 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1538 : TYPE(cp_logger_type), POINTER :: logger
1539 23020 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs, p_matrix
1540 23020 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
1541 : TYPE(dbcsr_type), POINTER :: mo_coeff, s_matrix
1542 : TYPE(dft_control_type), POINTER :: dft_control
1543 : TYPE(mp_para_env_type), POINTER :: para_env
1544 23020 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1545 : TYPE(qs_energy_type), POINTER :: energy
1546 23020 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1547 : TYPE(qs_ks_env_type), POINTER :: ks_env
1548 : TYPE(qs_rho_type), POINTER :: rho
1549 : TYPE(qs_scf_env_type), POINTER :: scf_env
1550 : TYPE(section_vals_type), POINTER :: scf_section
1551 : TYPE(xtb_atom_type), POINTER :: xtb_kind
1552 :
1553 23020 : CALL timeset(routineN, handle)
1554 :
1555 23020 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
1556 23020 : ks_matrix, rho, energy)
1557 23020 : CPASSERT(ASSOCIATED(qs_env))
1558 :
1559 23020 : logger => cp_get_default_logger()
1560 23020 : iounit = cp_logger_get_default_io_unit(logger)
1561 :
1562 : CALL get_qs_env(qs_env, &
1563 : dft_control=dft_control, &
1564 : atomic_kind_set=atomic_kind_set, &
1565 : qs_kind_set=qs_kind_set, &
1566 : matrix_h_kp=matrix_h, &
1567 : para_env=para_env, &
1568 : ks_env=ks_env, &
1569 : matrix_ks_kp=ks_matrix, &
1570 : rho=rho, &
1571 23020 : energy=energy)
1572 :
1573 23020 : IF (PRESENT(ext_ks_matrix)) THEN
1574 : ! remap pointer to allow for non-kpoint external ks matrix
1575 : ! ext_ks_matrix is used in linear response code
1576 16 : ns = SIZE(ext_ks_matrix)
1577 16 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
1578 : END IF
1579 :
1580 23020 : energy%hartree = 0.0_dp
1581 23020 : energy%qmmm_el = 0.0_dp
1582 23020 : energy%efield = 0.0_dp
1583 :
1584 23020 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
1585 23020 : nspins = dft_control%nspins
1586 23020 : nimg = dft_control%nimages
1587 23020 : CPASSERT(ASSOCIATED(matrix_h))
1588 23020 : CPASSERT(ASSOCIATED(rho))
1589 23020 : CPASSERT(SIZE(ks_matrix) > 0)
1590 :
1591 48192 : DO ispin = 1, nspins
1592 434174 : DO img = 1, nimg
1593 : ! copy the core matrix into the fock matrix
1594 411154 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
1595 : END DO
1596 : END DO
1597 :
1598 23020 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
1599 : dft_control%apply_efield_field) THEN
1600 : do_efield = .TRUE.
1601 : ELSE
1602 17128 : do_efield = .FALSE.
1603 : END IF
1604 :
1605 23020 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
1606 : ! Mulliken charges
1607 21138 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
1608 21138 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1609 21138 : natom = SIZE(particle_set)
1610 105690 : ALLOCATE (mcharge(natom), charges(natom, 5))
1611 1146428 : charges = 0.0_dp
1612 21138 : nkind = SIZE(atomic_kind_set)
1613 21138 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1614 84552 : ALLOCATE (aocg(nsgf, natom))
1615 1091628 : aocg = 0.0_dp
1616 21138 : IF (nimg > 1) THEN
1617 2644 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
1618 : ELSE
1619 18494 : p_matrix => matrix_p(:, 1)
1620 18494 : s_matrix => matrix_s(1, 1)%matrix
1621 18494 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1622 : END IF
1623 76960 : DO ikind = 1, nkind
1624 55822 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1625 55822 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1626 55822 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
1627 336702 : DO iatom = 1, na
1628 203920 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1629 1223520 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
1630 909932 : DO is = 1, natorb
1631 650190 : ns = lao(is) + 1
1632 854110 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
1633 : END DO
1634 : END DO
1635 : END DO
1636 21138 : DEALLOCATE (aocg)
1637 : ! charge mixing
1638 21138 : IF (dft_control%qs_control%do_ls_scf) THEN
1639 : !
1640 : ELSE
1641 20930 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
1642 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1643 20930 : charges, para_env, scf_env%iter_count)
1644 : END IF
1645 :
1646 246196 : DO iatom = 1, natom
1647 1246540 : mcharge(iatom) = SUM(charges(iatom, :))
1648 : END DO
1649 : END IF
1650 :
1651 23020 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1652 : CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
1653 21138 : calculate_forces, just_energy)
1654 : END IF
1655 :
1656 23020 : IF (do_efield) THEN
1657 5892 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
1658 : END IF
1659 :
1660 23020 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1661 21138 : IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
1662 20280 : pass_check = .TRUE.
1663 73404 : DO ikind = 1, nkind
1664 53124 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1665 53124 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1666 53124 : CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
1667 326124 : DO iatom = 1, na
1668 199596 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1669 199596 : achrg = mcharge(atom_a)
1670 252720 : IF (ABS(achrg) > chmax) THEN
1671 0 : IF (iounit > 0) THEN
1672 0 : WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
1673 0 : " Kind Atom=", ikind, atom_a, " Limit=", chmax, " Charge=", achrg
1674 : END IF
1675 : pass_check = .FALSE.
1676 : END IF
1677 : END DO
1678 : END DO
1679 20280 : IF (.NOT. pass_check) THEN
1680 : CALL cp_warn(__LOCATION__, "Atomic charges outside chemical range were detected."// &
1681 : " Switch-off CHECK_ATOMIC_CHARGES keyword in the &xTB section"// &
1682 0 : " if you want to force to continue the calculation.")
1683 0 : CPABORT("xTB Charges")
1684 : END IF
1685 : END IF
1686 : END IF
1687 :
1688 23020 : IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
1689 21138 : DEALLOCATE (mcharge, charges)
1690 : END IF
1691 :
1692 23020 : IF (qs_env%qmmm) THEN
1693 5146 : CPASSERT(SIZE(ks_matrix, 2) == 1)
1694 10292 : DO ispin = 1, nspins
1695 : ! If QM/MM sumup the 1el Hamiltonian
1696 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
1697 5146 : 1.0_dp, 1.0_dp)
1698 5146 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
1699 : ! Compute QM/MM Energy
1700 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
1701 5146 : matrix_p1(ispin)%matrix, qmmm_el)
1702 10292 : energy%qmmm_el = energy%qmmm_el + qmmm_el
1703 : END DO
1704 5146 : pc_ener = qs_env%ks_qmmm_env%pc_ener
1705 5146 : energy%qmmm_el = energy%qmmm_el + pc_ener
1706 : END IF
1707 :
1708 : energy%total = energy%core + energy%hartree + energy%efield + energy%qmmm_el + &
1709 23020 : energy%repulsive + energy%dispersion + energy%dftb3
1710 :
1711 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
1712 23020 : extension=".scfLog")
1713 23020 : IF (iounit > 0) THEN
1714 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
1715 0 : "Repulsive pair potential energy: ", energy%repulsive, &
1716 0 : "Zeroth order Hamiltonian energy: ", energy%core, &
1717 0 : "Charge fluctuation energy: ", energy%hartree, &
1718 0 : "London dispersion energy: ", energy%dispersion
1719 0 : IF (dft_control%qs_control%xtb_control%xb_interaction) &
1720 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
1721 0 : "Correction for halogen bonding: ", energy%xtb_xb_inter
1722 0 : IF (dft_control%qs_control%xtb_control%do_nonbonded) &
1723 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
1724 0 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
1725 0 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
1726 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
1727 0 : "Electric field interaction energy: ", energy%efield
1728 : END IF
1729 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
1730 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
1731 0 : IF (qs_env%qmmm) THEN
1732 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
1733 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
1734 : END IF
1735 : END IF
1736 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
1737 23020 : "PRINT%DETAILED_ENERGY")
1738 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
1739 23020 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
1740 7234 : CPASSERT(SIZE(ks_matrix, 2) == 1)
1741 : BLOCK
1742 7234 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1743 7234 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
1744 14516 : DO ispin = 1, SIZE(mo_derivs)
1745 7282 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
1746 7282 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
1747 0 : CPABORT("")
1748 : END IF
1749 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
1750 14516 : 0.0_dp, mo_derivs(ispin)%matrix)
1751 : END DO
1752 : END BLOCK
1753 : END IF
1754 :
1755 23020 : CALL timestop(handle)
1756 :
1757 46040 : END SUBROUTINE build_xtb_ks_matrix
1758 :
1759 : ! **************************************************************************************************
1760 : !> \brief Distributes the neighbor atom information to all processors
1761 : !>
1762 : !> \param neighbor_atoms ...
1763 : !> \param para_env ...
1764 : !> \par History
1765 : !> 1.2019 JGH
1766 : !> \version 1.1
1767 : ! **************************************************************************************************
1768 3184 : SUBROUTINE xb_neighbors(neighbor_atoms, para_env)
1769 : TYPE(neighbor_atoms_type), DIMENSION(:), &
1770 : INTENT(INOUT) :: neighbor_atoms
1771 : TYPE(mp_para_env_type), POINTER :: para_env
1772 :
1773 : INTEGER :: iatom, ikind, natom, nkind
1774 3184 : INTEGER, ALLOCATABLE, DIMENSION(:) :: matom
1775 3184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dmloc
1776 3184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
1777 :
1778 3184 : nkind = SIZE(neighbor_atoms)
1779 11012 : DO ikind = 1, nkind
1780 11012 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
1781 58 : natom = SIZE(neighbor_atoms(ikind)%rab)
1782 406 : ALLOCATE (dmloc(2*natom), matom(natom), coord(3, natom))
1783 174 : dmloc = 0.0_dp
1784 116 : DO iatom = 1, natom
1785 58 : dmloc(2*iatom - 1) = neighbor_atoms(ikind)%rab(iatom)
1786 116 : dmloc(2*iatom) = REAL(para_env%mepos, KIND=dp)
1787 : END DO
1788 58 : CALL para_env%minloc(dmloc)
1789 290 : coord = 0.0_dp
1790 116 : matom = 0
1791 116 : DO iatom = 1, natom
1792 58 : neighbor_atoms(ikind)%rab(iatom) = dmloc(2*iatom - 1)
1793 116 : IF (NINT(dmloc(2*iatom)) == para_env%mepos) THEN
1794 116 : coord(1:3, iatom) = neighbor_atoms(ikind)%coord(1:3, iatom)
1795 29 : matom(iatom) = neighbor_atoms(ikind)%katom(iatom)
1796 : END IF
1797 : END DO
1798 58 : CALL para_env%sum(coord)
1799 290 : neighbor_atoms(ikind)%coord(1:3, :) = coord(1:3, :)
1800 58 : CALL para_env%sum(matom)
1801 116 : neighbor_atoms(ikind)%katom(:) = matom(:)
1802 58 : DEALLOCATE (dmloc, matom, coord)
1803 : END IF
1804 : END DO
1805 :
1806 3184 : END SUBROUTINE xb_neighbors
1807 : ! **************************************************************************************************
1808 : !> \brief Computes a correction for nonbonded interactions based on a generic potential
1809 : !>
1810 : !> \param enonbonded energy contribution
1811 : !> \param force ...
1812 : !> \param qs_env ...
1813 : !> \param xtb_control ...
1814 : !> \param sab_xtb_nonbond ...
1815 : !> \param atomic_kind_set ...
1816 : !> \param calculate_forces ...
1817 : !> \param use_virial ...
1818 : !> \param virial ...
1819 : !> \param atprop ...
1820 : !> \param atom_of_kind ..
1821 : !> \par History
1822 : !> 12.2018 JGH
1823 : !> \version 1.1
1824 : ! **************************************************************************************************
1825 68 : SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1826 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1827 : REAL(dp), INTENT(INOUT) :: enonbonded
1828 : TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
1829 : POINTER :: force
1830 : TYPE(qs_environment_type), POINTER :: qs_env
1831 : TYPE(xtb_control_type), POINTER :: xtb_control
1832 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1833 : INTENT(IN), POINTER :: sab_xtb_nonbond
1834 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
1835 : POINTER :: atomic_kind_set
1836 : LOGICAL, INTENT(IN) :: calculate_forces, use_virial
1837 : TYPE(virial_type), INTENT(IN), POINTER :: virial
1838 : TYPE(atprop_type), INTENT(IN), POINTER :: atprop
1839 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind
1840 :
1841 : CHARACTER(len=*), PARAMETER :: routineN = 'nonbonded_correction'
1842 :
1843 : CHARACTER(LEN=default_string_length) :: def_error, this_error
1844 : INTEGER :: atom_i, atom_j, handle, iatom, ikind, &
1845 : jatom, jkind, kk, ntype
1846 : INTEGER, DIMENSION(3) :: cell
1847 : LOGICAL :: do_ewald
1848 : REAL(KIND=dp) :: dedf, dr, dx, energy_cutoff, err, fval, &
1849 : lerr, rcut
1850 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
1851 : TYPE(ewald_environment_type), POINTER :: ewald_env
1852 : TYPE(neighbor_list_iterator_p_type), &
1853 34 : DIMENSION(:), POINTER :: nl_iterator
1854 : TYPE(pair_potential_p_type), POINTER :: nonbonded
1855 : TYPE(pair_potential_pp_type), POINTER :: potparm
1856 : TYPE(pair_potential_single_type), POINTER :: pot
1857 : TYPE(section_vals_type), POINTER :: nonbonded_section
1858 :
1859 34 : CALL timeset(routineN, handle)
1860 :
1861 : NULLIFY (nonbonded)
1862 34 : NULLIFY (potparm)
1863 34 : NULLIFY (ewald_env)
1864 34 : nonbonded => xtb_control%nonbonded
1865 34 : do_ewald = xtb_control%do_ewald
1866 34 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
1867 :
1868 34 : ntype = SIZE(atomic_kind_set)
1869 34 : CALL pair_potential_pp_create(potparm, ntype)
1870 : !Assign input and potential info to potparm_nonbond
1871 34 : CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
1872 : !Initialize genetic potential
1873 34 : CALL init_genpot(potparm, ntype)
1874 :
1875 34 : NULLIFY (pot)
1876 34 : enonbonded = 0._dp
1877 34 : energy_cutoff = 0._dp
1878 :
1879 34 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond)
1880 227 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1881 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1882 193 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
1883 193 : pot => potparm%pot(ikind, jkind)%pot
1884 193 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
1885 193 : rcut = SQRT(pot%rcutsq)
1886 193 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
1887 25 : fval = 1.0_dp
1888 25 : IF (ikind == jkind) fval = 0.5_dp
1889 : ! splines not implemented
1890 25 : enonbonded = enonbonded + fval*ener_pot(pot, dr, energy_cutoff)
1891 25 : IF (atprop%energy) THEN
1892 16 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
1893 16 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
1894 : END IF
1895 : END IF
1896 :
1897 193 : IF (calculate_forces) THEN
1898 :
1899 93 : kk = SIZE(pot%type)
1900 93 : IF (kk /= 1) THEN
1901 0 : CALL cp_warn(__LOCATION__, "Generic potential with type > 1 not implemented.")
1902 0 : CPABORT("pot type")
1903 : END IF
1904 : ! rmin and rmax and rcut
1905 93 : IF ((pot%set(kk)%rmin /= not_initialized) .AND. (dr < pot%set(kk)%rmin)) CYCLE
1906 : ! An upper boundary for the potential definition was defined
1907 93 : IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) CYCLE
1908 : ! If within limits let's compute the potential...
1909 93 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
1910 :
1911 9 : NULLIFY (nonbonded_section)
1912 9 : nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED")
1913 9 : CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx)
1914 9 : CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr)
1915 :
1916 9 : dedf = fval*evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err)
1917 9 : IF (ABS(err) > lerr) THEN
1918 1 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1919 1 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1920 1 : CALL compress(this_error, .TRUE.)
1921 1 : CALL compress(def_error, .TRUE.)
1922 : CALL cp_warn(__LOCATION__, &
1923 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
1924 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
1925 1 : TRIM(def_error)//' .')
1926 : END IF
1927 :
1928 9 : atom_i = atom_of_kind(iatom)
1929 9 : atom_j = atom_of_kind(jatom)
1930 :
1931 36 : fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
1932 :
1933 36 : force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
1934 36 : force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:)
1935 :
1936 9 : IF (use_virial) THEN
1937 8 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
1938 8 : IF (atprop%stress) THEN
1939 8 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
1940 8 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij)
1941 : END IF
1942 : END IF
1943 :
1944 : END IF
1945 : END IF
1946 193 : NULLIFY (pot)
1947 : END DO
1948 34 : CALL neighbor_list_iterator_release(nl_iterator)
1949 34 : CALL finalizef()
1950 :
1951 34 : IF (ASSOCIATED(potparm)) THEN
1952 34 : CALL pair_potential_pp_release(potparm)
1953 : END IF
1954 :
1955 34 : CALL timestop(handle)
1956 :
1957 34 : END SUBROUTINE nonbonded_correction
1958 : ! **************************************************************************************************
1959 : !> \brief ...
1960 : !> \param atomic_kind_set ...
1961 : !> \param nonbonded ...
1962 : !> \param potparm ...
1963 : !> \param ewald_env ...
1964 : !> \param do_ewald ...
1965 : ! **************************************************************************************************
1966 34 : SUBROUTINE force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
1967 :
1968 : ! routine based on force_field_pack_nonbond
1969 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
1970 : POINTER :: atomic_kind_set
1971 : TYPE(pair_potential_p_type), INTENT(IN), POINTER :: nonbonded
1972 : TYPE(pair_potential_pp_type), INTENT(INOUT), &
1973 : POINTER :: potparm
1974 : TYPE(ewald_environment_type), INTENT(IN), POINTER :: ewald_env
1975 : LOGICAL, INTENT(IN) :: do_ewald
1976 :
1977 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
1978 : name_atm_b, name_atm_b_local
1979 : INTEGER :: ikind, ingp, iw, jkind
1980 : LOGICAL :: found
1981 : REAL(KIND=dp) :: ewald_rcut
1982 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1983 : TYPE(cp_logger_type), POINTER :: logger
1984 : TYPE(pair_potential_single_type), POINTER :: pot
1985 :
1986 34 : NULLIFY (pot, logger)
1987 :
1988 34 : logger => cp_get_default_logger()
1989 34 : iw = cp_logger_get_default_io_unit(logger)
1990 :
1991 170 : DO ikind = 1, SIZE(atomic_kind_set)
1992 136 : atomic_kind => atomic_kind_set(ikind)
1993 136 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
1994 510 : DO jkind = ikind, SIZE(atomic_kind_set)
1995 340 : atomic_kind => atomic_kind_set(jkind)
1996 340 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
1997 340 : found = .FALSE.
1998 340 : name_atm_a = name_atm_a_local
1999 340 : name_atm_b = name_atm_b_local
2000 340 : CALL uppercase(name_atm_a)
2001 340 : CALL uppercase(name_atm_b)
2002 340 : pot => potparm%pot(ikind, jkind)%pot
2003 340 : IF (ASSOCIATED(nonbonded)) THEN
2004 680 : DO ingp = 1, SIZE(nonbonded%pot)
2005 340 : IF ((TRIM(nonbonded%pot(ingp)%pot%at1) == "*") .OR. &
2006 : (TRIM(nonbonded%pot(ingp)%pot%at2) == "*")) CYCLE
2007 :
2008 : !IF (iw > 0) WRITE (iw, *) "TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), &
2009 : ! " with ", TRIM(nonbonded%pot(ingp)%pot%at1), &
2010 : ! TRIM(nonbonded%pot(ingp)%pot%at2)
2011 : IF ((((name_atm_a) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
2012 340 : ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. &
2013 : (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
2014 340 : ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN
2015 34 : CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot)
2016 : ! multiple potential not implemented, simply overwriting
2017 34 : IF (found) &
2018 : CALL cp_warn(__LOCATION__, &
2019 : "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// &
2020 0 : " and "//TRIM(name_atm_b)//" OVERWRITING! ")
2021 : !IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
2022 : found = .TRUE.
2023 : END IF
2024 : END DO
2025 : END IF
2026 476 : IF (.NOT. found) THEN
2027 306 : CALL pair_potential_single_clean(pot)
2028 : !IF (iw > 0) WRITE (iw, *) " NOTFOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
2029 : END IF
2030 : END DO !jkind
2031 : END DO !ikind
2032 :
2033 : ! Cutoff is defined always as the maximum between the FF and Ewald
2034 34 : IF (do_ewald) THEN
2035 16 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
2036 16 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
2037 : !IF (iw > 0) WRITE (iw, *) " RCUT ", SQRT(pot%rcutsq), ewald_rcut
2038 : END IF
2039 :
2040 34 : END SUBROUTINE force_field_pack_nonbond_pot_correction
2041 : ! **************************************************************************************************
2042 : !> \brief Computes the interaction term between Br/I/At and donor atoms
2043 : !>
2044 : !> \param exb ...
2045 : !> \param neighbor_atoms ...
2046 : !> \param atom_of_kind ...
2047 : !> \param kind_of ...
2048 : !> \param sab_xb ...
2049 : !> \param kx ...
2050 : !> \param kx2 ...
2051 : !> \param rcab ...
2052 : !> \param calculate_forces ...
2053 : !> \param use_virial ...
2054 : !> \param force ...
2055 : !> \param virial ...
2056 : !> \param atprop ...
2057 : !> \par History
2058 : !> 12.2018 JGH
2059 : !> \version 1.1
2060 : ! **************************************************************************************************
2061 3184 : SUBROUTINE xb_interaction(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
2062 : calculate_forces, use_virial, force, virial, atprop)
2063 : REAL(dp), INTENT(INOUT) :: exb
2064 : TYPE(neighbor_atoms_type), DIMENSION(:), &
2065 : INTENT(IN) :: neighbor_atoms
2066 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
2067 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2068 : POINTER :: sab_xb
2069 : REAL(dp), DIMENSION(:), INTENT(IN) :: kx
2070 : REAL(dp), INTENT(IN) :: kx2
2071 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rcab
2072 : LOGICAL, INTENT(IN) :: calculate_forces, use_virial
2073 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2074 : TYPE(virial_type), POINTER :: virial
2075 : TYPE(atprop_type), POINTER :: atprop
2076 :
2077 : INTEGER :: atom_a, atom_b, atom_c, iatom, ikind, &
2078 : jatom, jkind, katom, kkind
2079 : INTEGER, DIMENSION(3) :: cell
2080 : REAL(KIND=dp) :: alp, aterm, cosa, daterm, ddab, ddax, &
2081 : ddbx, ddr, ddr12, ddr6, deval, dr, &
2082 : drab, drax, drbx, eval, xy
2083 : REAL(KIND=dp), DIMENSION(3) :: fia, fij, fja, ria, rij, rja
2084 : TYPE(neighbor_list_iterator_p_type), &
2085 3184 : DIMENSION(:), POINTER :: nl_iterator
2086 :
2087 : ! exonent in angular term
2088 3184 : alp = 6.0_dp
2089 : ! loop over all atom pairs
2090 3184 : CALL neighbor_list_iterator_create(nl_iterator, sab_xb)
2091 3196 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2092 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2093 12 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
2094 : ! ikind, iatom : Halogen
2095 : ! jkind, jatom : Donor
2096 12 : atom_a = atom_of_kind(iatom)
2097 12 : katom = neighbor_atoms(ikind)%katom(atom_a)
2098 12 : IF (katom == 0) CYCLE
2099 12 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
2100 12 : ddr = rcab(ikind, jkind)/dr
2101 12 : ddr6 = ddr**6
2102 12 : ddr12 = ddr6*ddr6
2103 12 : eval = kx(ikind)*(ddr12 - kx2*ddr6)/(1.0_dp + ddr12)
2104 : ! angle
2105 48 : ria(1:3) = neighbor_atoms(ikind)%coord(1:3, atom_a)
2106 48 : rja(1:3) = rij(1:3) - ria(1:3)
2107 12 : drax = ria(1)**2 + ria(2)**2 + ria(3)**2
2108 12 : drbx = dr*dr
2109 12 : drab = rja(1)**2 + rja(2)**2 + rja(3)**2
2110 12 : xy = SQRT(drbx*drax)
2111 : ! cos angle B-X-A
2112 12 : cosa = (drbx + drax - drab)/xy
2113 12 : aterm = (0.5_dp - 0.25_dp*cosa)**alp
2114 : !
2115 12 : exb = exb + aterm*eval
2116 12 : IF (atprop%energy) THEN
2117 0 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*aterm*eval
2118 0 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*aterm*eval
2119 : END IF
2120 : !
2121 3196 : IF (calculate_forces) THEN
2122 6 : kkind = kind_of(katom)
2123 6 : atom_b = atom_of_kind(jatom)
2124 6 : atom_c = atom_of_kind(katom)
2125 : !
2126 6 : deval = 6.0_dp*kx(ikind)*ddr6*(kx2*ddr12 + 2.0_dp*ddr6 - kx2)/(1.0_dp + ddr12)**2
2127 6 : deval = -rcab(ikind, jkind)*deval/(dr*dr)/ddr
2128 24 : fij(1:3) = aterm*deval*rij(1:3)/dr
2129 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:)
2130 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:)
2131 6 : IF (use_virial) THEN
2132 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
2133 0 : IF (atprop%stress) THEN
2134 0 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
2135 0 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij)
2136 : END IF
2137 : END IF
2138 : !
2139 6 : fij(1:3) = 0.0_dp
2140 6 : fia(1:3) = 0.0_dp
2141 6 : fja(1:3) = 0.0_dp
2142 6 : daterm = -0.25_dp*alp*(0.5_dp - 0.25_dp*cosa)**(alp - 1.0_dp)
2143 6 : ddbx = 0.5_dp*(drab - drax + drbx)/xy/drbx
2144 6 : ddax = 0.5_dp*(drab + drax - drbx)/xy/drax
2145 6 : ddab = -1._dp/xy
2146 24 : fij(1:3) = 2.0_dp*daterm*ddbx*rij(1:3)*eval
2147 24 : fia(1:3) = 2.0_dp*daterm*ddax*ria(1:3)*eval
2148 24 : fja(1:3) = 2.0_dp*daterm*ddab*rja(1:3)*eval
2149 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:) - fia(:)
2150 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:) + fja(:)
2151 24 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fia(:) - fja(:)
2152 6 : IF (use_virial) THEN
2153 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
2154 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fia, ria)
2155 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fja, rja)
2156 0 : IF (atprop%stress) THEN
2157 0 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
2158 0 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij)
2159 0 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fia, ria)
2160 0 : CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fia, ria)
2161 0 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fja, rja)
2162 0 : CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fja, rja)
2163 : END IF
2164 : END IF
2165 : END IF
2166 : END DO
2167 3184 : CALL neighbor_list_iterator_release(nl_iterator)
2168 :
2169 3184 : END SUBROUTINE xb_interaction
2170 :
2171 : ! **************************************************************************************************
2172 : !> \brief ...
2173 : !> \param z ...
2174 : !> \return ...
2175 : ! **************************************************************************************************
2176 7866 : FUNCTION metal(z) RESULT(ismetal)
2177 : INTEGER :: z
2178 : LOGICAL :: ismetal
2179 :
2180 7866 : SELECT CASE (z)
2181 : CASE DEFAULT
2182 7816 : ismetal = .TRUE.
2183 : CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
2184 7866 : ismetal = .FALSE.
2185 : END SELECT
2186 :
2187 7866 : END FUNCTION metal
2188 :
2189 : ! **************************************************************************************************
2190 : !> \brief ...
2191 : !> \param z ...
2192 : !> \return ...
2193 : ! **************************************************************************************************
2194 7816 : FUNCTION early3d(z) RESULT(isearly3d)
2195 : INTEGER :: z
2196 : LOGICAL :: isearly3d
2197 :
2198 7816 : isearly3d = .FALSE.
2199 7816 : IF (z >= 21 .AND. z <= 24) isearly3d = .TRUE.
2200 :
2201 7816 : END FUNCTION early3d
2202 :
2203 0 : END MODULE xtb_matrices
2204 :
|
