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