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