Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 3996 : 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 3996 : 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 3996 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
111 : TYPE(atprop_type), POINTER :: atprop
112 : TYPE(neighbor_list_iterator_p_type), &
113 3996 : DIMENSION(:), POINTER :: nl_iterator
114 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
115 3996 : POINTER :: sab_xtb_pp
116 3996 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 3996 : 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 3996 : CALL timeset(routineN, handle)
122 :
123 3996 : 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 3996 : sab_xtb_pp=sab_xtb_pp)
130 3996 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
131 :
132 3996 : IF (calculate_forces) THEN
133 456 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
134 456 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
135 : END IF
136 :
137 3996 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
138 92820 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
139 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
140 88824 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
141 88824 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za, xtb_parameter=xtb_atom_a)
142 88824 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
143 88824 : IF (.NOT. defined) CYCLE
144 88824 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb, xtb_parameter=xtb_atom_b)
145 88824 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
146 88824 : IF (.NOT. defined) CYCLE
147 :
148 355296 : dr = SQRT(SUM(rij(:)**2))
149 : ! repulsive potential
150 92820 : IF (dr > 0.001_dp) THEN
151 :
152 : ! atomic parameters
153 72190 : CALL get_xtb_atom_param(xtb_atom_a, en=ena, alpha=alphaa, zneff=zneffa)
154 72190 : CALL get_xtb_atom_param(xtb_atom_b, en=enb, alpha=alphab, zneff=zneffb)
155 :
156 : ! scaling (not in papers! but in code)
157 72190 : den2 = (ena - enb)**2
158 72190 : den4 = den2*den2
159 72190 : sal = SQRT(alphaa*alphab)
160 72190 : ens = 1.0_dp + (0.01_dp*den2 + 0.01_dp*den4)*enscale
161 :
162 72190 : erepij = zneffa*zneffb/dr*EXP(-ens*sal*dr**kf)
163 72190 : erep = erep + erepij
164 72190 : 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 72190 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
169 7010 : derepij = -(1.0_dp/dr + ens*sal*kf*dr**(kf - 1.0_dp))*erepij
170 7010 : force_rr(1) = derepij*rij(1)/dr
171 7010 : force_rr(2) = derepij*rij(2)/dr
172 7010 : force_rr(3) = derepij*rij(3)/dr
173 7010 : atom_a = atom_of_kind(iatom)
174 7010 : atom_b = atom_of_kind(jatom)
175 28040 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
176 28040 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
177 7010 : IF (use_virial) THEN
178 3311 : f1 = 1.0_dp
179 3311 : IF (iatom == jatom) f1 = 0.5_dp
180 3311 : 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 3996 : CALL neighbor_list_iterator_release(nl_iterator)
187 :
188 3996 : CALL timestop(handle)
189 :
190 7992 : 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 810 : 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 810 : 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, fdika, fdikb, force_rr, rij, rik
226 810 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
227 : TYPE(atprop_type), POINTER :: atprop
228 : TYPE(neighbor_list_iterator_p_type), &
229 810 : DIMENSION(:), POINTER :: nl_iterator
230 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
231 810 : POINTER :: sab_xtb_pp
232 810 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
233 810 : 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 810 : CALL timeset(routineN, handle)
238 :
239 810 : 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 810 : 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 810 : kind_of=kind_of)
248 :
249 810 : IF (calculate_forces) THEN
250 28 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
251 28 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
252 : END IF
253 :
254 : ! SRB parameters
255 810 : ksrb = xtb_control%ksrb
256 810 : etasrb = xtb_control%esrb
257 810 : c1srb = xtb_control%c1srb*0.01_dp
258 810 : c2srb = xtb_control%c2srb*0.01_dp
259 810 : gscal = xtb_control%gscal
260 810 : shift = xtb_control%shift
261 :
262 810 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
263 21488 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
264 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
265 20678 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
266 20678 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
267 20678 : CALL get_xtb_atom_param(xtb_atom_a, z=za, electronegativity=enta, defined=defined)
268 20678 : IF (.NOT. defined) CYCLE
269 20678 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
270 20678 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, electronegativity=entb, defined=defined)
271 20678 : IF (.NOT. defined) CYCLE
272 :
273 82712 : dr = SQRT(SUM(rij(:)**2))
274 : ! short-ranged correction term
275 21488 : IF (dr > 0.001_dp) THEN
276 18490 : IF (za >= 5 .AND. za <= 9 .AND. zb >= 5 .AND. zb <= 9 .AND. za /= zb) THEN
277 263 : rra0 = r0srb(za) + cnfac(za)*cnumbers(iatom) + shift
278 263 : rrb0 = r0srb(zb) + cnfac(zb)*cnumbers(jatom) + shift
279 263 : den1 = ABS(ensrb(za) - ensrb(zb))
280 263 : dr0 = (rra0 + rrb0)*(1._dp - c1srb*den1 - c2srb*den1*den1)
281 263 : den2 = (enta - entb)**2
282 263 : esrbij = ksrb*EXP(-etasrb*(1._dp + gscal*den2)*(dr - dr0)**2)
283 263 : esrb = esrb + esrbij
284 263 : 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 263 : IF (calculate_forces) THEN
289 10 : desrbij = 2.0_dp*esrbij*(-etasrb*(1._dp + gscal*den2)*(dr - dr0))
290 10 : force_rr(1) = desrbij*rij(1)/dr
291 10 : force_rr(2) = desrbij*rij(2)/dr
292 10 : force_rr(3) = desrbij*rij(3)/dr
293 10 : atom_a = atom_of_kind(iatom)
294 10 : atom_b = atom_of_kind(jatom)
295 40 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - force_rr(:)
296 40 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + force_rr(:)
297 10 : IF (use_virial) THEN
298 0 : f1 = 1.0_dp
299 0 : IF (iatom == jatom) f1 = 0.5_dp
300 0 : CALL virial_pair_force(virial%pv_virial, -f1, force_rr, rij)
301 : END IF
302 : ! coordination number derivatives
303 : ! iatom
304 10 : fhua = -desrbij*cnfac(za)*(1._dp - c1srb*den1 - c2srb*den1*den1)
305 40 : DO i = 1, dcnum(iatom)%neighbors
306 30 : katom = dcnum(iatom)%nlist(i)
307 30 : kkind = kind_of(katom)
308 30 : atom_c = atom_of_kind(katom)
309 120 : rik = dcnum(iatom)%rik(:, i)
310 120 : drk = SQRT(SUM(rik(:)**2))
311 40 : IF (drk > 1.e-3_dp) THEN
312 120 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
313 120 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fdika(:)
314 120 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdika(:)
315 30 : IF (use_virial) THEN
316 0 : fdik = fdika + fdikb
317 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
318 : END IF
319 : END IF
320 : END DO
321 : ! jatom
322 10 : fhub = -desrbij*cnfac(zb)*(1._dp - c1srb*den1 - c2srb*den1*den1)
323 20 : DO i = 1, dcnum(jatom)%neighbors
324 10 : katom = dcnum(jatom)%nlist(i)
325 10 : kkind = kind_of(katom)
326 10 : atom_c = atom_of_kind(katom)
327 40 : rik = dcnum(jatom)%rik(:, i)
328 40 : drk = SQRT(SUM(rik(:)**2))
329 20 : IF (drk > 1.e-3_dp) THEN
330 40 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
331 40 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) - fdik(:)
332 40 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
333 10 : IF (use_virial) THEN
334 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
335 : END IF
336 : END IF
337 : END DO
338 : END IF
339 : END IF
340 : END IF
341 :
342 : END DO
343 810 : CALL neighbor_list_iterator_release(nl_iterator)
344 :
345 810 : CALL timestop(handle)
346 :
347 1620 : END SUBROUTINE srb_potential
348 :
349 : ! **************************************************************************************************
350 : !> \brief ...
351 : !> \param qs_kind_set ...
352 : !> \param ppradius ...
353 : !> \param eps_pair ...
354 : !> \param kfparam ...
355 : ! **************************************************************************************************
356 284 : SUBROUTINE xtb_pp_radius(qs_kind_set, ppradius, eps_pair, kfparam)
357 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
358 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: ppradius
359 : REAL(KIND=dp), INTENT(IN) :: eps_pair, kfparam
360 :
361 : INTEGER :: ikind, ir, jkind, nkind
362 : LOGICAL :: defa, defb
363 : REAL(KIND=dp) :: alphaa, alphab, erep, rab, rab0, rcova, &
364 : rcovb, saa, zneffa, zneffb
365 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
366 :
367 3120 : ppradius = 0.0_dp
368 284 : nkind = SIZE(ppradius, 1)
369 976 : DO ikind = 1, nkind
370 692 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
371 692 : CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova, alpha=alphaa, zneff=zneffa, defined=defa)
372 692 : IF (.NOT. defa) CYCLE
373 3082 : DO jkind = ikind, nkind
374 1416 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
375 1416 : CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb, alpha=alphab, zneff=zneffb, defined=defb)
376 1416 : IF (.NOT. defb) CYCLE
377 : rab = 0.0_dp
378 9064 : DO ir = 1, 24
379 9064 : rab = rab + 1.0_dp
380 9064 : saa = SQRT(alphaa*alphab)
381 9064 : erep = zneffa*zneffb/rab*EXP(-saa*rab**kfparam)
382 9064 : IF (erep < eps_pair) EXIT
383 : END DO
384 1414 : rab0 = rcova + rcovb
385 1414 : rab = MAX(rab, rab0 + 2.0_dp)
386 1414 : ppradius(ikind, jkind) = rab
387 3522 : ppradius(jkind, ikind) = ppradius(ikind, jkind)
388 : END DO
389 : END DO
390 :
391 284 : END SUBROUTINE xtb_pp_radius
392 :
393 : ! **************************************************************************************************
394 : !> \brief ...
395 : !> \param qs_env ...
396 : !> \param exb ...
397 : !> \param calculate_forces ...
398 : ! **************************************************************************************************
399 3186 : SUBROUTINE xb_interaction(qs_env, exb, calculate_forces)
400 :
401 : TYPE(qs_environment_type), POINTER :: qs_env
402 : REAL(KIND=dp), INTENT(INOUT) :: exb
403 : LOGICAL, INTENT(IN) :: calculate_forces
404 :
405 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xb_interaction'
406 :
407 : INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
408 : jatom, jkind, na, natom, nkind, zat
409 3186 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
410 : INTEGER, DIMENSION(3) :: cell
411 : LOGICAL :: defined, use_virial
412 : REAL(KIND=dp) :: dr, kx2, kxr, rcova, rcovb
413 3186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: kx
414 3186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcab
415 : REAL(KIND=dp), DIMENSION(3) :: rij
416 3186 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
417 : TYPE(atprop_type), POINTER :: atprop
418 : TYPE(dft_control_type), POINTER :: dft_control
419 : TYPE(mp_para_env_type), POINTER :: para_env
420 : TYPE(neighbor_atoms_type), ALLOCATABLE, &
421 3186 : DIMENSION(:) :: neighbor_atoms
422 : TYPE(neighbor_list_iterator_p_type), &
423 3186 : DIMENSION(:), POINTER :: nl_iterator
424 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
425 3186 : POINTER :: sab_xb, sab_xtb_pp
426 3186 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
427 3186 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
428 3186 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
429 : TYPE(virial_type), POINTER :: virial
430 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
431 : TYPE(xtb_control_type), POINTER :: xtb_control
432 :
433 3186 : CALL timeset(routineN, handle)
434 :
435 : CALL get_qs_env(qs_env=qs_env, &
436 : atomic_kind_set=atomic_kind_set, &
437 : qs_kind_set=qs_kind_set, &
438 : para_env=para_env, &
439 : atprop=atprop, &
440 : dft_control=dft_control, &
441 : sab_xb=sab_xb, &
442 3186 : sab_xtb_pp=sab_xtb_pp)
443 :
444 3186 : nkind = SIZE(atomic_kind_set)
445 3186 : xtb_control => dft_control%qs_control%xtb_control
446 :
447 : ! global parameters
448 3186 : kxr = xtb_control%kxr
449 3186 : kx2 = xtb_control%kx2
450 :
451 3186 : NULLIFY (particle_set)
452 3186 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
453 3186 : natom = SIZE(particle_set)
454 3186 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
455 3186 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
456 :
457 3186 : IF (calculate_forces) THEN
458 428 : CALL get_qs_env(qs_env=qs_env, virial=virial, force=force)
459 806 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
460 : END IF
461 :
462 : ! list of neighbor atoms for XB term
463 17392 : ALLOCATE (neighbor_atoms(nkind))
464 11020 : DO ikind = 1, nkind
465 7834 : NULLIFY (neighbor_atoms(ikind)%coord)
466 7834 : NULLIFY (neighbor_atoms(ikind)%rab)
467 7834 : NULLIFY (neighbor_atoms(ikind)%katom)
468 7834 : CALL get_atomic_kind(atomic_kind_set(ikind), z=zat, natom=na)
469 11020 : IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
470 174 : ALLOCATE (neighbor_atoms(ikind)%coord(3, na))
471 290 : neighbor_atoms(ikind)%coord(1:3, 1:na) = 0.0_dp
472 174 : ALLOCATE (neighbor_atoms(ikind)%rab(na))
473 116 : neighbor_atoms(ikind)%rab(1:na) = HUGE(0.0_dp)
474 174 : ALLOCATE (neighbor_atoms(ikind)%katom(na))
475 116 : neighbor_atoms(ikind)%katom(1:na) = 0
476 : END IF
477 : END DO
478 : ! kx parameters
479 9558 : ALLOCATE (kx(nkind))
480 11020 : DO ikind = 1, nkind
481 7834 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
482 11020 : CALL get_xtb_atom_param(xtb_atom_a, kx=kx(ikind))
483 : END DO
484 : !
485 12744 : ALLOCATE (rcab(nkind, nkind))
486 11020 : DO ikind = 1, nkind
487 7834 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
488 7834 : CALL get_xtb_atom_param(xtb_atom_a, rcov=rcova)
489 32766 : DO jkind = 1, nkind
490 21746 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
491 21746 : CALL get_xtb_atom_param(xtb_atom_b, rcov=rcovb)
492 29580 : rcab(ikind, jkind) = kxr*(rcova + rcovb)
493 : END DO
494 : END DO
495 :
496 3186 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_pp)
497 71332 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
498 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
499 68146 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
500 68146 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
501 68146 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
502 68146 : IF (.NOT. defined) CYCLE
503 68146 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
504 68146 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
505 68146 : IF (.NOT. defined) CYCLE
506 :
507 272584 : dr = SQRT(SUM(rij(:)**2))
508 :
509 : ! neighbor atom for XB term
510 71332 : IF (dr > 1.e-3_dp) THEN
511 53700 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
512 82 : atom_a = atom_of_kind(iatom)
513 82 : IF (dr < neighbor_atoms(ikind)%rab(atom_a)) THEN
514 41 : neighbor_atoms(ikind)%rab(atom_a) = dr
515 164 : neighbor_atoms(ikind)%coord(1:3, atom_a) = rij(1:3)
516 41 : neighbor_atoms(ikind)%katom(atom_a) = jatom
517 : END IF
518 : END IF
519 53700 : IF (ASSOCIATED(neighbor_atoms(jkind)%rab)) THEN
520 111 : atom_b = atom_of_kind(jatom)
521 111 : IF (dr < neighbor_atoms(jkind)%rab(atom_b)) THEN
522 74 : neighbor_atoms(jkind)%rab(atom_b) = dr
523 296 : neighbor_atoms(jkind)%coord(1:3, atom_b) = -rij(1:3)
524 74 : neighbor_atoms(jkind)%katom(atom_b) = iatom
525 : END IF
526 : END IF
527 : END IF
528 :
529 : END DO
530 3186 : CALL neighbor_list_iterator_release(nl_iterator)
531 :
532 3186 : exb = 0.0_dp
533 3186 : CALL xb_neighbors(neighbor_atoms, para_env)
534 : CALL xb_energy(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
535 3186 : calculate_forces, use_virial, force, virial, atprop)
536 :
537 11020 : DO ikind = 1, nkind
538 7834 : IF (ASSOCIATED(neighbor_atoms(ikind)%coord)) THEN
539 58 : DEALLOCATE (neighbor_atoms(ikind)%coord)
540 : END IF
541 7834 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
542 58 : DEALLOCATE (neighbor_atoms(ikind)%rab)
543 : END IF
544 11020 : IF (ASSOCIATED(neighbor_atoms(ikind)%katom)) THEN
545 58 : DEALLOCATE (neighbor_atoms(ikind)%katom)
546 : END IF
547 : END DO
548 3186 : DEALLOCATE (neighbor_atoms)
549 3186 : DEALLOCATE (kx, rcab)
550 :
551 3186 : CALL timestop(handle)
552 :
553 6372 : END SUBROUTINE xb_interaction
554 :
555 : ! **************************************************************************************************
556 : !> \brief Distributes the neighbor atom information to all processors
557 : !>
558 : !> \param neighbor_atoms ...
559 : !> \param para_env ...
560 : !> \par History
561 : !> 1.2019 JGH
562 : !> \version 1.1
563 : ! **************************************************************************************************
564 3186 : SUBROUTINE xb_neighbors(neighbor_atoms, para_env)
565 : TYPE(neighbor_atoms_type), DIMENSION(:), &
566 : INTENT(INOUT) :: neighbor_atoms
567 : TYPE(mp_para_env_type), POINTER :: para_env
568 :
569 : INTEGER :: iatom, ikind, natom, nkind
570 3186 : INTEGER, ALLOCATABLE, DIMENSION(:) :: matom
571 3186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dmloc
572 3186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
573 :
574 3186 : nkind = SIZE(neighbor_atoms)
575 11020 : DO ikind = 1, nkind
576 11020 : IF (ASSOCIATED(neighbor_atoms(ikind)%rab)) THEN
577 58 : natom = SIZE(neighbor_atoms(ikind)%rab)
578 406 : ALLOCATE (dmloc(2*natom), matom(natom), coord(3, natom))
579 174 : dmloc = 0.0_dp
580 116 : DO iatom = 1, natom
581 58 : dmloc(2*iatom - 1) = neighbor_atoms(ikind)%rab(iatom)
582 116 : dmloc(2*iatom) = REAL(para_env%mepos, KIND=dp)
583 : END DO
584 58 : CALL para_env%minloc(dmloc)
585 290 : coord = 0.0_dp
586 116 : matom = 0
587 116 : DO iatom = 1, natom
588 58 : neighbor_atoms(ikind)%rab(iatom) = dmloc(2*iatom - 1)
589 116 : IF (NINT(dmloc(2*iatom)) == para_env%mepos) THEN
590 116 : coord(1:3, iatom) = neighbor_atoms(ikind)%coord(1:3, iatom)
591 29 : matom(iatom) = neighbor_atoms(ikind)%katom(iatom)
592 : END IF
593 : END DO
594 58 : CALL para_env%sum(coord)
595 290 : neighbor_atoms(ikind)%coord(1:3, :) = coord(1:3, :)
596 58 : CALL para_env%sum(matom)
597 116 : neighbor_atoms(ikind)%katom(:) = matom(:)
598 58 : DEALLOCATE (dmloc, matom, coord)
599 : END IF
600 : END DO
601 :
602 3186 : END SUBROUTINE xb_neighbors
603 :
604 : ! **************************************************************************************************
605 : !> \brief Computes a correction for nonbonded interactions based on a generic potential
606 : !>
607 : !> \param enonbonded energy contribution
608 : !> \param force ...
609 : !> \param qs_env ...
610 : !> \param xtb_control ...
611 : !> \param sab_xtb_nonbond ...
612 : !> \param atomic_kind_set ...
613 : !> \param calculate_forces ...
614 : !> \param use_virial ...
615 : !> \param virial ...
616 : !> \param atprop ...
617 : !> \param atom_of_kind ..
618 : !> \par History
619 : !> 12.2018 JGH
620 : !> \version 1.1
621 : ! **************************************************************************************************
622 68 : SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
623 34 : atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
624 : REAL(dp), INTENT(INOUT) :: enonbonded
625 : TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
626 : POINTER :: force
627 : TYPE(qs_environment_type), POINTER :: qs_env
628 : TYPE(xtb_control_type), POINTER :: xtb_control
629 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
630 : INTENT(IN), POINTER :: sab_xtb_nonbond
631 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
632 : POINTER :: atomic_kind_set
633 : LOGICAL, INTENT(IN) :: calculate_forces, use_virial
634 : TYPE(virial_type), INTENT(IN), POINTER :: virial
635 : TYPE(atprop_type), INTENT(IN), POINTER :: atprop
636 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind
637 :
638 : CHARACTER(len=*), PARAMETER :: routineN = 'nonbonded_correction'
639 :
640 : CHARACTER(LEN=default_string_length) :: def_error, this_error
641 : INTEGER :: atom_i, atom_j, handle, iatom, ikind, &
642 : jatom, jkind, kk, ntype
643 : INTEGER, DIMENSION(3) :: cell
644 : LOGICAL :: do_ewald
645 : REAL(KIND=dp) :: dedf, dr, dx, energy_cutoff, err, fval, &
646 : lerr, rcut
647 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
648 : TYPE(ewald_environment_type), POINTER :: ewald_env
649 : TYPE(neighbor_list_iterator_p_type), &
650 34 : DIMENSION(:), POINTER :: nl_iterator
651 : TYPE(pair_potential_p_type), POINTER :: nonbonded
652 : TYPE(pair_potential_pp_type), POINTER :: potparm
653 : TYPE(pair_potential_single_type), POINTER :: pot
654 : TYPE(section_vals_type), POINTER :: nonbonded_section
655 :
656 34 : CALL timeset(routineN, handle)
657 :
658 : NULLIFY (nonbonded)
659 34 : NULLIFY (potparm)
660 34 : NULLIFY (ewald_env)
661 34 : nonbonded => xtb_control%nonbonded
662 34 : do_ewald = xtb_control%do_ewald
663 34 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
664 :
665 34 : ntype = SIZE(atomic_kind_set)
666 34 : CALL pair_potential_pp_create(potparm, ntype)
667 : !Assign input and potential info to potparm_nonbond
668 34 : CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
669 : !Initialize genetic potential
670 34 : CALL init_genpot(potparm, ntype)
671 :
672 34 : NULLIFY (pot)
673 34 : enonbonded = 0._dp
674 34 : energy_cutoff = 0._dp
675 :
676 34 : CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond)
677 227 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
678 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
679 193 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
680 193 : pot => potparm%pot(ikind, jkind)%pot
681 193 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
682 193 : rcut = SQRT(pot%rcutsq)
683 193 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
684 25 : fval = 1.0_dp
685 25 : IF (ikind == jkind) fval = 0.5_dp
686 : ! splines not implemented
687 25 : enonbonded = enonbonded + fval*ener_pot(pot, dr, energy_cutoff)
688 25 : IF (atprop%energy) THEN
689 16 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
690 16 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*fval*ener_pot(pot, dr, energy_cutoff)
691 : END IF
692 : END IF
693 :
694 193 : IF (calculate_forces) THEN
695 :
696 93 : kk = SIZE(pot%type)
697 93 : IF (kk /= 1) THEN
698 0 : CALL cp_warn(__LOCATION__, "Generic potential with type > 1 not implemented.")
699 0 : CPABORT("pot type")
700 : END IF
701 : ! rmin and rmax and rcut
702 93 : IF ((pot%set(kk)%rmin /= not_initialized) .AND. (dr < pot%set(kk)%rmin)) CYCLE
703 : ! An upper boundary for the potential definition was defined
704 93 : IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) CYCLE
705 : ! If within limits let's compute the potential...
706 93 : IF (dr <= rcut .AND. dr > 1.E-3_dp) THEN
707 :
708 9 : NULLIFY (nonbonded_section)
709 9 : nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED")
710 9 : CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx)
711 9 : CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr)
712 :
713 9 : dedf = fval*evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err)
714 9 : IF (ABS(err) > lerr) THEN
715 1 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
716 1 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
717 1 : CALL compress(this_error, .TRUE.)
718 1 : CALL compress(def_error, .TRUE.)
719 : CALL cp_warn(__LOCATION__, &
720 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
721 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
722 1 : TRIM(def_error)//' .')
723 : END IF
724 :
725 9 : atom_i = atom_of_kind(iatom)
726 9 : atom_j = atom_of_kind(jatom)
727 :
728 36 : fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
729 :
730 36 : force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
731 36 : force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:)
732 :
733 9 : IF (use_virial) THEN
734 8 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
735 : END IF
736 :
737 : END IF
738 : END IF
739 193 : NULLIFY (pot)
740 : END DO
741 34 : CALL neighbor_list_iterator_release(nl_iterator)
742 34 : CALL finalizef()
743 :
744 34 : IF (ASSOCIATED(potparm)) THEN
745 34 : CALL pair_potential_pp_release(potparm)
746 : END IF
747 :
748 34 : CALL timestop(handle)
749 :
750 34 : END SUBROUTINE nonbonded_correction
751 :
752 : ! **************************************************************************************************
753 : !> \brief ...
754 : !> \param atomic_kind_set ...
755 : !> \param nonbonded ...
756 : !> \param potparm ...
757 : !> \param ewald_env ...
758 : !> \param do_ewald ...
759 : ! **************************************************************************************************
760 34 : SUBROUTINE force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald)
761 :
762 : ! routine based on force_field_pack_nonbond
763 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
764 : POINTER :: atomic_kind_set
765 : TYPE(pair_potential_p_type), INTENT(IN), POINTER :: nonbonded
766 : TYPE(pair_potential_pp_type), INTENT(INOUT), &
767 : POINTER :: potparm
768 : TYPE(ewald_environment_type), INTENT(IN), POINTER :: ewald_env
769 : LOGICAL, INTENT(IN) :: do_ewald
770 :
771 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
772 : name_atm_b, name_atm_b_local
773 : INTEGER :: ikind, ingp, iw, jkind
774 : LOGICAL :: found
775 : REAL(KIND=dp) :: ewald_rcut
776 : TYPE(atomic_kind_type), POINTER :: atomic_kind
777 : TYPE(cp_logger_type), POINTER :: logger
778 : TYPE(pair_potential_single_type), POINTER :: pot
779 :
780 34 : NULLIFY (pot, logger)
781 :
782 34 : logger => cp_get_default_logger()
783 34 : iw = cp_logger_get_default_io_unit(logger)
784 :
785 170 : DO ikind = 1, SIZE(atomic_kind_set)
786 136 : atomic_kind => atomic_kind_set(ikind)
787 136 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
788 510 : DO jkind = ikind, SIZE(atomic_kind_set)
789 340 : atomic_kind => atomic_kind_set(jkind)
790 340 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
791 340 : found = .FALSE.
792 340 : name_atm_a = name_atm_a_local
793 340 : name_atm_b = name_atm_b_local
794 340 : CALL uppercase(name_atm_a)
795 340 : CALL uppercase(name_atm_b)
796 340 : pot => potparm%pot(ikind, jkind)%pot
797 340 : IF (ASSOCIATED(nonbonded)) THEN
798 680 : DO ingp = 1, SIZE(nonbonded%pot)
799 340 : IF ((TRIM(nonbonded%pot(ingp)%pot%at1) == "*") .OR. &
800 : (TRIM(nonbonded%pot(ingp)%pot%at2) == "*")) CYCLE
801 :
802 : !IF (iw > 0) WRITE (iw, *) "TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), &
803 : ! " with ", TRIM(nonbonded%pot(ingp)%pot%at1), &
804 : ! TRIM(nonbonded%pot(ingp)%pot%at2)
805 : IF ((((name_atm_a) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
806 340 : ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. &
807 : (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. &
808 340 : ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN
809 34 : CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot)
810 : ! multiple potential not implemented, simply overwriting
811 34 : IF (found) &
812 : CALL cp_warn(__LOCATION__, &
813 : "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// &
814 0 : " and "//TRIM(name_atm_b)//" OVERWRITING! ")
815 : !IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
816 : found = .TRUE.
817 : END IF
818 : END DO
819 : END IF
820 476 : IF (.NOT. found) THEN
821 306 : CALL pair_potential_single_clean(pot)
822 : !IF (iw > 0) WRITE (iw, *) " NOTFOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
823 : END IF
824 : END DO !jkind
825 : END DO !ikind
826 :
827 : ! Cutoff is defined always as the maximum between the FF and Ewald
828 34 : IF (do_ewald) THEN
829 16 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
830 16 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
831 : !IF (iw > 0) WRITE (iw, *) " RCUT ", SQRT(pot%rcutsq), ewald_rcut
832 : END IF
833 :
834 34 : END SUBROUTINE force_field_pack_nonbond_pot_correction
835 :
836 : ! **************************************************************************************************
837 : !> \brief Computes the interaction term between Br/I/At and donor atoms
838 : !>
839 : !> \param exb ...
840 : !> \param neighbor_atoms ...
841 : !> \param atom_of_kind ...
842 : !> \param kind_of ...
843 : !> \param sab_xb ...
844 : !> \param kx ...
845 : !> \param kx2 ...
846 : !> \param rcab ...
847 : !> \param calculate_forces ...
848 : !> \param use_virial ...
849 : !> \param force ...
850 : !> \param virial ...
851 : !> \param atprop ...
852 : !> \par History
853 : !> 12.2018 JGH
854 : !> \version 1.1
855 : ! **************************************************************************************************
856 3186 : SUBROUTINE xb_energy(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, &
857 : calculate_forces, use_virial, force, virial, atprop)
858 : REAL(dp), INTENT(INOUT) :: exb
859 : TYPE(neighbor_atoms_type), DIMENSION(:), &
860 : INTENT(IN) :: neighbor_atoms
861 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
862 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
863 : POINTER :: sab_xb
864 : REAL(dp), DIMENSION(:), INTENT(IN) :: kx
865 : REAL(dp), INTENT(IN) :: kx2
866 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rcab
867 : LOGICAL, INTENT(IN) :: calculate_forces, use_virial
868 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
869 : TYPE(virial_type), POINTER :: virial
870 : TYPE(atprop_type), POINTER :: atprop
871 :
872 : INTEGER :: atom_a, atom_b, atom_c, iatom, ikind, &
873 : jatom, jkind, katom, kkind
874 : INTEGER, DIMENSION(3) :: cell
875 : REAL(KIND=dp) :: alp, aterm, cosa, daterm, ddab, ddax, &
876 : ddbx, ddr, ddr12, ddr6, deval, dr, &
877 : drab, drax, drbx, eval, xy
878 : REAL(KIND=dp), DIMENSION(3) :: fia, fij, fja, ria, rij, rja
879 : TYPE(neighbor_list_iterator_p_type), &
880 3186 : DIMENSION(:), POINTER :: nl_iterator
881 :
882 : ! exonent in angular term
883 3186 : alp = 6.0_dp
884 : ! loop over all atom pairs
885 3186 : CALL neighbor_list_iterator_create(nl_iterator, sab_xb)
886 3198 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
887 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
888 12 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
889 : ! ikind, iatom : Halogen
890 : ! jkind, jatom : Donor
891 12 : atom_a = atom_of_kind(iatom)
892 12 : katom = neighbor_atoms(ikind)%katom(atom_a)
893 12 : IF (katom == 0) CYCLE
894 12 : dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2)
895 12 : ddr = rcab(ikind, jkind)/dr
896 12 : ddr6 = ddr**6
897 12 : ddr12 = ddr6*ddr6
898 12 : eval = kx(ikind)*(ddr12 - kx2*ddr6)/(1.0_dp + ddr12)
899 : ! angle
900 48 : ria(1:3) = neighbor_atoms(ikind)%coord(1:3, atom_a)
901 48 : rja(1:3) = rij(1:3) - ria(1:3)
902 12 : drax = ria(1)**2 + ria(2)**2 + ria(3)**2
903 12 : drbx = dr*dr
904 12 : drab = rja(1)**2 + rja(2)**2 + rja(3)**2
905 12 : xy = SQRT(drbx*drax)
906 : ! cos angle B-X-A
907 12 : cosa = (drbx + drax - drab)/xy
908 12 : aterm = (0.5_dp - 0.25_dp*cosa)**alp
909 : !
910 12 : exb = exb + aterm*eval
911 12 : IF (atprop%energy) THEN
912 0 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*aterm*eval
913 0 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*aterm*eval
914 : END IF
915 : !
916 3198 : IF (calculate_forces) THEN
917 6 : kkind = kind_of(katom)
918 6 : atom_b = atom_of_kind(jatom)
919 6 : atom_c = atom_of_kind(katom)
920 : !
921 6 : deval = 6.0_dp*kx(ikind)*ddr6*(kx2*ddr12 + 2.0_dp*ddr6 - kx2)/(1.0_dp + ddr12)**2
922 6 : deval = -rcab(ikind, jkind)*deval/(dr*dr)/ddr
923 24 : fij(1:3) = aterm*deval*rij(1:3)/dr
924 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:)
925 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:)
926 6 : IF (use_virial) THEN
927 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
928 : END IF
929 : !
930 6 : fij(1:3) = 0.0_dp
931 6 : fia(1:3) = 0.0_dp
932 6 : fja(1:3) = 0.0_dp
933 6 : daterm = -0.25_dp*alp*(0.5_dp - 0.25_dp*cosa)**(alp - 1.0_dp)
934 6 : ddbx = 0.5_dp*(drab - drax + drbx)/xy/drbx
935 6 : ddax = 0.5_dp*(drab + drax - drbx)/xy/drax
936 6 : ddab = -1._dp/xy
937 24 : fij(1:3) = 2.0_dp*daterm*ddbx*rij(1:3)*eval
938 24 : fia(1:3) = 2.0_dp*daterm*ddax*ria(1:3)*eval
939 24 : fja(1:3) = 2.0_dp*daterm*ddab*rja(1:3)*eval
940 24 : force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fij(:) - fia(:)
941 24 : force(jkind)%repulsive(:, atom_b) = force(jkind)%repulsive(:, atom_b) + fij(:) + fja(:)
942 24 : force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fia(:) - fja(:)
943 6 : IF (use_virial) THEN
944 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij)
945 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fia, ria)
946 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fja, rja)
947 : END IF
948 : END IF
949 : END DO
950 3186 : CALL neighbor_list_iterator_release(nl_iterator)
951 :
952 3186 : END SUBROUTINE xb_energy
953 :
954 0 : END MODULE xtb_potentials
|
