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