(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
15MODULE 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
125CONTAINS
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
2203END MODULE xtb_matrices
2204
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
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
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
mulliken::ao_charges
Definition mulliken.F:42
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
qs_overlap::create_sab_matrix
Definition qs_overlap.F:81
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::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_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_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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
block_p_types::block_p_type
Definition block_p_types.F:23
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_control_types::xtb_control_type
Definition cp_control_types.F:127
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:55
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pair_potential_types::pair_potential_p_type
Definition pair_potential_types.F:407
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:412
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:379
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_pairpot::dcnum_type
Definition qs_dispersion_pairpot.F:94
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:31
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:95
virial_types::virial_type
Definition virial_types.F:26
xtb_matrices::neighbor_atoms_type
Definition xtb_matrices.F:113
xtb_types::xtb_atom_type
Definition xtb_types.F:42