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 5560 : SUBROUTINE calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
78 5560 : 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 5560 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
93 : INTEGER, DIMENSION(3) :: cell_b, cell_c, ncell, periodic
94 5560 : 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 5560 : 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 5560 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
105 5560 : cnumfix, radd2
106 5560 : 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 5560 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
111 5560 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
112 : TYPE(atprop_type), POINTER :: atprop
113 : TYPE(cell_type), POINTER :: cell
114 5560 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
115 5560 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 : TYPE(neighbor_list_iterator_p_type), &
118 5560 : DIMENSION(:), POINTER :: nl_iterator
119 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
120 5560 : POINTER :: sab_vdw
121 5560 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
122 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a, disp_b, disp_c
123 5560 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
124 5560 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(virial_type), POINTER :: virial
126 :
127 5560 : CALL timeset(routineN, handle)
128 :
129 5560 : evdw = 0._dp
130 :
131 5560 : 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 5560 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, para_env=para_env, atprop=atprop)
135 :
136 5560 : debugall = dispersion_env%verbose
137 :
138 : ! atomic energy and stress arrays
139 5560 : atenergy = atprop%energy
140 5560 : 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 5560 : atex = .FALSE.
146 5560 : IF (PRESENT(atevdw)) THEN
147 2 : atex = .TRUE.
148 : END IF
149 :
150 5560 : NULLIFY (particle_set)
151 5560 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
152 5560 : natom = SIZE(particle_set)
153 :
154 5560 : NULLIFY (force)
155 5560 : CALL get_qs_env(qs_env=qs_env, force=force)
156 5560 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
157 5560 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
158 5560 : pv_virial_thread(:, :) = 0._dp
159 :
160 44480 : ALLOCATE (dodisp(nkind), exclude(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
161 18626 : DO ikind = 1, nkind
162 13066 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
163 13066 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
164 13066 : dodisp(ikind) = disp_a%defined
165 13066 : exclude(ikind) = ghost_a .OR. floating_a
166 13066 : atomnumber(ikind) = za
167 13066 : c6d2(ikind) = disp_a%c6
168 31692 : radd2(ikind) = disp_a%vdw_radii
169 : END DO
170 :
171 16680 : ALLOCATE (rcpbc(3, natom))
172 54090 : DO iatom = 1, natom
173 54090 : rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
174 : END DO
175 :
176 5560 : rcut = 2._dp*dispersion_env%rc_disp
177 5560 : rc2 = rcut*rcut
178 :
179 5560 : maxc = SIZE(dispersion_env%c6ab, 3)
180 5560 : max_elem = SIZE(dispersion_env%c6ab, 1)
181 5560 : alp6 = dispersion_env%alp
182 5560 : alp8 = alp6 + 2._dp
183 5560 : s6 = dispersion_env%s6
184 5560 : s8 = dispersion_env%s8
185 5560 : s9 = dispersion_env%s6
186 5560 : rs6 = dispersion_env%sr6
187 5560 : rs8 = 1._dp
188 5560 : a1 = dispersion_env%a1
189 5560 : a2 = dispersion_env%a2
190 5560 : eps_cn = dispersion_env%eps_cn
191 5560 : e6tot = 0._dp
192 5560 : e8tot = 0._dp
193 5560 : e9tot = 0._dp
194 5560 : esrb = 0._dp
195 5560 : domol = dispersion_env%domol
196 : ! molecule correction
197 5560 : kgc8 = dispersion_env%kgc8
198 5560 : 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 5560 : idmp = 0
209 5560 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
210 : idmp = 1
211 5288 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
212 5288 : idmp = 2
213 : END IF
214 : ! SRB parameters
215 5560 : ssrb = dispersion_env%srb_params(1)
216 5560 : gsrb = dispersion_env%srb_params(2)
217 5560 : t1srb = dispersion_env%srb_params(3)
218 5560 : t2srb = dispersion_env%srb_params(4)
219 :
220 5560 : 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 5560 : NULLIFY (particle_set)
244 5560 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
245 5560 : natom = SIZE(particle_set)
246 16680 : ALLOCATE (cnumbers(natom))
247 5560 : cnumbers = 0._dp
248 :
249 5560 : IF (calculate_forces .OR. debugall) THEN
250 12318 : ALLOCATE (dcnum(natom))
251 10690 : dcnum(:)%neighbors = 0
252 10690 : DO iatom = 1, natom
253 10690 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
254 : END DO
255 : ELSE
256 9492 : ALLOCATE (dcnum(1))
257 : END IF
258 :
259 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, exclude, atomnumber, &
260 5560 : calculate_forces, 1)
261 :
262 5560 : CALL para_env%sum(cnumbers)
263 :
264 : ! for parallel runs we have to update dcnum on all processors
265 5560 : IF (calculate_forces .OR. debugall) THEN
266 814 : CALL dcnum_distribute(dcnum, para_env)
267 814 : 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 5560 : nab = 0._dp
278 5560 : nabc = 0._dp
279 5560 : IF (dispersion_env%doabc) THEN
280 112 : rcc = 2._dp*dispersion_env%rc_disp
281 112 : CALL get_cell(cell=cell, periodic=periodic)
282 112 : sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
283 112 : sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
284 112 : sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
285 448 : ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
286 112 : 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 112 : IF (dispersion_env%c9cnst) THEN
294 66 : IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
295 198 : ALLOCATE (cnumfix(natom))
296 66 : cnumfix = 0._dp
297 : ! first use the default values
298 326 : DO iatom = 1, natom
299 260 : ikind = kind_of(iatom)
300 260 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
301 326 : cnumfix(iatom) = dispersion_env%cn(za)
302 : END DO
303 : ! now check for changes from default
304 66 : 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 66 : 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 66 : 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 5560 : sab_vdw => dispersion_env%sab_vdw
338 :
339 5560 : num_pe = 1
340 5560 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
341 :
342 5560 : mepos = 0
343 8888015 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
344 8882455 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
345 :
346 8882455 : IF (exclude(ikind) .OR. exclude(jkind)) CYCLE
347 :
348 8882398 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
349 :
350 8882164 : za = atomnumber(ikind)
351 8882164 : zb = atomnumber(jkind)
352 : ! vdW potential
353 35528656 : dr = SQRT(SUM(rij(:)**2))
354 8887724 : IF (dr <= rcut) THEN
355 8882164 : nab = nab + 1._dp
356 8882164 : fac = 1._dp
357 8882164 : IF (iatom == jatom) fac = 0.5_dp
358 8882164 : IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
359 8857596 : 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 8857364 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
367 8857364 : c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
368 8857364 : dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
369 8857364 : dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
370 8857364 : r6 = dr**6
371 8857364 : r8 = r6*dr*dr
372 8857364 : s8i = s8
373 8857364 : IF (domol) THEN
374 3568 : IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
375 1452 : s8i = kgc8
376 : END IF
377 : END IF
378 : ! damping
379 8857364 : IF (idmp == 1) THEN
380 : ! zero
381 3890544 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
382 3890544 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
383 3890544 : e6 = s6*fac*c6*fdab6/r6
384 3890544 : e8 = s8i*fac*c8*fdab8/r8
385 4966820 : ELSE IF (idmp == 2) THEN
386 : ! BJ
387 4966820 : r0ab = SQRT(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
388 4966820 : f0ab = a1*r0ab + a2
389 4966820 : fdab6 = 1.0_dp/(r6 + f0ab**6)
390 4966820 : fdab8 = 1.0_dp/(r8 + f0ab**8)
391 4966820 : e6 = s6*fac*c6*fdab6
392 4966820 : e8 = s8i*fac*c8*fdab8
393 : ELSE
394 0 : CPABORT("Unknown DFT-D3 damping function:")
395 : END IF
396 8857364 : 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 8857364 : evdw = evdw - e6 - e8
404 8857364 : e6tot = e6tot - e6
405 8857364 : e8tot = e8tot - e8
406 8857364 : 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 8857364 : 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 8857364 : IF (calculate_forces) THEN
415 : ! damping
416 3198770 : IF (idmp == 1) THEN
417 : ! zero
418 1999225 : de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
419 1999225 : de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
420 1199545 : ELSE IF (idmp == 2) THEN
421 : ! BJ
422 1199545 : de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
423 1199545 : de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
424 : ELSE
425 0 : CPABORT("Unknown DFT-D3 damping function:")
426 : END IF
427 12795080 : fdij(:) = (de6 + de8)*rij(:)/dr*fac
428 3198770 : 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 3198770 : atom_a = atom_of_kind(iatom)
432 3198770 : atom_b = atom_of_kind(jatom)
433 12795080 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
434 12795080 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
435 3198770 : IF (use_virial) THEN
436 2035489 : 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 3198770 : IF (idmp == 1) THEN
440 : ! zero
441 1999225 : dc6a = -s6*fac*dc6a*fdab6/r6
442 1999225 : dc6b = -s6*fac*dc6b*fdab6/r6
443 1999225 : dc8a = -s8i*fac*dc8a*fdab8/r8
444 1999225 : dc8b = -s8i*fac*dc8b*fdab8/r8
445 1199545 : ELSE IF (idmp == 2) THEN
446 : ! BJ
447 1199545 : dc6a = -s6*fac*dc6a*fdab6
448 1199545 : dc6b = -s6*fac*dc6b*fdab6
449 1199545 : dc8a = -s8i*fac*dc8a*fdab8
450 1199545 : dc8b = -s8i*fac*dc8b*fdab8
451 : ELSE
452 0 : CPABORT("Unknown DFT-D3 damping function:")
453 : END IF
454 43322651 : DO i = 1, dcnum(iatom)%neighbors
455 40123881 : katom = dcnum(iatom)%nlist(i)
456 40123881 : kkind = kind_of(katom)
457 160495524 : rik = dcnum(iatom)%rik(:, i)
458 160495524 : drk = SQRT(SUM(rik(:)**2))
459 160495524 : fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
460 40123881 : atom_c = atom_of_kind(katom)
461 160495524 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
462 160495524 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
463 43322651 : IF (use_virial) THEN
464 29100575 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
465 : END IF
466 : END DO
467 43321099 : DO i = 1, dcnum(jatom)%neighbors
468 40122329 : katom = dcnum(jatom)%nlist(i)
469 40122329 : kkind = kind_of(katom)
470 160489316 : rik = dcnum(jatom)%rik(:, i)
471 160489316 : drk = SQRT(SUM(rik(:)**2))
472 160489316 : fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
473 40122329 : atom_c = atom_of_kind(katom)
474 160489316 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
475 160489316 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
476 43321099 : IF (use_virial) THEN
477 29097827 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
478 : END IF
479 : END DO
480 : END IF
481 8857364 : IF (dispersion_env%doabc) THEN
482 16192 : CALL get_iterator_info(nl_iterator, cell=cell_b)
483 16192 : hashb = cellhash(cell_b, ncell)
484 25262 : is000 = (ALL(cell_b == 0))
485 259072 : rb0(:) = MATMUL(cell%hmat, cell_b)
486 16192 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
487 80960 : rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
488 64768 : r2ab = SUM((ra - rb)**2)
489 105328 : DO icx = -ncell(1), ncell(1)
490 665904 : DO icy = -ncell(2), ncell(2)
491 4632608 : DO icz = -ncell(3), ncell(3)
492 3982896 : cell_c(1) = icx
493 3982896 : cell_c(2) = icy
494 3982896 : cell_c(3) = icz
495 3982896 : hashc = cellhash(cell_c, ncell)
496 4648800 : IF (is000 .AND. (ALL(cell_c == 0))) THEN
497 : ! CASE 1: all atoms in (000), use only ordered triples
498 1206 : kstart = MAX(jatom + 1, iatom + 1)
499 1206 : fac0 = 1.0_dp
500 3981690 : 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 3948058 : IF (hashc == hashb) CYCLE
509 4582770 : 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 63246784 : rc0 = MATMUL(cell%hmat, cell_c)
515 15614318 : DO katom = kstart, natom
516 11100818 : kkind = kind_of(katom)
517 11100818 : IF (exclude(kkind) .OR. .NOT. dodisp(kkind)) CYCLE
518 44316664 : rc(:) = rcpbc(:, katom) + rc0(:)
519 44316664 : r2bc = SUM((rb - rc)**2)
520 11079166 : IF (r2bc >= rc2) CYCLE
521 9440716 : r2ca = SUM((rc - ra)**2)
522 2360179 : IF (r2ca >= rc2) CYCLE
523 1223808 : r2abc = r2ab*r2bc*r2ca
524 1223808 : IF (r2abc <= 0.001_dp) CYCLE
525 1223808 : 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 5173928 : IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
532 54030 : kkind = kind_of(katom)
533 54030 : atom_c = atom_of_kind(katom)
534 54030 : zc = atomnumber(kkind)
535 : ! avoid double counting!
536 54030 : fac = 1._dp
537 54030 : IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
538 54030 : IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
539 54030 : fac = fac*fac0
540 54030 : nabc = nabc + 1._dp
541 54030 : IF (dispersion_env%c9cnst) THEN
542 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
543 33752 : 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 33752 : 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 33752 : 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 20278 : 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 20278 : 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 20278 : cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
555 : END IF
556 54030 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
557 54030 : rabc = r2abc**(1._dp/6._dp)
558 : r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
559 54030 : 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 54030 : CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
563 54030 : s1 = r2ab + r2bc - r2ca
564 54030 : s2 = r2ab + r2ca - r2bc
565 54030 : s3 = r2ca + r2bc - r2ab
566 54030 : ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
567 :
568 54030 : e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
569 54030 : evdw = evdw - e9
570 54030 : e9tot = e9tot - e9
571 54030 : 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 54030 : 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 54030 : IF (calculate_forces) THEN
583 271500 : rab = rb - ra; rbc = rc - rb; rca = ra - rc
584 27150 : de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
585 27150 : de91 = -de91/3._dp*rabc + 3._dp*e9
586 27150 : dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
587 108600 : fdij(:) = de91*rab(:)/r2ab
588 108600 : fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
589 108600 : fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
590 108600 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
591 108600 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
592 27150 : IF (use_virial) THEN
593 16376 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
594 : END IF
595 108600 : fdij(:) = de91*rbc(:)/r2bc
596 108600 : fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
597 108600 : fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
598 108600 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
599 108600 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
600 27150 : IF (use_virial) THEN
601 16376 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
602 : END IF
603 108600 : fdij(:) = de91*rca(:)/r2ca
604 108600 : fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
605 108600 : fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
606 108600 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
607 108600 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
608 27150 : IF (use_virial) THEN
609 16376 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
610 : END IF
611 :
612 27150 : IF (.NOT. dispersion_env%c9cnst) THEN
613 : ! forces from the r-dependence of the coordination numbers
614 : ! atomic stress not implemented
615 :
616 13966 : de91 = 0.5_dp*e9/cc6ab
617 13966 : de921 = de91*dcc6aba
618 13966 : de922 = de91*dcc6abb
619 69354 : 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 69354 : IF (use_virial) THEN
631 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
632 : END IF
633 : END DO
634 69354 : 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 69354 : 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 13966 : de91 = 0.5_dp*e9/cc6bc
651 13966 : de921 = de91*dcc6bcb
652 13966 : de922 = de91*dcc6bcc
653 69354 : 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 69354 : IF (use_virial) THEN
665 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
666 : END IF
667 : END DO
668 69354 : 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 69354 : 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 13966 : de91 = 0.5_dp*e9/cc6ca
685 13966 : de921 = de91*dcc6cac
686 13966 : de922 = de91*dcc6caa
687 69354 : 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 69354 : IF (use_virial) THEN
699 55224 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
700 : END IF
701 : END DO
702 69354 : 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 69354 : 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 72280 : virial%pv_virial = virial%pv_virial + pv_virial_thread
733 :
734 5560 : CALL neighbor_list_iterator_release(nl_iterator)
735 :
736 5560 : elrc6 = 0._dp
737 5560 : elrc8 = 0._dp
738 5560 : elrc9 = 0._dp
739 : ! Long range correction (atomic contributions not implemented)
740 5560 : IF (dispersion_env%lrc) THEN
741 120 : ALLOCATE (cnkind(nkind))
742 40 : cnkind = 0._dp
743 : ! first use the default values
744 90 : DO ikind = 1, nkind
745 50 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
746 90 : cnkind(ikind) = dispersion_env%cn(za)
747 : END DO
748 : ! now check for changes from default
749 40 : IF (ASSOCIATED(dispersion_env%cnkind)) THEN
750 12 : DO i = 1, SIZE(dispersion_env%cnkind)
751 6 : ikind = dispersion_env%cnkind(i)%kind
752 12 : cnkind(ikind) = dispersion_env%cnkind(i)%cnum
753 : END DO
754 : END IF
755 90 : DO ikind = 1, nkind
756 50 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
757 50 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
758 50 : IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) CYCLE
759 204 : DO jkind = 1, nkind
760 66 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
761 66 : CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
762 66 : IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
763 64 : 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 58 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
770 58 : elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
771 58 : c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
772 58 : elrc8 = elrc8 - s8*twopi*REAL(na*nb, KIND=dp)*c8/(5._dp*rcut**5*cell%deth)
773 174 : IF (dispersion_env%doabc) THEN
774 144 : DO kkind = 1, nkind
775 86 : CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
776 86 : CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
777 86 : IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
778 84 : 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 82 : 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 82 : 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 82 : cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
789 82 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
790 226 : 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 40 : IF (use_virial) THEN
796 26 : IF (para_env%is_source()) THEN
797 52 : DO i = 1, 3
798 52 : virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
799 : END DO
800 : END IF
801 : END IF
802 40 : DEALLOCATE (cnkind)
803 : END IF
804 :
805 5560 : DEALLOCATE (cnumbers)
806 5560 : IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
807 66 : DEALLOCATE (cnumfix)
808 : END IF
809 5560 : IF (calculate_forces .OR. debugall) THEN
810 10690 : DO iatom = 1, natom
811 10690 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
812 : END DO
813 814 : DEALLOCATE (dcnum)
814 : ELSE
815 4746 : DEALLOCATE (dcnum)
816 : END IF
817 :
818 : ! set dispersion energy
819 :
820 5560 : CALL para_env%sum(e6tot)
821 5560 : CALL para_env%sum(e8tot)
822 5560 : CALL para_env%sum(e9tot)
823 5560 : CALL para_env%sum(evdw)
824 5560 : CALL para_env%sum(nab)
825 5560 : CALL para_env%sum(nabc)
826 5560 : e6tot = e6tot + elrc6
827 5560 : e8tot = e8tot + elrc8
828 5560 : e9tot = e9tot + elrc9
829 : ! For printing, we need all contributions
830 5560 : evdw = evdw + (elrc6 + elrc8 + elrc9)
831 5560 : IF (unit_nr > 0) THEN
832 45 : WRITE (unit_nr, "(A,F20.0)") " E6 vdW terms :", nab
833 45 : WRITE (unit_nr, *) " E6 vdW energy [au/kcal] :", e6tot, e6tot*kcalmol
834 45 : WRITE (unit_nr, *) " E8 vdW energy [au/kcal] :", e8tot, e8tot*kcalmol
835 45 : WRITE (unit_nr, *) " %E8 on total vdW energy :", e8tot/evdw*100._dp
836 45 : WRITE (unit_nr, "(A,F20.0)") " E9 vdW terms :", nabc
837 45 : WRITE (unit_nr, *) " E9 vdW energy [au/kcal] :", e9tot, e9tot*kcalmol
838 45 : WRITE (unit_nr, *) " %E9 on total vdW energy :", e9tot/evdw*100._dp
839 45 : IF (dispersion_env%lrc) THEN
840 1 : WRITE (unit_nr, *) " E LRC C6 [au/kcal] :", elrc6, elrc6*kcalmol
841 1 : WRITE (unit_nr, *) " E LRC C8 [au/kcal] :", elrc8, elrc8*kcalmol
842 1 : WRITE (unit_nr, *) " E LRC C9 [au/kcal] :", elrc9, elrc9*kcalmol
843 : END IF
844 : END IF
845 :
846 5560 : evdw = evdw/para_env%num_pe
847 :
848 5560 : DEALLOCATE (dodisp, exclude, atomnumber, rcpbc, radd2, c6d2)
849 :
850 5560 : IF (domol) THEN
851 2 : DEALLOCATE (atom2mol)
852 : END IF
853 :
854 5560 : CALL timestop(handle)
855 :
856 11120 : END SUBROUTINE calculate_dispersion_d3_pairpot
857 :
858 : ! **************************************************************************************************
859 : !> \brief ...
860 : !> \param c6ab ...
861 : !> \param maxci ...
862 : !> \param filename ...
863 : !> \param para_env ...
864 : ! **************************************************************************************************
865 466 : SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
866 :
867 : REAL(KIND=dp), DIMENSION(:, :, :, :, :) :: c6ab
868 : INTEGER, DIMENSION(:) :: maxci
869 : CHARACTER(LEN=*) :: filename
870 : TYPE(mp_para_env_type), POINTER :: para_env
871 :
872 : INTEGER :: funit, iadr, iat, jadr, jat, kk, nl, &
873 : nlines, nn, ref_stride
874 466 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pars
875 :
876 466 : IF (para_env%is_source()) THEN
877 : ! Read the DFT-D3 C6AB parameters from file "filename"
878 238 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
879 238 : READ (funit, *) nl, nlines
880 : END IF
881 466 : CALL para_env%bcast(nl)
882 466 : CALL para_env%bcast(nlines)
883 1398 : ALLOCATE (pars(nl))
884 466 : IF (para_env%is_source()) THEN
885 238 : READ (funit, *) pars(1:nl)
886 238 : CALL close_file(unit_number=funit)
887 : END IF
888 466 : CALL para_env%bcast(pars)
889 :
890 : ! Store C6AB coefficients in an array
891 733873576 : c6ab = -1._dp
892 48464 : maxci = 0
893 466 : ref_stride = get_ref_stride(pars, nlines)
894 466 : kk = 1
895 26698072 : DO nn = 1, nlines
896 26697606 : iat = NINT(pars(kk + 1))
897 26697606 : jat = NINT(pars(kk + 2))
898 26697606 : CALL limit(iat, jat, iadr, jadr, ref_stride)
899 26697606 : maxci(iat) = MAX(maxci(iat), iadr)
900 26697606 : maxci(jat) = MAX(maxci(jat), jadr)
901 26697606 : c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
902 26697606 : c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
903 26697606 : c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
904 :
905 26697606 : c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
906 26697606 : c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
907 26697606 : c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
908 26698072 : kk = (nn*5) + 1
909 : END DO
910 :
911 466 : DEALLOCATE (pars)
912 :
913 466 : END SUBROUTINE dftd3_c6_param
914 :
915 : ! **************************************************************************************************
916 : !> \brief Detects the reference-state encoding used in the DFT-D3 C6AB table.
917 : !> \param pars DFT-D3 C6AB parameter table.
918 : !> \param nlines Number of C6AB parameter lines.
919 : !> \return Reference-state stride.
920 : ! **************************************************************************************************
921 466 : FUNCTION get_ref_stride(pars, nlines) RESULT(ref_stride)
922 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pars
923 : INTEGER, INTENT(IN) :: nlines
924 : INTEGER :: ref_stride
925 :
926 : INTEGER :: kk, nn
927 :
928 466 : ref_stride = 100
929 466 : kk = 1
930 2496362 : DO nn = 1, nlines
931 2496362 : IF (NINT(pars(kk + 1)) > 1000 .OR. NINT(pars(kk + 2)) > 1000) THEN
932 : ref_stride = 1000
933 : EXIT
934 : END IF
935 2495896 : kk = (nn*5) + 1
936 : END DO
937 :
938 466 : END FUNCTION get_ref_stride
939 :
940 : ! **************************************************************************************************
941 : !> \brief Decodes element numbers and reference-state indices.
942 : !> \param iat Encoded first element number.
943 : !> \param jat Encoded second element number.
944 : !> \param iadr First reference-state index.
945 : !> \param jadr Second reference-state index.
946 : !> \param ref_stride Reference-state stride.
947 : ! **************************************************************************************************
948 26697606 : SUBROUTINE limit(iat, jat, iadr, jadr, ref_stride)
949 : INTEGER :: iat, jat, iadr, jadr, ref_stride
950 :
951 26697606 : iadr = 1
952 26697606 : jadr = 1
953 89391382 : DO WHILE (iat > ref_stride)
954 62693776 : iat = iat - ref_stride
955 89391382 : iadr = iadr + 1
956 : END DO
957 :
958 45992336 : DO WHILE (jat > ref_stride)
959 19294730 : jat = jat - ref_stride
960 19294730 : jadr = jadr + 1
961 : END DO
962 26697606 : END SUBROUTINE limit
963 :
964 : ! **************************************************************************************************
965 : !> \brief ...
966 : !> \param rab ...
967 : !> \param rcutab ...
968 : !> \param srn ...
969 : !> \param alpn ...
970 : !> \param rcut ...
971 : !> \param fdab ...
972 : !> \param dfdab ...
973 : ! **************************************************************************************************
974 7835118 : SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
975 :
976 : REAL(KIND=dp), INTENT(IN) :: rab, rcutab, srn, alpn, rcut
977 : REAL(KIND=dp), INTENT(OUT) :: fdab, dfdab
978 :
979 : REAL(KIND=dp) :: a, b, c, d, dd, dfab, dfcc, dz, fab, &
980 : fcc, rl, rr, ru, z, zz
981 :
982 7835118 : rl = rcut - 1._dp
983 7835118 : ru = rcut
984 7835118 : IF (rab >= ru) THEN
985 : fcc = 0._dp
986 : dfcc = 0._dp
987 7835118 : ELSEIF (rab <= rl) THEN
988 : fcc = 1._dp
989 : dfcc = 0._dp
990 : ELSE
991 712728 : z = rab*rab - rl*rl
992 712728 : dz = 2._dp*rab
993 712728 : zz = z*z*z
994 712728 : d = ru*ru - rl*rl
995 712728 : dd = d*d*d
996 712728 : a = -10._dp/dd
997 712728 : b = 15._dp/(dd*d)
998 712728 : c = -6._dp/(dd*d*d)
999 712728 : fcc = 1._dp + zz*(a + b*z + c*z*z)
1000 712728 : dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
1001 : END IF
1002 :
1003 7835118 : rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
1004 7835118 : fab = 1._dp/(1._dp + rr)
1005 7835118 : dfab = fab*fab*rr*alpn/rab
1006 7835118 : fdab = fab*fcc
1007 7835118 : dfdab = dfab*fcc + fab*dfcc
1008 :
1009 7835118 : END SUBROUTINE damping_d3
1010 :
1011 : ! **************************************************************************************************
1012 : !> \brief ...
1013 : !> \param maxc ...
1014 : !> \param max_elem ...
1015 : !> \param c6ab ...
1016 : !> \param mxc ...
1017 : !> \param iat ...
1018 : !> \param jat ...
1019 : !> \param nci ...
1020 : !> \param ncj ...
1021 : !> \param k3 ...
1022 : !> \param c6 ...
1023 : !> \param dc6a ...
1024 : !> \param dc6b ...
1025 : ! **************************************************************************************************
1026 9019758 : SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
1027 :
1028 : INTEGER, INTENT(IN) :: maxc, max_elem
1029 : REAL(KIND=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
1030 : INTEGER, INTENT(IN) :: mxc(max_elem), iat, jat
1031 : REAL(KIND=dp), INTENT(IN) :: nci, ncj, k3
1032 : REAL(KIND=dp), INTENT(OUT) :: c6, dc6a, dc6b
1033 :
1034 : INTEGER :: i, j
1035 : REAL(KIND=dp) :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
1036 : dwa, dwb, dza, dzb, r, rsave, rsum, &
1037 : tmp1
1038 :
1039 : ! the exponential is sensitive to numerics
1040 : ! when nci or ncj is much larger than cn1/cn2
1041 :
1042 9019758 : c6mem = -1.0e+99_dp
1043 9019758 : rsave = 1.0e+99_dp
1044 9019758 : rsum = 0.0_dp
1045 9019758 : csum = 0.0_dp
1046 9019758 : dza = 0.0_dp
1047 9019758 : dzb = 0.0_dp
1048 9019758 : dwa = 0.0_dp
1049 9019758 : dwb = 0.0_dp
1050 9019758 : c6 = 0.0_dp
1051 33117414 : DO i = 1, mxc(iat)
1052 103390894 : DO j = 1, mxc(jat)
1053 70273480 : c6 = c6ab(iat, jat, i, j, 1)
1054 94371136 : IF (c6 > 0.0_dp) THEN
1055 70273480 : cn1 = c6ab(iat, jat, i, j, 2)
1056 70273480 : cn2 = c6ab(iat, jat, i, j, 3)
1057 : ! distance
1058 70273480 : r = (cn1 - nci)**2 + (cn2 - ncj)**2
1059 70273480 : IF (r < rsave) THEN
1060 35061600 : rsave = r
1061 35061600 : c6mem = c6
1062 : END IF
1063 70273480 : tmp1 = EXP(k3*r)
1064 70273480 : dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
1065 70273480 : dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
1066 70273480 : rsum = rsum + tmp1
1067 70273480 : csum = csum + tmp1*c6
1068 70273480 : dza = dza + dtmpa*c6
1069 70273480 : dwa = dwa + dtmpa
1070 70273480 : dzb = dzb + dtmpb*c6
1071 70273480 : dwb = dwb + dtmpb
1072 : END IF
1073 : END DO
1074 : END DO
1075 :
1076 9019758 : IF (c6 == 0.0_dp) c6mem = 0.0_dp
1077 :
1078 9019758 : IF (rsum > 1.0e-66_dp) THEN
1079 9005390 : c6 = csum/rsum
1080 9005390 : dc6a = (dza - c6*dwa)/rsum
1081 9005390 : dc6b = (dzb - c6*dwb)/rsum
1082 : ELSE
1083 14368 : c6 = c6mem
1084 14368 : dc6a = 0._dp
1085 14368 : dc6b = 0._dp
1086 : END IF
1087 :
1088 9019758 : END SUBROUTINE getc6
1089 :
1090 : END MODULE qs_dispersion_d3
|
