Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 4274 : SUBROUTINE calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
78 4274 : 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 4274 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
93 : INTEGER, DIMENSION(3) :: cell_b, cell_c, ncell, periodic
94 4274 : 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 4274 : 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 4274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
105 4274 : cnumfix, radd2
106 4274 : 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 4274 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
111 4274 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
112 : TYPE(atprop_type), POINTER :: atprop
113 : TYPE(cell_type), POINTER :: cell
114 4274 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
115 4274 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 : TYPE(neighbor_list_iterator_p_type), &
118 4274 : DIMENSION(:), POINTER :: nl_iterator
119 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
120 4274 : POINTER :: sab_vdw
121 4274 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
122 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a, disp_b, disp_c
123 4274 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
124 4274 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(virial_type), POINTER :: virial
126 :
127 4274 : CALL timeset(routineN, handle)
128 :
129 4274 : evdw = 0._dp
130 :
131 4274 : 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 4274 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, para_env=para_env, atprop=atprop)
135 :
136 4274 : debugall = dispersion_env%verbose
137 :
138 : ! atomic energy and stress arrays
139 4274 : atenergy = atprop%energy
140 4274 : 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 4274 : atex = .FALSE.
146 4274 : IF (PRESENT(atevdw)) THEN
147 2 : atex = .TRUE.
148 : END IF
149 :
150 4274 : NULLIFY (particle_set)
151 4274 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
152 4274 : natom = SIZE(particle_set)
153 :
154 4274 : NULLIFY (force)
155 4274 : CALL get_qs_env(qs_env=qs_env, force=force)
156 4274 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
157 4274 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
158 4274 : pv_virial_thread(:, :) = 0._dp
159 :
160 34192 : ALLOCATE (dodisp(nkind), exclude(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
161 14730 : DO ikind = 1, nkind
162 10456 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
163 10456 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
164 10456 : dodisp(ikind) = disp_a%defined
165 10456 : exclude(ikind) = ghost_a .OR. floating_a
166 10456 : atomnumber(ikind) = za
167 10456 : c6d2(ikind) = disp_a%c6
168 25186 : radd2(ikind) = disp_a%vdw_radii
169 : END DO
170 :
171 12822 : ALLOCATE (rcpbc(3, natom))
172 46662 : DO iatom = 1, natom
173 46662 : rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
174 : END DO
175 :
176 4274 : rcut = 2._dp*dispersion_env%rc_disp
177 4274 : rc2 = rcut*rcut
178 :
179 4274 : maxc = SIZE(dispersion_env%c6ab, 3)
180 4274 : max_elem = SIZE(dispersion_env%c6ab, 1)
181 4274 : alp6 = dispersion_env%alp
182 4274 : alp8 = alp6 + 2._dp
183 4274 : s6 = dispersion_env%s6
184 4274 : s8 = dispersion_env%s8
185 4274 : s9 = dispersion_env%s6
186 4274 : rs6 = dispersion_env%sr6
187 4274 : rs8 = 1._dp
188 4274 : a1 = dispersion_env%a1
189 4274 : a2 = dispersion_env%a2
190 4274 : eps_cn = dispersion_env%eps_cn
191 4274 : e6tot = 0._dp
192 4274 : e8tot = 0._dp
193 4274 : e9tot = 0._dp
194 4274 : esrb = 0._dp
195 4274 : domol = dispersion_env%domol
196 : ! molecule correction
197 4274 : kgc8 = dispersion_env%kgc8
198 4274 : 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 4274 : idmp = 0
209 4274 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
210 : idmp = 1
211 4036 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
212 4036 : idmp = 2
213 : END IF
214 : ! SRB parameters
215 4274 : ssrb = dispersion_env%srb_params(1)
216 4274 : gsrb = dispersion_env%srb_params(2)
217 4274 : t1srb = dispersion_env%srb_params(3)
218 4274 : t2srb = dispersion_env%srb_params(4)
219 :
220 4274 : IF (unit_nr > 0) THEN
221 6 : WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
222 6 : WRITE (unit_nr, *) " Scaling parameter (s8) ", s8
223 6 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
224 4 : WRITE (unit_nr, *) " Zero Damping parameter (sr6)", rs6
225 4 : WRITE (unit_nr, *) " Zero Damping parameter (sr8)", rs8
226 2 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
227 2 : WRITE (unit_nr, *) " BJ Damping parameter (a1) ", a1
228 2 : WRITE (unit_nr, *) " BJ Damping parameter (a2) ", a2
229 : END IF
230 6 : WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
231 6 : IF (dispersion_env%lrc) THEN
232 1 : WRITE (unit_nr, *) " Apply a long range correction"
233 : END IF
234 6 : IF (dispersion_env%srb) THEN
235 0 : WRITE (unit_nr, *) " Apply a short range bond correction"
236 0 : WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
237 : END IF
238 6 : 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 4274 : NULLIFY (particle_set)
244 4274 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
245 4274 : natom = SIZE(particle_set)
246 12822 : ALLOCATE (cnumbers(natom))
247 46662 : cnumbers = 0._dp
248 :
249 4274 : IF (calculate_forces .OR. debugall) THEN
250 11044 : ALLOCATE (dcnum(natom))
251 9712 : dcnum(:)%neighbors = 0
252 9712 : DO iatom = 1, natom
253 9712 : ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
254 : END DO
255 : ELSE
256 7216 : ALLOCATE (dcnum(1))
257 : END IF
258 :
259 : CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, exclude, atomnumber, &
260 4274 : calculate_forces, 1)
261 :
262 4274 : CALL para_env%sum(cnumbers)
263 : ! for parallel runs we have to update dcnum on all processors
264 4274 : IF (calculate_forces .OR. debugall) THEN
265 666 : CALL dcnum_distribute(dcnum, para_env)
266 666 : IF (unit_nr > 0 .AND. SIZE(dcnum) > 0) THEN
267 2 : WRITE (unit_nr, *)
268 2 : WRITE (unit_nr, *) " ATOM Coordination Neighbors"
269 11 : DO i = 1, natom
270 11 : WRITE (unit_nr, "(I6,F20.10,I12)") i, cnumbers(i), dcnum(i)%neighbors
271 : END DO
272 2 : WRITE (unit_nr, *)
273 : END IF
274 : END IF
275 :
276 4274 : nab = 0._dp
277 4274 : nabc = 0._dp
278 4274 : IF (dispersion_env%doabc) THEN
279 104 : rcc = 2._dp*dispersion_env%rc_disp
280 104 : CALL get_cell(cell=cell, periodic=periodic)
281 104 : sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
282 104 : sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
283 104 : sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
284 416 : ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
285 104 : 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 104 : IF (dispersion_env%c9cnst) THEN
293 62 : IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
294 186 : ALLOCATE (cnumfix(natom))
295 306 : cnumfix = 0._dp
296 : ! first use the default values
297 306 : DO iatom = 1, natom
298 244 : ikind = kind_of(iatom)
299 244 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
300 306 : cnumfix(iatom) = dispersion_env%cn(za)
301 : END DO
302 : ! now check for changes from default
303 62 : 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 62 : 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 62 : 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 4274 : sab_vdw => dispersion_env%sab_vdw
337 :
338 4274 : num_pe = 1
339 4274 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
340 :
341 4274 : mepos = 0
342 8206895 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
343 8202621 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
344 :
345 8202621 : IF (exclude(ikind) .OR. exclude(jkind)) CYCLE
346 :
347 8202564 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
348 :
349 8202330 : za = atomnumber(ikind)
350 8202330 : zb = atomnumber(jkind)
351 : ! vdW potential
352 32809320 : dr = SQRT(SUM(rij(:)**2))
353 8206604 : IF (dr <= rcut) THEN
354 8202330 : nab = nab + 1._dp
355 8202330 : fac = 1._dp
356 8202330 : IF (iatom == jatom) fac = 0.5_dp
357 8202330 : IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
358 8180833 : 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 8180601 : cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
366 8180601 : c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
367 8180601 : dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
368 8180601 : dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
369 8180601 : r6 = dr**6
370 8180601 : r8 = r6*dr*dr
371 8180601 : s8i = s8
372 8180601 : IF (domol) THEN
373 3568 : IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
374 1452 : s8i = kgc8
375 : END IF
376 : END IF
377 : ! damping
378 8180601 : IF (idmp == 1) THEN
379 : ! zero
380 3857247 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
381 3857247 : CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
382 3857247 : e6 = s6*fac*c6*fdab6/r6
383 3857247 : e8 = s8i*fac*c8*fdab8/r8
384 4323354 : ELSE IF (idmp == 2) THEN
385 : ! BJ
386 4323354 : r0ab = SQRT(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
387 4323354 : f0ab = a1*r0ab + a2
388 4323354 : fdab6 = 1.0_dp/(r6 + f0ab**6)
389 4323354 : fdab8 = 1.0_dp/(r8 + f0ab**8)
390 4323354 : e6 = s6*fac*c6*fdab6
391 4323354 : e8 = s8i*fac*c8*fdab8
392 : ELSE
393 0 : CPABORT("Unknown DFT-D3 damping function:")
394 : END IF
395 8180601 : IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
396 65 : srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
397 65 : esrb = esrb + srbe
398 65 : evdw = evdw - srbe
399 : ELSE
400 : srbe = 0.0_dp
401 : END IF
402 8180601 : evdw = evdw - e6 - e8
403 8180601 : e6tot = e6tot - e6
404 8180601 : e8tot = e8tot - e8
405 8180601 : IF (atenergy) THEN
406 2765504 : atener(iatom) = atener(iatom) - 0.5_dp*(e6 + e8 + srbe)
407 2765504 : atener(jatom) = atener(jatom) - 0.5_dp*(e6 + e8 + srbe)
408 : END IF
409 8180601 : 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 8180601 : IF (calculate_forces) THEN
414 : ! damping
415 2980475 : IF (idmp == 1) THEN
416 : ! zero
417 1983165 : de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
418 1983165 : de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
419 997310 : ELSE IF (idmp == 2) THEN
420 : ! BJ
421 997310 : de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
422 997310 : de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
423 : ELSE
424 0 : CPABORT("Unknown DFT-D3 damping function:")
425 : END IF
426 11921900 : fdij(:) = (de6 + de8)*rij(:)/dr*fac
427 2980475 : IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
428 80 : fdij(:) = fdij(:) + srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
429 : END IF
430 2980475 : atom_a = atom_of_kind(iatom)
431 2980475 : atom_b = atom_of_kind(jatom)
432 11921900 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
433 11921900 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
434 2980475 : IF (use_virial) THEN
435 1821413 : 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 2980475 : IF (idmp == 1) THEN
439 : ! zero
440 1983165 : dc6a = -s6*fac*dc6a*fdab6/r6
441 1983165 : dc6b = -s6*fac*dc6b*fdab6/r6
442 1983165 : dc8a = -s8i*fac*dc8a*fdab8/r8
443 1983165 : dc8b = -s8i*fac*dc8b*fdab8/r8
444 997310 : ELSE IF (idmp == 2) THEN
445 : ! BJ
446 997310 : dc6a = -s6*fac*dc6a*fdab6
447 997310 : dc6b = -s6*fac*dc6b*fdab6
448 997310 : dc8a = -s8i*fac*dc8a*fdab8
449 997310 : dc8b = -s8i*fac*dc8b*fdab8
450 : ELSE
451 0 : CPABORT("Unknown DFT-D3 damping function:")
452 : END IF
453 32696183 : DO i = 1, dcnum(iatom)%neighbors
454 29715708 : katom = dcnum(iatom)%nlist(i)
455 29715708 : kkind = kind_of(katom)
456 118862832 : rik = dcnum(iatom)%rik(:, i)
457 118862832 : drk = SQRT(SUM(rik(:)**2))
458 118862832 : fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
459 29715708 : atom_c = atom_of_kind(katom)
460 118862832 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
461 118862832 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
462 32696183 : IF (use_virial) THEN
463 18703635 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
464 : END IF
465 : END DO
466 32694631 : DO i = 1, dcnum(jatom)%neighbors
467 29714156 : katom = dcnum(jatom)%nlist(i)
468 29714156 : kkind = kind_of(katom)
469 118856624 : rik = dcnum(jatom)%rik(:, i)
470 118856624 : drk = SQRT(SUM(rik(:)**2))
471 118856624 : fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
472 29714156 : atom_c = atom_of_kind(katom)
473 118856624 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
474 118856624 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
475 32694631 : IF (use_virial) THEN
476 18700887 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
477 : END IF
478 : END DO
479 : END IF
480 8180601 : IF (dispersion_env%doabc) THEN
481 16052 : CALL get_iterator_info(nl_iterator, cell=cell_b)
482 16052 : hashb = cellhash(cell_b, ncell)
483 24294 : is000 = (ALL(cell_b == 0))
484 256832 : rb0(:) = MATMUL(cell%hmat, cell_b)
485 16052 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
486 80260 : rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
487 64208 : r2ab = SUM((ra - rb)**2)
488 103248 : DO icx = -ncell(1), ncell(1)
489 626724 : DO icy = -ncell(2), ncell(2)
490 4092908 : DO icz = -ncell(3), ncell(3)
491 3482236 : cell_c(1) = icx
492 3482236 : cell_c(2) = icy
493 3482236 : cell_c(3) = icz
494 3482236 : hashc = cellhash(cell_c, ncell)
495 4108960 : IF (is000 .AND. (ALL(cell_c == 0))) THEN
496 : ! CASE 1: all atoms in (000), use only ordered triples
497 874 : kstart = MAX(jatom + 1, iatom + 1)
498 874 : fac0 = 1.0_dp
499 3481362 : 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 3446242 : IF (hashc == hashb) CYCLE
508 4042074 : 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 55230080 : rc0 = MATMUL(cell%hmat, cell_c)
514 14809501 : DO katom = kstart, natom
515 10834145 : kkind = kind_of(katom)
516 10834145 : IF (exclude(kkind) .OR. .NOT. dodisp(kkind)) CYCLE
517 43249972 : rc(:) = rcpbc(:, katom) + rc0(:)
518 43249972 : r2bc = SUM((rb - rc)**2)
519 10812493 : IF (r2bc >= rc2) CYCLE
520 9114324 : r2ca = SUM((rc - ra)**2)
521 2278581 : IF (r2ca >= rc2) CYCLE
522 1168819 : r2abc = r2ab*r2bc*r2ca
523 1168819 : IF (r2abc <= 0.001_dp) CYCLE
524 1168819 : 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 4617895 : IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
531 47775 : kkind = kind_of(katom)
532 47775 : atom_c = atom_of_kind(katom)
533 47775 : zc = atomnumber(kkind)
534 : ! avoid double counting!
535 47775 : fac = 1._dp
536 47775 : IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
537 47775 : IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
538 47775 : fac = fac*fac0
539 47775 : nabc = nabc + 1._dp
540 47775 : IF (dispersion_env%c9cnst) THEN
541 : CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
542 32520 : 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 32520 : 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 32520 : 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 15255 : 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 15255 : 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 15255 : cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
554 : END IF
555 47775 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
556 47775 : rabc = r2abc**(1._dp/6._dp)
557 : r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
558 47775 : 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 47775 : CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
562 47775 : s1 = r2ab + r2bc - r2ca
563 47775 : s2 = r2ab + r2ca - r2bc
564 47775 : s3 = r2ca + r2bc - r2ab
565 47775 : ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
566 :
567 47775 : e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
568 47775 : evdw = evdw - e9
569 47775 : e9tot = e9tot - e9
570 47775 : 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 47775 : 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 47775 : IF (calculate_forces) THEN
582 230390 : rab = rb - ra; rbc = rc - rb; rca = ra - rc
583 23039 : de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
584 23039 : de91 = -de91/3._dp*rabc + 3._dp*e9
585 23039 : dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
586 92156 : fdij(:) = de91*rab(:)/r2ab
587 92156 : fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
588 92156 : fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
589 92156 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
590 92156 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
591 23039 : IF (use_virial) THEN
592 15262 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
593 : END IF
594 92156 : fdij(:) = de91*rbc(:)/r2bc
595 92156 : fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
596 92156 : fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
597 92156 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
598 92156 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
599 23039 : IF (use_virial) THEN
600 15262 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
601 : END IF
602 92156 : fdij(:) = de91*rca(:)/r2ca
603 92156 : fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
604 92156 : fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
605 92156 : force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
606 92156 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
607 23039 : IF (use_virial) THEN
608 15262 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
609 : END IF
610 :
611 23039 : IF (.NOT. dispersion_env%c9cnst) THEN
612 : ! forces from the r-dependence of the coordination numbers
613 : ! atomic stress not implemented
614 :
615 11087 : de91 = 0.5_dp*e9/cc6ab
616 11087 : de921 = de91*dcc6aba
617 11087 : de922 = de91*dcc6abb
618 38781 : 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 38781 : IF (use_virial) THEN
630 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
631 : END IF
632 : END DO
633 38781 : 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 38781 : 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 11087 : de91 = 0.5_dp*e9/cc6bc
650 11087 : de921 = de91*dcc6bcb
651 11087 : de922 = de91*dcc6bcc
652 38781 : 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 38781 : IF (use_virial) THEN
664 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
665 : END IF
666 : END DO
667 38781 : 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 38781 : 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 11087 : de91 = 0.5_dp*e9/cc6ca
684 11087 : de921 = de91*dcc6cac
685 11087 : de922 = de91*dcc6caa
686 38781 : 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 38781 : IF (use_virial) THEN
698 27612 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
699 : END IF
700 : END DO
701 38781 : 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 38781 : 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 55562 : virial%pv_virial = virial%pv_virial + pv_virial_thread
732 :
733 4274 : CALL neighbor_list_iterator_release(nl_iterator)
734 :
735 4274 : elrc6 = 0._dp
736 4274 : elrc8 = 0._dp
737 4274 : elrc9 = 0._dp
738 : ! Long range correction (atomic contributions not implemented)
739 4274 : IF (dispersion_env%lrc) THEN
740 144 : ALLOCATE (cnkind(nkind))
741 106 : cnkind = 0._dp
742 : ! first use the default values
743 106 : DO ikind = 1, nkind
744 58 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
745 106 : cnkind(ikind) = dispersion_env%cn(za)
746 : END DO
747 : ! now check for changes from default
748 48 : 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 106 : DO ikind = 1, nkind
755 58 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
756 58 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
757 58 : IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) CYCLE
758 236 : DO jkind = 1, nkind
759 74 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
760 74 : CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
761 74 : IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
762 72 : 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 66 : cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
769 66 : elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
770 66 : c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
771 66 : elrc8 = elrc8 - s8*twopi*REAL(na*nb, KIND=dp)*c8/(5._dp*rcut**5*cell%deth)
772 198 : IF (dispersion_env%doabc) THEN
773 160 : DO kkind = 1, nkind
774 94 : CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
775 94 : CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
776 94 : IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
777 92 : 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 90 : 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 90 : 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 90 : cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
788 90 : c9 = -SQRT(cc6ab*cc6bc*cc6ca)
789 250 : 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 48 : IF (use_virial) THEN
795 34 : IF (para_env%is_source()) THEN
796 68 : DO i = 1, 3
797 68 : 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 48 : DEALLOCATE (cnkind)
802 : END IF
803 :
804 4274 : DEALLOCATE (cnumbers)
805 4274 : IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
806 62 : DEALLOCATE (cnumfix)
807 : END IF
808 4274 : IF (calculate_forces .OR. debugall) THEN
809 9712 : DO iatom = 1, natom
810 9712 : DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
811 : END DO
812 666 : DEALLOCATE (dcnum)
813 : ELSE
814 3608 : DEALLOCATE (dcnum)
815 : END IF
816 :
817 : ! set dispersion energy
818 4274 : IF (para_env%is_source()) THEN
819 2239 : evdw = evdw + (elrc6 + elrc8 + elrc9)
820 : END IF
821 :
822 4274 : CALL para_env%sum(e6tot)
823 4274 : CALL para_env%sum(e8tot)
824 4274 : CALL para_env%sum(e9tot)
825 4274 : CALL para_env%sum(nab)
826 4274 : CALL para_env%sum(nabc)
827 4274 : e6tot = e6tot + elrc6
828 4274 : e8tot = e8tot + elrc8
829 4274 : e9tot = e9tot + elrc9
830 4274 : IF (unit_nr > 0) THEN
831 6 : WRITE (unit_nr, "(A,F20.0)") " E6 vdW terms :", nab
832 6 : WRITE (unit_nr, *) " E6 vdW energy [au/kcal] :", e6tot, e6tot*kcalmol
833 6 : WRITE (unit_nr, *) " E8 vdW energy [au/kcal] :", e8tot, e8tot*kcalmol
834 6 : WRITE (unit_nr, *) " %E8 on total vdW energy :", e8tot/evdw*100._dp
835 6 : WRITE (unit_nr, "(A,F20.0)") " E9 vdW terms :", nabc
836 6 : WRITE (unit_nr, *) " E9 vdW energy [au/kcal] :", e9tot, e9tot*kcalmol
837 6 : WRITE (unit_nr, *) " %E9 on total vdW energy :", e9tot/evdw*100._dp
838 6 : 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 4274 : DEALLOCATE (dodisp, exclude, atomnumber, rcpbc, radd2, c6d2)
846 :
847 4274 : IF (domol) THEN
848 2 : DEALLOCATE (atom2mol)
849 : END IF
850 :
851 4274 : CALL timestop(handle)
852 :
853 8548 : END SUBROUTINE calculate_dispersion_d3_pairpot
854 :
855 : ! **************************************************************************************************
856 : !> \brief ...
857 : !> \param c6ab ...
858 : !> \param maxci ...
859 : !> \param filename ...
860 : !> \param para_env ...
861 : ! **************************************************************************************************
862 362 : SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
863 :
864 : REAL(KIND=dp), DIMENSION(:, :, :, :, :) :: c6ab
865 : INTEGER, DIMENSION(:) :: maxci
866 : CHARACTER(LEN=*) :: filename
867 : TYPE(mp_para_env_type), POINTER :: para_env
868 :
869 : INTEGER :: funit, iadr, iat, jadr, jat, kk, nl, &
870 : nlines, nn
871 362 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pars
872 :
873 362 : IF (para_env%is_source()) THEN
874 : ! Read the DFT-D3 C6AB parameters from file "filename"
875 186 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
876 186 : READ (funit, *) nl, nlines
877 : END IF
878 362 : CALL para_env%bcast(nl)
879 362 : CALL para_env%bcast(nlines)
880 1086 : ALLOCATE (pars(nl))
881 362 : IF (para_env%is_source()) THEN
882 186 : READ (funit, *) pars(1:nl)
883 186 : CALL close_file(unit_number=funit)
884 : END IF
885 362 : CALL para_env%bcast(pars)
886 :
887 : ! Store C6AB coefficients in an array
888 242483528 : c6ab = -1._dp
889 34390 : maxci = 0
890 362 : kk = 1
891 11723732 : DO nn = 1, nlines
892 11723370 : iat = NINT(pars(kk + 1))
893 11723370 : jat = NINT(pars(kk + 2))
894 11723370 : CALL limit(iat, jat, iadr, jadr)
895 11723370 : maxci(iat) = MAX(maxci(iat), iadr)
896 11723370 : maxci(jat) = MAX(maxci(jat), jadr)
897 11723370 : c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
898 11723370 : c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
899 11723370 : c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
900 :
901 11723370 : c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
902 11723370 : c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
903 11723370 : c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
904 11723732 : kk = (nn*5) + 1
905 : END DO
906 :
907 362 : DEALLOCATE (pars)
908 :
909 362 : END SUBROUTINE dftd3_c6_param
910 :
911 : ! **************************************************************************************************
912 : !> \brief ...
913 : !> \param iat ...
914 : !> \param jat ...
915 : !> \param iadr ...
916 : !> \param jadr ...
917 : ! **************************************************************************************************
918 11723370 : SUBROUTINE limit(iat, jat, iadr, jadr)
919 : INTEGER :: iat, jat, iadr, jadr
920 :
921 : INTEGER :: i
922 :
923 11723370 : iadr = 1
924 11723370 : jadr = 1
925 11723370 : i = 100
926 30160754 : DO WHILE (iat .GT. 100)
927 18437384 : iat = iat - 100
928 30160754 : iadr = iadr + 1
929 : END DO
930 :
931 17471206 : i = 100
932 17471206 : DO WHILE (jat .GT. 100)
933 5747836 : jat = jat - 100
934 5747836 : jadr = jadr + 1
935 : END DO
936 11723370 : END SUBROUTINE limit
937 :
938 : ! **************************************************************************************************
939 : !> \brief ...
940 : !> \param rab ...
941 : !> \param rcutab ...
942 : !> \param srn ...
943 : !> \param alpn ...
944 : !> \param rcut ...
945 : !> \param fdab ...
946 : !> \param dfdab ...
947 : ! **************************************************************************************************
948 7762269 : SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
949 :
950 : REAL(KIND=dp), INTENT(IN) :: rab, rcutab, srn, alpn, rcut
951 : REAL(KIND=dp), INTENT(OUT) :: fdab, dfdab
952 :
953 : REAL(KIND=dp) :: a, b, c, d, dd, dfab, dfcc, dz, fab, &
954 : fcc, rl, rr, ru, z, zz
955 :
956 7762269 : rl = rcut - 1._dp
957 7762269 : ru = rcut
958 7762269 : IF (rab >= ru) THEN
959 : fcc = 0._dp
960 : dfcc = 0._dp
961 7762269 : ELSEIF (rab <= rl) THEN
962 : fcc = 1._dp
963 : dfcc = 0._dp
964 : ELSE
965 709372 : z = rab*rab - rl*rl
966 709372 : dz = 2._dp*rab
967 709372 : zz = z*z*z
968 709372 : d = ru*ru - rl*rl
969 709372 : dd = d*d*d
970 709372 : a = -10._dp/dd
971 709372 : b = 15._dp/(dd*d)
972 709372 : c = -6._dp/(dd*d*d)
973 709372 : fcc = 1._dp + zz*(a + b*z + c*z*z)
974 709372 : dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
975 : END IF
976 :
977 7762269 : rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
978 7762269 : fab = 1._dp/(1._dp + rr)
979 7762269 : dfab = fab*fab*rr*alpn/rab
980 7762269 : fdab = fab*fcc
981 7762269 : dfdab = dfab*fcc + fab*dfcc
982 :
983 7762269 : END SUBROUTINE damping_d3
984 :
985 : ! **************************************************************************************************
986 : !> \brief ...
987 : !> \param maxc ...
988 : !> \param max_elem ...
989 : !> \param c6ab ...
990 : !> \param mxc ...
991 : !> \param iat ...
992 : !> \param jat ...
993 : !> \param nci ...
994 : !> \param ncj ...
995 : !> \param k3 ...
996 : !> \param c6 ...
997 : !> \param dc6a ...
998 : !> \param dc6b ...
999 : ! **************************************************************************************************
1000 8324262 : SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
1001 :
1002 : INTEGER, INTENT(IN) :: maxc, max_elem
1003 : REAL(KIND=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
1004 : INTEGER, INTENT(IN) :: mxc(max_elem), iat, jat
1005 : REAL(KIND=dp), INTENT(IN) :: nci, ncj, k3
1006 : REAL(KIND=dp), INTENT(OUT) :: c6, dc6a, dc6b
1007 :
1008 : INTEGER :: i, j
1009 : REAL(KIND=dp) :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
1010 : dwa, dwb, dza, dzb, r, rsave, rsum, &
1011 : tmp1
1012 :
1013 : ! the exponential is sensitive to numerics
1014 : ! when nci or ncj is much larger than cn1/cn2
1015 :
1016 8324262 : c6mem = -1.0e+99_dp
1017 8324262 : rsave = 1.0e+99_dp
1018 8324262 : rsum = 0.0_dp
1019 8324262 : csum = 0.0_dp
1020 8324262 : dza = 0.0_dp
1021 8324262 : dzb = 0.0_dp
1022 8324262 : dwa = 0.0_dp
1023 8324262 : dwb = 0.0_dp
1024 8324262 : c6 = 0.0_dp
1025 29992800 : DO i = 1, mxc(iat)
1026 91067074 : DO j = 1, mxc(jat)
1027 61074274 : c6 = c6ab(iat, jat, i, j, 1)
1028 82742812 : IF (c6 > 0.0_dp) THEN
1029 61074274 : cn1 = c6ab(iat, jat, i, j, 2)
1030 61074274 : cn2 = c6ab(iat, jat, i, j, 3)
1031 : ! distance
1032 61074274 : r = (cn1 - nci)**2 + (cn2 - ncj)**2
1033 61074274 : IF (r < rsave) THEN
1034 29409155 : rsave = r
1035 29409155 : c6mem = c6
1036 : END IF
1037 61074274 : tmp1 = EXP(k3*r)
1038 61074274 : dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
1039 61074274 : dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
1040 61074274 : rsum = rsum + tmp1
1041 61074274 : csum = csum + tmp1*c6
1042 61074274 : dza = dza + dtmpa*c6
1043 61074274 : dwa = dwa + dtmpa
1044 61074274 : dzb = dzb + dtmpb*c6
1045 61074274 : dwb = dwb + dtmpb
1046 : END IF
1047 : END DO
1048 : END DO
1049 :
1050 8324262 : IF (c6 == 0.0_dp) c6mem = 0.0_dp
1051 :
1052 8324262 : IF (rsum > 1.0e-66_dp) THEN
1053 8317078 : c6 = csum/rsum
1054 8317078 : dc6a = (dza - c6*dwa)/rsum
1055 8317078 : dc6b = (dzb - c6*dwb)/rsum
1056 : ELSE
1057 7184 : c6 = c6mem
1058 7184 : dc6a = 0._dp
1059 7184 : dc6b = 0._dp
1060 : END IF
1061 :
1062 8324262 : END SUBROUTINE getc6
1063 :
1064 : END MODULE qs_dispersion_d3
|
