(git:1f285aa)
xtb_matrices.F
Go to the documentation of this file.
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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
146  INTEGER, DIMENSION(25) :: laoa, laob, lval, naoa, naob
147  INTEGER, DIMENSION(3) :: cell
148  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
149  npgfb, nsgfa, nsgfb
150  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
151  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers, kx
160  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, huckel, kcnlk, owork, rcab
161  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
162  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  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
170  REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171  scon_a, scon_b, wblock, zeta, zetb
172  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173  TYPE(atprop_type), POINTER :: atprop
174  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
178  TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
179  TYPE(dft_control_type), POINTER :: dft_control
180  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  DIMENSION(:) :: neighbor_atoms
185  TYPE(neighbor_list_iterator_p_type), &
186  DIMENSION(:), POINTER :: nl_iterator
187  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
188  POINTER :: sab_orb, sab_xb, sab_xtb_nonbond
189  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
190  TYPE(qs_dispersion_type), POINTER :: dispersion_env
191  TYPE(qs_energy_type), POINTER :: energy
192  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
193  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  CALL timeset(routinen, handle)
203 
204  NULLIFY (logger, virial, atprop)
205  logger => cp_get_default_logger()
206 
207  NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
208  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  sab_orb=sab_orb)
220 
221  nkind = SIZE(atomic_kind_set)
222  xtb_control => dft_control%qs_control%xtb_control
223  xb_inter = xtb_control%xb_interaction
224  do_nonbonded = xtb_control%do_nonbonded
225  nimg = dft_control%nimages
226  nderivatives = 0
227  IF (calculate_forces) nderivatives = 1
228  IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
229  maxder = ncoset(nderivatives)
230 
231  IF (xb_inter) energy%xtb_xb_inter = 0.0_dp
232  IF (do_nonbonded) energy%xtb_nonbonded = 0.0_dp
233 
234  ! global parameters (Table 2 from Ref.)
235  ks = xtb_control%ks
236  kp = xtb_control%kp
237  kd = xtb_control%kd
238  ksp = xtb_control%ksp
239  k2sh = xtb_control%k2sh
240  kg = xtb_control%kg
241  kf = xtb_control%kf
242  kcns = xtb_control%kcns
243  kcnp = xtb_control%kcnp
244  kcnd = xtb_control%kcnd
245  ken = xtb_control%ken
246  kxr = xtb_control%kxr
247  kx2 = xtb_control%kx2
248 
249  NULLIFY (particle_set)
250  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
251  natom = SIZE(particle_set)
252  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
253  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
254 
255  IF (calculate_forces) THEN
256  NULLIFY (rho, force, matrix_w)
257  CALL get_qs_env(qs_env=qs_env, &
258  rho=rho, matrix_w_kp=matrix_w, &
259  virial=virial, force=force)
260  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
261 
262  IF (SIZE(matrix_p, 1) == 2) THEN
263  DO img = 1, nimg
264  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
265  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  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
268  END DO
269  END IF
270  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
271  END IF
272  ! atomic energy decomposition
273  IF (atprop%energy) THEN
274  CALL atprop_array_init(atprop%atecc, natom)
275  END IF
276 
277  NULLIFY (cell_to_index)
278  IF (nimg > 1) THEN
279  CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
280  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
281  END IF
282 
283  ! set up basis set lists
284  ALLOCATE (basis_set_list(nkind))
285  CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
286 
287  ! allocate overlap matrix
288  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  sab_orb, .true.)
291  CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
292 
293  ! initialize H matrix
294  CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
295  DO img = 1, nimg
296  ALLOCATE (matrix_h(1, img)%matrix)
297  CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
298  name="HAMILTONIAN MATRIX")
299  CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
300  END DO
301  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  ALLOCATE (cnumbers(natom))
307  cnumbers = 0._dp
308  IF (calculate_forces) THEN
309  ALLOCATE (dcnum(natom))
310  dcnum(:)%neighbors = 0
311  DO iatom = 1, natom
312  ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
313  END DO
314  ELSE
315  ALLOCATE (dcnum(1))
316  END IF
317 
318  ALLOCATE (ghost(nkind), floating(nkind), atomnumber(nkind))
319  DO ikind = 1, nkind
320  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
321  CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_a, floating=floating_a)
322  ghost(ikind) = ghost_a
323  floating(ikind) = floating_a
324  atomnumber(ikind) = za
325  END DO
326  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  calculate_forces, .false.)
329  DEALLOCATE (ghost, floating, atomnumber)
330  CALL para_env%sum(cnumbers)
331  IF (calculate_forces) THEN
332  CALL dcnum_distribute(dcnum, para_env)
333  END IF
334 
335  ! Calculate Huckel parameters
336  ! Eq 12
337  ! huckel(nshell,natom)
338  ALLOCATE (kcnlk(0:3, nkind))
339  DO ikind = 1, nkind
340  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
341  IF (metal(za)) THEN
342  kcnlk(0:3, ikind) = 0.0_dp
343  ELSEIF (early3d(za)) THEN
344  kcnlk(0, ikind) = kcns
345  kcnlk(1, ikind) = kcnp
346  kcnlk(2, ikind) = 0.005_dp
347  kcnlk(3, ikind) = 0.0_dp
348  ELSE
349  kcnlk(0, ikind) = kcns
350  kcnlk(1, ikind) = kcnp
351  kcnlk(2, ikind) = kcnd
352  kcnlk(3, ikind) = 0.0_dp
353  END IF
354  END DO
355  ALLOCATE (huckel(5, natom))
356  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
357  DO iatom = 1, natom
358  ikind = kind_of(iatom)
359  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
360  CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
361  huckel(:, iatom) = 0.0_dp
362  DO i = 1, nshell
363  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  kl(0) = ks
369  kl(1) = kp
370  kl(2) = kd
371  kl(3) = 0.0_dp
372  ALLOCATE (kijab(maxs, maxs, nkind, nkind))
373  kijab = 0.0_dp
374  DO ikind = 1, nkind
375  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
376  CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
377  IF (.NOT. defined .OR. natorb_a < 1) cycle
378  CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
379  DO jkind = 1, nkind
380  CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
381  CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
382  IF (.NOT. defined .OR. natorb_b < 1) cycle
383  CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
384  ! get Fen = (1+ken*deltaEN^2)
385  fen = 1.0_dp + ken*(ena - enb)**2
386  ! Kab
387  kab = xtb_set_kab(za, zb, xtb_control)
388  DO j = 1, natorb_b
389  lb = laob(j)
390  nb = naob(j)
391  DO i = 1, natorb_a
392  la = laoa(i)
393  na = naoa(i)
394  kia = kl(la)
395  kjb = kl(lb)
396  IF (zb == 1 .AND. nb == 2) kjb = k2sh
397  IF (za == 1 .AND. na == 2) kia = k2sh
398  IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
399  kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
400  ELSE
401  IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
402  kijab(i, j, ikind, jkind) = ksp*kab*fen
403  ELSE
404  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  IF (xb_inter) THEN
413  ! list of neighbor atoms for XB term
414  ALLOCATE (neighbor_atoms(nkind))
415  DO ikind = 1, nkind
416  NULLIFY (neighbor_atoms(ikind)%coord)
417  NULLIFY (neighbor_atoms(ikind)%rab)
418  NULLIFY (neighbor_atoms(ikind)%katom)
419  CALL get_atomic_kind(atomic_kind_set(ikind), z=zat, natom=na)
420  IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
421  ALLOCATE (neighbor_atoms(ikind)%coord(3, na))
422  neighbor_atoms(ikind)%coord(1:3, 1:na) = 0.0_dp
423  ALLOCATE (neighbor_atoms(ikind)%rab(na))
424  neighbor_atoms(ikind)%rab(1:na) = huge(0.0_dp)
425  ALLOCATE (neighbor_atoms(ikind)%katom(na))
426  neighbor_atoms(ikind)%katom(1:na) = 0
427  END IF
428  END DO
429  ! kx parameters
430  ALLOCATE (kx(nkind))
431  DO ikind = 1, nkind
432  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
433  CALL get_xtb_atom_param(xtb_atom_a, kx=kx(ikind))
434  END DO
435  !
436  ALLOCATE (rcab(nkind, nkind))
437  DO ikind = 1, nkind
438  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
439  CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova)
440  DO jkind = 1, nkind
441  CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
442  CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb)
443  rcab(ikind, jkind) = kxr*(rcova + rcovb)
444  END DO
445  END DO
446  !
447  CALL get_qs_env(qs_env=qs_env, sab_xb=sab_xb)
448  END IF
449 
450  ! initialize repulsion energy
451  erep = 0._dp
452 
453  ! loop over all atom pairs with a non-zero overlap (sab_orb)
454  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
455  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
456  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
457  iatom=iatom, jatom=jatom, r=rij, cell=cell)
458  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
459  CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
460  IF (.NOT. defined .OR. natorb_a < 1) cycle
461  CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
462  CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
463  IF (.NOT. defined .OR. natorb_b < 1) cycle
464 
465  dr = sqrt(sum(rij(:)**2))
466 
467  ! neighbor atom for XB term
468  IF (xb_inter .AND. (dr > 1.e-3_dp)) THEN
469  IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
470  atom_a = atom_of_kind(iatom)
471  IF (dr < neighbor_atoms(ikind)%rab(atom_a)) THEN
472  neighbor_atoms(ikind)%rab(atom_a) = dr
473  neighbor_atoms(ikind)%coord(1:3, atom_a) = rij(1:3)
474  neighbor_atoms(ikind)%katom(atom_a) = jatom
475  END IF
476  END IF
477  IF (ASSOCIATED(neighbor_atoms(jkind)%rab)) THEN
478  atom_b = atom_of_kind(jatom)
479  IF (dr < neighbor_atoms(jkind)%rab(atom_b)) THEN
480  neighbor_atoms(jkind)%rab(atom_b) = dr
481  neighbor_atoms(jkind)%coord(1:3, atom_b) = -rij(1:3)
482  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  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  kappa=kappab, hen=henb)
494 
495  IF (nimg == 1) THEN
496  ic = 1
497  ELSE
498  ic = cell_to_index(cell(1), cell(2), cell(3))
499  cpassert(ic > 0)
500  END IF
501 
502  icol = max(iatom, jatom)
503  irow = min(iatom, jatom)
504  NULLIFY (sblock, fblock)
505  CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
506  row=irow, col=icol, block=sblock, found=found)
507  cpassert(found)
508  CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
509  row=irow, col=icol, block=fblock, found=found)
510  cpassert(found)
511 
512  IF (calculate_forces) THEN
513  NULLIFY (pblock)
514  CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
515  row=irow, col=icol, block=pblock, found=found)
516  cpassert(ASSOCIATED(pblock))
517  NULLIFY (wblock)
518  CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
519  row=irow, col=icol, block=wblock, found=found)
520  cpassert(ASSOCIATED(wblock))
521  DO i = 2, 4
522  NULLIFY (dsblocks(i)%block)
523  CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
524  row=irow, col=icol, block=dsblocks(i)%block, found=found)
525  cpassert(found)
526  END DO
527  END IF
528 
529  ! overlap
530  basis_set_a => basis_set_list(ikind)%gto_basis_set
531  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
532  basis_set_b => basis_set_list(jkind)%gto_basis_set
533  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
534  atom_a = atom_of_kind(iatom)
535  atom_b = atom_of_kind(jatom)
536  ! basis ikind
537  first_sgfa => basis_set_a%first_sgf
538  la_max => basis_set_a%lmax
539  la_min => basis_set_a%lmin
540  npgfa => basis_set_a%npgf
541  nseta = basis_set_a%nset
542  nsgfa => basis_set_a%nsgf_set
543  rpgfa => basis_set_a%pgf_radius
544  set_radius_a => basis_set_a%set_radius
545  scon_a => basis_set_a%scon
546  zeta => basis_set_a%zet
547  ! basis jkind
548  first_sgfb => basis_set_b%first_sgf
549  lb_max => basis_set_b%lmax
550  lb_min => basis_set_b%lmin
551  npgfb => basis_set_b%npgf
552  nsetb = basis_set_b%nset
553  nsgfb => basis_set_b%nsgf_set
554  rpgfb => basis_set_b%pgf_radius
555  set_radius_b => basis_set_b%set_radius
556  scon_b => basis_set_b%scon
557  zetb => basis_set_b%zet
558 
559  ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
560  ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
561  ALLOCATE (sint(natorb_a, natorb_b, maxder))
562  sint = 0.0_dp
563 
564  DO iset = 1, nseta
565  ncoa = npgfa(iset)*ncoset(la_max(iset))
566  n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
567  sgfa = first_sgfa(1, iset)
568  DO jset = 1, nsetb
569  IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
570  ncob = npgfb(jset)*ncoset(lb_max(jset))
571  n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
572  sgfb = first_sgfb(1, jset)
573  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  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  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  cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
585  CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
586  IF (calculate_forces) THEN
587  DO i = 2, 4
588  CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
589  cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
590  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  IF (calculate_forces) THEN
597  DO i = 1, 3
598  IF (iatom <= jatom) THEN
599  force_ab(i) = sum(sint(:, :, i + 1)*wblock(:, :))
600  ELSE
601  force_ab(i) = sum(sint(:, :, i + 1)*transpose(wblock(:, :)))
602  END IF
603  END DO
604  f1 = 2.0_dp
605  force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
606  force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
607  IF (use_virial .AND. dr > 1.e-3_dp) THEN
608  IF (iatom == jatom) f1 = 1.0_dp
609  CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
610  IF (atprop%stress) THEN
611  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_ab, rij)
612  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  IF (iatom <= jatom) THEN
618  sblock(:, :) = sblock(:, :) + sint(:, :, 1)
619  ELSE
620  sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
621  END IF
622  IF (calculate_forces) THEN
623  DO i = 2, 4
624  IF (iatom <= jatom) THEN
625  dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
626  ELSE
627  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  rcovab = rcova + rcovb
634  rrab = sqrt(dr/rcovab)
635  DO i = 1, nsa
636  pia(i) = 1._dp + kpolya(i)*rrab
637  END DO
638  DO i = 1, nsb
639  pib(i) = 1._dp + kpolyb(i)*rrab
640  END DO
641  IF (calculate_forces) THEN
642  IF (dr > 1.e-6_dp) THEN
643  drx = 0.5_dp/rrab/rcovab
644  ELSE
645  drx = 0.0_dp
646  END IF
647  dpia(1:nsa) = drx*kpolya(1:nsa)
648  dpib(1:nsb) = drx*kpolyb(1:nsb)
649  END IF
650 
651  ! diagonal block
652  diagblock = .false.
653  IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .true.
654  !
655  ! Eq. 10
656  !
657  IF (diagblock) THEN
658  DO i = 1, natorb_a
659  na = naoa(i)
660  fblock(i, i) = fblock(i, i) + huckel(na, iatom)
661  END DO
662  ELSE
663  DO j = 1, natorb_b
664  nb = naob(j)
665  DO i = 1, natorb_a
666  na = naoa(i)
667  hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
668  IF (iatom <= jatom) THEN
669  fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
670  ELSE
671  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  IF (calculate_forces) THEN
677  f0 = 1.0_dp
678  IF (irow == iatom) f0 = -1.0_dp
679  ! Derivative wrt coordination number
680  fhua = 0.0_dp
681  fhub = 0.0_dp
682  fhud = 0.0_dp
683  IF (diagblock) THEN
684  DO i = 1, natorb_a
685  la = laoa(i)
686  na = naoa(i)
687  fhud = fhud + pblock(i, i)*kcnlk(la, ikind)*hena(na)
688  END DO
689  ELSE
690  DO j = 1, natorb_b
691  lb = laob(j)
692  nb = naob(j)
693  DO i = 1, natorb_a
694  la = laoa(i)
695  na = naoa(i)
696  hij = 0.5_dp*pia(na)*pib(nb)
697  IF (iatom <= jatom) THEN
698  fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(la, ikind)*hena(na)
699  fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(lb, jkind)*henb(nb)
700  ELSE
701  fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(la, ikind)*hena(na)
702  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  IF (iatom /= jatom) THEN
707  fhua = 2.0_dp*fhua
708  fhub = 2.0_dp*fhub
709  END IF
710  END IF
711  ! iatom
712  atom_a = atom_of_kind(iatom)
713  DO i = 1, dcnum(iatom)%neighbors
714  katom = dcnum(iatom)%nlist(i)
715  kkind = kind_of(katom)
716  atom_c = atom_of_kind(katom)
717  rik = dcnum(iatom)%rik(:, i)
718  drk = sqrt(sum(rik(:)**2))
719  IF (drk > 1.e-3_dp) THEN
720  fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
721  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
722  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
723  fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
724  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
725  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
726  IF (use_virial) THEN
727  fdik = fdika + fdikb
728  CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
729  IF (atprop%stress) THEN
730  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
731  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  atom_b = atom_of_kind(jatom)
738  DO i = 1, dcnum(jatom)%neighbors
739  katom = dcnum(jatom)%nlist(i)
740  kkind = kind_of(katom)
741  atom_c = atom_of_kind(katom)
742  rik = dcnum(jatom)%rik(:, i)
743  drk = sqrt(sum(rik(:)**2))
744  IF (drk > 1.e-3_dp) THEN
745  fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
746  force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
747  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
748  IF (use_virial) THEN
749  CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
750  IF (atprop%stress) THEN
751  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
752  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  IF (diagblock) THEN
758  force_ab = 0._dp
759  ELSE
760  ! force from R dendent Huckel element
761  n1 = SIZE(fblock, 1)
762  n2 = SIZE(fblock, 2)
763  ALLOCATE (dfblock(n1, n2))
764  dfblock = 0.0_dp
765  DO j = 1, natorb_b
766  lb = laob(j)
767  nb = naob(j)
768  DO i = 1, natorb_a
769  la = laoa(i)
770  na = naoa(i)
771  dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
772  IF (iatom <= jatom) THEN
773  dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
774  ELSE
775  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  dfp = f0*sum(dfblock(:, :)*pblock(:, :))
780  DO ir = 1, 3
781  foab = 2.0_dp*dfp*rij(ir)/dr
782  ! force from overlap matrix contribution to H
783  DO j = 1, natorb_b
784  lb = laob(j)
785  nb = naob(j)
786  DO i = 1, natorb_a
787  la = laoa(i)
788  na = naoa(i)
789  hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
790  IF (iatom <= jatom) THEN
791  foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
792  ELSE
793  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  force_ab(ir) = foab
798  END DO
799  DEALLOCATE (dfblock)
800  END IF
801  END IF
802 
803  IF (calculate_forces) THEN
804  atom_a = atom_of_kind(iatom)
805  atom_b = atom_of_kind(jatom)
806  IF (irow == iatom) force_ab = -force_ab
807  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
808  force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
809  IF (use_virial) THEN
810  f1 = 1.0_dp
811  IF (iatom == jatom) f1 = 0.5_dp
812  CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
813  IF (atprop%stress) THEN
814  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_ab, rij)
815  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  DEALLOCATE (oint, owork, sint)
821 
822  ! repulsive potential
823  IF (dr > 0.001_dp) THEN
824  erepij = zneffa*zneffb/dr*exp(-sqrt(alphaa*alphab)*dr**kf)
825  erep = erep + erepij
826  IF (atprop%energy) THEN
827  atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
828  atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
829  END IF
830  IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
831  derepij = -(1.0_dp/dr + sqrt(alphaa*alphab)*kf*dr**(kf - 1.0_dp))*erepij
832  force_rr(1) = derepij*rij(1)/dr
833  force_rr(2) = derepij*rij(2)/dr
834  force_rr(3) = derepij*rij(3)/dr
835  atom_a = atom_of_kind(iatom)
836  atom_b = atom_of_kind(jatom)
837  force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
838  force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
839  IF (use_virial) THEN
840  f1 = 1.0_dp
841  IF (iatom == jatom) f1 = 0.5_dp
842  CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
843  IF (atprop%stress) THEN
844  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*f1, force_rr, rij)
845  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  CALL neighbor_list_iterator_release(nl_iterator)
853 
854  DO i = 1, SIZE(matrix_h, 1)
855  DO img = 1, nimg
856  CALL dbcsr_finalize(matrix_h(i, img)%matrix)
857  CALL dbcsr_finalize(matrix_s(i, img)%matrix)
858  END DO
859  END DO
860 
861  exb = 0.0_dp
862  IF (xb_inter) THEN
863  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  calculate_forces, use_virial, force, virial, atprop)
866  END IF
867 
868  enonbonded = 0.0_dp
869  IF (do_nonbonded) THEN
870  ! nonbonded interactions
871  NULLIFY (sab_xtb_nonbond)
872  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  atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
875  END IF
876  ! set repulsive energy
877  erep = erep + exb + enonbonded
878  IF (xb_inter) THEN
879  CALL para_env%sum(exb)
880  energy%xtb_xb_inter = exb
881  END IF
882  IF (do_nonbonded) THEN
883  CALL para_env%sum(enonbonded)
884  energy%xtb_nonbonded = enonbonded
885  END IF
886  CALL para_env%sum(erep)
887  energy%repulsive = erep
888 
889  ! deallocate coordination numbers
890  DEALLOCATE (cnumbers)
891  IF (calculate_forces) THEN
892  DO iatom = 1, natom
893  DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
894  END DO
895  END IF
896  DEALLOCATE (dcnum)
897 
898  ! deallocate Huckel parameters
899  DEALLOCATE (huckel)
900  ! deallocate KAB parameters
901  DEALLOCATE (kijab)
902 
903  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
904  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  extension=".Log")
908  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
909  after = min(max(after, 1), 16)
910  DO img = 1, nimg
911  CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
912  output_unit=iw, omit_headers=omit_headers)
913  END DO
914  CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
915  END IF
916 
917  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  extension=".Log")
921  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
922  after = min(max(after, 1), 16)
923  DO img = 1, nimg
924  CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
925  output_unit=iw, omit_headers=omit_headers)
926  IF (btest(cp_print_key_should_output(logger%iter_info, &
927  qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
928  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  output_unit=iw, omit_headers=omit_headers)
931  END DO
932  END IF
933  END DO
934  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  IF (.NOT. calculate_forces) THEN
939  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  extension=".Log")
943  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
944  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
945  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
946  CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
947  CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
948  END IF
949  END IF
950 
951  DEALLOCATE (basis_set_list)
952  IF (calculate_forces) THEN
953  IF (SIZE(matrix_p, 1) == 2) THEN
954  DO img = 1, nimg
955  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
956  beta_scalar=-1.0_dp)
957  END DO
958  END IF
959  END IF
960 
961  IF (xb_inter) THEN
962  DO ikind = 1, nkind
963  IF (ASSOCIATED(neighbor_atoms(ikind)%coord)) THEN
964  DEALLOCATE (neighbor_atoms(ikind)%coord)
965  END IF
966  IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
967  DEALLOCATE (neighbor_atoms(ikind)%rab)
968  END IF
969  IF (ASSOCIATED(neighbor_atoms(ikind)%katom)) THEN
970  DEALLOCATE (neighbor_atoms(ikind)%katom)
971  END IF
972  END DO
973  DEALLOCATE (neighbor_atoms)
974  DEALLOCATE (kx, rcab)
975  END IF
976 
977  CALL timestop(handle)
978 
979  END SUBROUTINE build_xtb_matrices
980 
981 ! **************************************************************************************************
982 !> \brief ...
983 !> \param qs_env ...
984 !> \param p_matrix ...
985 ! **************************************************************************************************
986  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
998  INTEGER, DIMENSION(25) :: laoa, laob, lval, naoa, naob
999  INTEGER, DIMENSION(3) :: cell
1000  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1001  npgfb, nsgfa, nsgfb
1002  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1003  LOGICAL :: defined, diagblock, floating_a, found, &
1004  ghost_a, use_virial
1005  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
1010  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, huckel, kcnlk, owork
1011  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
1012  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  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1018  REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
1019  scon_a, scon_b, zeta, zetb
1020  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1021  TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
1022  TYPE(cp_logger_type), POINTER :: logger
1023  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1024  TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
1025  TYPE(dft_control_type), POINTER :: dft_control
1026  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  DIMENSION(:), POINTER :: nl_iterator
1031  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1032  POINTER :: sab_orb
1033  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1034  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1035  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1036  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  CALL timeset(routinen, handle)
1042 
1043  NULLIFY (logger)
1044  logger => cp_get_default_logger()
1045 
1046  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  sab_orb=sab_orb)
1054 
1055  nkind = SIZE(atomic_kind_set)
1056  xtb_control => dft_control%qs_control%xtb_control
1057  nimg = dft_control%nimages
1058  nderivatives = 1
1059  maxder = ncoset(nderivatives)
1060 
1061  ! global parameters (Table 2 from Ref.)
1062  ks = xtb_control%ks
1063  kp = xtb_control%kp
1064  kd = xtb_control%kd
1065  ksp = xtb_control%ksp
1066  k2sh = xtb_control%k2sh
1067  kg = xtb_control%kg
1068  kf = xtb_control%kf
1069  kcns = xtb_control%kcns
1070  kcnp = xtb_control%kcnp
1071  kcnd = xtb_control%kcnd
1072  ken = xtb_control%ken
1073  kxr = xtb_control%kxr
1074  kx2 = xtb_control%kx2
1075 
1076  NULLIFY (particle_set)
1077  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
1078  natom = SIZE(particle_set)
1079  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1080  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
1081 
1082  NULLIFY (force)
1083  CALL get_qs_env(qs_env=qs_env, force=force)
1084  use_virial = .false.
1085  cpassert(nimg == 1)
1086 
1087  ! set up basis set lists
1088  ALLOCATE (basis_set_list(nkind))
1089  CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1090 
1091  ! allocate overlap matrix
1092  CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
1093  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  sab_orb, .true.)
1096  ! initialize H matrix
1097  CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
1098  DO img = 1, nimg
1099  ALLOCATE (matrix_h(1, img)%matrix)
1100  CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
1101  name="HAMILTONIAN MATRIX")
1102  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  ALLOCATE (cnumbers(natom))
1109  cnumbers = 0._dp
1110  ALLOCATE (dcnum(natom))
1111  dcnum(:)%neighbors = 0
1112  DO iatom = 1, natom
1113  ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
1114  END DO
1115 
1116  ALLOCATE (ghost(nkind), floating(nkind), atomnumber(nkind))
1117  DO ikind = 1, nkind
1118  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1119  CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_a, floating=floating_a)
1120  ghost(ikind) = ghost_a
1121  floating(ikind) = floating_a
1122  atomnumber(ikind) = za
1123  END DO
1124  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  .true., .false.)
1127  DEALLOCATE (ghost, floating, atomnumber)
1128  CALL para_env%sum(cnumbers)
1129  CALL dcnum_distribute(dcnum, para_env)
1130 
1131  ! Calculate Huckel parameters
1132  ! Eq 12
1133  ! huckel(nshell,natom)
1134  ALLOCATE (kcnlk(0:3, nkind))
1135  DO ikind = 1, nkind
1136  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1137  IF (metal(za)) THEN
1138  kcnlk(0:3, ikind) = 0.0_dp
1139  ELSEIF (early3d(za)) THEN
1140  kcnlk(0, ikind) = kcns
1141  kcnlk(1, ikind) = kcnp
1142  kcnlk(2, ikind) = 0.005_dp
1143  kcnlk(3, ikind) = 0.0_dp
1144  ELSE
1145  kcnlk(0, ikind) = kcns
1146  kcnlk(1, ikind) = kcnp
1147  kcnlk(2, ikind) = kcnd
1148  kcnlk(3, ikind) = 0.0_dp
1149  END IF
1150  END DO
1151  ALLOCATE (huckel(5, natom))
1152  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1153  DO iatom = 1, natom
1154  ikind = kind_of(iatom)
1155  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1156  CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
1157  huckel(:, iatom) = 0.0_dp
1158  DO i = 1, nshell
1159  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  kl(0) = ks
1165  kl(1) = kp
1166  kl(2) = kd
1167  kl(3) = 0.0_dp
1168  ALLOCATE (kijab(maxs, maxs, nkind, nkind))
1169  kijab = 0.0_dp
1170  DO ikind = 1, nkind
1171  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1172  CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
1173  IF (.NOT. defined .OR. natorb_a < 1) cycle
1174  CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
1175  DO jkind = 1, nkind
1176  CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1177  CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
1178  IF (.NOT. defined .OR. natorb_b < 1) cycle
1179  CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
1180  ! get Fen = (1+ken*deltaEN^2)
1181  fen = 1.0_dp + ken*(ena - enb)**2
1182  ! Kab
1183  kab = xtb_set_kab(za, zb, xtb_control)
1184  DO j = 1, natorb_b
1185  lb = laob(j)
1186  nb = naob(j)
1187  DO i = 1, natorb_a
1188  la = laoa(i)
1189  na = naoa(i)
1190  kia = kl(la)
1191  kjb = kl(lb)
1192  IF (zb == 1 .AND. nb == 2) kjb = k2sh
1193  IF (za == 1 .AND. na == 2) kia = k2sh
1194  IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
1195  kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
1196  ELSE
1197  IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
1198  kijab(i, j, ikind, jkind) = ksp*kab*fen
1199  ELSE
1200  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  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1210  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1211  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1212  iatom=iatom, jatom=jatom, r=rij, cell=cell)
1213  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1214  CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
1215  IF (.NOT. defined .OR. natorb_a < 1) cycle
1216  CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1217  CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
1218  IF (.NOT. defined .OR. natorb_b < 1) cycle
1219 
1220  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  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  kappa=kappab, hen=henb)
1229 
1230  ic = 1
1231 
1232  icol = max(iatom, jatom)
1233  irow = min(iatom, jatom)
1234  NULLIFY (sblock, fblock)
1235  CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1236  row=irow, col=icol, block=sblock, found=found)
1237  cpassert(found)
1238  CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
1239  row=irow, col=icol, block=fblock, found=found)
1240  cpassert(found)
1241 
1242  NULLIFY (pblock)
1243  CALL dbcsr_get_block_p(matrix=p_matrix, &
1244  row=irow, col=icol, block=pblock, found=found)
1245  cpassert(ASSOCIATED(pblock))
1246  DO i = 2, 4
1247  NULLIFY (dsblocks(i)%block)
1248  CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
1249  row=irow, col=icol, block=dsblocks(i)%block, found=found)
1250  cpassert(found)
1251  END DO
1252 
1253  ! overlap
1254  basis_set_a => basis_set_list(ikind)%gto_basis_set
1255  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1256  basis_set_b => basis_set_list(jkind)%gto_basis_set
1257  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1258  atom_a = atom_of_kind(iatom)
1259  atom_b = atom_of_kind(jatom)
1260  ! basis ikind
1261  first_sgfa => basis_set_a%first_sgf
1262  la_max => basis_set_a%lmax
1263  la_min => basis_set_a%lmin
1264  npgfa => basis_set_a%npgf
1265  nseta = basis_set_a%nset
1266  nsgfa => basis_set_a%nsgf_set
1267  rpgfa => basis_set_a%pgf_radius
1268  set_radius_a => basis_set_a%set_radius
1269  scon_a => basis_set_a%scon
1270  zeta => basis_set_a%zet
1271  ! basis jkind
1272  first_sgfb => basis_set_b%first_sgf
1273  lb_max => basis_set_b%lmax
1274  lb_min => basis_set_b%lmin
1275  npgfb => basis_set_b%npgf
1276  nsetb = basis_set_b%nset
1277  nsgfb => basis_set_b%nsgf_set
1278  rpgfb => basis_set_b%pgf_radius
1279  set_radius_b => basis_set_b%set_radius
1280  scon_b => basis_set_b%scon
1281  zetb => basis_set_b%zet
1282 
1283  ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1284  ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1285  ALLOCATE (sint(natorb_a, natorb_b, maxder))
1286  sint = 0.0_dp
1287 
1288  DO iset = 1, nseta
1289  ncoa = npgfa(iset)*ncoset(la_max(iset))
1290  n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1291  sgfa = first_sgfa(1, iset)
1292  DO jset = 1, nsetb
1293  IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
1294  ncob = npgfb(jset)*ncoset(lb_max(jset))
1295  n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1296  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  rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1300  ! Contraction
1301  DO i = 1, 4
1302  CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1303  cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1304  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  IF (iatom <= jatom) THEN
1310  sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1311  ELSE
1312  sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
1313  END IF
1314  DO i = 2, 4
1315  IF (iatom <= jatom) THEN
1316  dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1317  ELSE
1318  dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - transpose(sint(:, :, i))
1319  END IF
1320  END DO
1321 
1322  ! Calculate Pi = Pia * Pib (Eq. 11)
1323  rcovab = rcova + rcovb
1324  rrab = sqrt(dr/rcovab)
1325  DO i = 1, nsa
1326  pia(i) = 1._dp + kpolya(i)*rrab
1327  END DO
1328  DO i = 1, nsb
1329  pib(i) = 1._dp + kpolyb(i)*rrab
1330  END DO
1331  IF (dr > 1.e-6_dp) THEN
1332  drx = 0.5_dp/rrab/rcovab
1333  ELSE
1334  drx = 0.0_dp
1335  END IF
1336  dpia(1:nsa) = drx*kpolya(1:nsa)
1337  dpib(1:nsb) = drx*kpolyb(1:nsb)
1338 
1339  ! diagonal block
1340  diagblock = .false.
1341  IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .true.
1342  !
1343  ! Eq. 10
1344  !
1345  IF (diagblock) THEN
1346  DO i = 1, natorb_a
1347  na = naoa(i)
1348  fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1349  END DO
1350  ELSE
1351  DO j = 1, natorb_b
1352  nb = naob(j)
1353  DO i = 1, natorb_a
1354  na = naoa(i)
1355  hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1356  IF (iatom <= jatom) THEN
1357  fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1358  ELSE
1359  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  f0 = 1.0_dp
1366  IF (irow == iatom) f0 = -1.0_dp
1367  ! Derivative wrt coordination number
1368  fhua = 0.0_dp
1369  fhub = 0.0_dp
1370  fhud = 0.0_dp
1371  IF (diagblock) THEN
1372  DO i = 1, natorb_a
1373  la = laoa(i)
1374  na = naoa(i)
1375  fhud = fhud + pblock(i, i)*kcnlk(la, ikind)*hena(na)
1376  END DO
1377  ELSE
1378  DO j = 1, natorb_b
1379  lb = laob(j)
1380  nb = naob(j)
1381  DO i = 1, natorb_a
1382  la = laoa(i)
1383  na = naoa(i)
1384  hij = 0.5_dp*pia(na)*pib(nb)
1385  IF (iatom <= jatom) THEN
1386  fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(la, ikind)*hena(na)
1387  fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*kcnlk(lb, jkind)*henb(nb)
1388  ELSE
1389  fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*kcnlk(la, ikind)*hena(na)
1390  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  IF (iatom /= jatom) THEN
1395  fhua = 2.0_dp*fhua
1396  fhub = 2.0_dp*fhub
1397  END IF
1398  END IF
1399  ! iatom
1400  atom_a = atom_of_kind(iatom)
1401  DO i = 1, dcnum(iatom)%neighbors
1402  katom = dcnum(iatom)%nlist(i)
1403  kkind = kind_of(katom)
1404  atom_c = atom_of_kind(katom)
1405  rik = dcnum(iatom)%rik(:, i)
1406  drk = sqrt(sum(rik(:)**2))
1407  IF (drk > 1.e-3_dp) THEN
1408  fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1409  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1410  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1411  fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1412  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1413  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1414  END IF
1415  END DO
1416  ! jatom
1417  atom_b = atom_of_kind(jatom)
1418  DO i = 1, dcnum(jatom)%neighbors
1419  katom = dcnum(jatom)%nlist(i)
1420  kkind = kind_of(katom)
1421  atom_c = atom_of_kind(katom)
1422  rik = dcnum(jatom)%rik(:, i)
1423  drk = sqrt(sum(rik(:)**2))
1424  IF (drk > 1.e-3_dp) THEN
1425  fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1426  force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1427  force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1428  END IF
1429  END DO
1430  IF (diagblock) THEN
1431  force_ab = 0._dp
1432  ELSE
1433  ! force from R dendent Huckel element
1434  n1 = SIZE(fblock, 1)
1435  n2 = SIZE(fblock, 2)
1436  ALLOCATE (dfblock(n1, n2))
1437  dfblock = 0.0_dp
1438  DO j = 1, natorb_b
1439  lb = laob(j)
1440  nb = naob(j)
1441  DO i = 1, natorb_a
1442  la = laoa(i)
1443  na = naoa(i)
1444  dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1445  IF (iatom <= jatom) THEN
1446  dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1447  ELSE
1448  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  dfp = f0*sum(dfblock(:, :)*pblock(:, :))
1453  DO ir = 1, 3
1454  foab = 2.0_dp*dfp*rij(ir)/dr
1455  ! force from overlap matrix contribution to H
1456  DO j = 1, natorb_b
1457  lb = laob(j)
1458  nb = naob(j)
1459  DO i = 1, natorb_a
1460  la = laoa(i)
1461  na = naoa(i)
1462  hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1463  IF (iatom <= jatom) THEN
1464  foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1465  ELSE
1466  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  force_ab(ir) = foab
1471  END DO
1472  DEALLOCATE (dfblock)
1473  END IF
1474 
1475  atom_a = atom_of_kind(iatom)
1476  atom_b = atom_of_kind(jatom)
1477  IF (irow == iatom) force_ab = -force_ab
1478  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1479  force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1480 
1481  DEALLOCATE (oint, owork, sint)
1482 
1483  END DO
1484  CALL neighbor_list_iterator_release(nl_iterator)
1485 
1486  DO i = 1, SIZE(matrix_h, 1)
1487  DO img = 1, nimg
1488  CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1489  CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1490  END DO
1491  END DO
1492  CALL dbcsr_deallocate_matrix_set(matrix_s)
1493  CALL dbcsr_deallocate_matrix_set(matrix_h)
1494 
1495  ! deallocate coordination numbers
1496  DEALLOCATE (cnumbers)
1497  DO iatom = 1, natom
1498  DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
1499  END DO
1500  DEALLOCATE (dcnum)
1501 
1502  ! deallocate Huckel parameters
1503  DEALLOCATE (huckel)
1504  ! deallocate KAB parameters
1505  DEALLOCATE (kijab)
1506 
1507  DEALLOCATE (basis_set_list)
1508 
1509  CALL timestop(handle)
1510 
1511  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
1536  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, charges
1537  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1538  TYPE(cp_logger_type), POINTER :: logger
1539  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs, p_matrix
1540  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1545  TYPE(qs_energy_type), POINTER :: energy
1546  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  CALL timeset(routinen, handle)
1554 
1555  NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
1556  ks_matrix, rho, energy)
1557  cpassert(ASSOCIATED(qs_env))
1558 
1559  logger => cp_get_default_logger()
1560  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  energy=energy)
1572 
1573  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  ns = SIZE(ext_ks_matrix)
1577  ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
1578  END IF
1579 
1580  energy%hartree = 0.0_dp
1581  energy%qmmm_el = 0.0_dp
1582  energy%efield = 0.0_dp
1583 
1584  scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
1585  nspins = dft_control%nspins
1586  nimg = dft_control%nimages
1587  cpassert(ASSOCIATED(matrix_h))
1588  cpassert(ASSOCIATED(rho))
1589  cpassert(SIZE(ks_matrix) > 0)
1590 
1591  DO ispin = 1, nspins
1592  DO img = 1, nimg
1593  ! copy the core matrix into the fock matrix
1594  CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
1595  END DO
1596  END DO
1597 
1598  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  do_efield = .false.
1603  END IF
1604 
1605  IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
1606  ! Mulliken charges
1607  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
1608  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1609  natom = SIZE(particle_set)
1610  ALLOCATE (mcharge(natom), charges(natom, 5))
1611  charges = 0.0_dp
1612  nkind = SIZE(atomic_kind_set)
1613  CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1614  ALLOCATE (aocg(nsgf, natom))
1615  aocg = 0.0_dp
1616  IF (nimg > 1) THEN
1617  CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
1618  ELSE
1619  p_matrix => matrix_p(:, 1)
1620  s_matrix => matrix_s(1, 1)%matrix
1621  CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1622  END IF
1623  DO ikind = 1, nkind
1624  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1625  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1626  CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
1627  DO iatom = 1, na
1628  atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1629  charges(atom_a, :) = real(occ(:), kind=dp)
1630  DO is = 1, natorb
1631  ns = lao(is) + 1
1632  charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
1633  END DO
1634  END DO
1635  END DO
1636  DEALLOCATE (aocg)
1637  ! charge mixing
1638  IF (dft_control%qs_control%do_ls_scf) THEN
1639  !
1640  ELSE
1641  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  charges, para_env, scf_env%iter_count)
1644  END IF
1645 
1646  DO iatom = 1, natom
1647  mcharge(iatom) = sum(charges(iatom, :))
1648  END DO
1649  END IF
1650 
1651  IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1652  CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
1653  calculate_forces, just_energy)
1654  END IF
1655 
1656  IF (do_efield) THEN
1657  CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
1658  END IF
1659 
1660  IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1661  IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
1662  pass_check = .true.
1663  DO ikind = 1, nkind
1664  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1665  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1666  CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
1667  DO iatom = 1, na
1668  atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1669  achrg = mcharge(atom_a)
1670  IF (abs(achrg) > chmax) THEN
1671  IF (iounit > 0) THEN
1672  WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
1673  " 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  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  " if you want to force to continue the calculation.")
1683  cpabort("xTB Charges")
1684  END IF
1685  END IF
1686  END IF
1687 
1688  IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
1689  DEALLOCATE (mcharge, charges)
1690  END IF
1691 
1692  IF (qs_env%qmmm) THEN
1693  cpassert(SIZE(ks_matrix, 2) == 1)
1694  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  1.0_dp, 1.0_dp)
1698  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  matrix_p1(ispin)%matrix, qmmm_el)
1702  energy%qmmm_el = energy%qmmm_el + qmmm_el
1703  END DO
1704  pc_ener = qs_env%ks_qmmm_env%pc_ener
1705  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  energy%repulsive + energy%dispersion + energy%dftb3
1710 
1711  iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
1712  extension=".scfLog")
1713  IF (iounit > 0) THEN
1714  WRITE (unit=iounit, fmt="(/,(T9,A,T60,F20.10))") &
1715  "Repulsive pair potential energy: ", energy%repulsive, &
1716  "Zeroth order Hamiltonian energy: ", energy%core, &
1717  "Charge fluctuation energy: ", energy%hartree, &
1718  "London dispersion energy: ", energy%dispersion
1719  IF (dft_control%qs_control%xtb_control%xb_interaction) &
1720  WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
1721  "Correction for halogen bonding: ", energy%xtb_xb_inter
1722  IF (dft_control%qs_control%xtb_control%do_nonbonded) &
1723  WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
1724  "Correction for nonbonded interactions: ", energy%xtb_nonbonded
1725  IF (abs(energy%efield) > 1.e-12_dp) THEN
1726  WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
1727  "Electric field interaction energy: ", energy%efield
1728  END IF
1729  WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
1730  "DFTB3 3rd Order Energy Correction ", energy%dftb3
1731  IF (qs_env%qmmm) THEN
1732  WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
1733  "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  "PRINT%DETAILED_ENERGY")
1738  ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
1739  IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
1740  cpassert(SIZE(ks_matrix, 2) == 1)
1741  block
1742  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1743  CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
1744  DO ispin = 1, SIZE(mo_derivs)
1745  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
1746  IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
1747  cpabort("")
1748  END IF
1749  CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
1750  0.0_dp, mo_derivs(ispin)%matrix)
1751  END DO
1752  END block
1753  END IF
1754 
1755  CALL timestop(handle)
1756 
1757  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: matom
1775  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dmloc
1776  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
1777 
1778  nkind = SIZE(neighbor_atoms)
1779  DO ikind = 1, nkind
1780  IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
1781  natom = SIZE(neighbor_atoms(ikind)%rab)
1782  ALLOCATE (dmloc(2*natom), matom(natom), coord(3, natom))
1783  dmloc = 0.0_dp
1784  DO iatom = 1, natom
1785  dmloc(2*iatom - 1) = neighbor_atoms(ikind)%rab(iatom)
1786  dmloc(2*iatom) = real(para_env%mepos, kind=dp)
1787  END DO
1788  CALL para_env%minloc(dmloc)
1789  coord = 0.0_dp
1790  matom = 0
1791  DO iatom = 1, natom
1792  neighbor_atoms(ikind)%rab(iatom) = dmloc(2*iatom - 1)
1793  IF (nint(dmloc(2*iatom)) == para_env%mepos) THEN
1794  coord(1:3, iatom) = neighbor_atoms(ikind)%coord(1:3, iatom)
1795  matom(iatom) = neighbor_atoms(ikind)%katom(iatom)
1796  END IF
1797  END DO
1798  CALL para_env%sum(coord)
1799  neighbor_atoms(ikind)%coord(1:3, :) = coord(1:3, :)
1800  CALL para_env%sum(matom)
1801  neighbor_atoms(ikind)%katom(:) = matom(:)
1802  DEALLOCATE (dmloc, matom, coord)
1803  END IF
1804  END DO
1805 
1806  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  SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1826  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  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  CALL timeset(routinen, handle)
1860 
1861  NULLIFY (nonbonded)
1862  NULLIFY (potparm)
1863  NULLIFY (ewald_env)
1864  nonbonded => xtb_control%nonbonded
1865  do_ewald = xtb_control%do_ewald
1866  CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
1867 
1868  ntype = SIZE(atomic_kind_set)
1869  CALL pair_potential_pp_create(potparm, ntype)
1870  !Assign input and potential info to potparm_nonbond
1871  CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
1872  !Initialize genetic potential
1873  CALL init_genpot(potparm, ntype)
1874 
1875  NULLIFY (pot)
1876  enonbonded = 0._dp
1877  energy_cutoff = 0._dp
1878 
1879  CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond)
1880  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1881  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1882  iatom=iatom, jatom=jatom, r=rij, cell=cell)
1883  pot => potparm%pot(ikind, jkind)%pot
1884  dr = sqrt(rij(1)**2 + rij(2)**2 + rij(3)**2)
1885  rcut = sqrt(pot%rcutsq)
1886  IF (dr <= rcut .AND. dr > 1.e-3_dp) THEN
1887  fval = 1.0_dp
1888  IF (ikind == jkind) fval = 0.5_dp
1889  ! splines not implemented
1890  enonbonded = enonbonded + fval*ener_pot(pot, dr, energy_cutoff)
1891  IF (atprop%energy) THEN
1892  atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
1893  atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
1894  END IF
1895  END IF
1896 
1897  IF (calculate_forces) THEN
1898 
1899  kk = SIZE(pot%type)
1900  IF (kk /= 1) THEN
1901  CALL cp_warn(__location__, "Generic potential with type > 1 not implemented.")
1902  cpabort("pot type")
1903  END IF
1904  ! rmin and rmax and rcut
1905  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  IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) cycle
1908  ! If within limits let's compute the potential...
1909  IF (dr <= rcut .AND. dr > 1.e-3_dp) THEN
1910 
1911  NULLIFY (nonbonded_section)
1912  nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED")
1913  CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx)
1914  CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr)
1915 
1916  dedf = fval*evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err)
1917  IF (abs(err) > lerr) THEN
1918  WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1919  WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1920  CALL compress(this_error, .true.)
1921  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  trim(def_error)//' .')
1926  END IF
1927 
1928  atom_i = atom_of_kind(iatom)
1929  atom_j = atom_of_kind(jatom)
1930 
1931  fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
1932 
1933  force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
1934  force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:)
1935 
1936  IF (use_virial) THEN
1937  CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
1938  IF (atprop%stress) THEN
1939  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
1940  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  NULLIFY (pot)
1947  END DO
1948  CALL neighbor_list_iterator_release(nl_iterator)
1949  CALL finalizef()
1950 
1951  IF (ASSOCIATED(potparm)) THEN
1952  CALL pair_potential_pp_release(potparm)
1953  END IF
1954 
1955  CALL timestop(handle)
1956 
1957  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  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  NULLIFY (pot, logger)
1987 
1988  logger => cp_get_default_logger()
1989  iw = cp_logger_get_default_io_unit(logger)
1990 
1991  DO ikind = 1, SIZE(atomic_kind_set)
1992  atomic_kind => atomic_kind_set(ikind)
1993  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
1994  DO jkind = ikind, SIZE(atomic_kind_set)
1995  atomic_kind => atomic_kind_set(jkind)
1996  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
1997  found = .false.
1998  name_atm_a = name_atm_a_local
1999  name_atm_b = name_atm_b_local
2000  CALL uppercase(name_atm_a)
2001  CALL uppercase(name_atm_b)
2002  pot => potparm%pot(ikind, jkind)%pot
2003  IF (ASSOCIATED(nonbonded)) THEN
2004  DO ingp = 1, SIZE(nonbonded%pot)
2005  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  ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. &
2013  (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
2014  ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN
2015  CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot)
2016  ! multiple potential not implemented, simply overwriting
2017  IF (found) &
2018  CALL cp_warn(__location__, &
2019  "Multiple NONBONDED declaration: "//trim(name_atm_a)// &
2020  " 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  IF (.NOT. found) THEN
2027  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  IF (do_ewald) THEN
2035  CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
2036  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  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  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  DIMENSION(:), POINTER :: nl_iterator
2086 
2087  ! exonent in angular term
2088  alp = 6.0_dp
2089  ! loop over all atom pairs
2090  CALL neighbor_list_iterator_create(nl_iterator, sab_xb)
2091  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2092  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2093  iatom=iatom, jatom=jatom, r=rij, cell=cell)
2094  ! ikind, iatom : Halogen
2095  ! jkind, jatom : Donor
2096  atom_a = atom_of_kind(iatom)
2097  katom = neighbor_atoms(ikind)%katom(atom_a)
2098  IF (katom == 0) cycle
2099  dr = sqrt(rij(1)**2 + rij(2)**2 + rij(3)**2)
2100  ddr = rcab(ikind, jkind)/dr
2101  ddr6 = ddr**6
2102  ddr12 = ddr6*ddr6
2103  eval = kx(ikind)*(ddr12 - kx2*ddr6)/(1.0_dp + ddr12)
2104  ! angle
2105  ria(1:3) = neighbor_atoms(ikind)%coord(1:3, atom_a)
2106  rja(1:3) = rij(1:3) - ria(1:3)
2107  drax = ria(1)**2 + ria(2)**2 + ria(3)**2
2108  drbx = dr*dr
2109  drab = rja(1)**2 + rja(2)**2 + rja(3)**2
2110  xy = sqrt(drbx*drax)
2111  ! cos angle B-X-A
2112  cosa = (drbx + drax - drab)/xy
2113  aterm = (0.5_dp - 0.25_dp*cosa)**alp
2114  !
2115  exb = exb + aterm*eval
2116  IF (atprop%energy) THEN
2117  atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*aterm*eval
2118  atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*aterm*eval
2119  END IF
2120  !
2121  IF (calculate_forces) THEN
2122  kkind = kind_of(katom)
2123  atom_b = atom_of_kind(jatom)
2124  atom_c = atom_of_kind(katom)
2125  !
2126  deval = 6.0_dp*kx(ikind)*ddr6*(kx2*ddr12 + 2.0_dp*ddr6 - kx2)/(1.0_dp + ddr12)**2
2127  deval = -rcab(ikind, jkind)*deval/(dr*dr)/ddr
2128  fij(1:3) = aterm*deval*rij(1:3)/dr
2129  force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:)
2130  force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:)
2131  IF (use_virial) THEN
2132  CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
2133  IF (atprop%stress) THEN
2134  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
2135  CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij)
2136  END IF
2137  END IF
2138  !
2139  fij(1:3) = 0.0_dp
2140  fia(1:3) = 0.0_dp
2141  fja(1:3) = 0.0_dp
2142  daterm = -0.25_dp*alp*(0.5_dp - 0.25_dp*cosa)**(alp - 1.0_dp)
2143  ddbx = 0.5_dp*(drab - drax + drbx)/xy/drbx
2144  ddax = 0.5_dp*(drab + drax - drbx)/xy/drax
2145  ddab = -1._dp/xy
2146  fij(1:3) = 2.0_dp*daterm*ddbx*rij(1:3)*eval
2147  fia(1:3) = 2.0_dp*daterm*ddax*ria(1:3)*eval
2148  fja(1:3) = 2.0_dp*daterm*ddab*rja(1:3)*eval
2149  force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:) - fia(:)
2150  force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:) + fja(:)
2151  force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fia(:) - fja(:)
2152  IF (use_virial) THEN
2153  CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
2154  CALL virial_pair_force(virial%pv_virial, -1._dp, fia, ria)
2155  CALL virial_pair_force(virial%pv_virial, -1._dp, fja, rja)
2156  IF (atprop%stress) THEN
2157  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij)
2158  CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij)
2159  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fia, ria)
2160  CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fia, ria)
2161  CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fja, rja)
2162  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  CALL neighbor_list_iterator_release(nl_iterator)
2168 
2169  END SUBROUTINE xb_interaction
2170 
2171 ! **************************************************************************************************
2172 !> \brief ...
2173 !> \param z ...
2174 !> \return ...
2175 ! **************************************************************************************************
2176  FUNCTION metal(z) RESULT(ismetal)
2177  INTEGER :: z
2178  LOGICAL :: ismetal
2179 
2180  SELECT CASE (z)
2181  CASE DEFAULT
2182  ismetal = .true.
2183  CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
2184  ismetal = .false.
2185  END SELECT
2186 
2187  END FUNCTION metal
2188 
2189 ! **************************************************************************************************
2190 !> \brief ...
2191 !> \param z ...
2192 !> \return ...
2193 ! **************************************************************************************************
2194  FUNCTION early3d(z) RESULT(isearly3d)
2195  INTEGER :: z
2196  LOGICAL :: isearly3d
2197 
2198  isearly3d = .false.
2199  IF (z >= 21 .AND. z <= 24) isearly3d = .true.
2200 
2201  END FUNCTION early3d
2202 
2203 END MODULE xtb_matrices
2204 
uppercase
subroutine uppercase(string)
...
Definition: dumpdcd.F:1376
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition: ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition: ai_overlap.F:680
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition: atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition: atprop_types.F:98
basis_set_types
Definition: basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition: block_p_types.F:14
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_output
DBCSR output in CP2K.
Definition: cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition: cp_dbcsr_output.F:163
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
efield_tb_methods
Calculation of electric field contributions in TB.
Definition: efield_tb_methods.F:12
efield_tb_methods::efield_tb_matrix
subroutine, public efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Definition: efield_tb_methods.F:76
ewald_environment_types
Definition: ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition: ewald_environment_types.F:123
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition: fparser.F:17
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition: fparser.F:976
fparser::finalizef
subroutine, public finalizef()
...
Definition: fparser.F:101
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition: mulliken.F:13
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
pair_potential_types
Definition: pair_potential_types.F:14
pair_potential_types::pair_potential_pp_release
subroutine, public pair_potential_pp_release(potparm)
Release Data-structure that constains potential parameters.
Definition: pair_potential_types.F:969
pair_potential_types::pair_potential_single_copy
subroutine, public pair_potential_single_copy(potparm_source, potparm_dest)
Copy two potential parameter type.
Definition: pair_potential_types.F:700
pair_potential_types::pair_potential_single_clean
subroutine, public pair_potential_single_clean(potparm)
Cleans the potential parameter type.
Definition: pair_potential_types.F:651
pair_potential_types::not_initialized
real(kind=dp), parameter, public not_initialized
Definition: pair_potential_types.F:88
pair_potential_types::pair_potential_pp_create
subroutine, public pair_potential_pp_create(potparm, nkinds)
Data-structure that constains potential parameters.
Definition: pair_potential_types.F:938
pair_potential_util
Definition: pair_potential_util.F:14
pair_potential_util::ener_pot
real(kind=dp) function, public ener_pot(pot, r, energy_cutoff)
Evaluates the nonbond potential energy for the implemented FF kinds.
Definition: pair_potential_util.F:47
pair_potential
Definition: pair_potential.F:16
pair_potential::init_genpot
subroutine, public init_genpot(potparm, ntype)
Initialize genpot.
Definition: pair_potential.F:78
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_charge_mixing
Definition: qs_charge_mixing.F:9
qs_charge_mixing::charge_mixing
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count)
Driver for the charge mixing, calls the proper routine given the requested method.
Definition: qs_charge_mixing.F:42
qs_condnum
Calculation of overlap matrix condition numbers.
Definition: qs_condnum.F:13
qs_condnum::overlap_condnum
subroutine, public overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
Calculation of the overlap matrix Condition Number.
Definition: qs_condnum.F:64
qs_dispersion_pairpot
Calculation of dispersion using pair potentials.
Definition: qs_dispersion_pairpot.F:12
qs_dispersion_pairpot::dcnum_distribute
subroutine, public dcnum_distribute(dcnum, para_env)
...
Definition: qs_dispersion_pairpot.F:3369
qs_dispersion_pairpot::d3_cnumber
subroutine, public d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, calculate_forces, debugall)
...
Definition: qs_dispersion_pairpot.F:3449
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition: qs_dispersion_types.F:12
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition: qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition: qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition: string_utilities.F:3190
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13
xtb_coulomb
Calculation of Coulomb contributions in xTB.
Definition: xtb_coulomb.F:12
xtb_coulomb::build_xtb_coulomb
subroutine, public build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, calculate_forces, just_energy)
...
Definition: xtb_coulomb.F:103
xtb_matrices
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition: xtb_matrices.F:15
xtb_matrices::build_xtb_matrices
subroutine, public build_xtb_matrices(qs_env, para_env, calculate_forces)
...
Definition: xtb_matrices.F:134
xtb_matrices::xtb_hab_force
subroutine, public xtb_hab_force(qs_env, p_matrix)
...
Definition: xtb_matrices.F:987
xtb_matrices::build_xtb_ks_matrix
subroutine, public build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Definition: xtb_matrices.F:1521
xtb_parameters
Read xTB parameters.
Definition: xtb_parameters.F:12
xtb_parameters::xtb_set_kab
real(kind=dp) function, public xtb_set_kab(za, zb, xtb_control)
...
Definition: xtb_parameters.F:557
xtb_types
Definition of the xTB parameter types.
Definition: xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175
xtb_matrices::neighbor_atoms_type
Definition: xtb_matrices.F:113