Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief ...
87 : !> \param qs_env ...
88 : !> \param erep ...
89 : !> \param kf ...
90 : !> \param enscale ...
91 : !> \param calculate_forces ...
92 : ! **************************************************************************************************
93 5084 : 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 5084 : 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 5084 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
111 : TYPE(atprop_type), POINTER :: atprop
112 : TYPE(neighbor_list_iterator_p_type), &
113 5084 : DIMENSION(:), POINTER :: nl_iterator
114 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
115 5084 : POINTER :: sab_xtb_pp
116 5084 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 5084 : 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 5084 : CALL timeset(routineN, handle)
122 :
123 5084 : 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 5084 : sab_xtb_pp=sab_xtb_pp)
130 5084 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
131 :
132 5084 : IF (calculate_forces) THEN
133 486 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
134 486 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
135 : END IF
136 :
137 5084 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
138 99614 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
139 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
140 94530 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
141 94530 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za, xtb_parameter=xtb_atom_a)
142 94530 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
143 94530 : IF (.NOT. defined) CYCLE
144 94530 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb, xtb_parameter=xtb_atom_b)
145 94530 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
146 94530 : IF (.NOT. defined) CYCLE
147 :
148 378120 : dr = SQRT(SUM(rij(:)**2))
149 : ! repulsive potential
150 99614 : IF (dr > 0.001_dp) THEN
151 :
152 : ! atomic parameters
153 75705 : CALL get_xtb_atom_param(xtb_atom_a, en=ena, alpha=alphaa, zneff=zneffa)
154 75705 : CALL get_xtb_atom_param(xtb_atom_b, en=enb, alpha=alphab, zneff=zneffb)
155 :
156 : ! scaling (not in papers! but in code)
157 75705 : den2 = (ena - enb)**2
158 75705 : den4 = den2*den2
159 75705 : sal = SQRT(alphaa*alphab)
160 75705 : ens = 1.0_dp + (0.01_dp*den2 + 0.01_dp*den4)*enscale
161 :
162 75705 : erepij = zneffa*zneffb/dr*EXP(-ens*sal*dr**kf)
163 75705 : erep = erep + erepij
164 75705 : IF (atprop%energy) THEN
165 3458 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
166 3458 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
167 : END IF
168 75705 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
169 7100 : derepij = -(1.0_dp/dr + ens*sal*kf*dr**(kf - 1.0_dp))*erepij
170 7100 : force_rr(1) = derepij*rij(1)/dr
171 7100 : force_rr(2) = derepij*rij(2)/dr
172 7100 : force_rr(3) = derepij*rij(3)/dr
173 7100 : atom_a = atom_of_kind(iatom)
174 7100 : atom_b = atom_of_kind(jatom)
175 28400 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
176 28400 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
177 7100 : IF (use_virial) THEN
178 3329 : f1 = 1.0_dp
179 3329 : IF (iatom == jatom) f1 = 0.5_dp
180 3329 : 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 5084 : CALL neighbor_list_iterator_release(nl_iterator)
187 :
188 5084 : CALL timestop(handle)
189 :
190 10168 : 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 1646 : 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 1646 : 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 1646 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
227 : TYPE(atprop_type), POINTER :: atprop
228 : TYPE(neighbor_list_iterator_p_type), &
229 1646 : DIMENSION(:), POINTER :: nl_iterator
230 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
231 1646 : POINTER :: sab_xtb_pp
232 1646 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
233 1646 : 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 1646 : CALL timeset(routineN, handle)
238 :
239 1646 : 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 1646 : 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 1646 : kind_of=kind_of)
248 :
249 1646 : IF (calculate_forces) THEN
250 42 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
251 42 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
252 : END IF
253 :
254 : ! SRB parameters
255 1646 : ksrb = xtb_control%ksrb
256 1646 : etasrb = xtb_control%esrb
257 1646 : c1srb = xtb_control%c1srb*0.01_dp
258 1646 : c2srb = xtb_control%c2srb*0.01_dp
259 1646 : gscal = xtb_control%gscal
260 1646 : shift = xtb_control%shift
261 :
262 1646 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
263 26770 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
264 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
265 25124 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
266 25124 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
267 25124 : CALL get_xtb_atom_param(xtb_atom_a, z=za, electronegativity=enta, defined=defined)
268 25124 : IF (.NOT. defined) CYCLE
269 25124 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
270 25124 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, electronegativity=entb, defined=defined)
271 25124 : IF (.NOT. defined) CYCLE
272 :
273 100496 : dr = SQRT(SUM(rij(:)**2))
274 : ! short-ranged correction term
275 26770 : IF (dr > 0.001_dp) THEN
276 21249 : IF (za >= 5 .AND. za <= 9 .AND. zb >= 5 .AND. zb <= 9 .AND. za /= zb) THEN
277 539 : rra0 = r0srb(za) + cnfac(za)*cnumbers(iatom) + shift
278 539 : rrb0 = r0srb(zb) + cnfac(zb)*cnumbers(jatom) + shift
279 539 : den1 = ABS(ensrb(za) - ensrb(zb))
280 539 : dr0 = (rra0 + rrb0)*(1._dp - c1srb*den1 - c2srb*den1*den1)
281 539 : den2 = (enta - entb)**2
282 539 : esrbij = ksrb*EXP(-etasrb*(1._dp + gscal*den2)*(dr - dr0)**2)
283 539 : esrb = esrb + esrbij
284 539 : IF (atprop%energy) THEN
285 0 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*esrbij
286 0 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*esrbij
287 : END IF
288 539 : IF (calculate_forces) THEN
289 17 : desrbij = 2.0_dp*esrbij*(-etasrb*(1._dp + gscal*den2)*(dr - dr0))
290 17 : force_rr(1) = desrbij*rij(1)/dr
291 17 : force_rr(2) = desrbij*rij(2)/dr
292 17 : force_rr(3) = desrbij*rij(3)/dr
293 17 : atom_a = atom_of_kind(iatom)
294 17 : atom_b = atom_of_kind(jatom)
295 68 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
296 68 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
297 17 : IF (use_virial) THEN
298 1 : f1 = 1.0_dp
299 1 : IF (iatom == jatom) f1 = 0.5_dp
300 1 : CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
301 : END IF
302 : ! coordination number derivatives
303 : ! iatom
304 17 : fhua = -desrbij*cnfac(za)*(1._dp - c1srb*den1 - c2srb*den1*den1)
305 68 : DO i = 1, dcnum(iatom)%neighbors
306 51 : katom = dcnum(iatom)%nlist(i)
307 51 : kkind = kind_of(katom)
308 51 : atom_c = atom_of_kind(katom)
309 204 : rik = dcnum(iatom)%rik(:, i)
310 204 : drk = SQRT(SUM(rik(:)**2))
311 68 : IF (drk > 1.e-3_dp) THEN
312 204 : fdik(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
313 204 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fdik(:)
314 204 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
315 51 : IF (use_virial) THEN
316 3 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
317 : END IF
318 : END IF
319 : END DO
320 : ! jatom
321 17 : fhub = -desrbij*cnfac(zb)*(1._dp - c1srb*den1 - c2srb*den1*den1)
322 34 : DO i = 1, dcnum(jatom)%neighbors
323 17 : katom = dcnum(jatom)%nlist(i)
324 17 : kkind = kind_of(katom)
325 17 : atom_c = atom_of_kind(katom)
326 68 : rik = dcnum(jatom)%rik(:, i)
327 68 : drk = SQRT(SUM(rik(:)**2))
328 34 : IF (drk > 1.e-3_dp) THEN
329 68 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
330 68 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) - fdik(:)
331 68 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
332 17 : IF (use_virial) THEN
333 1 : 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 1646 : CALL neighbor_list_iterator_release(nl_iterator)
343 :
344 1646 : CALL timestop(handle)
345 :
346 3292 : 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 940 : 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 8316 : ppradius = 0.0_dp
367 940 : nkind = SIZE(ppradius, 1)
368 3030 : DO ikind = 1, nkind
369 2090 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
370 2090 : CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova, alpha=alphaa, zneff=zneffa, defined=defa)
371 2090 : IF (.NOT. defa) CYCLE
372 8804 : DO jkind = ikind, nkind
373 3686 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
374 3686 : CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb, alpha=alphab, zneff=zneffb, defined=defb)
375 3686 : IF (.NOT. defb) CYCLE
376 : rab = 0.0_dp
377 25952 : DO ir = 1, 24
378 25952 : rab = rab + 1.0_dp
379 25952 : saa = SQRT(alphaa*alphab)
380 25952 : erep = zneffa*zneffb/rab*EXP(-saa*rab**kfparam)
381 25952 : IF (erep < eps_pair) EXIT
382 : END DO
383 3684 : rab0 = rcova + rcovb
384 3684 : rab = MAX(rab, rab0 + 2.0_dp)
385 3684 : ppradius(ikind, jkind) = rab
386 9460 : ppradius(jkind, ikind) = ppradius(ikind, jkind)
387 : END DO
388 : END DO
389 :
390 940 : END SUBROUTINE xtb_pp_radius
391 :
392 : ! **************************************************************************************************
393 : !> \brief ...
394 : !> \param qs_env ...
395 : !> \param exb ...
396 : !> \param calculate_forces ...
397 : ! **************************************************************************************************
398 3438 : 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 3438 : 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 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: kx
413 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcab
414 : REAL(KIND=dp), DIMENSION(3) :: rij
415 3438 : 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 3438 : DIMENSION(:) :: neighbor_atoms
421 : TYPE(neighbor_list_iterator_p_type), &
422 3438 : DIMENSION(:), POINTER :: nl_iterator
423 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
424 3438 : POINTER :: sab_xb, sab_xtb_pp
425 3438 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426 3438 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
427 3438 : 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 3438 : 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 3438 : sab_xtb_pp=sab_xtb_pp)
442 :
443 3438 : nkind = SIZE(atomic_kind_set)
444 3438 : xtb_control => dft_control%qs_control%xtb_control
445 :
446 : ! global parameters
447 3438 : kxr = xtb_control%kxr
448 3438 : kx2 = xtb_control%kx2
449 :
450 3438 : NULLIFY (particle_set)
451 3438 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
452 3438 : natom = SIZE(particle_set)
453 3438 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
454 3438 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
455 :
456 3438 : IF (calculate_forces) THEN
457 444 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
458 834 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
459 : END IF
460 :
461 : ! list of neighbor atoms for XB term
462 18904 : ALLOCATE (neighbor_atoms(nkind))
463 12028 : DO ikind = 1, nkind
464 8590 : NULLIFY (neighbor_atoms(ikind)%coord)
465 8590 : NULLIFY (neighbor_atoms(ikind)%rab)
466 8590 : NULLIFY (neighbor_atoms(ikind)%katom)
467 8590 : CALL get_atomic_kind(atomic_kind_set(ikind), z=zat, natom=na)
468 12028 : IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
469 174 : ALLOCATE (neighbor_atoms(ikind)%coord(3, na))
470 290 : neighbor_atoms(ikind)%coord(1:3, 1:na) = 0.0_dp
471 174 : ALLOCATE (neighbor_atoms(ikind)%rab(na))
472 116 : neighbor_atoms(ikind)%rab(1:na) = HUGE(0.0_dp)
473 174 : ALLOCATE (neighbor_atoms(ikind)%katom(na))
474 116 : neighbor_atoms(ikind)%katom(1:na) = 0
475 : END IF
476 : END DO
477 : ! kx parameters
478 10314 : ALLOCATE (kx(nkind))
479 12028 : DO ikind = 1, nkind
480 8590 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
481 12028 : CALL get_xtb_atom_param(xtb_atom_a, kx=kx(ikind))
482 : END DO
483 : !
484 13752 : ALLOCATE (rcab(nkind, nkind))
485 12028 : DO ikind = 1, nkind
486 8590 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
487 8590 : CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova)
488 36042 : DO jkind = 1, nkind
489 24014 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
490 24014 : CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb)
491 32604 : rcab(ikind, jkind) = kxr*(rcova + rcovb)
492 : END DO
493 : END DO
494 :
495 3438 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
496 72844 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
497 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
498 69406 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
499 69406 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
500 69406 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
501 69406 : IF (.NOT. defined) CYCLE
502 69406 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
503 69406 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
504 69406 : IF (.NOT. defined) CYCLE
505 :
506 277624 : dr = SQRT(SUM(rij(:)**2))
507 :
508 : ! neighbor atom for XB term
509 72844 : IF (dr > 1.e-3_dp) THEN
510 54456 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
511 82 : atom_a = atom_of_kind(iatom)
512 82 : IF (dr < neighbor_atoms(ikind)%rab(atom_a)) THEN
513 41 : neighbor_atoms(ikind)%rab(atom_a) = dr
514 164 : neighbor_atoms(ikind)%coord(1:3, atom_a) = rij(1:3)
515 41 : neighbor_atoms(ikind)%katom(atom_a) = jatom
516 : END IF
517 : END IF
518 54456 : IF (ASSOCIATED(neighbor_atoms(jkind)%rab)) THEN
519 111 : atom_b = atom_of_kind(jatom)
520 111 : IF (dr < neighbor_atoms(jkind)%rab(atom_b)) THEN
521 74 : neighbor_atoms(jkind)%rab(atom_b) = dr
522 296 : neighbor_atoms(jkind)%coord(1:3, atom_b) = -rij(1:3)
523 74 : neighbor_atoms(jkind)%katom(atom_b) = iatom
524 : END IF
525 : END IF
526 : END IF
527 :
528 : END DO
529 3438 : CALL neighbor_list_iterator_release(nl_iterator)
530 :
531 3438 : exb = 0.0_dp
532 3438 : 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 3438 : calculate_forces, use_virial, force, virial, atprop)
535 :
536 12028 : DO ikind = 1, nkind
537 8590 : IF (ASSOCIATED(neighbor_atoms(ikind)%coord)) THEN
538 58 : DEALLOCATE (neighbor_atoms(ikind)%coord)
539 : END IF
540 8590 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
541 58 : DEALLOCATE (neighbor_atoms(ikind)%rab)
542 : END IF
543 12028 : IF (ASSOCIATED(neighbor_atoms(ikind)%katom)) THEN
544 58 : DEALLOCATE (neighbor_atoms(ikind)%katom)
545 : END IF
546 : END DO
547 3438 : DEALLOCATE (neighbor_atoms)
548 3438 : DEALLOCATE (kx, rcab)
549 :
550 3438 : CALL timestop(handle)
551 :
552 6876 : 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 3438 : 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 3438 : INTEGER, ALLOCATABLE, DIMENSION(:) :: matom
570 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dmloc
571 3438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
572 :
573 3438 : nkind = SIZE(neighbor_atoms)
574 12028 : DO ikind = 1, nkind
575 12028 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
576 58 : natom = SIZE(neighbor_atoms(ikind)%rab)
577 406 : ALLOCATE (dmloc(2*natom), matom(natom), coord(3, natom))
578 174 : dmloc = 0.0_dp
579 116 : DO iatom = 1, natom
580 58 : dmloc(2*iatom - 1) = neighbor_atoms(ikind)%rab(iatom)
581 116 : dmloc(2*iatom) = REAL(para_env%mepos, KIND=dp)
582 : END DO
583 58 : CALL para_env%minloc(dmloc)
584 290 : coord = 0.0_dp
585 116 : matom = 0
586 116 : DO iatom = 1, natom
587 58 : neighbor_atoms(ikind)%rab(iatom) = dmloc(2*iatom - 1)
588 116 : IF (NINT(dmloc(2*iatom)) == para_env%mepos) THEN
589 116 : coord(1:3, iatom) = neighbor_atoms(ikind)%coord(1:3, iatom)
590 29 : matom(iatom) = neighbor_atoms(ikind)%katom(iatom)
591 : END IF
592 : END DO
593 58 : CALL para_env%sum(coord)
594 290 : neighbor_atoms(ikind)%coord(1:3, :) = coord(1:3, :)
595 58 : CALL para_env%sum(matom)
596 116 : neighbor_atoms(ikind)%katom(:) = matom(:)
597 58 : DEALLOCATE (dmloc, matom, coord)
598 : END IF
599 : END DO
600 :
601 3438 : 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 68 : SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
622 34 : 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 34 : 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 34 : CALL timeset(routineN, handle)
656 :
657 : NULLIFY (nonbonded)
658 34 : NULLIFY (potparm)
659 34 : NULLIFY (ewald_env)
660 34 : nonbonded => xtb_control%nonbonded
661 34 : do_ewald = xtb_control%do_ewald
662 34 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
663 :
664 34 : ntype = SIZE(atomic_kind_set)
665 34 : CALL pair_potential_pp_create(potparm, ntype)
666 : !Assign input and potential info to potparm_nonbond
667 34 : CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
668 : !Initialize genetic potential
669 34 : CALL init_genpot(potparm, ntype)
670 :
671 34 : NULLIFY (pot)
672 34 : enonbonded = 0._dp
673 34 : energy_cutoff = 0._dp
674 :
675 34 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond)
676 227 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
677 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
678 193 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
679 193 : pot => potparm%pot(ikind, jkind)%pot
680 193 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
681 193 : rcut = SQRT(pot%rcutsq)
682 193 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
683 25 : fval = 1.0_dp
684 25 : IF (ikind == jkind) fval = 0.5_dp
685 : ! splines not implemented
686 25 : enonbonded = enonbonded + fval*ener_pot(pot, dr, energy_cutoff)
687 25 : IF (atprop%energy) THEN
688 16 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
689 16 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
690 : END IF
691 : END IF
692 :
693 193 : IF (calculate_forces) THEN
694 :
695 93 : kk = SIZE(pot%type)
696 93 : IF (kk /= 1) THEN
697 0 : CALL cp_warn(__LOCATION__, "Generic potential with type > 1 not implemented.")
698 0 : CPABORT("pot type")
699 : END IF
700 : ! rmin and rmax and rcut
701 93 : 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 93 : IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) CYCLE
704 : ! If within limits let's compute the potential...
705 93 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
706 :
707 9 : NULLIFY (nonbonded_section)
708 9 : nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED")
709 9 : CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx)
710 9 : CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr)
711 :
712 9 : dedf = fval*evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err)
713 9 : IF (ABS(err) > lerr) THEN
714 1 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
715 1 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
716 1 : CALL compress(this_error, .TRUE.)
717 1 : 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 1 : TRIM(def_error)//' .')
722 : END IF
723 :
724 9 : atom_i = atom_of_kind(iatom)
725 9 : atom_j = atom_of_kind(jatom)
726 :
727 36 : fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
728 :
729 36 : force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
730 36 : force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:)
731 :
732 9 : IF (use_virial) THEN
733 8 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
734 : END IF
735 :
736 : END IF
737 : END IF
738 193 : NULLIFY (pot)
739 : END DO
740 34 : CALL neighbor_list_iterator_release(nl_iterator)
741 34 : CALL finalizef()
742 :
743 34 : IF (ASSOCIATED(potparm)) THEN
744 34 : CALL pair_potential_pp_release(potparm)
745 : END IF
746 :
747 34 : CALL timestop(handle)
748 :
749 34 : 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 34 : 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 34 : NULLIFY (pot, logger)
780 :
781 34 : logger => cp_get_default_logger()
782 34 : iw = cp_logger_get_default_io_unit(logger)
783 :
784 170 : DO ikind = 1, SIZE(atomic_kind_set)
785 136 : atomic_kind => atomic_kind_set(ikind)
786 136 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
787 510 : DO jkind = ikind, SIZE(atomic_kind_set)
788 340 : atomic_kind => atomic_kind_set(jkind)
789 340 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
790 340 : found = .FALSE.
791 340 : name_atm_a = name_atm_a_local
792 340 : name_atm_b = name_atm_b_local
793 340 : CALL uppercase(name_atm_a)
794 340 : CALL uppercase(name_atm_b)
795 340 : pot => potparm%pot(ikind, jkind)%pot
796 340 : IF (ASSOCIATED(nonbonded)) THEN
797 680 : DO ingp = 1, SIZE(nonbonded%pot)
798 340 : 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 340 : ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. &
806 : (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
807 340 : ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN
808 34 : CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot)
809 : ! multiple potential not implemented, simply overwriting
810 34 : IF (found) &
811 : CALL cp_warn(__LOCATION__, &
812 : "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// &
813 0 : " 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 476 : IF (.NOT. found) THEN
820 306 : 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 34 : IF (do_ewald) THEN
828 16 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
829 16 : 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 34 : 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 3438 : 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 3438 : DIMENSION(:), POINTER :: nl_iterator
880 :
881 : ! exonent in angular term
882 3438 : alp = 6.0_dp
883 : ! loop over all atom pairs
884 3438 : CALL neighbor_list_iterator_create(nl_iterator, sab_xb)
885 3450 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
886 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
887 12 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
888 : ! ikind, iatom : Halogen
889 : ! jkind, jatom : Donor
890 12 : atom_a = atom_of_kind(iatom)
891 12 : katom = neighbor_atoms(ikind)%katom(atom_a)
892 12 : IF (katom == 0) CYCLE
893 12 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
894 12 : ddr = rcab(ikind, jkind)/dr
895 12 : ddr6 = ddr**6
896 12 : ddr12 = ddr6*ddr6
897 12 : eval = kx(ikind)*(ddr12 - kx2*ddr6)/(1.0_dp + ddr12)
898 : ! angle
899 48 : ria(1:3) = neighbor_atoms(ikind)%coord(1:3, atom_a)
900 48 : rja(1:3) = rij(1:3) - ria(1:3)
901 12 : drax = ria(1)**2 + ria(2)**2 + ria(3)**2
902 12 : drbx = dr*dr
903 12 : drab = rja(1)**2 + rja(2)**2 + rja(3)**2
904 12 : xy = SQRT(drbx*drax)
905 : ! cos angle B-X-A
906 12 : cosa = (drbx + drax - drab)/xy
907 12 : aterm = (0.5_dp - 0.25_dp*cosa)**alp
908 : !
909 12 : exb = exb + aterm*eval
910 12 : IF (atprop%energy) THEN
911 0 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*aterm*eval
912 0 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*aterm*eval
913 : END IF
914 : !
915 3450 : IF (calculate_forces) THEN
916 6 : kkind = kind_of(katom)
917 6 : atom_b = atom_of_kind(jatom)
918 6 : atom_c = atom_of_kind(katom)
919 : !
920 6 : deval = 6.0_dp*kx(ikind)*ddr6*(kx2*ddr12 + 2.0_dp*ddr6 - kx2)/(1.0_dp + ddr12)**2
921 6 : deval = -rcab(ikind, jkind)*deval/(dr*dr)/ddr
922 24 : fij(1:3) = aterm*deval*rij(1:3)/dr
923 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:)
924 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:)
925 6 : IF (use_virial) THEN
926 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
927 : END IF
928 : !
929 6 : fij(1:3) = 0.0_dp
930 6 : fia(1:3) = 0.0_dp
931 6 : fja(1:3) = 0.0_dp
932 6 : daterm = -0.25_dp*alp*(0.5_dp - 0.25_dp*cosa)**(alp - 1.0_dp)
933 6 : ddbx = 0.5_dp*(drab - drax + drbx)/xy/drbx
934 6 : ddax = 0.5_dp*(drab + drax - drbx)/xy/drax
935 6 : ddab = -1._dp/xy
936 24 : fij(1:3) = 2.0_dp*daterm*ddbx*rij(1:3)*eval
937 24 : fia(1:3) = 2.0_dp*daterm*ddax*ria(1:3)*eval
938 24 : fja(1:3) = 2.0_dp*daterm*ddab*rja(1:3)*eval
939 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:) - fia(:)
940 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:) + fja(:)
941 24 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fia(:) - fja(:)
942 6 : IF (use_virial) THEN
943 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
944 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fia, ria)
945 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fja, rja)
946 : END IF
947 : END IF
948 : END DO
949 3438 : CALL neighbor_list_iterator_release(nl_iterator)
950 :
951 3438 : END SUBROUTINE xb_energy
952 :
953 0 : END MODULE xtb_potentials
|
