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 dispersion using pair potentials
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dispersion_pairpot
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 bibliography, ONLY: &
20 : Caldeweyher2017, Caldeweyher2019, Caldeweyher2020, Goerigk2017, Wittmann2024, &
21 : cite_reference, grimme2006, grimme2010, grimme2011
22 : USE cell_types, ONLY: cell_type
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_io_unit,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
27 : cp_print_key_unit_nr
28 : USE cp_parser_methods, ONLY: parser_get_next_line,&
29 : parser_get_object
30 : USE cp_parser_types, ONLY: cp_parser_type,&
31 : parser_create,&
32 : parser_release
33 : USE eeq_input, ONLY: read_eeq_param
34 : USE input_constants, ONLY: vdw_pairpot_dftd2,&
35 : vdw_pairpot_dftd3,&
36 : vdw_pairpot_dftd3bj,&
37 : vdw_pairpot_dftd4,&
38 : xc_vdw_fun_pairpot
39 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
40 : section_vals_type,&
41 : section_vals_val_get
42 : USE kinds, ONLY: default_path_length,&
43 : default_string_length,&
44 : dp
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE physcon, ONLY: bohr,&
47 : kcalmol,&
48 : kjmol
49 : USE qs_dispersion_cnum, ONLY: get_cn_radius,&
50 : setcn,&
51 : seten,&
52 : setr0ab,&
53 : setrcov
54 : USE qs_dispersion_d2, ONLY: calculate_dispersion_d2_pairpot,&
55 : dftd2_param
56 : USE qs_dispersion_d3, ONLY: calculate_dispersion_d3_pairpot,&
57 : dftd3_c6_param
58 : USE qs_dispersion_d4, ONLY: calculate_dispersion_d4_pairpot
59 : USE qs_dispersion_s_dftd3, ONLY: dftd3_param_from_library
60 : USE qs_dispersion_types, ONLY: dftd2_pp,&
61 : dftd3_pp,&
62 : dftd4_pp,&
63 : qs_atom_dispersion_type,&
64 : qs_dispersion_type
65 : USE qs_environment_types, ONLY: get_qs_env,&
66 : qs_environment_type
67 : USE qs_force_types, ONLY: qs_force_type
68 : USE qs_kind_types, ONLY: get_qs_kind,&
69 : qs_kind_type,&
70 : set_qs_kind
71 : USE virial_types, ONLY: virial_type
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : PRIVATE
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_pairpot'
79 :
80 : PUBLIC :: qs_dispersion_pairpot_init, calculate_dispersion_pairpot
81 :
82 : ! **************************************************************************************************
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief ...
88 : !> \param atomic_kind_set ...
89 : !> \param qs_kind_set ...
90 : !> \param dispersion_env ...
91 : !> \param pp_section ...
92 : !> \param para_env ...
93 : ! **************************************************************************************************
94 1274 : SUBROUTINE qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
95 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
97 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
98 : TYPE(section_vals_type), OPTIONAL, POINTER :: pp_section
99 : TYPE(mp_para_env_type), POINTER :: para_env
100 :
101 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_dispersion_pairpot_init'
102 :
103 : CHARACTER(LEN=2) :: symbol
104 : CHARACTER(LEN=default_path_length) :: filename
105 : CHARACTER(LEN=default_string_length) :: aname, error_msg
106 : CHARACTER(LEN=default_string_length), &
107 1274 : DIMENSION(:), POINTER :: tmpstringlist
108 : INTEGER :: elem, handle, i, ikind, j, max_elem, &
109 : maxc, n_rep, nkind, nl, vdw_pp_type, &
110 : vdw_type
111 1274 : INTEGER, DIMENSION(:), POINTER :: exlist
112 : LOGICAL :: at_end, explicit, found, is_available
113 : REAL(KIND=dp) :: dum
114 : TYPE(qs_atom_dispersion_type), POINTER :: disp
115 : TYPE(section_vals_type), POINTER :: eeq_section
116 :
117 1274 : CALL timeset(routineN, handle)
118 :
119 1274 : nkind = SIZE(atomic_kind_set)
120 :
121 1274 : vdw_type = dispersion_env%type
122 1274 : SELECT CASE (vdw_type)
123 : CASE DEFAULT
124 : ! do nothing
125 : CASE (xc_vdw_fun_pairpot)
126 : ! setup information on pair potentials
127 1274 : vdw_pp_type = dispersion_env%type
128 1274 : SELECT CASE (dispersion_env%pp_type)
129 : CASE DEFAULT
130 : ! do nothing
131 : CASE (vdw_pairpot_dftd2)
132 36 : CALL cite_reference(Grimme2006)
133 102 : DO ikind = 1, nkind
134 66 : CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
135 66 : ALLOCATE (disp)
136 66 : disp%type = dftd2_pp
137 : ! get filename of parameter file
138 66 : filename = dispersion_env%parameter_file_name
139 : ! check for local parameters
140 66 : found = .FALSE.
141 66 : IF (PRESENT(pp_section)) THEN
142 60 : CALL section_vals_val_get(pp_section, "ATOMPARM", n_rep_val=n_rep)
143 60 : DO i = 1, n_rep
144 : CALL section_vals_val_get(pp_section, "ATOMPARM", i_rep_val=i, &
145 0 : c_vals=tmpstringlist)
146 60 : IF (TRIM(tmpstringlist(1)) == TRIM(symbol)) THEN
147 : ! we assume the parameters are in atomic units!
148 0 : READ (tmpstringlist(2), *) disp%c6
149 0 : READ (tmpstringlist(3), *) disp%vdw_radii
150 0 : found = .TRUE.
151 0 : EXIT
152 : END IF
153 : END DO
154 : END IF
155 66 : IF (.NOT. found) THEN
156 : ! check for internal parameters
157 66 : CALL dftd2_param(elem, disp%c6, disp%vdw_radii, found)
158 : END IF
159 66 : IF (.NOT. found) THEN
160 : ! check on file
161 0 : INQUIRE (FILE=filename, EXIST=is_available)
162 0 : IF (is_available) THEN
163 0 : BLOCK
164 : TYPE(cp_parser_type) :: parser
165 0 : CALL parser_create(parser, filename, para_env=para_env)
166 : DO
167 : at_end = .FALSE.
168 0 : CALL parser_get_next_line(parser, 1, at_end)
169 0 : IF (at_end) EXIT
170 0 : CALL parser_get_object(parser, aname)
171 0 : IF (TRIM(aname) == TRIM(symbol)) THEN
172 0 : CALL parser_get_object(parser, disp%c6)
173 : ! we have to change the units J*nm^6*mol^-1 -> Hartree*Bohr^6
174 0 : disp%c6 = disp%c6*1000._dp*bohr**6/kjmol
175 0 : CALL parser_get_object(parser, disp%vdw_radii)
176 0 : disp%vdw_radii = disp%vdw_radii*bohr
177 0 : found = .TRUE.
178 0 : EXIT
179 : END IF
180 : END DO
181 0 : CALL parser_release(parser)
182 : END BLOCK
183 : END IF
184 : END IF
185 66 : IF (found) THEN
186 66 : disp%defined = .TRUE.
187 : ELSE
188 0 : disp%defined = .FALSE.
189 : END IF
190 : ! Check if the parameter is defined
191 66 : IF (.NOT. disp%defined) THEN
192 : CALL cp_abort(__LOCATION__, &
193 : "Dispersion parameters for element ("//TRIM(symbol)//") are not defined! "// &
194 : "Please provide a valid set of parameters through the input section or "// &
195 0 : "through an external file! ")
196 : END IF
197 168 : CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
198 : END DO
199 : CASE (vdw_pairpot_dftd3, vdw_pairpot_dftd3bj)
200 : !DFT-D3 Method initial setup
201 556 : CALL cite_reference(Grimme2010)
202 556 : CALL cite_reference(Grimme2011)
203 556 : CALL cite_reference(Goerigk2017)
204 556 : CALL cite_reference(Wittmann2024)
205 556 : max_elem = 103
206 556 : maxc = 7
207 556 : dispersion_env%max_elem = max_elem
208 556 : dispersion_env%maxc = maxc
209 556 : ALLOCATE (dispersion_env%maxci(max_elem))
210 556 : ALLOCATE (dispersion_env%c6ab(max_elem, max_elem, maxc, maxc, 3))
211 556 : ALLOCATE (dispersion_env%r0ab(max_elem, max_elem))
212 556 : ALLOCATE (dispersion_env%rcov(max_elem))
213 556 : ALLOCATE (dispersion_env%eneg(max_elem))
214 556 : ALLOCATE (dispersion_env%r2r4(max_elem))
215 556 : ALLOCATE (dispersion_env%cn(max_elem))
216 :
217 556 : IF (dispersion_env%d3_reference_code) THEN
218 : CALL dftd3_param_from_library(dispersion_env%c6ab, dispersion_env%maxci, &
219 : dispersion_env%r0ab, dispersion_env%rcov, &
220 : dispersion_env%r2r4, &
221 : dispersion_env%pp_type, dispersion_env%ref_functional, &
222 : dispersion_env%s6, dispersion_env%s8, &
223 : dispersion_env%a1, dispersion_env%a2, &
224 : dispersion_env%sr6, para_env, error=error_msg, &
225 52 : calc_scaling=.NOT. dispersion_env%d3_scaling_explicit)
226 52 : IF (error_msg /= "") THEN
227 0 : CALL cp_abort(__LOCATION__, error_msg)
228 : END IF
229 : ELSE
230 504 : filename = dispersion_env%parameter_file_name
231 504 : CALL dftd3_c6_param(dispersion_env%c6ab, dispersion_env%maxci, filename, para_env)
232 504 : CALL setrcov(dispersion_env%rcov)
233 504 : CALL setr0ab(dispersion_env%r0ab, dispersion_env%rcov, dispersion_env%r2r4)
234 : END IF
235 : ! Electronegativity
236 556 : CALL seten(dispersion_env%eneg)
237 : ! the default coordination numbers
238 556 : CALL setcn(dispersion_env%cn)
239 : ! scale r4/r2 values of the atoms by sqrt(Z)
240 : ! sqrt is also globally close to optimum
241 : ! together with the factor 1/2 this yield reasonable
242 : ! c8 for he, ne and ar. for larger Z, C8 becomes too large
243 : ! which effectively mimics higher R^n terms neglected due
244 : ! to stability reasons
245 556 : IF (.NOT. dispersion_env%d3_reference_code) THEN
246 52416 : DO i = 1, max_elem
247 51912 : dum = 0.5_dp*dispersion_env%r2r4(i)*REAL(i, dp)**0.5_dp
248 : ! store it as sqrt because the geom. av. is taken
249 52416 : dispersion_env%r2r4(i) = SQRT(dum)
250 : END DO
251 : END IF
252 : ! parameters
253 556 : dispersion_env%k1 = 16.0_dp
254 556 : dispersion_env%k2 = 4._dp/3._dp
255 : ! reasonable choices are between 3 and 5
256 : ! this gives smoth curves with maxima around the integer values
257 : ! k3=3 give for CN=0 a slightly smaller value than computed
258 : ! for the free atom. This also yields to larger CN for atoms
259 : ! in larger molecules but with the same chem. environment
260 : ! which is physically not right
261 : ! values >5 might lead to bumps in the potential
262 556 : dispersion_env%k3 = -4._dp
263 556 : IF (.NOT. dispersion_env%d3_reference_code) THEN
264 52416 : dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
265 : END IF
266 : ! alpha default parameter
267 556 : dispersion_env%alp = 14._dp
268 : !
269 1778 : DO ikind = 1, nkind
270 1222 : CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
271 1222 : ALLOCATE (disp)
272 1222 : disp%type = dftd3_pp
273 1222 : IF (elem <= max_elem) THEN
274 1222 : disp%defined = .TRUE.
275 : ELSE
276 : disp%defined = .FALSE.
277 : END IF
278 1222 : IF (.NOT. disp%defined) THEN
279 : CALL cp_abort(__LOCATION__, &
280 : "Dispersion parameters for element ("//TRIM(symbol)//") are not defined! "// &
281 : "Please provide a valid set of parameters through the input section or "// &
282 0 : "through an external file! ")
283 : END IF
284 3000 : CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
285 : END DO
286 :
287 556 : IF (PRESENT(pp_section)) THEN
288 : ! Check for coordination numbers
289 154 : CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", n_rep_val=n_rep)
290 154 : IF (n_rep > 0) THEN
291 24 : ALLOCATE (dispersion_env%cnkind(n_rep))
292 12 : DO i = 1, n_rep
293 : CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", i_rep_val=i, &
294 6 : c_vals=tmpstringlist)
295 6 : READ (tmpstringlist(1), *) dispersion_env%cnkind(i)%cnum
296 12 : READ (tmpstringlist(2), *) dispersion_env%cnkind(i)%kind
297 : END DO
298 : END IF
299 154 : CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", n_rep_val=n_rep)
300 154 : IF (n_rep > 0) THEN
301 10 : ALLOCATE (dispersion_env%cnlist(n_rep))
302 6 : DO i = 1, n_rep
303 : CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", i_rep_val=i, &
304 4 : c_vals=tmpstringlist)
305 4 : nl = SIZE(tmpstringlist)
306 12 : ALLOCATE (dispersion_env%cnlist(i)%atom(nl - 1))
307 4 : dispersion_env%cnlist(i)%natom = nl - 1
308 4 : READ (tmpstringlist(1), *) dispersion_env%cnlist(i)%cnum
309 14 : DO j = 1, nl - 1
310 12 : READ (tmpstringlist(j + 1), *) dispersion_env%cnlist(i)%atom(j)
311 : END DO
312 : END DO
313 : END IF
314 : ! Check for exclusion lists
315 154 : CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", explicit=explicit)
316 154 : IF (explicit) THEN
317 2 : CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", i_vals=exlist)
318 4 : DO j = 1, SIZE(exlist)
319 2 : ikind = exlist(j)
320 2 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp)
321 4 : disp%defined = .FALSE.
322 : END DO
323 : END IF
324 154 : CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", n_rep_val=n_rep)
325 154 : dispersion_env%nd3_exclude_pair = n_rep
326 154 : IF (n_rep > 0) THEN
327 6 : ALLOCATE (dispersion_env%d3_exclude_pair(n_rep, 2))
328 6 : DO i = 1, n_rep
329 : CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", i_rep_val=i, &
330 4 : i_vals=exlist)
331 22 : dispersion_env%d3_exclude_pair(i, :) = exlist
332 : END DO
333 : END IF
334 : END IF
335 : CASE (vdw_pairpot_dftd4)
336 : !most checks are done by the library
337 682 : CALL cite_reference(Caldeweyher2017)
338 682 : CALL cite_reference(Caldeweyher2019)
339 682 : CALL cite_reference(Caldeweyher2020)
340 2122 : DO ikind = 1, nkind
341 1440 : CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
342 1440 : ALLOCATE (disp)
343 1440 : disp%type = dftd4_pp
344 1440 : disp%defined = .TRUE.
345 3562 : CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
346 : END DO
347 : ! maybe needed in cnumber calculations
348 682 : max_elem = 103
349 682 : maxc = 7
350 682 : dispersion_env%max_elem = max_elem
351 682 : dispersion_env%maxc = maxc
352 682 : ALLOCATE (dispersion_env%maxci(max_elem))
353 682 : ALLOCATE (dispersion_env%rcov(max_elem))
354 682 : ALLOCATE (dispersion_env%eneg(max_elem))
355 682 : ALLOCATE (dispersion_env%cn(max_elem))
356 : ! the default covalent radii
357 682 : CALL setrcov(dispersion_env%rcov)
358 : ! the default coordination numbers
359 682 : CALL setcn(dispersion_env%cn)
360 : ! Electronegativity
361 682 : CALL seten(dispersion_env%eneg)
362 : ! parameters
363 682 : dispersion_env%k1 = 16.0_dp
364 682 : dispersion_env%k2 = 4._dp/3._dp
365 682 : dispersion_env%k3 = -4._dp
366 70928 : dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
367 682 : dispersion_env%alp = 14._dp
368 : !
369 682 : dispersion_env%cnfun = 3
370 682 : IF (dispersion_env%rc_cn < 0.0_dp) THEN
371 682 : dispersion_env%rc_cn = get_cn_radius(dispersion_env)
372 : END IF
373 1956 : IF (PRESENT(pp_section)) THEN
374 42 : eeq_section => section_vals_get_subs_vals(pp_section, "EEQ")
375 42 : CALL read_eeq_param(eeq_section, dispersion_env%eeq_sparam)
376 : END IF
377 : END SELECT
378 : END SELECT
379 :
380 1274 : CALL timestop(handle)
381 :
382 1274 : END SUBROUTINE qs_dispersion_pairpot_init
383 :
384 : ! **************************************************************************************************
385 : !> \brief ...
386 : !> \param qs_env ...
387 : !> \param dispersion_env ...
388 : !> \param energy ...
389 : !> \param calculate_forces ...
390 : !> \param atevdw ...
391 : ! **************************************************************************************************
392 24642 : SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
393 :
394 : TYPE(qs_environment_type), POINTER :: qs_env
395 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
396 : REAL(KIND=dp), INTENT(INOUT) :: energy
397 : LOGICAL, INTENT(IN) :: calculate_forces
398 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atevdw
399 :
400 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_pairpot'
401 :
402 : INTEGER :: atom_a, handle, iatom, ikind, iw, natom, &
403 : nkind, unit_nr
404 24642 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
405 : LOGICAL :: atenergy, atex, debugall, use_virial
406 : REAL(KIND=dp) :: evdw, gnorm
407 24642 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atomic_energy
408 : REAL(KIND=dp), DIMENSION(3) :: fdij
409 : REAL(KIND=dp), DIMENSION(3, 3) :: dvirial, pv_loc, pv_virial_thread
410 24642 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
411 : TYPE(atprop_type), POINTER :: atprop
412 : TYPE(cell_type), POINTER :: cell
413 : TYPE(cp_logger_type), POINTER :: logger
414 : TYPE(mp_para_env_type), POINTER :: para_env
415 24642 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
416 : TYPE(virial_type), POINTER :: virial
417 :
418 24642 : energy = 0._dp
419 : ! make valgrind happy
420 24642 : use_virial = .FALSE.
421 :
422 24642 : IF (dispersion_env%type /= xc_vdw_fun_pairpot) THEN
423 : RETURN
424 : END IF
425 :
426 6746 : CALL timeset(routineN, handle)
427 :
428 6746 : NULLIFY (atomic_kind_set)
429 :
430 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
431 6746 : cell=cell, virial=virial, para_env=para_env, atprop=atprop)
432 :
433 6746 : debugall = dispersion_env%verbose
434 :
435 6746 : NULLIFY (logger)
436 6746 : logger => cp_get_default_logger()
437 6746 : IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
438 : unit_nr = cp_print_key_unit_nr(logger, dispersion_env%dftd_section, "PRINT_DFTD", &
439 378 : extension=".dftd")
440 : ELSE
441 6368 : unit_nr = -1
442 : END IF
443 :
444 : ! atomic energy and stress arrays
445 6746 : atenergy = atprop%energy
446 : ! external atomic energy
447 6746 : atex = .FALSE.
448 6746 : IF (PRESENT(atevdw)) THEN
449 2 : atex = .TRUE.
450 : END IF
451 :
452 6746 : IF (unit_nr > 0) THEN
453 68 : WRITE (unit_nr, *)
454 68 : WRITE (unit_nr, *) " Pair potential vdW calculation"
455 68 : IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
456 7 : WRITE (unit_nr, *) " Dispersion potential type: DFT-D2"
457 7 : WRITE (unit_nr, *) " Scaling parameter (s6) ", dispersion_env%scaling
458 7 : WRITE (unit_nr, *) " Exponential prefactor ", dispersion_env%exp_pre
459 61 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
460 20 : WRITE (unit_nr, *) " Dispersion potential type: DFT-D3"
461 41 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
462 25 : WRITE (unit_nr, *) " Dispersion potential type: DFT-D3(BJ)"
463 16 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
464 16 : WRITE (unit_nr, *) " Dispersion potential type: DFT-D4"
465 : END IF
466 : END IF
467 :
468 6746 : CALL get_qs_env(qs_env=qs_env, force=force)
469 6746 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
470 714 : IF (use_virial .AND. debugall) THEN
471 104 : dvirial = virial%pv_virial
472 : END IF
473 6746 : IF (use_virial) THEN
474 9282 : pv_loc = virial%pv_virial
475 : END IF
476 :
477 6746 : evdw = 0._dp
478 : pv_virial_thread(:, :) = 0._dp
479 :
480 6746 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
481 :
482 6746 : IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
483 144 : CALL calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
484 6674 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
485 : dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
486 : CALL calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
487 11526 : unit_nr, atevdw)
488 910 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
489 910 : IF (dispersion_env%lrc) THEN
490 0 : CPABORT("Long range correction with DFTD4 not implemented")
491 : END IF
492 910 : IF (dispersion_env%srb) THEN
493 0 : CPABORT("Short range bond correction with DFTD4 not implemented")
494 : END IF
495 910 : IF (dispersion_env%domol) THEN
496 0 : CPABORT("Molecular approximation with DFTD4 not implemented")
497 : END IF
498 : !
499 910 : iw = -1
500 910 : IF (dispersion_env%verbose) iw = cp_logger_get_default_io_unit(logger)
501 : !
502 910 : IF (atenergy .OR. atex) THEN
503 0 : ALLOCATE (atomic_energy(natom))
504 : CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
505 0 : iw, atomic_energy=atomic_energy)
506 : ELSE
507 910 : CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw)
508 : END IF
509 : !
510 910 : IF (atex) THEN
511 0 : atevdw(1:natom) = atomic_energy(1:natom)
512 : END IF
513 910 : IF (atenergy) THEN
514 0 : CALL atprop_array_init(atprop%atevdw, natom)
515 0 : atprop%atevdw(1:natom) = atomic_energy(1:natom)
516 : END IF
517 910 : IF (atenergy .OR. atex) THEN
518 0 : DEALLOCATE (atomic_energy)
519 : END IF
520 : END IF
521 :
522 : ! set dispersion energy
523 6746 : CALL para_env%sum(evdw)
524 6746 : energy = evdw
525 6746 : IF (unit_nr > 0) THEN
526 68 : WRITE (unit_nr, *) " Total vdW energy [au] :", evdw
527 68 : WRITE (unit_nr, *) " Total vdW energy [kcal] :", evdw*kcalmol
528 68 : WRITE (unit_nr, *)
529 : END IF
530 6746 : IF (calculate_forces .AND. debugall) THEN
531 42 : IF (unit_nr > 0) THEN
532 1 : WRITE (unit_nr, *) " Dispersion Forces "
533 1 : WRITE (unit_nr, *) " Atom Kind Forces "
534 : END IF
535 42 : gnorm = 0._dp
536 226 : DO iatom = 1, natom
537 184 : ikind = kind_of(iatom)
538 184 : atom_a = atom_of_kind(iatom)
539 736 : fdij(1:3) = force(ikind)%dispersion(:, atom_a)
540 184 : CALL para_env%sum(fdij)
541 736 : gnorm = gnorm + SUM(ABS(fdij))
542 226 : IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
543 : END DO
544 42 : IF (unit_nr > 0) THEN
545 1 : WRITE (unit_nr, *)
546 1 : WRITE (unit_nr, *) "|G| = ", gnorm
547 1 : WRITE (unit_nr, *)
548 : END IF
549 42 : IF (use_virial) THEN
550 78 : dvirial = virial%pv_virial - dvirial
551 6 : CALL para_env%sum(dvirial)
552 6 : IF (unit_nr > 0) THEN
553 0 : WRITE (unit_nr, *) "Stress Tensor (dispersion)"
554 0 : WRITE (unit_nr, "(3G20.12)") dvirial
555 0 : WRITE (unit_nr, *) " Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
556 0 : WRITE (unit_nr, *)
557 : END IF
558 : END IF
559 : END IF
560 :
561 6746 : IF (calculate_forces .AND. use_virial) THEN
562 3588 : virial%pv_vdw = virial%pv_vdw + (virial%pv_virial - pv_loc)
563 : END IF
564 :
565 6746 : IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
566 378 : CALL cp_print_key_finished_output(unit_nr, logger, dispersion_env%dftd_section, "PRINT_DFTD")
567 : END IF
568 :
569 6746 : CALL timestop(handle)
570 :
571 31388 : END SUBROUTINE calculate_dispersion_pairpot
572 :
573 : END MODULE qs_dispersion_pairpot
|
