(git:d18deda)
Loading...
Searching...
No Matches
xtb_potentials.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief xTB (repulsive) pair potentials
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_potentials
16 USE atomic_kind_types, ONLY: atomic_kind_type,&
17 get_atomic_kind,&
18 get_atomic_kind_set
19 USE atprop_types, ONLY: atprop_type
20 USE cp_control_types, ONLY: dft_control_type,&
21 xtb_control_type
22 USE cp_log_handling, ONLY: cp_get_default_logger,&
23 cp_logger_get_default_io_unit,&
24 cp_logger_type,&
25 cp_to_string
26 USE ewald_environment_types, ONLY: ewald_env_get,&
27 ewald_environment_type
28 USE fparser, ONLY: evalfd,&
29 finalizef
30 USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 section_vals_type,&
32 section_vals_val_get
33 USE kinds, ONLY: default_string_length,&
34 dp
35 USE message_passing, ONLY: mp_para_env_type
36 USE pair_potential, ONLY: init_genpot
37 USE pair_potential_types, ONLY: not_initialized,&
38 pair_potential_p_type,&
39 pair_potential_pp_create,&
40 pair_potential_pp_release,&
41 pair_potential_pp_type,&
42 pair_potential_single_clean,&
43 pair_potential_single_copy,&
44 pair_potential_single_type
45 USE pair_potential_util, ONLY: ener_pot
46 USE particle_types, ONLY: particle_type
47 USE qs_dispersion_cnum, ONLY: dcnum_type
48 USE qs_environment_types, ONLY: get_qs_env,&
49 qs_environment_type
50 USE qs_force_types, ONLY: qs_force_type
51 USE qs_kind_types, ONLY: get_qs_kind,&
52 qs_kind_type
53 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
54 neighbor_list_iterate,&
55 neighbor_list_iterator_create,&
56 neighbor_list_iterator_p_type,&
57 neighbor_list_iterator_release,&
58 neighbor_list_set_p_type
59 USE string_utilities, ONLY: compress,&
60 uppercase
61 USE virial_methods, ONLY: virial_pair_force
62 USE virial_types, ONLY: virial_type
63 USE xtb_types, ONLY: get_xtb_atom_param,&
64 xtb_atom_type
65#include "./base/base_uses.f90"
66
67 IMPLICIT NONE
68
69 TYPE neighbor_atoms_type
70 REAL(kind=dp), DIMENSION(:, :), POINTER :: coord => null()
71 REAL(kind=dp), DIMENSION(:), POINTER :: rab => null()
72 INTEGER, DIMENSION(:), POINTER :: katom => null()
73 END TYPE neighbor_atoms_type
74
75 PRIVATE
76
77 CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'xtb_potentials'
78
79 PUBLIC :: xtb_pp_radius, repulsive_potential, srb_potential
80 PUBLIC :: nonbonded_correction, xb_interaction
81 PUBLIC :: neighbor_atoms_type
82
83CONTAINS
84
85! **************************************************************************************************
86!> \brief ...
87!> \param qs_env ...
88!> \param erep ...
89!> \param kf ...
90!> \param enscale ...
91!> \param calculate_forces ...
92! **************************************************************************************************
93 SUBROUTINE repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
94 TYPE(qs_environment_type), POINTER :: qs_env
95 REAL(dp), INTENT(INOUT) :: erep
96 REAL(dp), INTENT(IN) :: kf, enscale
97 LOGICAL, INTENT(IN) :: calculate_forces
98
99 CHARACTER(len=*), PARAMETER :: routinen = 'repulsive_potential'
100
101 INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
102 jatom, jkind, za, zb
103 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
104 INTEGER, DIMENSION(3) :: cell
105 LOGICAL :: defined, use_virial
106 REAL(kind=dp) :: alphaa, alphab, den2, den4, derepij, dr, &
107 ena, enb, ens, erepij, f1, sal, &
108 zneffa, zneffb
109 REAL(kind=dp), DIMENSION(3) :: force_rr, rij
110 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
111 TYPE(atprop_type), POINTER :: atprop
112 TYPE(neighbor_list_iterator_p_type), &
113 DIMENSION(:), POINTER :: nl_iterator
114 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
115 POINTER :: sab_xtb_pp
116 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118 TYPE(virial_type), POINTER :: virial
119 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
120
121 CALL timeset(routinen, handle)
122
123 erep = 0._dp
124
125 CALL get_qs_env(qs_env=qs_env, &
126 atomic_kind_set=atomic_kind_set, &
127 qs_kind_set=qs_kind_set, &
128 atprop=atprop, &
129 sab_xtb_pp=sab_xtb_pp)
130 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
131
132 IF (calculate_forces) THEN
133 CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
134 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
135 END IF
136
137 CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
138 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
139 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
140 iatom=iatom, jatom=jatom, r=rij, cell=cell)
141 CALL get_qs_kind(qs_kind_set(ikind), zatom=za, xtb_parameter=xtb_atom_a)
142 CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
143 IF (.NOT. defined) cycle
144 CALL get_qs_kind(qs_kind_set(jkind), zatom=zb, xtb_parameter=xtb_atom_b)
145 CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
146 IF (.NOT. defined) cycle
147
148 dr = sqrt(sum(rij(:)**2))
149 ! repulsive potential
150 IF (dr > 0.001_dp) THEN
151
152 ! atomic parameters
153 CALL get_xtb_atom_param(xtb_atom_a, en=ena, alpha=alphaa, zneff=zneffa)
154 CALL get_xtb_atom_param(xtb_atom_b, en=enb, alpha=alphab, zneff=zneffb)
155
156 ! scaling (not in papers! but in code)
157 den2 = (ena - enb)**2
158 den4 = den2*den2
159 sal = sqrt(alphaa*alphab)
160 ens = 1.0_dp + (0.01_dp*den2 + 0.01_dp*den4)*enscale
161
162 erepij = zneffa*zneffb/dr*exp(-ens*sal*dr**kf)
163 erep = erep + erepij
164 IF (atprop%energy) THEN
165 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
166 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
167 END IF
168 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
169 derepij = -(1.0_dp/dr + ens*sal*kf*dr**(kf - 1.0_dp))*erepij
170 force_rr(1) = derepij*rij(1)/dr
171 force_rr(2) = derepij*rij(2)/dr
172 force_rr(3) = derepij*rij(3)/dr
173 atom_a = atom_of_kind(iatom)
174 atom_b = atom_of_kind(jatom)
175 force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
176 force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
177 IF (use_virial) THEN
178 f1 = 1.0_dp
179 IF (iatom == jatom) f1 = 0.5_dp
180 CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
181 END IF
182 END IF
183 END IF
184
185 END DO
186 CALL neighbor_list_iterator_release(nl_iterator)
187
188 CALL timestop(handle)
189
190 END SUBROUTINE repulsive_potential
191
192! **************************************************************************************************
193!> \brief ...
194!> \param qs_env ...
195!> \param esrb ...
196!> \param calculate_forces ...
197!> \param xtb_control ...
198!> \param cnumbers ...
199!> \param dcnum ...
200! **************************************************************************************************
201 SUBROUTINE srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
202 TYPE(qs_environment_type), POINTER :: qs_env
203 REAL(dp), INTENT(INOUT) :: esrb
204 LOGICAL, INTENT(IN) :: calculate_forces
205 TYPE(xtb_control_type), POINTER :: xtb_control
206 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: cnumbers
207 TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
208
209 CHARACTER(len=*), PARAMETER :: routinen = 'srb_potential'
210 REAL(kind=dp), DIMENSION(5:9), PARAMETER :: &
211 cnfac = (/0.05646607_dp, 0.10514203_dp, 0.09753494_dp, 0.30470380_dp, 0.23261783_dp/), &
212 ensrb = (/2.20568300_dp, 2.49640820_dp, 2.81007174_dp, 4.51078438_dp, 4.67476223_dp/), &
213 r0srb = (/1.35974851_dp, 0.98310699_dp, 0.98423007_dp, 0.76716063_dp, 1.06139799_dp/)
214
215 INTEGER :: atom_a, atom_b, atom_c, handle, i, &
216 iatom, ikind, jatom, jkind, katom, &
217 kkind, za, zb
218 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
219 INTEGER, DIMENSION(3) :: cell
220 LOGICAL :: defined, use_virial
221 REAL(kind=dp) :: c1srb, c2srb, den1, den2, desrbij, dr, &
222 dr0, drk, enta, entb, esrbij, etasrb, &
223 f1, fhua, fhub, gscal, ksrb, rra0, &
224 rrb0, shift
225 REAL(kind=dp), DIMENSION(3) :: fdik, force_rr, rij, rik
226 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
227 TYPE(atprop_type), POINTER :: atprop
228 TYPE(neighbor_list_iterator_p_type), &
229 DIMENSION(:), POINTER :: nl_iterator
230 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
231 POINTER :: sab_xtb_pp
232 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
233 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
234 TYPE(virial_type), POINTER :: virial
235 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
236
237 CALL timeset(routinen, handle)
238
239 esrb = 0._dp
240
241 CALL get_qs_env(qs_env=qs_env, &
242 atomic_kind_set=atomic_kind_set, &
243 qs_kind_set=qs_kind_set, &
244 atprop=atprop, &
245 sab_xtb_pp=sab_xtb_pp)
246 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
247 kind_of=kind_of)
248
249 IF (calculate_forces) THEN
250 CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
251 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
252 END IF
253
254 ! SRB parameters
255 ksrb = xtb_control%ksrb
256 etasrb = xtb_control%esrb
257 c1srb = xtb_control%c1srb*0.01_dp
258 c2srb = xtb_control%c2srb*0.01_dp
259 gscal = xtb_control%gscal
260 shift = xtb_control%shift
261
262 CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
263 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
264 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
265 iatom=iatom, jatom=jatom, r=rij, cell=cell)
266 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
267 CALL get_xtb_atom_param(xtb_atom_a, z=za, electronegativity=enta, defined=defined)
268 IF (.NOT. defined) cycle
269 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
270 CALL get_xtb_atom_param(xtb_atom_b, z=zb, electronegativity=entb, defined=defined)
271 IF (.NOT. defined) cycle
272
273 dr = sqrt(sum(rij(:)**2))
274 ! short-ranged correction term
275 IF (dr > 0.001_dp) THEN
276 IF (za >= 5 .AND. za <= 9 .AND. zb >= 5 .AND. zb <= 9 .AND. za /= zb) THEN
277 rra0 = r0srb(za) + cnfac(za)*cnumbers(iatom) + shift
278 rrb0 = r0srb(zb) + cnfac(zb)*cnumbers(jatom) + shift
279 den1 = abs(ensrb(za) - ensrb(zb))
280 dr0 = (rra0 + rrb0)*(1._dp - c1srb*den1 - c2srb*den1*den1)
281 den2 = (enta - entb)**2
282 esrbij = ksrb*exp(-etasrb*(1._dp + gscal*den2)*(dr - dr0)**2)
283 esrb = esrb + esrbij
284 IF (atprop%energy) THEN
285 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*esrbij
286 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*esrbij
287 END IF
288 IF (calculate_forces) THEN
289 desrbij = 2.0_dp*esrbij*(-etasrb*(1._dp + gscal*den2)*(dr - dr0))
290 force_rr(1) = desrbij*rij(1)/dr
291 force_rr(2) = desrbij*rij(2)/dr
292 force_rr(3) = desrbij*rij(3)/dr
293 atom_a = atom_of_kind(iatom)
294 atom_b = atom_of_kind(jatom)
295 force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
296 force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
297 IF (use_virial) THEN
298 f1 = 1.0_dp
299 IF (iatom == jatom) f1 = 0.5_dp
300 CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
301 END IF
302 ! coordination number derivatives
303 ! iatom
304 fhua = -desrbij*cnfac(za)*(1._dp - c1srb*den1 - c2srb*den1*den1)
305 DO i = 1, dcnum(iatom)%neighbors
306 katom = dcnum(iatom)%nlist(i)
307 kkind = kind_of(katom)
308 atom_c = atom_of_kind(katom)
309 rik = dcnum(iatom)%rik(:, i)
310 drk = sqrt(sum(rik(:)**2))
311 IF (drk > 1.e-3_dp) THEN
312 fdik(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
313 force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fdik(:)
314 force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
315 IF (use_virial) THEN
316 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
317 END IF
318 END IF
319 END DO
320 ! jatom
321 fhub = -desrbij*cnfac(zb)*(1._dp - c1srb*den1 - c2srb*den1*den1)
322 DO i = 1, dcnum(jatom)%neighbors
323 katom = dcnum(jatom)%nlist(i)
324 kkind = kind_of(katom)
325 atom_c = atom_of_kind(katom)
326 rik = dcnum(jatom)%rik(:, i)
327 drk = sqrt(sum(rik(:)**2))
328 IF (drk > 1.e-3_dp) THEN
329 fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
330 force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) - fdik(:)
331 force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
332 IF (use_virial) THEN
333 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
334 END IF
335 END IF
336 END DO
337 END IF
338 END IF
339 END IF
340
341 END DO
342 CALL neighbor_list_iterator_release(nl_iterator)
343
344 CALL timestop(handle)
345
346 END SUBROUTINE srb_potential
347
348! **************************************************************************************************
349!> \brief ...
350!> \param qs_kind_set ...
351!> \param ppradius ...
352!> \param eps_pair ...
353!> \param kfparam ...
354! **************************************************************************************************
355 SUBROUTINE xtb_pp_radius(qs_kind_set, ppradius, eps_pair, kfparam)
356 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
357 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: ppradius
358 REAL(kind=dp), INTENT(IN) :: eps_pair, kfparam
359
360 INTEGER :: ikind, ir, jkind, nkind
361 LOGICAL :: defa, defb
362 REAL(kind=dp) :: alphaa, alphab, erep, rab, rab0, rcova, &
363 rcovb, saa, zneffa, zneffb
364 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
365
366 ppradius = 0.0_dp
367 nkind = SIZE(ppradius, 1)
368 DO ikind = 1, nkind
369 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
370 CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova, alpha=alphaa, zneff=zneffa, defined=defa)
371 IF (.NOT. defa) cycle
372 DO jkind = ikind, nkind
373 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
374 CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb, alpha=alphab, zneff=zneffb, defined=defb)
375 IF (.NOT. defb) cycle
376 rab = 0.0_dp
377 DO ir = 1, 24
378 rab = rab + 1.0_dp
379 saa = sqrt(alphaa*alphab)
380 erep = zneffa*zneffb/rab*exp(-saa*rab**kfparam)
381 IF (erep < eps_pair) EXIT
382 END DO
383 rab0 = rcova + rcovb
384 rab = max(rab, rab0 + 2.0_dp)
385 ppradius(ikind, jkind) = rab
386 ppradius(jkind, ikind) = ppradius(ikind, jkind)
387 END DO
388 END DO
389
390 END SUBROUTINE xtb_pp_radius
391
392! **************************************************************************************************
393!> \brief ...
394!> \param qs_env ...
395!> \param exb ...
396!> \param calculate_forces ...
397! **************************************************************************************************
398 SUBROUTINE xb_interaction(qs_env, exb, calculate_forces)
399
400 TYPE(qs_environment_type), POINTER :: qs_env
401 REAL(kind=dp), INTENT(INOUT) :: exb
402 LOGICAL, INTENT(IN) :: calculate_forces
403
404 CHARACTER(LEN=*), PARAMETER :: routinen = 'xb_interaction'
405
406 INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
407 jatom, jkind, na, natom, nkind, zat
408 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
409 INTEGER, DIMENSION(3) :: cell
410 LOGICAL :: defined, use_virial
411 REAL(kind=dp) :: dr, kx2, kxr, rcova, rcovb
412 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kx
413 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rcab
414 REAL(kind=dp), DIMENSION(3) :: rij
415 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
416 TYPE(atprop_type), POINTER :: atprop
417 TYPE(dft_control_type), POINTER :: dft_control
418 TYPE(mp_para_env_type), POINTER :: para_env
419 TYPE(neighbor_atoms_type), ALLOCATABLE, &
420 DIMENSION(:) :: neighbor_atoms
421 TYPE(neighbor_list_iterator_p_type), &
422 DIMENSION(:), POINTER :: nl_iterator
423 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
424 POINTER :: sab_xb, sab_xtb_pp
425 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
427 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
428 TYPE(virial_type), POINTER :: virial
429 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
430 TYPE(xtb_control_type), POINTER :: xtb_control
431
432 CALL timeset(routinen, handle)
433
434 CALL get_qs_env(qs_env=qs_env, &
435 atomic_kind_set=atomic_kind_set, &
436 qs_kind_set=qs_kind_set, &
437 para_env=para_env, &
438 atprop=atprop, &
439 dft_control=dft_control, &
440 sab_xb=sab_xb, &
441 sab_xtb_pp=sab_xtb_pp)
442
443 nkind = SIZE(atomic_kind_set)
444 xtb_control => dft_control%qs_control%xtb_control
445
446 ! global parameters
447 kxr = xtb_control%kxr
448 kx2 = xtb_control%kx2
449
450 NULLIFY (particle_set)
451 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
452 natom = SIZE(particle_set)
453 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
454 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
455
456 IF (calculate_forces) THEN
457 CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
458 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
459 END IF
460
461 ! list of neighbor atoms for XB term
462 ALLOCATE (neighbor_atoms(nkind))
463 DO ikind = 1, nkind
464 NULLIFY (neighbor_atoms(ikind)%coord)
465 NULLIFY (neighbor_atoms(ikind)%rab)
466 NULLIFY (neighbor_atoms(ikind)%katom)
467 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat, natom=na)
468 IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
469 ALLOCATE (neighbor_atoms(ikind)%coord(3, na))
470 neighbor_atoms(ikind)%coord(1:3, 1:na) = 0.0_dp
471 ALLOCATE (neighbor_atoms(ikind)%rab(na))
472 neighbor_atoms(ikind)%rab(1:na) = huge(0.0_dp)
473 ALLOCATE (neighbor_atoms(ikind)%katom(na))
474 neighbor_atoms(ikind)%katom(1:na) = 0
475 END IF
476 END DO
477 ! kx parameters
478 ALLOCATE (kx(nkind))
479 DO ikind = 1, nkind
480 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
481 CALL get_xtb_atom_param(xtb_atom_a, kx=kx(ikind))
482 END DO
483 !
484 ALLOCATE (rcab(nkind, nkind))
485 DO ikind = 1, nkind
486 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
487 CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova)
488 DO jkind = 1, nkind
489 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
490 CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb)
491 rcab(ikind, jkind) = kxr*(rcova + rcovb)
492 END DO
493 END DO
494
495 CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
496 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
497 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
498 iatom=iatom, jatom=jatom, r=rij, cell=cell)
499 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
500 CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
501 IF (.NOT. defined) cycle
502 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
503 CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
504 IF (.NOT. defined) cycle
505
506 dr = sqrt(sum(rij(:)**2))
507
508 ! neighbor atom for XB term
509 IF (dr > 1.e-3_dp) THEN
510 IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
511 atom_a = atom_of_kind(iatom)
512 IF (dr < neighbor_atoms(ikind)%rab(atom_a)) THEN
513 neighbor_atoms(ikind)%rab(atom_a) = dr
514 neighbor_atoms(ikind)%coord(1:3, atom_a) = rij(1:3)
515 neighbor_atoms(ikind)%katom(atom_a) = jatom
516 END IF
517 END IF
518 IF (ASSOCIATED(neighbor_atoms(jkind)%rab)) THEN
519 atom_b = atom_of_kind(jatom)
520 IF (dr < neighbor_atoms(jkind)%rab(atom_b)) THEN
521 neighbor_atoms(jkind)%rab(atom_b) = dr
522 neighbor_atoms(jkind)%coord(1:3, atom_b) = -rij(1:3)
523 neighbor_atoms(jkind)%katom(atom_b) = iatom
524 END IF
525 END IF
526 END IF
527
528 END DO
529 CALL neighbor_list_iterator_release(nl_iterator)
530
531 exb = 0.0_dp
532 CALL xb_neighbors(neighbor_atoms, para_env)
533 CALL xb_energy(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
534 calculate_forces, use_virial, force, virial, atprop)
535
536 DO ikind = 1, nkind
537 IF (ASSOCIATED(neighbor_atoms(ikind)%coord)) THEN
538 DEALLOCATE (neighbor_atoms(ikind)%coord)
539 END IF
540 IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
541 DEALLOCATE (neighbor_atoms(ikind)%rab)
542 END IF
543 IF (ASSOCIATED(neighbor_atoms(ikind)%katom)) THEN
544 DEALLOCATE (neighbor_atoms(ikind)%katom)
545 END IF
546 END DO
547 DEALLOCATE (neighbor_atoms)
548 DEALLOCATE (kx, rcab)
549
550 CALL timestop(handle)
551
552 END SUBROUTINE xb_interaction
553
554! **************************************************************************************************
555!> \brief Distributes the neighbor atom information to all processors
556!>
557!> \param neighbor_atoms ...
558!> \param para_env ...
559!> \par History
560!> 1.2019 JGH
561!> \version 1.1
562! **************************************************************************************************
563 SUBROUTINE xb_neighbors(neighbor_atoms, para_env)
564 TYPE(neighbor_atoms_type), DIMENSION(:), &
565 INTENT(INOUT) :: neighbor_atoms
566 TYPE(mp_para_env_type), POINTER :: para_env
567
568 INTEGER :: iatom, ikind, natom, nkind
569 INTEGER, ALLOCATABLE, DIMENSION(:) :: matom
570 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dmloc
571 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
572
573 nkind = SIZE(neighbor_atoms)
574 DO ikind = 1, nkind
575 IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
576 natom = SIZE(neighbor_atoms(ikind)%rab)
577 ALLOCATE (dmloc(2*natom), matom(natom), coord(3, natom))
578 dmloc = 0.0_dp
579 DO iatom = 1, natom
580 dmloc(2*iatom - 1) = neighbor_atoms(ikind)%rab(iatom)
581 dmloc(2*iatom) = real(para_env%mepos, kind=dp)
582 END DO
583 CALL para_env%minloc(dmloc)
584 coord = 0.0_dp
585 matom = 0
586 DO iatom = 1, natom
587 neighbor_atoms(ikind)%rab(iatom) = dmloc(2*iatom - 1)
588 IF (nint(dmloc(2*iatom)) == para_env%mepos) THEN
589 coord(1:3, iatom) = neighbor_atoms(ikind)%coord(1:3, iatom)
590 matom(iatom) = neighbor_atoms(ikind)%katom(iatom)
591 END IF
592 END DO
593 CALL para_env%sum(coord)
594 neighbor_atoms(ikind)%coord(1:3, :) = coord(1:3, :)
595 CALL para_env%sum(matom)
596 neighbor_atoms(ikind)%katom(:) = matom(:)
597 DEALLOCATE (dmloc, matom, coord)
598 END IF
599 END DO
600
601 END SUBROUTINE xb_neighbors
602
603! **************************************************************************************************
604!> \brief Computes a correction for nonbonded interactions based on a generic potential
605!>
606!> \param enonbonded energy contribution
607!> \param force ...
608!> \param qs_env ...
609!> \param xtb_control ...
610!> \param sab_xtb_nonbond ...
611!> \param atomic_kind_set ...
612!> \param calculate_forces ...
613!> \param use_virial ...
614!> \param virial ...
615!> \param atprop ...
616!> \param atom_of_kind ..
617!> \par History
618!> 12.2018 JGH
619!> \version 1.1
620! **************************************************************************************************
621 SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
622 atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
623 REAL(dp), INTENT(INOUT) :: enonbonded
624 TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
625 POINTER :: force
626 TYPE(qs_environment_type), POINTER :: qs_env
627 TYPE(xtb_control_type), POINTER :: xtb_control
628 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
629 INTENT(IN), POINTER :: sab_xtb_nonbond
630 TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
631 POINTER :: atomic_kind_set
632 LOGICAL, INTENT(IN) :: calculate_forces, use_virial
633 TYPE(virial_type), INTENT(IN), POINTER :: virial
634 TYPE(atprop_type), INTENT(IN), POINTER :: atprop
635 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind
636
637 CHARACTER(len=*), PARAMETER :: routinen = 'nonbonded_correction'
638
639 CHARACTER(LEN=default_string_length) :: def_error, this_error
640 INTEGER :: atom_i, atom_j, handle, iatom, ikind, &
641 jatom, jkind, kk, ntype
642 INTEGER, DIMENSION(3) :: cell
643 LOGICAL :: do_ewald
644 REAL(kind=dp) :: dedf, dr, dx, energy_cutoff, err, fval, &
645 lerr, rcut
646 REAL(kind=dp), DIMENSION(3) :: fij, rij
647 TYPE(ewald_environment_type), POINTER :: ewald_env
648 TYPE(neighbor_list_iterator_p_type), &
649 DIMENSION(:), POINTER :: nl_iterator
650 TYPE(pair_potential_p_type), POINTER :: nonbonded
651 TYPE(pair_potential_pp_type), POINTER :: potparm
652 TYPE(pair_potential_single_type), POINTER :: pot
653 TYPE(section_vals_type), POINTER :: nonbonded_section
654
655 CALL timeset(routinen, handle)
656
657 NULLIFY (nonbonded)
658 NULLIFY (potparm)
659 NULLIFY (ewald_env)
660 nonbonded => xtb_control%nonbonded
661 do_ewald = xtb_control%do_ewald
662 CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
663
664 ntype = SIZE(atomic_kind_set)
665 CALL pair_potential_pp_create(potparm, ntype)
666 !Assign input and potential info to potparm_nonbond
667 CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
668 !Initialize genetic potential
669 CALL init_genpot(potparm, ntype)
670
671 NULLIFY (pot)
672 enonbonded = 0._dp
673 energy_cutoff = 0._dp
674
675 CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond)
676 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
677 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
678 iatom=iatom, jatom=jatom, r=rij, cell=cell)
679 pot => potparm%pot(ikind, jkind)%pot
680 dr = sqrt(rij(1)**2 + rij(2)**2 + rij(3)**2)
681 rcut = sqrt(pot%rcutsq)
682 IF (dr <= rcut .AND. dr > 1.e-3_dp) THEN
683 fval = 1.0_dp
684 IF (ikind == jkind) fval = 0.5_dp
685 ! splines not implemented
686 enonbonded = enonbonded + fval*ener_pot(pot, dr, energy_cutoff)
687 IF (atprop%energy) THEN
688 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
689 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
690 END IF
691 END IF
692
693 IF (calculate_forces) THEN
694
695 kk = SIZE(pot%type)
696 IF (kk /= 1) THEN
697 CALL cp_warn(__location__, "Generic potential with type > 1 not implemented.")
698 cpabort("pot type")
699 END IF
700 ! rmin and rmax and rcut
701 IF ((pot%set(kk)%rmin /= not_initialized) .AND. (dr < pot%set(kk)%rmin)) cycle
702 ! An upper boundary for the potential definition was defined
703 IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) cycle
704 ! If within limits let's compute the potential...
705 IF (dr <= rcut .AND. dr > 1.e-3_dp) THEN
706
707 NULLIFY (nonbonded_section)
708 nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED")
709 CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx)
710 CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr)
711
712 dedf = fval*evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err)
713 IF (abs(err) > lerr) THEN
714 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
715 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
716 CALL compress(this_error, .true.)
717 CALL compress(def_error, .true.)
718 CALL cp_warn(__location__, &
719 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
720 ' Error '//trim(this_error)//' in computing numerical derivatives larger then'// &
721 trim(def_error)//' .')
722 END IF
723
724 atom_i = atom_of_kind(iatom)
725 atom_j = atom_of_kind(jatom)
726
727 fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
728
729 force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
730 force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:)
731
732 IF (use_virial) THEN
733 CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
734 END IF
735
736 END IF
737 END IF
738 NULLIFY (pot)
739 END DO
740 CALL neighbor_list_iterator_release(nl_iterator)
741 CALL finalizef()
742
743 IF (ASSOCIATED(potparm)) THEN
744 CALL pair_potential_pp_release(potparm)
745 END IF
746
747 CALL timestop(handle)
748
749 END SUBROUTINE nonbonded_correction
750
751! **************************************************************************************************
752!> \brief ...
753!> \param atomic_kind_set ...
754!> \param nonbonded ...
755!> \param potparm ...
756!> \param ewald_env ...
757!> \param do_ewald ...
758! **************************************************************************************************
759 SUBROUTINE force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
760
761 ! routine based on force_field_pack_nonbond
762 TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
763 POINTER :: atomic_kind_set
764 TYPE(pair_potential_p_type), INTENT(IN), POINTER :: nonbonded
765 TYPE(pair_potential_pp_type), INTENT(INOUT), &
766 POINTER :: potparm
767 TYPE(ewald_environment_type), INTENT(IN), POINTER :: ewald_env
768 LOGICAL, INTENT(IN) :: do_ewald
769
770 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
771 name_atm_b, name_atm_b_local
772 INTEGER :: ikind, ingp, iw, jkind
773 LOGICAL :: found
774 REAL(kind=dp) :: ewald_rcut
775 TYPE(atomic_kind_type), POINTER :: atomic_kind
776 TYPE(cp_logger_type), POINTER :: logger
777 TYPE(pair_potential_single_type), POINTER :: pot
778
779 NULLIFY (pot, logger)
780
781 logger => cp_get_default_logger()
782 iw = cp_logger_get_default_io_unit(logger)
783
784 DO ikind = 1, SIZE(atomic_kind_set)
785 atomic_kind => atomic_kind_set(ikind)
786 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
787 DO jkind = ikind, SIZE(atomic_kind_set)
788 atomic_kind => atomic_kind_set(jkind)
789 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
790 found = .false.
791 name_atm_a = name_atm_a_local
792 name_atm_b = name_atm_b_local
793 CALL uppercase(name_atm_a)
794 CALL uppercase(name_atm_b)
795 pot => potparm%pot(ikind, jkind)%pot
796 IF (ASSOCIATED(nonbonded)) THEN
797 DO ingp = 1, SIZE(nonbonded%pot)
798 IF ((trim(nonbonded%pot(ingp)%pot%at1) == "*") .OR. &
799 (trim(nonbonded%pot(ingp)%pot%at2) == "*")) cycle
800
801 !IF (iw > 0) WRITE (iw, *) "TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), &
802 ! " with ", TRIM(nonbonded%pot(ingp)%pot%at1), &
803 ! TRIM(nonbonded%pot(ingp)%pot%at2)
804 IF ((((name_atm_a) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
805 ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. &
806 (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
807 ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN
808 CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot)
809 ! multiple potential not implemented, simply overwriting
810 IF (found) &
811 CALL cp_warn(__location__, &
812 "Multiple NONBONDED declaration: "//trim(name_atm_a)// &
813 " and "//trim(name_atm_b)//" OVERWRITING! ")
814 !IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
815 found = .true.
816 END IF
817 END DO
818 END IF
819 IF (.NOT. found) THEN
820 CALL pair_potential_single_clean(pot)
821 !IF (iw > 0) WRITE (iw, *) " NOTFOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
822 END IF
823 END DO !jkind
824 END DO !ikind
825
826 ! Cutoff is defined always as the maximum between the FF and Ewald
827 IF (do_ewald) THEN
828 CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
829 pot%rcutsq = max(pot%rcutsq, ewald_rcut*ewald_rcut)
830 !IF (iw > 0) WRITE (iw, *) " RCUT ", SQRT(pot%rcutsq), ewald_rcut
831 END IF
832
833 END SUBROUTINE force_field_pack_nonbond_pot_correction
834
835! **************************************************************************************************
836!> \brief Computes the interaction term between Br/I/At and donor atoms
837!>
838!> \param exb ...
839!> \param neighbor_atoms ...
840!> \param atom_of_kind ...
841!> \param kind_of ...
842!> \param sab_xb ...
843!> \param kx ...
844!> \param kx2 ...
845!> \param rcab ...
846!> \param calculate_forces ...
847!> \param use_virial ...
848!> \param force ...
849!> \param virial ...
850!> \param atprop ...
851!> \par History
852!> 12.2018 JGH
853!> \version 1.1
854! **************************************************************************************************
855 SUBROUTINE xb_energy(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
856 calculate_forces, use_virial, force, virial, atprop)
857 REAL(dp), INTENT(INOUT) :: exb
858 TYPE(neighbor_atoms_type), DIMENSION(:), &
859 INTENT(IN) :: neighbor_atoms
860 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
861 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
862 POINTER :: sab_xb
863 REAL(dp), DIMENSION(:), INTENT(IN) :: kx
864 REAL(dp), INTENT(IN) :: kx2
865 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rcab
866 LOGICAL, INTENT(IN) :: calculate_forces, use_virial
867 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
868 TYPE(virial_type), POINTER :: virial
869 TYPE(atprop_type), POINTER :: atprop
870
871 INTEGER :: atom_a, atom_b, atom_c, iatom, ikind, &
872 jatom, jkind, katom, kkind
873 INTEGER, DIMENSION(3) :: cell
874 REAL(kind=dp) :: alp, aterm, cosa, daterm, ddab, ddax, &
875 ddbx, ddr, ddr12, ddr6, deval, dr, &
876 drab, drax, drbx, eval, xy
877 REAL(kind=dp), DIMENSION(3) :: fia, fij, fja, ria, rij, rja
878 TYPE(neighbor_list_iterator_p_type), &
879 DIMENSION(:), POINTER :: nl_iterator
880
881 ! exonent in angular term
882 alp = 6.0_dp
883 ! loop over all atom pairs
884 CALL neighbor_list_iterator_create(nl_iterator, sab_xb)
885 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
886 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
887 iatom=iatom, jatom=jatom, r=rij, cell=cell)
888 ! ikind, iatom : Halogen
889 ! jkind, jatom : Donor
890 atom_a = atom_of_kind(iatom)
891 katom = neighbor_atoms(ikind)%katom(atom_a)
892 IF (katom == 0) cycle
893 dr = sqrt(rij(1)**2 + rij(2)**2 + rij(3)**2)
894 ddr = rcab(ikind, jkind)/dr
895 ddr6 = ddr**6
896 ddr12 = ddr6*ddr6
897 eval = kx(ikind)*(ddr12 - kx2*ddr6)/(1.0_dp + ddr12)
898 ! angle
899 ria(1:3) = neighbor_atoms(ikind)%coord(1:3, atom_a)
900 rja(1:3) = rij(1:3) - ria(1:3)
901 drax = ria(1)**2 + ria(2)**2 + ria(3)**2
902 drbx = dr*dr
903 drab = rja(1)**2 + rja(2)**2 + rja(3)**2
904 xy = sqrt(drbx*drax)
905 ! cos angle B-X-A
906 cosa = (drbx + drax - drab)/xy
907 aterm = (0.5_dp - 0.25_dp*cosa)**alp
908 !
909 exb = exb + aterm*eval
910 IF (atprop%energy) THEN
911 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*aterm*eval
912 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*aterm*eval
913 END IF
914 !
915 IF (calculate_forces) THEN
916 kkind = kind_of(katom)
917 atom_b = atom_of_kind(jatom)
918 atom_c = atom_of_kind(katom)
919 !
920 deval = 6.0_dp*kx(ikind)*ddr6*(kx2*ddr12 + 2.0_dp*ddr6 - kx2)/(1.0_dp + ddr12)**2
921 deval = -rcab(ikind, jkind)*deval/(dr*dr)/ddr
922 fij(1:3) = aterm*deval*rij(1:3)/dr
923 force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:)
924 force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:)
925 IF (use_virial) THEN
926 CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
927 END IF
928 !
929 fij(1:3) = 0.0_dp
930 fia(1:3) = 0.0_dp
931 fja(1:3) = 0.0_dp
932 daterm = -0.25_dp*alp*(0.5_dp - 0.25_dp*cosa)**(alp - 1.0_dp)
933 ddbx = 0.5_dp*(drab - drax + drbx)/xy/drbx
934 ddax = 0.5_dp*(drab + drax - drbx)/xy/drax
935 ddab = -1._dp/xy
936 fij(1:3) = 2.0_dp*daterm*ddbx*rij(1:3)*eval
937 fia(1:3) = 2.0_dp*daterm*ddax*ria(1:3)*eval
938 fja(1:3) = 2.0_dp*daterm*ddab*rja(1:3)*eval
939 force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:) - fia(:)
940 force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:) + fja(:)
941 force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fia(:) - fja(:)
942 IF (use_virial) THEN
943 CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
944 CALL virial_pair_force(virial%pv_virial, -1._dp, fia, ria)
945 CALL virial_pair_force(virial%pv_virial, -1._dp, fja, rja)
946 END IF
947 END IF
948 END DO
949 CALL neighbor_list_iterator_release(nl_iterator)
950
951 END SUBROUTINE xb_energy
952
953END MODULE xtb_potentials
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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
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_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
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:124
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:731
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:1047
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
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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:973
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:704
pair_potential_types::pair_potential_single_clean
subroutine, public pair_potential_single_clean(potparm)
Cleans the potential parameter type.
Definition pair_potential_types.F:655
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:942
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_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_force_types
Definition qs_force_types.F:13
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, zatom, 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_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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
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
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:3191
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
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:142
virial_types
Definition virial_types.F:13
xtb_potentials
xTB (repulsive) pair potentials Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
Definition xtb_potentials.F:15
xtb_potentials::nonbonded_correction
subroutine, public nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
Computes a correction for nonbonded interactions based on a generic potential.
Definition xtb_potentials.F:623
xtb_potentials::srb_potential
subroutine, public srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
...
Definition xtb_potentials.F:202
xtb_potentials::xb_interaction
subroutine, public xb_interaction(qs_env, exb, calculate_forces)
...
Definition xtb_potentials.F:399
xtb_potentials::xtb_pp_radius
subroutine, public xtb_pp_radius(qs_kind_set, ppradius, eps_pair, kfparam)
...
Definition xtb_potentials.F:356
xtb_potentials::repulsive_potential
subroutine, public repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
...
Definition xtb_potentials.F:94
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, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
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
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::xtb_control_type
Definition cp_control_types.F:138
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:56
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
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:411
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:416
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:383
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
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_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
virial_types::virial_type
Definition virial_types.F:26
xtb_potentials::neighbor_atoms_type
Definition xtb_potentials.F:69
xtb_types::xtb_atom_type
Definition xtb_types.F:42