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 5490 : SUBROUTINE calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
78 5490 : 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 5490 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
93 : INTEGER, DIMENSION(3) :: cell_b, cell_c, ncell, periodic
94 5490 : 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 5490 : 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 5490 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
105 5490 : cnumfix, radd2
106 5490 : 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 5490 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
111 5490 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
112 : TYPE(atprop_type), POINTER :: atprop
113 : TYPE(cell_type), POINTER :: cell
114 5490 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
115 5490 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 : TYPE(neighbor_list_iterator_p_type), &
118 5490 : DIMENSION(:), POINTER :: nl_iterator
119 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
120 5490 : POINTER :: sab_vdw
121 5490 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
122 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a, disp_b, disp_c
123 5490 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
124 5490 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(virial_type), POINTER :: virial
126 :
127 5490 : CALL timeset(routineN, handle)
128 :
129 5490 : evdw = 0._dp
130 :
131 5490 : 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 5490 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, para_env=para_env, atprop=atprop)
135 :
136 5490 : debugall = dispersion_env%verbose
137 :
138 : ! atomic energy and stress arrays
139 5490 : atenergy = atprop%energy
140 5490 : 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 5490 : atex = .FALSE.
146 5490 : IF (PRESENT(atevdw)) THEN
147 2 : atex = .TRUE.
148 : END IF
149 :
150 5490 : NULLIFY (particle_set)
151 5490 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
152 5490 : natom = SIZE(particle_set)
153 :
154 5490 : NULLIFY (force)
155 5490 : CALL get_qs_env(qs_env=qs_env, force=force)
156 5490 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
157 5490 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
158 5490 : pv_virial_thread(:, :) = 0._dp
159 :
160 43920 : ALLOCATE (dodisp(nkind), exclude(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
161 18446 : DO ikind = 1, nkind
162 12956 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
163 12956 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
164 12956 : dodisp(ikind) = disp_a%defined
165 12956 : exclude(ikind) = ghost_a .OR. floating_a
166 12956 : atomnumber(ikind) = za
167 12956 : c6d2(ikind) = disp_a%c6
168 31402 : radd2(ikind) = disp_a%vdw_radii
169 : END DO
170 :
171 16470 : ALLOCATE (rcpbc(3, natom))
172 53642 : DO iatom = 1, natom
173 53642 : rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
174 : END DO
175 :
176 5490 : rcut = 2._dp*dispersion_env%rc_disp
177 5490 : rc2 = rcut*rcut
178 :
179 5490 : maxc = SIZE(dispersion_env%c6ab, 3)
180 5490 : max_elem = SIZE(dispersion_env%c6ab, 1)
181 5490 : alp6 = dispersion_env%alp
182 5490 : alp8 = alp6 + 2._dp
183 5490 : s6 = dispersion_env%s6
184 5490 : s8 = dispersion_env%s8
185 5490 : s9 = dispersion_env%s6
186 5490 : rs6 = dispersion_env%sr6
187 5490 : rs8 = 1._dp
188 5490 : a1 = dispersion_env%a1
189 5490 : a2 = dispersion_env%a2
190 5490 : eps_cn = dispersion_env%eps_cn
191 5490 : e6tot = 0._dp
192 5490 : e8tot = 0._dp
193 5490 : e9tot = 0._dp
194 5490 : esrb = 0._dp
195 5490 : domol = dispersion_env%domol
196 : ! molecule correction
197 5490 : kgc8 = dispersion_env%kgc8
198 5490 : 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 5490 : idmp = 0
209 5490 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
210 : idmp = 1
211 5238 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
212 5238 : idmp = 2
213 : END IF
214 : ! SRB parameters
215 5490 : ssrb = dispersion_env%srb_params(1)
216 5490 : gsrb = dispersion_env%srb_params(2)
217 5490 : t1srb = dispersion_env%srb_params(3)
218 5490 : t2srb = dispersion_env%srb_params(4)
219 :
220 5490 : IF (unit_nr > 0) THEN
221 25 : WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
222 25 : WRITE (unit_nr, *) " Scaling parameter (s8) ", s8
223 25 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
224 12 : WRITE (unit_nr, *) " Zero Damping parameter (sr6)", rs6
225 12 : WRITE (unit_nr, *) " Zero Damping parameter (sr8)", rs8
226 13 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
227 13 : WRITE (unit_nr, *) " BJ Damping parameter (a1) ", a1
228 13 : WRITE (unit_nr, *) " BJ Damping parameter (a2) ", a2
229 : END IF
230 25 : WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
231 25 : IF (dispersion_env%lrc) THEN
232 1 : WRITE (unit_nr, *) " Apply a long range correction"
233 : END IF
234 25 : IF (dispersion_env%srb) THEN
235 2 : WRITE (unit_nr, *) " Apply a short range bond correction"
236 2 : WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
237 : END IF
238 25 : 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 5490 : NULLIFY (particle_set)
244 5490 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
245 5490 : natom = SIZE(particle_set)
246 16470 : ALLOCATE (cnumbers(natom))
247 5490 : cnumbers = 0._dp
248 :
249 5490 : IF (calculate_forces .OR. debugall) THEN
250 12114 : ALLOCATE (dcnum(natom))
251 10546 : dcnum(:)%neighbors = 0
252 10546 : DO iatom = 1, natom
253 10546 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
254 : END DO
255 : ELSE
256 9412 : ALLOCATE (dcnum(1))
257 : END IF
258 :
259 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, exclude, atomnumber, &
260 5490 : calculate_forces, 1)
261 :
262 5490 : CALL para_env%sum(cnumbers)
263 : ! for parallel runs we have to update dcnum on all processors
264 5490 : IF (calculate_forces .OR. debugall) THEN
265 784 : CALL dcnum_distribute(dcnum, para_env)
266 784 : IF (unit_nr > 0 .AND. SIZE(dcnum) > 0) THEN
267 3 : WRITE (unit_nr, *)
268 3 : WRITE (unit_nr, *) " ATOM Coordination Neighbors"
269 15 : DO i = 1, natom
270 15 : WRITE (unit_nr, "(I6,F20.10,I12)") i, cnumbers(i), dcnum(i)%neighbors
271 : END DO
272 3 : WRITE (unit_nr, *)
273 : END IF
274 : END IF
275 :
276 5490 : nab = 0._dp
277 5490 : nabc = 0._dp
278 5490 : IF (dispersion_env%doabc) THEN
279 80 : rcc = 2._dp*dispersion_env%rc_disp
280 80 : CALL get_cell(cell=cell, periodic=periodic)
281 80 : sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
282 80 : sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
283 80 : sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
284 320 : ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
285 80 : IF (unit_nr > 0) THEN
286 3 : WRITE (unit_nr, *) " Calculate C9 Terms"
287 3 : WRITE (unit_nr, "(A,T20,I3,A,I3)") " Search in cells ", -ncell(1), ":", ncell(1)
288 3 : WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(2), ":", ncell(2)
289 3 : WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(3), ":", ncell(3)
290 3 : WRITE (unit_nr, *)
291 : END IF
292 80 : IF (dispersion_env%c9cnst) THEN
293 46 : IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
294 138 : ALLOCATE (cnumfix(natom))
295 46 : cnumfix = 0._dp
296 : ! first use the default values
297 226 : DO iatom = 1, natom
298 180 : ikind = kind_of(iatom)
299 180 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
300 226 : cnumfix(iatom) = dispersion_env%cn(za)
301 : END DO
302 : ! now check for changes from default
303 46 : IF (ASSOCIATED(dispersion_env%cnkind)) THEN
304 12 : DO i = 1, SIZE(dispersion_env%cnkind)
305 6 : ikind = dispersion_env%cnkind(i)%kind
306 6 : cnum = dispersion_env%cnkind(i)%cnum
307 6 : CPASSERT(ikind <= nkind)
308 6 : CPASSERT(ikind > 0)
309 6 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, atom_list=atom_list)
310 32 : DO ia = 1, na
311 20 : iatom = atom_list(ia)
312 26 : cnumfix(iatom) = cnum
313 : END DO
314 : END DO
315 : END IF
316 46 : IF (ASSOCIATED(dispersion_env%cnlist)) THEN
317 6 : DO i = 1, SIZE(dispersion_env%cnlist)
318 14 : DO ilist = 1, dispersion_env%cnlist(i)%natom
319 8 : iatom = dispersion_env%cnlist(i)%atom(ilist)
320 12 : cnumfix(iatom) = dispersion_env%cnlist(i)%cnum
321 : END DO
322 : END DO
323 : END IF
324 46 : IF (unit_nr > 0) THEN
325 5 : DO i = 1, natom
326 5 : IF (ABS(cnumbers(i) - cnumfix(i)) > 0.5_dp) THEN
327 0 : WRITE (unit_nr, "(A,T20,A,I6,T41,2F10.3)") " Difference in CN ", "Atom:", &
328 0 : i, cnumbers(i), cnumfix(i)
329 : END IF
330 : END DO
331 1 : WRITE (unit_nr, *)
332 : END IF
333 : END IF
334 : END IF
335 :
336 5490 : sab_vdw => dispersion_env%sab_vdw
337 :
338 5490 : num_pe = 1
339 5490 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
340 :
341 5490 : mepos = 0
342 8868018 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
343 8862528 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
344 :
345 8862528 : IF (exclude(ikind) .OR. exclude(jkind)) CYCLE
346 :
347 8862471 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
348 :
349 8862237 : za = atomnumber(ikind)
350 8862237 : zb = atomnumber(jkind)
351 : ! vdW potential
352 35448948 : dr = SQRT(SUM(rij(:)**2))
353 8867727 : IF (dr <= rcut) THEN
354 8862237 : nab = nab + 1._dp
355 8862237 : fac = 1._dp
356 8862237 : IF (iatom == jatom) fac = 0.5_dp
357 8862237 : IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
358 8837858 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
359 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
360 320 : exclude=exclude_pair)
361 320 : IF (exclude_pair) CYCLE
362 : END IF
363 : ! C6 coefficient
364 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
365 8837626 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
366 8837626 : c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
367 8837626 : dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
368 8837626 : dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
369 8837626 : r6 = dr**6
370 8837626 : r8 = r6*dr*dr
371 8837626 : s8i = s8
372 8837626 : IF (domol) THEN
373 3568 : IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
374 1452 : s8i = kgc8
375 : END IF
376 : END IF
377 : ! damping
378 8837626 : IF (idmp == 1) THEN
379 : ! zero
380 3874798 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
381 3874798 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
382 3874798 : e6 = s6*fac*c6*fdab6/r6
383 3874798 : e8 = s8i*fac*c8*fdab8/r8
384 4962828 : ELSE IF (idmp == 2) THEN
385 : ! BJ
386 4962828 : r0ab = SQRT(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
387 4962828 : f0ab = a1*r0ab + a2
388 4962828 : fdab6 = 1.0_dp/(r6 + f0ab**6)
389 4962828 : fdab8 = 1.0_dp/(r8 + f0ab**8)
390 4962828 : e6 = s6*fac*c6*fdab6
391 4962828 : e8 = s8i*fac*c8*fdab8
392 : ELSE
393 0 : CPABORT("Unknown DFT-D3 damping function:")
394 : END IF
395 8837626 : IF (dispersion_env%srb .AND. dr < 30.0d0) THEN
396 69 : srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
397 69 : esrb = esrb + srbe
398 69 : evdw = evdw - srbe
399 : ELSE
400 : srbe = 0.0_dp
401 : END IF
402 8837626 : evdw = evdw - e6 - e8
403 8837626 : e6tot = e6tot - e6
404 8837626 : e8tot = e8tot - e8
405 8837626 : IF (atenergy) THEN
406 2765506 : atener(iatom) = atener(iatom) - 0.5_dp*(e6 + e8 + srbe)
407 2765506 : atener(jatom) = atener(jatom) - 0.5_dp*(e6 + e8 + srbe)
408 : END IF
409 8837626 : IF (atex) THEN
410 850 : atevdw(iatom) = atevdw(iatom) - 0.5_dp*(e6 + e8 + srbe)
411 850 : atevdw(jatom) = atevdw(jatom) - 0.5_dp*(e6 + e8 + srbe)
412 : END IF
413 8837626 : IF (calculate_forces) THEN
414 : ! damping
415 3195111 : IF (idmp == 1) THEN
416 : ! zero
417 1998786 : de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
418 1998786 : de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
419 1196325 : ELSE IF (idmp == 2) THEN
420 : ! BJ
421 1196325 : de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
422 1196325 : de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
423 : ELSE
424 0 : CPABORT("Unknown DFT-D3 damping function:")
425 : END IF
426 12780444 : fdij(:) = (de6 + de8)*rij(:)/dr*fac
427 3195111 : IF (dispersion_env%srb .AND. dr < 30.0d0) THEN
428 92 : fdij(:) = fdij(:) + srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
429 : END IF
430 3195111 : atom_a = atom_of_kind(iatom)
431 3195111 : atom_b = atom_of_kind(jatom)
432 12780444 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
433 12780444 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
434 3195111 : IF (use_virial) THEN
435 2033773 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
436 : END IF
437 : ! forces from the r-dependence of the coordination numbers
438 3195111 : IF (idmp == 1) THEN
439 : ! zero
440 1998786 : dc6a = -s6*fac*dc6a*fdab6/r6
441 1998786 : dc6b = -s6*fac*dc6b*fdab6/r6
442 1998786 : dc8a = -s8i*fac*dc8a*fdab8/r8
443 1998786 : dc8b = -s8i*fac*dc8b*fdab8/r8
444 1196325 : ELSE IF (idmp == 2) THEN
445 : ! BJ
446 1196325 : dc6a = -s6*fac*dc6a*fdab6
447 1196325 : dc6b = -s6*fac*dc6b*fdab6
448 1196325 : dc8a = -s8i*fac*dc8a*fdab8
449 1196325 : dc8b = -s8i*fac*dc8b*fdab8
450 : ELSE
451 0 : CPABORT("Unknown DFT-D3 damping function:")
452 : END IF
453 43316290 : DO i = 1, dcnum(iatom)%neighbors
454 40121179 : katom = dcnum(iatom)%nlist(i)
455 40121179 : kkind = kind_of(katom)
456 160484716 : rik = dcnum(iatom)%rik(:, i)
457 160484716 : drk = SQRT(SUM(rik(:)**2))
458 160484716 : fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
459 40121179 : atom_c = atom_of_kind(katom)
460 160484716 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
461 160484716 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
462 43316290 : IF (use_virial) THEN
463 29097959 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
464 : END IF
465 : END DO
466 43314738 : DO i = 1, dcnum(jatom)%neighbors
467 40119627 : katom = dcnum(jatom)%nlist(i)
468 40119627 : kkind = kind_of(katom)
469 160478508 : rik = dcnum(jatom)%rik(:, i)
470 160478508 : drk = SQRT(SUM(rik(:)**2))
471 160478508 : fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
472 40119627 : atom_c = atom_of_kind(katom)
473 160478508 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
474 160478508 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
475 43314738 : IF (use_virial) THEN
476 29095211 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
477 : END IF
478 : END DO
479 : END IF
480 8837626 : IF (dispersion_env%doabc) THEN
481 12212 : CALL get_iterator_info(nl_iterator, cell=cell_b)
482 12212 : hashb = cellhash(cell_b, ncell)
483 18774 : is000 = (ALL(cell_b == 0))
484 195392 : rb0(:) = MATMUL(cell%hmat, cell_b)
485 12212 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
486 61060 : rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
487 48848 : r2ab = SUM((ra - rb)**2)
488 80208 : DO icx = -ncell(1), ncell(1)
489 507684 : DO icy = -ncell(2), ncell(2)
490 3497708 : DO icz = -ncell(3), ncell(3)
491 3002236 : cell_c(1) = icx
492 3002236 : cell_c(2) = icy
493 3002236 : cell_c(3) = icz
494 3002236 : hashc = cellhash(cell_c, ncell)
495 3509920 : IF (is000 .AND. (ALL(cell_c == 0))) THEN
496 : ! CASE 1: all atoms in (000), use only ordered triples
497 802 : kstart = MAX(jatom + 1, iatom + 1)
498 802 : fac0 = 1.0_dp
499 3001434 : ELSE IF (is000) THEN
500 : ! CASE 2: AB in (000), C in other cell
501 : ! This case covers also all instances with BC in same
502 : ! cell not (000)
503 : kstart = 1
504 : fac0 = 1.0_dp
505 : ELSE
506 : ! These are case 2 again, cycle
507 2975242 : IF (hashc == hashb) CYCLE
508 3459498 : IF (ALL(cell_c == 0)) CYCLE
509 : ! CASE 3: A in (000) and B and C in different cells
510 : kstart = 1
511 : fac0 = 1.0_dp/3.0_dp
512 : END IF
513 47670656 : rc0 = MATMUL(cell%hmat, cell_c)
514 12351421 : DO katom = kstart, natom
515 8944529 : kkind = kind_of(katom)
516 8944529 : IF (exclude(kkind) .OR. .NOT. dodisp(kkind)) CYCLE
517 35691508 : rc(:) = rcpbc(:, katom) + rc0(:)
518 35691508 : r2bc = SUM((rb - rc)**2)
519 8922877 : IF (r2bc >= rc2) CYCLE
520 7408212 : r2ca = SUM((rc - ra)**2)
521 1852053 : IF (r2ca >= rc2) CYCLE
522 951043 : r2abc = r2ab*r2bc*r2ca
523 951043 : IF (r2abc <= 0.001_dp) CYCLE
524 951043 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
525 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
526 5066 : kkind, exclude_pair)
527 5066 : IF (exclude_pair) CYCLE
528 : END IF
529 : ! this is an empirical scaling
530 3927655 : IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
531 42663 : kkind = kind_of(katom)
532 42663 : atom_c = atom_of_kind(katom)
533 42663 : zc = atomnumber(kkind)
534 : ! avoid double counting!
535 42663 : fac = 1._dp
536 42663 : IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
537 42663 : IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
538 42663 : fac = fac*fac0
539 42663 : nabc = nabc + 1._dp
540 42663 : IF (dispersion_env%c9cnst) THEN
541 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
542 29152 : cnumfix(iatom), cnumfix(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
543 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
544 29152 : cnumfix(jatom), cnumfix(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
545 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
546 29152 : cnumfix(katom), cnumfix(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
547 : ELSE
548 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
549 13511 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
550 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
551 13511 : cnumbers(jatom), cnumbers(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
552 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
553 13511 : cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
554 : END IF
555 42663 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
556 42663 : rabc = r2abc**(1._dp/6._dp)
557 : r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
558 42663 : dispersion_env%r0ab(zc, za))**(1._dp/3._dp)
559 : ! bug fixed 3.10.2017
560 : ! correct value from alp6=14 to 16 as used in original paper
561 42663 : CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
562 42663 : s1 = r2ab + r2bc - r2ca
563 42663 : s2 = r2ab + r2ca - r2bc
564 42663 : s3 = r2ca + r2bc - r2ab
565 42663 : ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
566 :
567 42663 : e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
568 42663 : evdw = evdw - e9
569 42663 : e9tot = e9tot - e9
570 42663 : IF (atenergy) THEN
571 20568 : atener(iatom) = atener(iatom) - e9/3._dp
572 20568 : atener(jatom) = atener(jatom) - e9/3._dp
573 20568 : atener(katom) = atener(katom) - e9/3._dp
574 : END IF
575 42663 : IF (atex) THEN
576 10284 : atevdw(iatom) = atevdw(iatom) - e9/3._dp
577 10284 : atevdw(jatom) = atevdw(jatom) - e9/3._dp
578 10284 : atevdw(katom) = atevdw(katom) - e9/3._dp
579 : END IF
580 :
581 42663 : IF (calculate_forces) THEN
582 179270 : rab = rb - ra; rbc = rc - rb; rca = ra - rc
583 17927 : de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
584 17927 : de91 = -de91/3._dp*rabc + 3._dp*e9
585 17927 : dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
586 71708 : fdij(:) = de91*rab(:)/r2ab
587 71708 : fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
588 71708 : fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
589 71708 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
590 71708 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
591 17927 : IF (use_virial) THEN
592 10150 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
593 : END IF
594 71708 : fdij(:) = de91*rbc(:)/r2bc
595 71708 : fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
596 71708 : fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
597 71708 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
598 71708 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
599 17927 : IF (use_virial) THEN
600 10150 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
601 : END IF
602 71708 : fdij(:) = de91*rca(:)/r2ca
603 71708 : fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
604 71708 : fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
605 71708 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
606 71708 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
607 17927 : IF (use_virial) THEN
608 10150 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
609 : END IF
610 :
611 17927 : IF (.NOT. dispersion_env%c9cnst) THEN
612 : ! forces from the r-dependence of the coordination numbers
613 : ! atomic stress not implemented
614 :
615 9343 : de91 = 0.5_dp*e9/cc6ab
616 9343 : de921 = de91*dcc6aba
617 9343 : de922 = de91*dcc6abb
618 37037 : DO i = 1, dcnum(iatom)%neighbors
619 27694 : latom = dcnum(iatom)%nlist(i)
620 27694 : lkind = kind_of(latom)
621 27694 : rik(1) = dcnum(iatom)%rik(1, i)
622 27694 : rik(2) = dcnum(iatom)%rik(2, i)
623 27694 : rik(3) = dcnum(iatom)%rik(3, i)
624 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
625 110776 : fdik(:) = -de921*dcnum(iatom)%dvals(i)*rik(:)/drk
626 27694 : atom_d = atom_of_kind(latom)
627 110776 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
628 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
629 37037 : IF (use_virial) THEN
630 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
631 : END IF
632 : END DO
633 37037 : DO i = 1, dcnum(jatom)%neighbors
634 27694 : latom = dcnum(jatom)%nlist(i)
635 27694 : lkind = kind_of(latom)
636 27694 : rik(1) = dcnum(jatom)%rik(1, i)
637 27694 : rik(2) = dcnum(jatom)%rik(2, i)
638 27694 : rik(3) = dcnum(jatom)%rik(3, i)
639 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
640 110776 : fdik(:) = -de922*dcnum(jatom)%dvals(i)*rik(:)/drk
641 27694 : atom_d = atom_of_kind(latom)
642 110776 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
643 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
644 37037 : IF (use_virial) THEN
645 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
646 : END IF
647 : END DO
648 :
649 9343 : de91 = 0.5_dp*e9/cc6bc
650 9343 : de921 = de91*dcc6bcb
651 9343 : de922 = de91*dcc6bcc
652 37037 : DO i = 1, dcnum(jatom)%neighbors
653 27694 : latom = dcnum(jatom)%nlist(i)
654 27694 : lkind = kind_of(latom)
655 27694 : rik(1) = dcnum(jatom)%rik(1, i)
656 27694 : rik(2) = dcnum(jatom)%rik(2, i)
657 27694 : rik(3) = dcnum(jatom)%rik(3, i)
658 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
659 110776 : fdik(:) = -de921*dcnum(jatom)%dvals(i)*rik(:)/drk
660 27694 : atom_d = atom_of_kind(latom)
661 110776 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
662 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
663 37037 : IF (use_virial) THEN
664 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
665 : END IF
666 : END DO
667 37037 : DO i = 1, dcnum(katom)%neighbors
668 27694 : latom = dcnum(katom)%nlist(i)
669 27694 : lkind = kind_of(latom)
670 27694 : rik(1) = dcnum(katom)%rik(1, i)
671 27694 : rik(2) = dcnum(katom)%rik(2, i)
672 27694 : rik(3) = dcnum(katom)%rik(3, i)
673 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
674 110776 : fdik(:) = -de922*dcnum(katom)%dvals(i)*rik(:)/drk
675 27694 : atom_d = atom_of_kind(latom)
676 110776 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
677 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
678 37037 : IF (use_virial) THEN
679 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
680 : END IF
681 : END DO
682 :
683 9343 : de91 = 0.5_dp*e9/cc6ca
684 9343 : de921 = de91*dcc6cac
685 9343 : de922 = de91*dcc6caa
686 37037 : DO i = 1, dcnum(katom)%neighbors
687 27694 : latom = dcnum(katom)%nlist(i)
688 27694 : lkind = kind_of(latom)
689 27694 : rik(1) = dcnum(katom)%rik(1, i)
690 27694 : rik(2) = dcnum(katom)%rik(2, i)
691 27694 : rik(3) = dcnum(katom)%rik(3, i)
692 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
693 110776 : fdik(:) = -de921*dcnum(katom)%dvals(i)*rik(:)/drk
694 27694 : atom_d = atom_of_kind(latom)
695 110776 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
696 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
697 37037 : IF (use_virial) THEN
698 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
699 : END IF
700 : END DO
701 37037 : DO i = 1, dcnum(iatom)%neighbors
702 27694 : latom = dcnum(iatom)%nlist(i)
703 27694 : lkind = kind_of(latom)
704 27694 : rik(1) = dcnum(iatom)%rik(1, i)
705 27694 : rik(2) = dcnum(iatom)%rik(2, i)
706 27694 : rik(3) = dcnum(iatom)%rik(3, i)
707 27694 : drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
708 110776 : fdik(:) = -de922*dcnum(iatom)%dvals(i)*rik(:)/drk
709 27694 : atom_d = atom_of_kind(latom)
710 110776 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
711 110776 : force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
712 37037 : IF (use_virial) THEN
713 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
714 : END IF
715 : END DO
716 : END IF
717 :
718 : END IF
719 :
720 : END IF
721 : END DO
722 : END DO
723 : END DO
724 : END DO
725 :
726 : END IF
727 : END IF
728 : END IF
729 : END DO
730 :
731 71370 : virial%pv_virial = virial%pv_virial + pv_virial_thread
732 :
733 5490 : CALL neighbor_list_iterator_release(nl_iterator)
734 :
735 5490 : elrc6 = 0._dp
736 5490 : elrc8 = 0._dp
737 5490 : elrc9 = 0._dp
738 : ! Long range correction (atomic contributions not implemented)
739 5490 : IF (dispersion_env%lrc) THEN
740 96 : ALLOCATE (cnkind(nkind))
741 32 : cnkind = 0._dp
742 : ! first use the default values
743 74 : DO ikind = 1, nkind
744 42 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
745 74 : cnkind(ikind) = dispersion_env%cn(za)
746 : END DO
747 : ! now check for changes from default
748 32 : IF (ASSOCIATED(dispersion_env%cnkind)) THEN
749 12 : DO i = 1, SIZE(dispersion_env%cnkind)
750 6 : ikind = dispersion_env%cnkind(i)%kind
751 12 : cnkind(ikind) = dispersion_env%cnkind(i)%cnum
752 : END DO
753 : END IF
754 74 : DO ikind = 1, nkind
755 42 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
756 42 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
757 42 : IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) CYCLE
758 172 : DO jkind = 1, nkind
759 58 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
760 58 : CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
761 58 : IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
762 56 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
763 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
764 8 : exclude=exclude_pair)
765 8 : IF (exclude_pair) CYCLE
766 : END IF
767 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
768 50 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
769 50 : elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
770 50 : c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
771 50 : elrc8 = elrc8 - s8*twopi*REAL(na*nb, KIND=dp)*c8/(5._dp*rcut**5*cell%deth)
772 150 : IF (dispersion_env%doabc) THEN
773 128 : DO kkind = 1, nkind
774 78 : CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
775 78 : CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
776 78 : IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
777 76 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
778 : CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, kkind, &
779 4 : exclude_pair)
780 4 : IF (exclude_pair) CYCLE
781 : END IF
782 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
783 74 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
784 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
785 74 : cnkind(kkind), cnkind(ikind), dispersion_env%k3, cc6ca, dcc6aba, dcc6abb)
786 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
787 74 : cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
788 74 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
789 202 : elrc9 = elrc9 - s9*64._dp*twopi*REAL(na*nb*nc, KIND=dp)*c9/(6._dp*rcut**3*cell%deth**2)
790 : END DO
791 : END IF
792 : END DO
793 : END DO
794 32 : IF (use_virial) THEN
795 18 : IF (para_env%is_source()) THEN
796 36 : DO i = 1, 3
797 36 : virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
798 : END DO
799 : END IF
800 : END IF
801 32 : DEALLOCATE (cnkind)
802 : END IF
803 :
804 5490 : DEALLOCATE (cnumbers)
805 5490 : IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
806 46 : DEALLOCATE (cnumfix)
807 : END IF
808 5490 : IF (calculate_forces .OR. debugall) THEN
809 10546 : DO iatom = 1, natom
810 10546 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
811 : END DO
812 784 : DEALLOCATE (dcnum)
813 : ELSE
814 4706 : DEALLOCATE (dcnum)
815 : END IF
816 :
817 : ! set dispersion energy
818 :
819 5490 : CALL para_env%sum(e6tot)
820 5490 : CALL para_env%sum(e8tot)
821 5490 : CALL para_env%sum(e9tot)
822 5490 : CALL para_env%sum(evdw)
823 5490 : CALL para_env%sum(nab)
824 5490 : CALL para_env%sum(nabc)
825 5490 : e6tot = e6tot + elrc6
826 5490 : e8tot = e8tot + elrc8
827 5490 : e9tot = e9tot + elrc9
828 : ! For printing, we need all contributions
829 5490 : evdw = evdw + (elrc6 + elrc8 + elrc9)
830 5490 : IF (unit_nr > 0) THEN
831 25 : WRITE (unit_nr, "(A,F20.0)") " E6 vdW terms :", nab
832 25 : WRITE (unit_nr, *) " E6 vdW energy [au/kcal] :", e6tot, e6tot*kcalmol
833 25 : WRITE (unit_nr, *) " E8 vdW energy [au/kcal] :", e8tot, e8tot*kcalmol
834 25 : WRITE (unit_nr, *) " %E8 on total vdW energy :", e8tot/evdw*100._dp
835 25 : WRITE (unit_nr, "(A,F20.0)") " E9 vdW terms :", nabc
836 25 : WRITE (unit_nr, *) " E9 vdW energy [au/kcal] :", e9tot, e9tot*kcalmol
837 25 : WRITE (unit_nr, *) " %E9 on total vdW energy :", e9tot/evdw*100._dp
838 25 : IF (dispersion_env%lrc) THEN
839 1 : WRITE (unit_nr, *) " E LRC C6 [au/kcal] :", elrc6, elrc6*kcalmol
840 1 : WRITE (unit_nr, *) " E LRC C8 [au/kcal] :", elrc8, elrc8*kcalmol
841 1 : WRITE (unit_nr, *) " E LRC C9 [au/kcal] :", elrc9, elrc9*kcalmol
842 : END IF
843 : END IF
844 :
845 5490 : evdw = evdw/para_env%num_pe
846 :
847 5490 : DEALLOCATE (dodisp, exclude, atomnumber, rcpbc, radd2, c6d2)
848 :
849 5490 : IF (domol) THEN
850 2 : DEALLOCATE (atom2mol)
851 : END IF
852 :
853 5490 : CALL timestop(handle)
854 :
855 10980 : END SUBROUTINE calculate_dispersion_d3_pairpot
856 :
857 : ! **************************************************************************************************
858 : !> \brief ...
859 : !> \param c6ab ...
860 : !> \param maxci ...
861 : !> \param filename ...
862 : !> \param para_env ...
863 : ! **************************************************************************************************
864 466 : SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
865 :
866 : REAL(KIND=dp), DIMENSION(:, :, :, :, :) :: c6ab
867 : INTEGER, DIMENSION(:) :: maxci
868 : CHARACTER(LEN=*) :: filename
869 : TYPE(mp_para_env_type), POINTER :: para_env
870 :
871 : INTEGER :: funit, iadr, iat, jadr, jat, kk, nl, &
872 : nlines, nn, ref_stride
873 466 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pars
874 :
875 466 : IF (para_env%is_source()) THEN
876 : ! Read the DFT-D3 C6AB parameters from file "filename"
877 238 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
878 238 : READ (funit, *) nl, nlines
879 : END IF
880 466 : CALL para_env%bcast(nl)
881 466 : CALL para_env%bcast(nlines)
882 1398 : ALLOCATE (pars(nl))
883 466 : IF (para_env%is_source()) THEN
884 238 : READ (funit, *) pars(1:nl)
885 238 : CALL close_file(unit_number=funit)
886 : END IF
887 466 : CALL para_env%bcast(pars)
888 :
889 : ! Store C6AB coefficients in an array
890 733873576 : c6ab = -1._dp
891 48464 : maxci = 0
892 466 : ref_stride = get_ref_stride(pars, nlines)
893 466 : kk = 1
894 26698072 : DO nn = 1, nlines
895 26697606 : iat = NINT(pars(kk + 1))
896 26697606 : jat = NINT(pars(kk + 2))
897 26697606 : CALL limit(iat, jat, iadr, jadr, ref_stride)
898 26697606 : maxci(iat) = MAX(maxci(iat), iadr)
899 26697606 : maxci(jat) = MAX(maxci(jat), jadr)
900 26697606 : c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
901 26697606 : c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
902 26697606 : c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
903 :
904 26697606 : c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
905 26697606 : c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
906 26697606 : c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
907 26698072 : kk = (nn*5) + 1
908 : END DO
909 :
910 466 : DEALLOCATE (pars)
911 :
912 466 : END SUBROUTINE dftd3_c6_param
913 :
914 : ! **************************************************************************************************
915 : !> \brief Detects the reference-state encoding used in the DFT-D3 C6AB table.
916 : !> \param pars DFT-D3 C6AB parameter table.
917 : !> \param nlines Number of C6AB parameter lines.
918 : !> \return Reference-state stride.
919 : ! **************************************************************************************************
920 466 : FUNCTION get_ref_stride(pars, nlines) RESULT(ref_stride)
921 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pars
922 : INTEGER, INTENT(IN) :: nlines
923 : INTEGER :: ref_stride
924 :
925 : INTEGER :: kk, nn
926 :
927 466 : ref_stride = 100
928 466 : kk = 1
929 2496362 : DO nn = 1, nlines
930 2496362 : IF (NINT(pars(kk + 1)) > 1000 .OR. NINT(pars(kk + 2)) > 1000) THEN
931 : ref_stride = 1000
932 : EXIT
933 : END IF
934 2495896 : kk = (nn*5) + 1
935 : END DO
936 :
937 466 : END FUNCTION get_ref_stride
938 :
939 : ! **************************************************************************************************
940 : !> \brief Decodes element numbers and reference-state indices.
941 : !> \param iat Encoded first element number.
942 : !> \param jat Encoded second element number.
943 : !> \param iadr First reference-state index.
944 : !> \param jadr Second reference-state index.
945 : !> \param ref_stride Reference-state stride.
946 : ! **************************************************************************************************
947 26697606 : SUBROUTINE limit(iat, jat, iadr, jadr, ref_stride)
948 : INTEGER :: iat, jat, iadr, jadr, ref_stride
949 :
950 26697606 : iadr = 1
951 26697606 : jadr = 1
952 89391382 : DO WHILE (iat > ref_stride)
953 62693776 : iat = iat - ref_stride
954 89391382 : iadr = iadr + 1
955 : END DO
956 :
957 45992336 : DO WHILE (jat > ref_stride)
958 19294730 : jat = jat - ref_stride
959 19294730 : jadr = jadr + 1
960 : END DO
961 26697606 : END SUBROUTINE limit
962 :
963 : ! **************************************************************************************************
964 : !> \brief ...
965 : !> \param rab ...
966 : !> \param rcutab ...
967 : !> \param srn ...
968 : !> \param alpn ...
969 : !> \param rcut ...
970 : !> \param fdab ...
971 : !> \param dfdab ...
972 : ! **************************************************************************************************
973 7792259 : SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
974 :
975 : REAL(KIND=dp), INTENT(IN) :: rab, rcutab, srn, alpn, rcut
976 : REAL(KIND=dp), INTENT(OUT) :: fdab, dfdab
977 :
978 : REAL(KIND=dp) :: a, b, c, d, dd, dfab, dfcc, dz, fab, &
979 : fcc, rl, rr, ru, z, zz
980 :
981 7792259 : rl = rcut - 1._dp
982 7792259 : ru = rcut
983 7792259 : IF (rab >= ru) THEN
984 : fcc = 0._dp
985 : dfcc = 0._dp
986 7792259 : ELSEIF (rab <= rl) THEN
987 : fcc = 1._dp
988 : dfcc = 0._dp
989 : ELSE
990 711180 : z = rab*rab - rl*rl
991 711180 : dz = 2._dp*rab
992 711180 : zz = z*z*z
993 711180 : d = ru*ru - rl*rl
994 711180 : dd = d*d*d
995 711180 : a = -10._dp/dd
996 711180 : b = 15._dp/(dd*d)
997 711180 : c = -6._dp/(dd*d*d)
998 711180 : fcc = 1._dp + zz*(a + b*z + c*z*z)
999 711180 : dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
1000 : END IF
1001 :
1002 7792259 : rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
1003 7792259 : fab = 1._dp/(1._dp + rr)
1004 7792259 : dfab = fab*fab*rr*alpn/rab
1005 7792259 : fdab = fab*fcc
1006 7792259 : dfdab = dfab*fcc + fab*dfcc
1007 :
1008 7792259 : END SUBROUTINE damping_d3
1009 :
1010 : ! **************************************************************************************************
1011 : !> \brief ...
1012 : !> \param maxc ...
1013 : !> \param max_elem ...
1014 : !> \param c6ab ...
1015 : !> \param mxc ...
1016 : !> \param iat ...
1017 : !> \param jat ...
1018 : !> \param nci ...
1019 : !> \param ncj ...
1020 : !> \param k3 ...
1021 : !> \param c6 ...
1022 : !> \param dc6a ...
1023 : !> \param dc6b ...
1024 : ! **************************************************************************************************
1025 8965887 : SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
1026 :
1027 : INTEGER, INTENT(IN) :: maxc, max_elem
1028 : REAL(KIND=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
1029 : INTEGER, INTENT(IN) :: mxc(max_elem), iat, jat
1030 : REAL(KIND=dp), INTENT(IN) :: nci, ncj, k3
1031 : REAL(KIND=dp), INTENT(OUT) :: c6, dc6a, dc6b
1032 :
1033 : INTEGER :: i, j
1034 : REAL(KIND=dp) :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
1035 : dwa, dwb, dza, dzb, r, rsave, rsum, &
1036 : tmp1
1037 :
1038 : ! the exponential is sensitive to numerics
1039 : ! when nci or ncj is much larger than cn1/cn2
1040 :
1041 8965887 : c6mem = -1.0e+99_dp
1042 8965887 : rsave = 1.0e+99_dp
1043 8965887 : rsum = 0.0_dp
1044 8965887 : csum = 0.0_dp
1045 8965887 : dza = 0.0_dp
1046 8965887 : dzb = 0.0_dp
1047 8965887 : dwa = 0.0_dp
1048 8965887 : dwb = 0.0_dp
1049 8965887 : c6 = 0.0_dp
1050 32972158 : DO i = 1, mxc(iat)
1051 103044297 : DO j = 1, mxc(jat)
1052 70072139 : c6 = c6ab(iat, jat, i, j, 1)
1053 94078410 : IF (c6 > 0.0_dp) THEN
1054 70072139 : cn1 = c6ab(iat, jat, i, j, 2)
1055 70072139 : cn2 = c6ab(iat, jat, i, j, 3)
1056 : ! distance
1057 70072139 : r = (cn1 - nci)**2 + (cn2 - ncj)**2
1058 70072139 : IF (r < rsave) THEN
1059 34956258 : rsave = r
1060 34956258 : c6mem = c6
1061 : END IF
1062 70072139 : tmp1 = EXP(k3*r)
1063 70072139 : dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
1064 70072139 : dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
1065 70072139 : rsum = rsum + tmp1
1066 70072139 : csum = csum + tmp1*c6
1067 70072139 : dza = dza + dtmpa*c6
1068 70072139 : dwa = dwa + dtmpa
1069 70072139 : dzb = dzb + dtmpb*c6
1070 70072139 : dwb = dwb + dtmpb
1071 : END IF
1072 : END DO
1073 : END DO
1074 :
1075 8965887 : IF (c6 == 0.0_dp) c6mem = 0.0_dp
1076 :
1077 8965887 : IF (rsum > 1.0e-66_dp) THEN
1078 8951519 : c6 = csum/rsum
1079 8951519 : dc6a = (dza - c6*dwa)/rsum
1080 8951519 : dc6b = (dzb - c6*dwb)/rsum
1081 : ELSE
1082 14368 : c6 = c6mem
1083 14368 : dc6a = 0._dp
1084 14368 : dc6b = 0._dp
1085 : END IF
1086 :
1087 8965887 : END SUBROUTINE getc6
1088 :
1089 : END MODULE qs_dispersion_d3
|
