Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of D3 dispersion
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dispersion_d3
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind,&
16 : get_atomic_kind_set
17 : USE atprop_types, ONLY: atprop_array_init,&
18 : atprop_type
19 : USE cell_types, ONLY: cell_type,&
20 : get_cell,&
21 : pbc,&
22 : plane_distance
23 : USE cp_files, ONLY: close_file,&
24 : open_file
25 : USE input_constants, ONLY: vdw_pairpot_dftd3,&
26 : vdw_pairpot_dftd3bj
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: twopi
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE molecule_types, ONLY: molecule_type
31 : USE particle_types, ONLY: particle_type
32 : USE physcon, ONLY: kcalmol
33 : USE qs_dispersion_cnum, ONLY: d3_cnumber,&
34 : dcnum_distribute,&
35 : dcnum_type,&
36 : exclude_d3_kind_pair
37 : USE qs_dispersion_types, ONLY: dftd3_pp,&
38 : qs_atom_dispersion_type,&
39 : qs_dispersion_type
40 : USE qs_dispersion_utils, ONLY: cellhash
41 : USE qs_environment_types, ONLY: get_qs_env,&
42 : qs_environment_type
43 : USE qs_force_types, ONLY: qs_force_type
44 : USE qs_kind_types, ONLY: get_qs_kind,&
45 : qs_kind_type
46 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
47 : neighbor_list_iterate,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type
52 : USE virial_methods, ONLY: virial_pair_force
53 : USE virial_types, ONLY: virial_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d3'
61 :
62 : PUBLIC :: calculate_dispersion_d3_pairpot, dftd3_c6_param
63 :
64 : ! **************************************************************************************************
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param qs_env ...
71 : !> \param dispersion_env ...
72 : !> \param evdw ...
73 : !> \param calculate_forces ...
74 : !> \param unit_nr ...
75 : !> \param atevdw ...
76 : ! **************************************************************************************************
77 5822 : SUBROUTINE calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
78 5822 : unit_nr, atevdw)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
82 : REAL(KIND=dp), INTENT(OUT) :: evdw
83 : LOGICAL, INTENT(IN) :: calculate_forces
84 : INTEGER, INTENT(IN) :: unit_nr
85 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atevdw
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d3_pairpot'
88 :
89 : INTEGER :: atom_a, atom_b, atom_c, atom_d, handle, hashb, hashc, i, ia, iat, iatom, icx, &
90 : icy, icz, idmp, ikind, ilist, imol, jatom, jkind, katom, kkind, kstart, latom, lkind, &
91 : max_elem, maxc, mepos, na, natom, nb, nc, nkind, num_pe, za, zb, zc
92 5822 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
93 : INTEGER, DIMENSION(3) :: cell_b, cell_c, ncell, periodic
94 5822 : INTEGER, DIMENSION(:), POINTER :: atom_list
95 : LOGICAL :: atenergy, atex, debugall, domol, exclude_pair, floating_a, floating_b, &
96 : floating_c, ghost_a, ghost_b, ghost_c, is000, use_virial
97 5822 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: dodisp, exclude
98 : REAL(KIND=dp) :: a1, a2, alp6, alp8, ang, c6, c8, c9, cc6ab, cc6bc, cc6ca, cnum, dc6a, dc6b, &
99 : dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, de6, de8, de91, de921, &
100 : de922, dea, dfdab6, dfdab8, dfdabc, dr, drk, e6, e6tot, e8, e8tot, e9, e9tot, elrc6, &
101 : elrc8, elrc9, eps_cn, esrb, f0ab, fac, fac0, fdab6, fdab8, fdabc, gsrb, kgc8, nab, nabc, &
102 : r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, rabc, rc2, rcc, rcut, rs6, rs8, s1, s2, s3, &
103 : s6, s8, s8i, s9, srbe, ssrb, t1srb, t2srb
104 5822 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
105 5822 : cnumfix, radd2
106 5822 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
107 : REAL(KIND=dp), DIMENSION(3) :: fdij, fdik, ra, rab, rb, rb0, rbc, rc, &
108 : rc0, rca, rij, rik, sab_max
109 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread
110 5822 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
111 5822 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
112 : TYPE(atprop_type), POINTER :: atprop
113 : TYPE(cell_type), POINTER :: cell
114 5822 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
115 5822 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 : TYPE(neighbor_list_iterator_p_type), &
118 5822 : DIMENSION(:), POINTER :: nl_iterator
119 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
120 5822 : POINTER :: sab_vdw
121 5822 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
122 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a, disp_b, disp_c
123 5822 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
124 5822 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(virial_type), POINTER :: virial
126 :
127 5822 : CALL timeset(routineN, handle)
128 :
129 5822 : evdw = 0._dp
130 :
131 5822 : NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
132 :
133 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
134 5822 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, para_env=para_env, atprop=atprop)
135 :
136 5822 : debugall = dispersion_env%verbose
137 :
138 : ! atomic energy and stress arrays
139 5822 : atenergy = atprop%energy
140 5822 : IF (atenergy) THEN
141 88 : CALL atprop_array_init(atprop%atevdw, natom)
142 88 : atener => atprop%atevdw
143 : END IF
144 : ! external atomic energy
145 5822 : atex = .FALSE.
146 5822 : IF (PRESENT(atevdw)) THEN
147 2 : atex = .TRUE.
148 : END IF
149 :
150 5822 : NULLIFY (particle_set)
151 5822 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
152 5822 : natom = SIZE(particle_set)
153 :
154 5822 : NULLIFY (force)
155 5822 : CALL get_qs_env(qs_env=qs_env, force=force)
156 5822 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
157 5822 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
158 5822 : pv_virial_thread(:, :) = 0._dp
159 :
160 46576 : ALLOCATE (dodisp(nkind), exclude(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
161 19498 : DO ikind = 1, nkind
162 13676 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
163 13676 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
164 13676 : dodisp(ikind) = disp_a%defined
165 13676 : exclude(ikind) = ghost_a .OR. floating_a
166 13676 : atomnumber(ikind) = za
167 13676 : c6d2(ikind) = disp_a%c6
168 33174 : radd2(ikind) = disp_a%vdw_radii
169 : END DO
170 :
171 17466 : ALLOCATE (rcpbc(3, natom))
172 55522 : DO iatom = 1, natom
173 55522 : rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
174 : END DO
175 :
176 5822 : rcut = 2._dp*dispersion_env%rc_disp
177 5822 : rc2 = rcut*rcut
178 :
179 5822 : maxc = SIZE(dispersion_env%c6ab, 3)
180 5822 : max_elem = SIZE(dispersion_env%c6ab, 1)
181 5822 : alp6 = dispersion_env%alp
182 5822 : alp8 = alp6 + 2._dp
183 5822 : s6 = dispersion_env%s6
184 5822 : s8 = dispersion_env%s8
185 5822 : s9 = dispersion_env%s6
186 5822 : rs6 = dispersion_env%sr6
187 5822 : rs8 = dispersion_env%sr8
188 5822 : a1 = dispersion_env%a1
189 5822 : a2 = dispersion_env%a2
190 5822 : eps_cn = dispersion_env%eps_cn
191 5822 : e6tot = 0._dp
192 5822 : e8tot = 0._dp
193 5822 : e9tot = 0._dp
194 5822 : esrb = 0._dp
195 5822 : domol = dispersion_env%domol
196 : ! molecule correction
197 5822 : kgc8 = dispersion_env%kgc8
198 5822 : IF (domol) THEN
199 2 : CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set)
200 6 : ALLOCATE (atom2mol(natom))
201 6 : DO imol = 1, SIZE(molecule_set)
202 16 : DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom
203 14 : atom2mol(iat) = imol
204 : END DO
205 : END DO
206 : END IF
207 : ! damping type
208 5822 : idmp = 0
209 5822 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
210 : idmp = 1
211 5542 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
212 5542 : idmp = 2
213 : END IF
214 : ! SRB parameters
215 5822 : ssrb = dispersion_env%srb_params(1)
216 5822 : gsrb = dispersion_env%srb_params(2)
217 5822 : t1srb = dispersion_env%srb_params(3)
218 5822 : t2srb = dispersion_env%srb_params(4)
219 :
220 5822 : IF (unit_nr > 0) THEN
221 45 : WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
222 45 : WRITE (unit_nr, *) " Scaling parameter (s8) ", s8
223 45 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
224 20 : WRITE (unit_nr, *) " Zero Damping parameter (sr6)", rs6
225 20 : WRITE (unit_nr, *) " Zero Damping parameter (sr8)", rs8
226 25 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
227 25 : WRITE (unit_nr, *) " BJ Damping parameter (a1) ", a1
228 25 : WRITE (unit_nr, *) " BJ Damping parameter (a2) ", a2
229 : END IF
230 45 : WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
231 45 : IF (dispersion_env%lrc) THEN
232 1 : WRITE (unit_nr, *) " Apply a long range correction"
233 : END IF
234 45 : IF (dispersion_env%srb) THEN
235 3 : WRITE (unit_nr, *) " Apply a short range bond correction"
236 3 : WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
237 : END IF
238 45 : IF (domol) THEN
239 1 : WRITE (unit_nr, *) " Inter-molecule scaling parameter (s8) ", kgc8
240 : END IF
241 : END IF
242 : ! Calculate coordination numbers
243 5822 : NULLIFY (particle_set)
244 5822 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
245 5822 : natom = SIZE(particle_set)
246 17466 : ALLOCATE (cnumbers(natom))
247 5822 : cnumbers = 0._dp
248 :
249 5822 : IF (calculate_forces .OR. debugall) THEN
250 12662 : ALLOCATE (dcnum(natom))
251 10926 : dcnum(:)%neighbors = 0
252 10926 : DO iatom = 1, natom
253 10926 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
254 : END DO
255 : ELSE
256 9908 : ALLOCATE (dcnum(1))
257 : END IF
258 :
259 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, exclude, atomnumber, &
260 5822 : calculate_forces, 1)
261 :
262 5822 : CALL para_env%sum(cnumbers)
263 :
264 : ! for parallel runs we have to update dcnum on all processors
265 5822 : IF (calculate_forces .OR. debugall) THEN
266 868 : CALL dcnum_distribute(dcnum, para_env)
267 868 : IF (unit_nr > 0 .AND. SIZE(dcnum) > 0) THEN
268 4 : WRITE (unit_nr, *)
269 4 : WRITE (unit_nr, *) " ATOM Coordination Neighbors"
270 19 : DO i = 1, natom
271 19 : WRITE (unit_nr, "(I6,F20.10,I12)") i, cnumbers(i), dcnum(i)%neighbors
272 : END DO
273 4 : WRITE (unit_nr, *)
274 : END IF
275 : END IF
276 :
277 5822 : nab = 0._dp
278 5822 : nabc = 0._dp
279 5822 : IF (dispersion_env%doabc) THEN
280 120 : rcc = 2._dp*dispersion_env%rc_disp
281 120 : CALL get_cell(cell=cell, periodic=periodic)
282 120 : sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
283 120 : sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
284 120 : sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
285 480 : ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
286 120 : IF (unit_nr > 0) THEN
287 4 : WRITE (unit_nr, *) " Calculate C9 Terms"
288 4 : WRITE (unit_nr, "(A,T20,I3,A,I3)") " Search in cells ", -ncell(1), ":", ncell(1)
289 4 : WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(2), ":", ncell(2)
290 4 : WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(3), ":", ncell(3)
291 4 : WRITE (unit_nr, *)
292 : END IF
293 120 : IF (dispersion_env%c9cnst) THEN
294 72 : IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
295 216 : ALLOCATE (cnumfix(natom))
296 72 : cnumfix = 0._dp
297 : ! first use the default values
298 356 : DO iatom = 1, natom
299 284 : ikind = kind_of(iatom)
300 284 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
301 356 : cnumfix(iatom) = dispersion_env%cn(za)
302 : END DO
303 : ! now check for changes from default
304 72 : IF (ASSOCIATED(dispersion_env%cnkind)) THEN
305 12 : DO i = 1, SIZE(dispersion_env%cnkind)
306 6 : ikind = dispersion_env%cnkind(i)%kind
307 6 : cnum = dispersion_env%cnkind(i)%cnum
308 6 : CPASSERT(ikind <= nkind)
309 6 : CPASSERT(ikind > 0)
310 6 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, atom_list=atom_list)
311 32 : DO ia = 1, na
312 20 : iatom = atom_list(ia)
313 26 : cnumfix(iatom) = cnum
314 : END DO
315 : END DO
316 : END IF
317 72 : IF (ASSOCIATED(dispersion_env%cnlist)) THEN
318 6 : DO i = 1, SIZE(dispersion_env%cnlist)
319 14 : DO ilist = 1, dispersion_env%cnlist(i)%natom
320 8 : iatom = dispersion_env%cnlist(i)%atom(ilist)
321 12 : cnumfix(iatom) = dispersion_env%cnlist(i)%cnum
322 : END DO
323 : END DO
324 : END IF
325 72 : IF (unit_nr > 0) THEN
326 5 : DO i = 1, natom
327 5 : IF (ABS(cnumbers(i) - cnumfix(i)) > 0.5_dp) THEN
328 0 : WRITE (unit_nr, "(A,T20,A,I6,T41,2F10.3)") " Difference in CN ", "Atom:", &
329 0 : i, cnumbers(i), cnumfix(i)
330 : END IF
331 : END DO
332 1 : WRITE (unit_nr, *)
333 : END IF
334 : END IF
335 : END IF
336 :
337 5822 : sab_vdw => dispersion_env%sab_vdw
338 :
339 5822 : num_pe = 1
340 5822 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
341 :
342 5822 : mepos = 0
343 9028093 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
344 9022271 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
345 :
346 9022271 : IF (exclude(ikind) .OR. exclude(jkind)) CYCLE
347 :
348 9022214 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
349 :
350 9021980 : za = atomnumber(ikind)
351 9021980 : zb = atomnumber(jkind)
352 : ! vdW potential
353 36087920 : dr = SQRT(SUM(rij(:)**2))
354 9027802 : IF (dr <= rcut) THEN
355 9021980 : nab = nab + 1._dp
356 9021980 : fac = 1._dp
357 9021980 : IF (iatom == jatom) fac = 0.5_dp
358 9021980 : IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
359 8996827 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
360 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
361 320 : exclude=exclude_pair)
362 320 : IF (exclude_pair) CYCLE
363 : END IF
364 : ! C6 coefficient
365 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
366 8996595 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
367 8996595 : c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
368 8996595 : dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
369 8996595 : dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
370 8996595 : r6 = dr**6
371 8996595 : r8 = r6*dr*dr
372 8996595 : s8i = s8
373 8996595 : IF (domol) THEN
374 3568 : IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
375 1452 : s8i = kgc8
376 : END IF
377 : END IF
378 : ! damping
379 8996595 : IF (idmp == 1) THEN
380 : ! zero
381 3892811 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
382 3892811 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
383 3892811 : e6 = s6*fac*c6*fdab6/r6
384 3892811 : e8 = s8i*fac*c8*fdab8/r8
385 5103784 : ELSE IF (idmp == 2) THEN
386 : ! BJ
387 5103784 : r0ab = SQRT(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
388 5103784 : f0ab = a1*r0ab + a2
389 5103784 : fdab6 = 1.0_dp/(r6 + f0ab**6)
390 5103784 : fdab8 = 1.0_dp/(r8 + f0ab**8)
391 5103784 : e6 = s6*fac*c6*fdab6
392 5103784 : e8 = s8i*fac*c8*fdab8
393 : ELSE
394 0 : CPABORT("Unknown DFT-D3 damping function:")
395 : END IF
396 8996595 : IF (dispersion_env%srb .AND. dr < 30.0d0) THEN
397 135 : srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
398 135 : esrb = esrb + srbe
399 135 : evdw = evdw - srbe
400 : ELSE
401 : srbe = 0.0_dp
402 : END IF
403 8996595 : evdw = evdw - e6 - e8
404 8996595 : e6tot = e6tot - e6
405 8996595 : e8tot = e8tot - e8
406 8996595 : IF (atenergy) THEN
407 2765506 : atener(iatom) = atener(iatom) - 0.5_dp*(e6 + e8 + srbe)
408 2765506 : atener(jatom) = atener(jatom) - 0.5_dp*(e6 + e8 + srbe)
409 : END IF
410 8996595 : IF (atex) THEN
411 850 : atevdw(iatom) = atevdw(iatom) - 0.5_dp*(e6 + e8 + srbe)
412 850 : atevdw(jatom) = atevdw(jatom) - 0.5_dp*(e6 + e8 + srbe)
413 : END IF
414 8996595 : IF (calculate_forces) THEN
415 : ! damping
416 3222753 : IF (idmp == 1) THEN
417 : ! zero
418 2003084 : de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
419 2003084 : de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
420 1219669 : ELSE IF (idmp == 2) THEN
421 : ! BJ
422 1219669 : de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
423 1219669 : de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
424 : ELSE
425 0 : CPABORT("Unknown DFT-D3 damping function:")
426 : END IF
427 12891012 : fdij(:) = (de6 + de8)*rij(:)/dr*fac
428 3222753 : IF (dispersion_env%srb .AND. dr < 30.0d0) THEN
429 172 : fdij(:) = fdij(:) + srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
430 : END IF
431 3222753 : atom_a = atom_of_kind(iatom)
432 3222753 : atom_b = atom_of_kind(jatom)
433 12891012 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
434 12891012 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
435 3222753 : IF (use_virial) THEN
436 2054743 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
437 : END IF
438 : ! forces from the r-dependence of the coordination numbers
439 3222753 : IF (idmp == 1) THEN
440 : ! zero
441 2003084 : dc6a = -s6*fac*dc6a*fdab6/r6
442 2003084 : dc6b = -s6*fac*dc6b*fdab6/r6
443 2003084 : dc8a = -s8i*fac*dc8a*fdab8/r8
444 2003084 : dc8b = -s8i*fac*dc8b*fdab8/r8
445 1219669 : ELSE IF (idmp == 2) THEN
446 : ! BJ
447 1219669 : dc6a = -s6*fac*dc6a*fdab6
448 1219669 : dc6b = -s6*fac*dc6b*fdab6
449 1219669 : dc8a = -s8i*fac*dc8a*fdab8
450 1219669 : dc8b = -s8i*fac*dc8b*fdab8
451 : ELSE
452 0 : CPABORT("Unknown DFT-D3 damping function:")
453 : END IF
454 43851792 : DO i = 1, dcnum(iatom)%neighbors
455 40629039 : katom = dcnum(iatom)%nlist(i)
456 40629039 : kkind = kind_of(katom)
457 162516156 : rik = dcnum(iatom)%rik(:, i)
458 162516156 : drk = SQRT(SUM(rik(:)**2))
459 162516156 : fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
460 40629039 : atom_c = atom_of_kind(katom)
461 162516156 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
462 162516156 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
463 43851792 : IF (use_virial) THEN
464 29564842 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
465 : END IF
466 : END DO
467 43848080 : DO i = 1, dcnum(jatom)%neighbors
468 40625327 : katom = dcnum(jatom)%nlist(i)
469 40625327 : kkind = kind_of(katom)
470 162501308 : rik = dcnum(jatom)%rik(:, i)
471 162501308 : drk = SQRT(SUM(rik(:)**2))
472 162501308 : fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
473 40625327 : atom_c = atom_of_kind(katom)
474 162501308 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
475 162501308 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
476 43848080 : IF (use_virial) THEN
477 29561014 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
478 : END IF
479 : END DO
480 : END IF
481 8996595 : IF (dispersion_env%doabc) THEN
482 17472 : CALL get_iterator_info(nl_iterator, cell=cell_b)
483 17472 : hashb = cellhash(cell_b, ncell)
484 27102 : is000 = (ALL(cell_b == 0))
485 279552 : rb0(:) = MATMUL(cell%hmat, cell_b)
486 17472 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
487 87360 : rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
488 69888 : r2ab = SUM((ra - rb)**2)
489 113008 : DO icx = -ncell(1), ncell(1)
490 705584 : DO icy = -ncell(2), ncell(2)
491 4831008 : DO icz = -ncell(3), ncell(3)
492 4142896 : cell_c(1) = icx
493 4142896 : cell_c(2) = icy
494 4142896 : cell_c(3) = icz
495 4142896 : hashc = cellhash(cell_c, ncell)
496 4848480 : IF (is000 .AND. (ALL(cell_c == 0))) THEN
497 : ! CASE 1: all atoms in (000), use only ordered triples
498 1230 : kstart = MAX(jatom + 1, iatom + 1)
499 1230 : fac0 = 1.0_dp
500 4141666 : ELSE IF (is000) THEN
501 : ! CASE 2: AB in (000), C in other cell
502 : ! This case covers also all instances with BC in same
503 : ! cell not (000)
504 : kstart = 1
505 : fac0 = 1.0_dp
506 : ELSE
507 : ! These are case 2 again, cycle
508 4105058 : IF (hashc == hashb) CYCLE
509 4776962 : IF (ALL(cell_c == 0)) CYCLE
510 : ! CASE 3: A in (000) and B and C in different cells
511 : kstart = 1
512 : fac0 = 1.0_dp/3.0_dp
513 : END IF
514 65766592 : rc0 = MATMUL(cell%hmat, cell_c)
515 16433678 : DO katom = kstart, natom
516 11730690 : kkind = kind_of(katom)
517 11730690 : IF (exclude(kkind) .OR. .NOT. dodisp(kkind)) CYCLE
518 46836152 : rc(:) = rcpbc(:, katom) + rc0(:)
519 46836152 : r2bc = SUM((rb - rc)**2)
520 11709038 : IF (r2bc >= rc2) CYCLE
521 10009420 : r2ca = SUM((rc - ra)**2)
522 2502355 : IF (r2ca >= rc2) CYCLE
523 1296400 : r2abc = r2ab*r2bc*r2ca
524 1296400 : IF (r2abc <= 0.001_dp) CYCLE
525 1296400 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
526 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
527 5066 : kkind, exclude_pair)
528 5066 : IF (exclude_pair) CYCLE
529 : END IF
530 : ! this is an empirical scaling
531 5404008 : IF (r2abc <= 0.01_dp*rc2*rc2*rc2) THEN
532 55654 : kkind = kind_of(katom)
533 55654 : atom_c = atom_of_kind(katom)
534 55654 : zc = atomnumber(kkind)
535 : ! avoid double counting!
536 55654 : fac = 1._dp
537 55654 : IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
538 55654 : IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
539 55654 : fac = fac*fac0
540 55654 : nabc = nabc + 1._dp
541 55654 : IF (dispersion_env%c9cnst) THEN
542 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
543 34940 : cnumfix(iatom), cnumfix(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
544 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
545 34940 : cnumfix(jatom), cnumfix(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
546 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
547 34940 : cnumfix(katom), cnumfix(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
548 : ELSE
549 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
550 20714 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
551 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
552 20714 : cnumbers(jatom), cnumbers(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
553 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
554 20714 : cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
555 : END IF
556 55654 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
557 55654 : rabc = r2abc**(1._dp/6._dp)
558 : r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
559 55654 : dispersion_env%r0ab(zc, za))**(1._dp/3._dp)
560 : ! bug fixed 3.10.2017
561 : ! correct value from alp6=14 to 16 as used in original paper
562 55654 : CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
563 55654 : s1 = r2ab + r2bc - r2ca
564 55654 : s2 = r2ab + r2ca - r2bc
565 55654 : s3 = r2ca + r2bc - r2ab
566 55654 : ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
567 :
568 55654 : e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
569 55654 : evdw = evdw - e9
570 55654 : e9tot = e9tot - e9
571 55654 : IF (atenergy) THEN
572 20568 : atener(iatom) = atener(iatom) - e9/3._dp
573 20568 : atener(jatom) = atener(jatom) - e9/3._dp
574 20568 : atener(katom) = atener(katom) - e9/3._dp
575 : END IF
576 55654 : IF (atex) THEN
577 10284 : atevdw(iatom) = atevdw(iatom) - e9/3._dp
578 10284 : atevdw(jatom) = atevdw(jatom) - e9/3._dp
579 10284 : atevdw(katom) = atevdw(katom) - e9/3._dp
580 : END IF
581 :
582 55654 : IF (calculate_forces) THEN
583 287740 : rab = rb - ra; rbc = rc - rb; rca = ra - rc
584 28774 : de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
585 28774 : de91 = -de91/3._dp*rabc + 3._dp*e9
586 28774 : dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
587 115096 : fdij(:) = de91*rab(:)/r2ab
588 115096 : fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
589 115096 : fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
590 115096 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
591 115096 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
592 28774 : IF (use_virial) THEN
593 18000 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
594 : END IF
595 115096 : fdij(:) = de91*rbc(:)/r2bc
596 115096 : fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
597 115096 : fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
598 115096 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
599 115096 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
600 28774 : IF (use_virial) THEN
601 18000 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
602 : END IF
603 115096 : fdij(:) = de91*rca(:)/r2ca
604 115096 : fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
605 115096 : fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
606 115096 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
607 115096 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
608 28774 : IF (use_virial) THEN
609 18000 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
610 : END IF
611 :
612 28774 : IF (.NOT. dispersion_env%c9cnst) THEN
613 : ! forces from the r-dependence of the coordination numbers
614 : ! atomic stress not implemented
615 :
616 14402 : de91 = 0.5_dp*e9/cc6ab
617 14402 : de921 = de91*dcc6aba
618 14402 : de922 = de91*dcc6abb
619 69790 : DO i = 1, dcnum(iatom)%neighbors
620 55388 : latom = dcnum(iatom)%nlist(i)
621 55388 : lkind = kind_of(latom)
622 55388 : rik(1) = dcnum(iatom)%rik(1, i)
623 55388 : rik(2) = dcnum(iatom)%rik(2, i)
624 55388 : rik(3) = dcnum(iatom)%rik(3, i)
625 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
626 221552 : fdik(:) = -de921*dcnum(iatom)%dvals(i)*rik(:)/drk
627 55388 : atom_d = atom_of_kind(latom)
628 221552 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
629 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
630 69790 : IF (use_virial) THEN
631 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
632 : END IF
633 : END DO
634 69790 : DO i = 1, dcnum(jatom)%neighbors
635 55388 : latom = dcnum(jatom)%nlist(i)
636 55388 : lkind = kind_of(latom)
637 55388 : rik(1) = dcnum(jatom)%rik(1, i)
638 55388 : rik(2) = dcnum(jatom)%rik(2, i)
639 55388 : rik(3) = dcnum(jatom)%rik(3, i)
640 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
641 221552 : fdik(:) = -de922*dcnum(jatom)%dvals(i)*rik(:)/drk
642 55388 : atom_d = atom_of_kind(latom)
643 221552 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
644 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
645 69790 : IF (use_virial) THEN
646 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
647 : END IF
648 : END DO
649 :
650 14402 : de91 = 0.5_dp*e9/cc6bc
651 14402 : de921 = de91*dcc6bcb
652 14402 : de922 = de91*dcc6bcc
653 69790 : DO i = 1, dcnum(jatom)%neighbors
654 55388 : latom = dcnum(jatom)%nlist(i)
655 55388 : lkind = kind_of(latom)
656 55388 : rik(1) = dcnum(jatom)%rik(1, i)
657 55388 : rik(2) = dcnum(jatom)%rik(2, i)
658 55388 : rik(3) = dcnum(jatom)%rik(3, i)
659 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
660 221552 : fdik(:) = -de921*dcnum(jatom)%dvals(i)*rik(:)/drk
661 55388 : atom_d = atom_of_kind(latom)
662 221552 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
663 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
664 69790 : IF (use_virial) THEN
665 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
666 : END IF
667 : END DO
668 69790 : DO i = 1, dcnum(katom)%neighbors
669 55388 : latom = dcnum(katom)%nlist(i)
670 55388 : lkind = kind_of(latom)
671 55388 : rik(1) = dcnum(katom)%rik(1, i)
672 55388 : rik(2) = dcnum(katom)%rik(2, i)
673 55388 : rik(3) = dcnum(katom)%rik(3, i)
674 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
675 221552 : fdik(:) = -de922*dcnum(katom)%dvals(i)*rik(:)/drk
676 55388 : atom_d = atom_of_kind(latom)
677 221552 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
678 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
679 69790 : IF (use_virial) THEN
680 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
681 : END IF
682 : END DO
683 :
684 14402 : de91 = 0.5_dp*e9/cc6ca
685 14402 : de921 = de91*dcc6cac
686 14402 : de922 = de91*dcc6caa
687 69790 : DO i = 1, dcnum(katom)%neighbors
688 55388 : latom = dcnum(katom)%nlist(i)
689 55388 : lkind = kind_of(latom)
690 55388 : rik(1) = dcnum(katom)%rik(1, i)
691 55388 : rik(2) = dcnum(katom)%rik(2, i)
692 55388 : rik(3) = dcnum(katom)%rik(3, i)
693 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
694 221552 : fdik(:) = -de921*dcnum(katom)%dvals(i)*rik(:)/drk
695 55388 : atom_d = atom_of_kind(latom)
696 221552 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
697 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
698 69790 : IF (use_virial) THEN
699 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
700 : END IF
701 : END DO
702 69790 : DO i = 1, dcnum(iatom)%neighbors
703 55388 : latom = dcnum(iatom)%nlist(i)
704 55388 : lkind = kind_of(latom)
705 55388 : rik(1) = dcnum(iatom)%rik(1, i)
706 55388 : rik(2) = dcnum(iatom)%rik(2, i)
707 55388 : rik(3) = dcnum(iatom)%rik(3, i)
708 55388 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
709 221552 : fdik(:) = -de922*dcnum(iatom)%dvals(i)*rik(:)/drk
710 55388 : atom_d = atom_of_kind(latom)
711 221552 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
712 221552 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
713 69790 : IF (use_virial) THEN
714 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
715 : END IF
716 : END DO
717 : END IF
718 :
719 : END IF
720 :
721 : END IF
722 : END DO
723 : END DO
724 : END DO
725 : END DO
726 :
727 : END IF
728 : END IF
729 : END IF
730 : END DO
731 :
732 75686 : virial%pv_virial = virial%pv_virial + pv_virial_thread
733 :
734 5822 : CALL neighbor_list_iterator_release(nl_iterator)
735 :
736 5822 : elrc6 = 0._dp
737 5822 : elrc8 = 0._dp
738 5822 : elrc9 = 0._dp
739 : ! Long range correction (atomic contributions not implemented)
740 5822 : IF (dispersion_env%lrc) THEN
741 144 : ALLOCATE (cnkind(nkind))
742 48 : cnkind = 0._dp
743 : ! first use the default values
744 106 : DO ikind = 1, nkind
745 58 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
746 106 : cnkind(ikind) = dispersion_env%cn(za)
747 : END DO
748 : ! now check for changes from default
749 48 : IF (ASSOCIATED(dispersion_env%cnkind)) THEN
750 16 : DO i = 1, SIZE(dispersion_env%cnkind)
751 8 : ikind = dispersion_env%cnkind(i)%kind
752 16 : cnkind(ikind) = dispersion_env%cnkind(i)%cnum
753 : END DO
754 : END IF
755 106 : DO ikind = 1, nkind
756 58 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
757 58 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
758 58 : IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) CYCLE
759 236 : DO jkind = 1, nkind
760 74 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
761 74 : CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
762 74 : IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
763 72 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
764 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
765 8 : exclude=exclude_pair)
766 8 : IF (exclude_pair) CYCLE
767 : END IF
768 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
769 66 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
770 66 : elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
771 66 : c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
772 66 : elrc8 = elrc8 - s8*twopi*REAL(na*nb, KIND=dp)*c8/(5._dp*rcut**5*cell%deth)
773 198 : IF (dispersion_env%doabc) THEN
774 156 : DO kkind = 1, nkind
775 92 : CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
776 92 : CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
777 92 : IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
778 90 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
779 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, kkind, &
780 4 : exclude_pair)
781 4 : IF (exclude_pair) CYCLE
782 : END IF
783 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
784 88 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
785 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
786 88 : cnkind(kkind), cnkind(ikind), dispersion_env%k3, cc6ca, dcc6aba, dcc6abb)
787 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
788 88 : cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
789 88 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
790 244 : elrc9 = elrc9 - s9*64._dp*twopi*REAL(na*nb*nc, KIND=dp)*c9/(6._dp*rcut**3*cell%deth**2)
791 : END DO
792 : END IF
793 : END DO
794 : END DO
795 48 : IF (idmp == 1) THEN
796 : ! Zero damping: the switch (cutoff_switch_d3) removes part of the pairs in [rcut - 1 bohr, rcut] from
797 : ! the explicit sum. Add that part with g(r) = 1, as GROMACS does for switched potentials. The factors
798 : ! depend only on rcut, so forces are unchanged and the tail pressure stays E_lrc/V = -dE_lrc/dV.
799 26 : elrc6 = elrc6*(1._dp + switched_shell_fraction(rcut, 6))
800 26 : elrc8 = elrc8*(1._dp + switched_shell_fraction(rcut, 8))
801 : END IF
802 48 : IF (use_virial) THEN
803 : ! Pressure of the pairs not contained in the explicit sum (g(r) = 1 beyond the cutoff).
804 : ! BJ damping: the pair sum is truncated sharply at rcut and its virial misses the pairs crossing rcut,
805 : ! so the tail pressure is (n/3) E_lrc,n / V, i.e. 2*E6/V and 8/3*E8/V (as pair_modify tail in LAMMPS
806 : ! or DispCorr in GROMACS).
807 : ! Zero damping: damping_d3 switches the pairs off smoothly over [rcut - 1 bohr, rcut]; the virial of
808 : ! the switching force already contains the cutoff contribution and the tail pressure is E_lrc,n / V.
809 : ! The three-body term (2*elrc9) has not been checked.
810 34 : IF (para_env%is_source()) THEN
811 68 : DO i = 1, 3
812 68 : IF (idmp == 2) THEN
813 33 : virial%pv_virial(i, i) = virial%pv_virial(i, i) + (2._dp*elrc6 + 8._dp/3._dp*elrc8 + 2._dp*elrc9)
814 : ELSE
815 18 : virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
816 : END IF
817 : END DO
818 : END IF
819 : END IF
820 48 : DEALLOCATE (cnkind)
821 : END IF
822 :
823 5822 : DEALLOCATE (cnumbers)
824 5822 : IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
825 72 : DEALLOCATE (cnumfix)
826 : END IF
827 5822 : IF (calculate_forces .OR. debugall) THEN
828 10926 : DO iatom = 1, natom
829 10926 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
830 : END DO
831 868 : DEALLOCATE (dcnum)
832 : ELSE
833 4954 : DEALLOCATE (dcnum)
834 : END IF
835 :
836 : ! set dispersion energy
837 :
838 5822 : CALL para_env%sum(e6tot)
839 5822 : CALL para_env%sum(e8tot)
840 5822 : CALL para_env%sum(e9tot)
841 5822 : CALL para_env%sum(evdw)
842 5822 : CALL para_env%sum(nab)
843 5822 : CALL para_env%sum(nabc)
844 5822 : e6tot = e6tot + elrc6
845 5822 : e8tot = e8tot + elrc8
846 5822 : e9tot = e9tot + elrc9
847 : ! For printing, we need all contributions
848 5822 : evdw = evdw + (elrc6 + elrc8 + elrc9)
849 5822 : IF (unit_nr > 0) THEN
850 45 : WRITE (unit_nr, "(A,F20.0)") " E6 vdW terms :", nab
851 45 : WRITE (unit_nr, *) " E6 vdW energy [au/kcal] :", e6tot, e6tot*kcalmol
852 45 : WRITE (unit_nr, *) " E8 vdW energy [au/kcal] :", e8tot, e8tot*kcalmol
853 45 : WRITE (unit_nr, *) " %E8 on total vdW energy :", e8tot/evdw*100._dp
854 45 : WRITE (unit_nr, "(A,F20.0)") " E9 vdW terms :", nabc
855 45 : WRITE (unit_nr, *) " E9 vdW energy [au/kcal] :", e9tot, e9tot*kcalmol
856 45 : WRITE (unit_nr, *) " %E9 on total vdW energy :", e9tot/evdw*100._dp
857 45 : IF (dispersion_env%lrc) THEN
858 1 : WRITE (unit_nr, *) " E LRC C6 [au/kcal] :", elrc6, elrc6*kcalmol
859 1 : WRITE (unit_nr, *) " E LRC C8 [au/kcal] :", elrc8, elrc8*kcalmol
860 1 : WRITE (unit_nr, *) " E LRC C9 [au/kcal] :", elrc9, elrc9*kcalmol
861 : END IF
862 : END IF
863 :
864 5822 : evdw = evdw/para_env%num_pe
865 :
866 5822 : DEALLOCATE (dodisp, exclude, atomnumber, rcpbc, radd2, c6d2)
867 :
868 5822 : IF (domol) THEN
869 2 : DEALLOCATE (atom2mol)
870 : END IF
871 :
872 5822 : CALL timestop(handle)
873 :
874 11644 : END SUBROUTINE calculate_dispersion_d3_pairpot
875 :
876 : ! **************************************************************************************************
877 : !> \brief ...
878 : !> \param c6ab ...
879 : !> \param maxci ...
880 : !> \param filename ...
881 : !> \param para_env ...
882 : ! **************************************************************************************************
883 532 : SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
884 :
885 : REAL(KIND=dp), DIMENSION(:, :, :, :, :) :: c6ab
886 : INTEGER, DIMENSION(:) :: maxci
887 : CHARACTER(LEN=*) :: filename
888 : TYPE(mp_para_env_type), POINTER :: para_env
889 :
890 : INTEGER :: funit, iadr, iat, jadr, jat, kk, nl, &
891 : nlines, nn, ref_stride
892 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pars
893 :
894 532 : IF (para_env%is_source()) THEN
895 : ! Read the DFT-D3 C6AB parameters from file "filename"
896 271 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
897 271 : READ (funit, *) nl, nlines
898 : END IF
899 532 : CALL para_env%bcast(nl)
900 532 : CALL para_env%bcast(nlines)
901 1596 : ALLOCATE (pars(nl))
902 532 : IF (para_env%is_source()) THEN
903 271 : READ (funit, *) pars(1:nl)
904 271 : CALL close_file(unit_number=funit)
905 : END IF
906 532 : CALL para_env%bcast(pars)
907 :
908 : ! Store C6AB coefficients in an array
909 837812752 : c6ab = -1._dp
910 55328 : maxci = 0
911 532 : ref_stride = get_ref_stride(pars, nlines)
912 532 : kk = 1
913 30479344 : DO nn = 1, nlines
914 30478812 : iat = NINT(pars(kk + 1))
915 30478812 : jat = NINT(pars(kk + 2))
916 30478812 : CALL limit(iat, jat, iadr, jadr, ref_stride)
917 30478812 : maxci(iat) = MAX(maxci(iat), iadr)
918 30478812 : maxci(jat) = MAX(maxci(jat), jadr)
919 30478812 : c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
920 30478812 : c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
921 30478812 : c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
922 :
923 30478812 : c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
924 30478812 : c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
925 30478812 : c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
926 30479344 : kk = (nn*5) + 1
927 : END DO
928 :
929 532 : DEALLOCATE (pars)
930 :
931 532 : END SUBROUTINE dftd3_c6_param
932 :
933 : ! **************************************************************************************************
934 : !> \brief Detects the reference-state encoding used in the DFT-D3 C6AB table.
935 : !> \param pars DFT-D3 C6AB parameter table.
936 : !> \param nlines Number of C6AB parameter lines.
937 : !> \return Reference-state stride.
938 : ! **************************************************************************************************
939 532 : FUNCTION get_ref_stride(pars, nlines) RESULT(ref_stride)
940 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pars
941 : INTEGER, INTENT(IN) :: nlines
942 : INTEGER :: ref_stride
943 :
944 : INTEGER :: kk, nn
945 :
946 532 : ref_stride = 100
947 532 : kk = 1
948 2849924 : DO nn = 1, nlines
949 2849924 : IF (NINT(pars(kk + 1)) > 1000 .OR. NINT(pars(kk + 2)) > 1000) THEN
950 : ref_stride = 1000
951 : EXIT
952 : END IF
953 2849392 : kk = (nn*5) + 1
954 : END DO
955 :
956 532 : END FUNCTION get_ref_stride
957 :
958 : ! **************************************************************************************************
959 : !> \brief Decodes element numbers and reference-state indices.
960 : !> \param iat Encoded first element number.
961 : !> \param jat Encoded second element number.
962 : !> \param iadr First reference-state index.
963 : !> \param jadr Second reference-state index.
964 : !> \param ref_stride Reference-state stride.
965 : ! **************************************************************************************************
966 30478812 : SUBROUTINE limit(iat, jat, iadr, jadr, ref_stride)
967 : INTEGER :: iat, jat, iadr, jadr, ref_stride
968 :
969 30478812 : iadr = 1
970 30478812 : jadr = 1
971 102051964 : DO WHILE (iat > ref_stride)
972 71573152 : iat = iat - ref_stride
973 102051964 : iadr = iadr + 1
974 : END DO
975 :
976 52506272 : DO WHILE (jat > ref_stride)
977 22027460 : jat = jat - ref_stride
978 22027460 : jadr = jadr + 1
979 : END DO
980 30478812 : END SUBROUTINE limit
981 :
982 : ! **************************************************************************************************
983 : !> \brief ...
984 : !> \param rab ...
985 : !> \param rcutab ...
986 : !> \param srn ...
987 : !> \param alpn ...
988 : !> \param rcut ...
989 : !> \param fdab ...
990 : !> \param dfdab ...
991 : ! **************************************************************************************************
992 7841276 : SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
993 :
994 : REAL(KIND=dp), INTENT(IN) :: rab, rcutab, srn, alpn, rcut
995 : REAL(KIND=dp), INTENT(OUT) :: fdab, dfdab
996 :
997 : REAL(KIND=dp) :: dfab, dfcc, fab, fcc, rr
998 :
999 7841276 : CALL cutoff_switch_d3(rab, rcut, fcc, dfcc)
1000 :
1001 7841276 : rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
1002 7841276 : fab = 1._dp/(1._dp + rr)
1003 7841276 : dfab = fab*fab*rr*alpn/rab
1004 7841276 : fdab = fab*fcc
1005 7841276 : dfdab = dfab*fcc + fab*dfcc
1006 :
1007 7841276 : END SUBROUTINE damping_d3
1008 :
1009 : ! **************************************************************************************************
1010 : !> \brief Switching function of zero damping, from 1 at rcut - 1 bohr to 0 at rcut
1011 : !> \param rab distance
1012 : !> \param rcut cutoff
1013 : !> \param fcc switching function
1014 : !> \param dfcc its derivative with respect to rab
1015 : ! **************************************************************************************************
1016 7842108 : SUBROUTINE cutoff_switch_d3(rab, rcut, fcc, dfcc)
1017 :
1018 : REAL(KIND=dp), INTENT(IN) :: rab, rcut
1019 : REAL(KIND=dp), INTENT(OUT) :: fcc, dfcc
1020 :
1021 : REAL(KIND=dp) :: a, b, c, d, dd, dz, rl, ru, z, zz
1022 :
1023 7842108 : rl = rcut - 1._dp
1024 7842108 : ru = rcut
1025 7842108 : IF (rab >= ru) THEN
1026 0 : fcc = 0._dp
1027 0 : dfcc = 0._dp
1028 7842108 : ELSE IF (rab <= rl) THEN
1029 7128132 : fcc = 1._dp
1030 7128132 : dfcc = 0._dp
1031 : ELSE
1032 713976 : z = rab*rab - rl*rl
1033 713976 : dz = 2._dp*rab
1034 713976 : zz = z*z*z
1035 713976 : d = ru*ru - rl*rl
1036 713976 : dd = d*d*d
1037 713976 : a = -10._dp/dd
1038 713976 : b = 15._dp/(dd*d)
1039 713976 : c = -6._dp/(dd*d*d)
1040 713976 : fcc = 1._dp + zz*(a + b*z + c*z*z)
1041 713976 : dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
1042 : END IF
1043 :
1044 7842108 : END SUBROUTINE cutoff_switch_d3
1045 :
1046 : ! **************************************************************************************************
1047 : !> \brief Long-range correction of an r^-n term carried by the pairs that the zero-damping switch removes
1048 : !> from [rcut - 1 bohr, rcut], relative to the correction beyond rcut (g(r) = 1):
1049 : !> (n-3) rcut^(n-3) int_{rcut-1}^{rcut} r^(2-n) (1 - fcc(r)) dr
1050 : !> evaluated with a 16-point Gauss-Legendre rule (xgl: positive nodes on [-1, 1], wgl: their weights)
1051 : !> \param rcut cutoff
1052 : !> \param n power of the term
1053 : !> \return relative correction
1054 : ! **************************************************************************************************
1055 52 : FUNCTION switched_shell_fraction(rcut, n) RESULT(frac)
1056 :
1057 : REAL(KIND=dp), INTENT(IN) :: rcut
1058 : INTEGER, INTENT(IN) :: n
1059 : REAL(KIND=dp) :: frac
1060 :
1061 : REAL(KIND=dp), DIMENSION(8), PARAMETER :: wgl = [0.18945061045506849629_dp, &
1062 : 0.18260341504492358887_dp, 0.16915651939500253819_dp, 0.14959598881657673208_dp, &
1063 : 0.12462897125553387205_dp, 0.09515851168249278481_dp, 0.062253523938647892863_dp, &
1064 : 0.027152459411754094852_dp], xgl = [0.095012509837637440185_dp, 0.28160355077925891323_dp,&
1065 : 0.45801677765722738634_dp, 0.61787624440264374845_dp, 0.7554044083550030339_dp, &
1066 : 0.86563120238783174388_dp, 0.94457502307323257608_dp, 0.9894009349916499326_dp]
1067 :
1068 : INTEGER :: i, j
1069 : REAL(KIND=dp) :: dfcc, fcc, r
1070 :
1071 52 : frac = 0._dp
1072 468 : DO i = 1, 8
1073 52 : DO j = -1, 1, 2
1074 832 : r = rcut - 0.5_dp + 0.5_dp*REAL(j, KIND=dp)*xgl(i)
1075 832 : CALL cutoff_switch_d3(r, rcut, fcc, dfcc)
1076 832 : frac = frac + wgl(i)*r**(2 - n)*(1._dp - fcc)
1077 : END DO
1078 : END DO
1079 52 : frac = 0.5_dp*REAL(n - 3, KIND=dp)*rcut**(n - 3)*frac
1080 :
1081 52 : END FUNCTION switched_shell_fraction
1082 :
1083 : ! **************************************************************************************************
1084 : !> \brief ...
1085 : !> \param maxc ...
1086 : !> \param max_elem ...
1087 : !> \param c6ab ...
1088 : !> \param mxc ...
1089 : !> \param iat ...
1090 : !> \param jat ...
1091 : !> \param nci ...
1092 : !> \param ncj ...
1093 : !> \param k3 ...
1094 : !> \param c6 ...
1095 : !> \param dc6a ...
1096 : !> \param dc6b ...
1097 : ! **************************************************************************************************
1098 9163887 : SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
1099 :
1100 : INTEGER, INTENT(IN) :: maxc, max_elem
1101 : REAL(KIND=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
1102 : INTEGER, INTENT(IN) :: mxc(max_elem), iat, jat
1103 : REAL(KIND=dp), INTENT(IN) :: nci, ncj, k3
1104 : REAL(KIND=dp), INTENT(OUT) :: c6, dc6a, dc6b
1105 :
1106 : INTEGER :: i, j
1107 : REAL(KIND=dp) :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
1108 : dwa, dwb, dza, dzb, r, rsave, rsum, &
1109 : tmp1
1110 :
1111 : ! the exponential is sensitive to numerics
1112 : ! when nci or ncj is much larger than cn1/cn2
1113 :
1114 9163887 : c6mem = -1.0e+99_dp
1115 9163887 : rsave = 1.0e+99_dp
1116 9163887 : rsum = 0.0_dp
1117 9163887 : csum = 0.0_dp
1118 9163887 : dza = 0.0_dp
1119 9163887 : dzb = 0.0_dp
1120 9163887 : dwa = 0.0_dp
1121 9163887 : dwb = 0.0_dp
1122 9163887 : c6 = 0.0_dp
1123 33690384 : DO i = 1, mxc(iat)
1124 105276093 : DO j = 1, mxc(jat)
1125 71585709 : c6 = c6ab(iat, jat, i, j, 1)
1126 96112206 : IF (c6 > 0.0_dp) THEN
1127 71585709 : cn1 = c6ab(iat, jat, i, j, 2)
1128 71585709 : cn2 = c6ab(iat, jat, i, j, 3)
1129 : ! distance
1130 71585709 : r = (cn1 - nci)**2 + (cn2 - ncj)**2
1131 71585709 : IF (r < rsave) THEN
1132 35608282 : rsave = r
1133 35608282 : c6mem = c6
1134 : END IF
1135 71585709 : tmp1 = EXP(k3*r)
1136 71585709 : dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
1137 71585709 : dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
1138 71585709 : rsum = rsum + tmp1
1139 71585709 : csum = csum + tmp1*c6
1140 71585709 : dza = dza + dtmpa*c6
1141 71585709 : dwa = dwa + dtmpa
1142 71585709 : dzb = dzb + dtmpb*c6
1143 71585709 : dwb = dwb + dtmpb
1144 : END IF
1145 : END DO
1146 : END DO
1147 :
1148 9163887 : IF (c6 == 0.0_dp) c6mem = 0.0_dp
1149 :
1150 9163887 : IF (rsum > 1.0e-66_dp) THEN
1151 9149519 : c6 = csum/rsum
1152 9149519 : dc6a = (dza - c6*dwa)/rsum
1153 9149519 : dc6b = (dzb - c6*dwb)/rsum
1154 : ELSE
1155 14368 : c6 = c6mem
1156 14368 : dc6a = 0._dp
1157 14368 : dc6b = 0._dp
1158 : END IF
1159 :
1160 9163887 : END SUBROUTINE getc6
1161 :
1162 : END MODULE qs_dispersion_d3
|
