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