(git:f2099e5)
Loading...
Searching...
No Matches
qs_dispersion_d4.F
Go to the documentation of this file.
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 Johann Pototschnig
11! **************************************************************************************************
12MODULE qs_dispersion_d4
13 USE atomic_kind_types, ONLY: atomic_kind_type, &
14 get_atomic_kind, &
15 get_atomic_kind_set
16 USE distribution_1d_types, ONLY: distribution_1d_type
17 USE eeq_method, ONLY: eeq_charges, eeq_forces
18 USE machine, ONLY: m_flush, &
19 m_walltime
20 USE cell_types, ONLY: cell_type, &
21 plane_distance, &
22 pbc, &
23 get_cell
24 USE qs_environment_types, ONLY: get_qs_env, &
25 qs_environment_type
26 USE qs_force_types, ONLY: qs_force_type
27 USE qs_kind_types, ONLY: get_qs_kind, &
28 qs_kind_type, &
29 set_qs_kind
30 USE qs_neighbor_list_types, ONLY: get_iterator_info, &
31 neighbor_list_iterate, &
32 neighbor_list_iterator_create, &
33 neighbor_list_iterator_p_type, &
34 neighbor_list_iterator_release, &
35 neighbor_list_set_p_type
36 USE virial_methods, ONLY: virial_pair_force
37 USE virial_types, ONLY: virial_type
38 USE kinds, ONLY: dp
39 USE particle_types, ONLY: particle_type
40 USE qs_dispersion_types, ONLY: qs_dispersion_type
41 USE qs_dispersion_utils, ONLY: cellhash
42 USE qs_dispersion_cnum, ONLY: cnumber_init, dcnum_type, cnumber_release
43 USE message_passing, ONLY: mp_para_env_type
44
45#if defined(__DFTD4)
46!&<
47 USE dftd4, ONLY: d4_model, &
48 damping_param, &
49 get_dispersion, &
50 get_rational_damping, &
51 new, &
52 new_d4_model, &
53 realspace_cutoff, &
54 structure_type, &
55 rational_damping_param, &
56 get_coordination_number, &
57 get_lattice_points
58 USE multicharge, ONLY: get_charges
59 USE mctc_env, ONLY: error_type
60!&>
61#endif
62#include "./base/base_uses.f90"
63
64 IMPLICIT NONE
65
66 PRIVATE
67
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
69
70 PUBLIC :: calculate_dispersion_d4_pairpot
71
72! **************************************************************************************************
73
74CONTAINS
75
76#if defined(__DFTD4)
77! **************************************************************************************************
78!> \brief ...
79!> \param qs_env ...
80!> \param dispersion_env ...
81!> \param evdw ...
82!> \param calculate_forces ...
83!> \param iw ...
84!> \param atomic_energy ...
85! **************************************************************************************************
86 SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, &
87 atomic_energy)
88 TYPE(qs_environment_type), POINTER :: qs_env
89 TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
90 REAL(KIND=dp), INTENT(INOUT) :: evdw
91 LOGICAL, INTENT(IN) :: calculate_forces
92 INTEGER, INTENT(IN) :: iw
93 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atomic_energy
94
95 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
96
97 INTEGER :: atoma, cnfun, enshift, handle, i, iatom, &
98 ifull, ikind, mref, natom, natom_full, &
99 ncoup, nghost
100 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_map, atom_map_back, atom_of_kind, &
101 atomtype, kind_of, species, t_atomtype
102 INTEGER, DIMENSION(3) :: periodic
103 LOGICAL :: debug, grad, ifloating, ighost, &
104 use_virial
105 LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_ghost, exclude_ghost
106 LOGICAL, DIMENSION(3) :: lperiod
107 REAL(KIND=dp) :: ed2, ed3, ev1, ev2, ev3, ev4, pd2, pd3, &
108 ta, tb, tc, td, te, ts
109 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cn, cn_red, cnd, dedcn, dedq, edcn, edq, &
110 enerd2, enerd3, energies, energies3, &
111 q_red, qd
112 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ga, gradient, t_xyz, tvec, xyz
113 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gdeb, gwdcn, gwdq, gwvec
114 REAL(KIND=dp), DIMENSION(3, 3) :: sigma, stress
115 REAL(KIND=dp), DIMENSION(3, 3, 4) :: sdeb
116 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
117 TYPE(cell_type), POINTER :: cell
118 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
119 TYPE(mp_para_env_type), POINTER :: para_env
120 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
121 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
122 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
123 TYPE(virial_type), POINTER :: virial
124
125 CLASS(damping_param), ALLOCATABLE :: param
126 TYPE(d4_model) :: disp
127 TYPE(structure_type) :: mol
128 TYPE(realspace_cutoff) :: cutoff
129
130 TYPE(error_type), ALLOCATABLE :: error
131
132 CALL timeset(routinen, handle)
133
134 debug = dispersion_env%d4_debug
135
136 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
137 cell=cell, force=force, virial=virial, para_env=para_env)
138 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
139
140 !get information about particles
141 natom_full = SIZE(particle_set)
142 nghost = 0
143 ALLOCATE (t_xyz(3, natom_full), t_atomtype(natom_full), a_ghost(natom_full))
144 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
145 DO iatom = 1, natom_full
146 t_xyz(:, iatom) = particle_set(iatom)%r(:)
147 ikind = kind_of(iatom)
148 CALL get_qs_kind(qs_kind_set(ikind), zatom=t_atomtype(iatom), ghost=ighost, floating=ifloating)
149 a_ghost(iatom) = ighost .OR. ifloating
150 IF (a_ghost(iatom)) nghost = nghost + 1
151 END DO
152
153 natom = natom_full - nghost
154 ! Build atom mapping: full index -> reduced index (0 for ghost)
155 ALLOCATE (atom_map(natom_full), atom_map_back(natom))
156 atom_map = 0
157 iatom = 0
158 ALLOCATE (xyz(3, natom), atomtype(natom))
159 DO i = 1, natom_full
160 IF (.NOT. a_ghost(i)) THEN
161 iatom = iatom + 1
162 atom_map(i) = iatom
163 atom_map_back(iatom) = i
164 xyz(:, iatom) = t_xyz(:, i)
165 atomtype(iatom) = t_atomtype(i)
166 END IF
167 END DO
168 DEALLOCATE (a_ghost, t_xyz, t_atomtype)
169
170 !get information about cell / lattice
171 CALL get_cell(cell=cell, periodic=periodic)
172 lperiod(1) = periodic(1) == 1
173 lperiod(2) = periodic(2) == 1
174 lperiod(3) = periodic(3) == 1
175 ! enforce en shift method 1 (original/molecular)
176 ! method 2 from paper on PBC seems not to work
177 enshift = 1
178 !IF (ALL(periodic == 0)) enshift = 1
179
180 !prepare for the call to the dispersion function
181 CALL new(mol, atomtype, xyz, lattice=cell%hmat, periodic=lperiod)
182 CALL new_d4_model(error, disp, mol)
183 IF (ALLOCATED(error)) THEN
184 cpabort(error%message)
185 END IF
186
187 ! Number of coupling
188 ncoup = disp%ncoup
189
190 ! Build species mapping: full atom index -> D4 species ID (0 for ghost)
191 ALLOCATE (species(natom_full))
192 species = 0
193 DO i = 1, natom
194 species(atom_map_back(i)) = mol%id(i)
195 END DO
196
197 ! Build per-kind exclusion mask for EEQ (ghost/floating kinds excluded)
198 ALLOCATE (exclude_ghost(SIZE(qs_kind_set)))
199 exclude_ghost = .false.
200 DO i = 1, SIZE(qs_kind_set)
201 CALL get_qs_kind(qs_kind_set(i), ghost=ighost, floating=ifloating)
202 exclude_ghost(i) = ighost .OR. ifloating
203 END DO
204
205 IF (dispersion_env%ref_functional == "none") THEN
206 CALL get_rational_damping("pbe", param, s9=0.0_dp)
207 IF (.NOT. ALLOCATED(param)) THEN
208 cpabort("D4: Failed to get rational damping parameters for default functional")
209 END IF
210 SELECT TYPE (param)
211 TYPE is (rational_damping_param)
212 param%s6 = dispersion_env%s6
213 param%s8 = dispersion_env%s8
214 param%a1 = dispersion_env%a1
215 param%a2 = dispersion_env%a2
216 param%alp = dispersion_env%alp
217 END SELECT
218 ELSE
219 CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
220 IF (.NOT. ALLOCATED(param)) THEN
221 cpabort("D4: Unknown reference functional '"//trim(dispersion_env%ref_functional)//"'")
222 END IF
223 SELECT TYPE (param)
224 TYPE is (rational_damping_param)
225 dispersion_env%s6 = param%s6
226 dispersion_env%s8 = param%s8
227 dispersion_env%a1 = param%a1
228 dispersion_env%a2 = param%a2
229 dispersion_env%alp = param%alp
230 END SELECT
231 END IF
232
233 ! Coordination number cutoff
234 cutoff%cn = dispersion_env%rc_cn
235 ! Two-body interaction cutoff
236 cutoff%disp2 = dispersion_env%rc_d4*2._dp
237 ! Three-body interaction cutoff
238 cutoff%disp3 = dispersion_env%rc_disp*2._dp
239 IF (cutoff%disp3 > cutoff%disp2) THEN
240 cpabort("D4: Three-body cutoff should be smaller than two-body cutoff")
241 END IF
242 cutoff%width2 = dispersion_env%d4_cutoff_width
243 cutoff%width3 = dispersion_env%d4_3b_cutoff_width
244 IF (cutoff%width2 < 0.0_dp .OR. cutoff%width2 >= cutoff%disp2) THEN
245 cpabort("D4: Two-body cutoff width must be non-negative and smaller than the cutoff")
246 END IF
247 IF (cutoff%width3 < 0.0_dp .OR. cutoff%width3 >= cutoff%disp3) THEN
248 cpabort("D4: Three-body cutoff width must be non-negative and smaller than the cutoff")
249 END IF
250
251 IF (calculate_forces) THEN
252 grad = .true.
253 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
254 ELSE
255 grad = .false.
256 use_virial = .false.
257 END IF
258
259 IF (dispersion_env%d4_reference_code) THEN
260
261 !> Wrapper to handle the evaluation of dispersion energy and derivatives
262 IF (.NOT. dispersion_env%doabc) THEN
263 cpwarn("Using D4_REFERENCE_CODE enforces calculation of C9 term.")
264 END IF
265 IF (grad) THEN
266 ALLOCATE (gradient(3, natom))
267 CALL get_dispersion(mol, disp, param, cutoff, evdw, gradient, stress)
268 IF (calculate_forces) THEN
269 IF (use_virial) THEN
270 virial%pv_virial = virial%pv_virial - stress/para_env%num_pe
271 END IF
272 DO iatom = 1, natom
273 ifull = atom_map_back(iatom)
274 ikind = kind_of(ifull)
275 atoma = atom_of_kind(ifull)
276 force(ikind)%dispersion(:, atoma) = &
277 force(ikind)%dispersion(:, atoma) + gradient(:, iatom)/para_env%num_pe
278 END DO
279 END IF
280 DEALLOCATE (gradient)
281 ELSE
282 CALL get_dispersion(mol, disp, param, cutoff, evdw)
283 END IF
284 !dispersion energy is computed by every MPI process
285 evdw = evdw/para_env%num_pe
286 IF (dispersion_env%ext_charges) dispersion_env%dcharges = 0.0_dp
287 IF (PRESENT(atomic_energy)) THEN
288 cpwarn("Atomic energies not available for D4 reference code")
289 atomic_energy = 0.0_dp
290 END IF
291
292 ELSE
293
294 IF (iw > 0) THEN
295 WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
296 WRITE (iw, fmt="(T32,A)") "DEBUG D4 DISPERSION"
297 WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
298 WRITE (iw, '(A,T71,A10)') " DEBUG D4| Reference functional ", trim(dispersion_env%ref_functional)
299 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s6) ", dispersion_env%s6
300 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s8) ", dispersion_env%s8
301 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a1) ", dispersion_env%a1
302 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a2) ", dispersion_env%a2
303 WRITE (iw, '(A,T71,E10.4)') " DEBUG D4| Cutoff value coordination numbers ", dispersion_env%eps_cn
304 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius coordination numbers ", dispersion_env%rc_cn
305 WRITE (iw, '(A,T71,I10)') " DEBUG D4| Coordination number function type ", dispersion_env%cnfun
306 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 2-body terms [bohr]", 2._dp*dispersion_env%rc_d4
307 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 3-body terms [bohr]", 2._dp*dispersion_env%rc_disp
308 END IF
309
310 td = 0.0_dp
311 IF (debug .AND. iw > 0) THEN
312 ts = m_walltime()
313 CALL refd4_debug(param, disp, mol, cutoff, grad, dispersion_env%doabc, &
314 enerd2, enerd3, cnd, qd, edcn, edq, gdeb, sdeb)
315 te = m_walltime()
316 td = te - ts
317 END IF
318
319 tc = 0.0_dp
320 ts = m_walltime()
321
322 mref = maxval(disp%ref)
323 ! Coordination numbers (full-size from qs_env; ghosts excluded internally)
324 cnfun = dispersion_env%cnfun
325 CALL cnumber_init(qs_env, cn, dcnum, cnfun, grad)
326 ! cn has size natom_full; ghost entries are 0
327
328 ! Filter CN to reduced space for D4 model
329 ALLOCATE (cn_red(natom))
330 DO i = 1, natom
331 cn_red(i) = cn(atom_map_back(i))
332 END DO
333 IF (debug .AND. iw > 0) THEN
334 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (max)", maxval(abs(cn_red - cnd))
335 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (ave)", sum(abs(cn_red - cnd))/natom
336 END IF
337
338 ! EEQ charges
339 ! Use CP2K's MPI-parallel EEQ solver with ghost exclusion.
340 ! Ghost kinds get huge hardness + zero coupling -> q_ghost = 0 exactly,
341 ! without influencing real atom charges. Preserves full MPI parallelism.
342 IF (dispersion_env%ext_charges) THEN
343 ALLOCATE (q_red(natom))
344 q_red(1:natom) = dispersion_env%charges(1:natom)
345 ELSE
346 ! eeq_charges writes natom_full entries (from qs_env), so allocate full-size
347 ALLOCATE (q_red(natom_full))
348 CALL eeq_charges(qs_env, q_red, dispersion_env%eeq_sparam, 2, enshift, &
349 exclude=exclude_ghost, cn_max=8.0_dp)
350 ! Filter to reduced size for D4 model
351 block
352 REAL(KIND=dp), ALLOCATABLE :: q_tmp(:)
353 ALLOCATE (q_tmp(natom))
354 DO i = 1, natom
355 q_tmp(i) = q_red(atom_map_back(i))
356 END DO
357 DEALLOCATE (q_red)
358 CALL move_alloc(q_tmp, q_red)
359 END block
360 END IF
361 IF (debug .AND. iw > 0) THEN
362 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (max)", maxval(abs(q_red - qd))
363 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (ave)", sum(abs(q_red - qd))/natom
364 END IF
365 ! Weights for C6 calculation (reduced space)
366 ALLOCATE (gwvec(mref, natom, ncoup))
367 IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
368 CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
369
370 ! Energies and derivatives (full-size for CP2K infrastructure compatibility)
371 ALLOCATE (energies(natom_full))
372 energies(:) = 0.0_dp
373 IF (grad) THEN
374 ALLOCATE (gradient(3, natom_full), ga(3, natom_full))
375 ALLOCATE (dedcn(natom_full), dedq(natom_full))
376 dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
377 ga(:, :) = 0.0_dp
378 sigma(:, :) = 0.0_dp
379 END IF
380 CALL dispersion_2b(dispersion_env, cutoff%disp2, disp%r4r2, &
381 gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
382 energies, dedcn, dedq, grad, ga, sigma, &
383 atom_map, species)
384 IF (grad) THEN
385 gradient(1:3, 1:natom_full) = ga(1:3, 1:natom_full)
386 stress = sigma
387 IF (debug) THEN
388 CALL para_env%sum(ga)
389 CALL para_env%sum(sigma)
390 IF (iw > 0) THEN
391 CALL gerror(ga, gdeb(:, :, 1), ev1, ev2, ev3, ev4)
392 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [2B]", ev1, ev2, " %"
393 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [2B]", ev3, ev4, " %"
394 IF (use_virial) THEN
395 CALL serror(sigma, sdeb(:, :, 1), ev1, ev2)
396 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [2B]", ev1, ev2, " %"
397 END IF
398 END IF
399 END IF
400 END IF
401 ! no contribution from dispersion_3b as q=0 (but q is changed!)
402 ! so we calculate this here
403 IF (grad) THEN
404 IF (dispersion_env%ext_charges) THEN
405 dispersion_env%dcharges = dedq
406 ELSE
407 CALL para_env%sum(dedq)
408 ! Reconstruct full-sized charges for eeq_forces
409 block
410 REAL(KIND=dp), ALLOCATABLE :: q_full(:)
411 ALLOCATE (q_full(natom_full))
412 q_full = 0.0_dp
413 DO i = 1, natom
414 q_full(atom_map_back(i)) = q_red(i)
415 END DO
416 ga(:, :) = 0.0_dp
417 sigma = 0.0_dp
418 CALL eeq_forces(qs_env, q_full, dedq, ga, sigma, dispersion_env%eeq_sparam, &
419 2, enshift, response_only=.true., exclude=exclude_ghost, &
420 cn_max=8.0_dp)
421 DEALLOCATE (q_full)
422 END block
423 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
424 stress = stress + sigma
425 IF (debug) THEN
426 CALL para_env%sum(ga)
427 CALL para_env%sum(sigma)
428 IF (iw > 0) THEN
429 CALL verror(dedq, edq, ev1, ev2)
430 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdq", ev1, ev2, " %"
431 CALL gerror(ga, gdeb(:, :, 2), ev1, ev2, ev3, ev4)
432 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdq]", ev1, ev2, " %"
433 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdq]", ev3, ev4, " %"
434 IF (use_virial) THEN
435 CALL serror(sigma, sdeb(:, :, 2), ev1, ev2)
436 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdq]", ev1, ev2, " %"
437 END IF
438 END IF
439 END IF
440 END IF
441 END IF
442
443 IF (dispersion_env%doabc) THEN
444 ALLOCATE (energies3(natom_full))
445 energies3(:) = 0.0_dp
446 q_red(:) = 0.0_dp
447 ! i.e. dc6dq = dEdq = 0
448 CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
449 !
450 IF (grad) THEN
451 gwdq = 0.0_dp
452 ga(:, :) = 0.0_dp
453 sigma = 0.0_dp
454 END IF
455 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, tvec)
456 CALL dispersion_3b(qs_env, dispersion_env, tvec, cutoff%disp3, disp%r4r2, &
457 gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
458 energies3, dedcn, dedq, grad, ga, sigma, &
459 atom_map, species)
460 IF (grad) THEN
461 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
462 stress = stress + sigma
463 IF (debug) THEN
464 CALL para_env%sum(ga)
465 CALL para_env%sum(sigma)
466 IF (iw > 0) THEN
467 CALL gerror(ga, gdeb(:, :, 3), ev1, ev2, ev3, ev4)
468 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [3B]", ev1, ev2, " %"
469 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [3B]", ev3, ev4, " %"
470 IF (use_virial) THEN
471 CALL serror(sigma, sdeb(:, :, 3), ev1, ev2)
472 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [3B]", ev1, ev2, " %"
473 END IF
474 END IF
475 END IF
476 END IF
477 END IF
478
479 IF (grad) THEN
480 CALL para_env%sum(dedcn)
481 ga(:, :) = 0.0_dp
482 sigma = 0.0_dp
483 CALL dedcn_force(qs_env, dedcn, dcnum, ga, sigma)
484 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
485 stress = stress + sigma
486 IF (debug) THEN
487 CALL para_env%sum(ga)
488 CALL para_env%sum(sigma)
489 IF (iw > 0) THEN
490 CALL verror(dedcn, edcn, ev1, ev2)
491 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdcn", ev1, ev2, " %"
492 CALL gerror(ga, gdeb(:, :, 4), ev1, ev2, ev3, ev4)
493 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdcn]", ev1, ev2, " %"
494 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdcn]", ev3, ev4, " %"
495 IF (use_virial) THEN
496 CALL serror(sigma, sdeb(:, :, 4), ev1, ev2)
497 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdcn]", ev1, ev2, " %"
498 END IF
499 END IF
500 END IF
501 END IF
502 DEALLOCATE (q_red, cn_red)
503 CALL cnumber_release(cn, dcnum, grad)
504 te = m_walltime()
505 tc = tc + te - ts
506
507 IF (debug) THEN
508 ta = sum(energies)
509 CALL para_env%sum(ta)
510 IF (iw > 0) THEN
511 tb = sum(enerd2)
512 ed2 = ta - tb
513 pd2 = abs(ed2)/abs(tb)*100._dp
514 WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 2-body", ed2, pd2, " %"
515 END IF
516 IF (dispersion_env%doabc) THEN
517 ta = sum(energies3)
518 CALL para_env%sum(ta)
519 IF (iw > 0) THEN
520 tb = sum(enerd3)
521 ed3 = ta - tb
522 pd3 = abs(ed3)/abs(tb)*100._dp
523 WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 3-body", ed3, pd3, " %"
524 END IF
525 END IF
526 IF (iw > 0) THEN
527 WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for reference code [s]", td
528 WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for production code [s]", tc
529 END IF
530 END IF
531
532 IF (dispersion_env%doabc) THEN
533 energies(:) = energies(:) + energies3(:)
534 END IF
535 evdw = sum(energies)
536 IF (PRESENT(atomic_energy)) THEN
537 atomic_energy(1:natom_full) = energies(1:natom_full)
538 END IF
539
540 IF (use_virial .AND. calculate_forces) THEN
541 virial%pv_virial = virial%pv_virial - stress
542 END IF
543 IF (calculate_forces) THEN
544 DO iatom = 1, natom_full
545 IF (atom_map(iatom) == 0) cycle
546 ikind = kind_of(iatom)
547 atoma = atom_of_kind(iatom)
548 force(ikind)%dispersion(:, atoma) = &
549 force(ikind)%dispersion(:, atoma) + gradient(:, iatom)
550 END DO
551 END IF
552
553 DEALLOCATE (energies)
554 IF (dispersion_env%doabc) DEALLOCATE (energies3)
555 IF (grad) THEN
556 DEALLOCATE (gradient, ga)
557 END IF
558
559 END IF
560
561 DEALLOCATE (xyz, atomtype, atom_map, atom_map_back, species, exclude_ghost)
562
563 CALL timestop(handle)
564
565 END SUBROUTINE calculate_dispersion_d4_pairpot
566
567! **************************************************************************************************
568!> \brief ...
569!> \param param ...
570!> \param disp ...
571!> \param mol ...
572!> \param cutoff ...
573!> \param grad ...
574!> \param doabc ...
575!> \param enerd2 ...
576!> \param enerd3 ...
577!> \param cnd ...
578!> \param qd ...
579!> \param dEdcn ...
580!> \param dEdq ...
581!> \param gradient ...
582!> \param stress ...
583! **************************************************************************************************
584 SUBROUTINE refd4_debug(param, disp, mol, cutoff, grad, doabc, &
585 enerd2, enerd3, cnd, qd, dEdcn, dEdq, gradient, stress)
586 CLASS(damping_param) :: param
587 TYPE(d4_model) :: disp
588 TYPE(structure_type) :: mol
589 TYPE(realspace_cutoff) :: cutoff
590 TYPE(error_type), ALLOCATABLE :: error
591 LOGICAL, INTENT(IN) :: grad, doabc
592 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: enerd2, enerd3, cnd, qd, dedcn, dedq
593 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gradient
594 REAL(KIND=dp), DIMENSION(3, 3, 4) :: stress
595
596 INTEGER :: mref, natom, i, ncoup
597 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q, qq
598 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: lattr, c6, dc6dcn, dc6dq
599 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cndr, cndl, qdr, qdl, gwdcn, gwdq, gwvec
600
601 mref = maxval(disp%ref)
602 natom = mol%nat
603 ncoup = disp%ncoup
604
605 ! Coordination numbers
606 ALLOCATE (cnd(natom))
607 IF (grad) ALLOCATE (cndr(3, natom, natom), cndl(3, 3, natom))
608 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
609 CALL get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, disp%en, &
610 cnd, cndr, cndl)
611 ! EEQ charges
612 ALLOCATE (qd(natom))
613 IF (grad) ALLOCATE (qdr(3, natom, natom), qdl(3, 3, natom))
614 CALL get_charges(disp%mchrg, mol, error, qd, qdr, qdl)
615 IF (ALLOCATED(error)) THEN
616 cpabort(error%message)
617 END IF
618 ! C6 interpolation
619 ALLOCATE (gwvec(mref, natom, ncoup))
620 IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
621 CALL disp%weight_references(mol, cnd, qd, gwvec, gwdcn, gwdq)
622 ALLOCATE (c6(natom, natom))
623 IF (grad) ALLOCATE (dc6dcn(natom, natom), dc6dq(natom, natom))
624 CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
625 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp2, lattr)
626 !
627 IF (grad) THEN
628 ALLOCATE (gradient(3, natom, 4))
629 gradient = 0.0_dp
630 stress = 0.0_dp
631 END IF
632 !
633 ALLOCATE (enerd2(natom))
634 enerd2(:) = 0.0_dp
635 IF (grad) THEN
636 ALLOCATE (dedcn(natom), dedq(natom))
637 dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
638 END IF
639 CALL param%get_dispersion2(mol, lattr, cutoff%disp2, cutoff%width2, disp%r4r2, c6, &
640 dc6dcn, dc6dq, enerd2, dedcn, dedq, gradient(:, :, 1), &
641 stress(:, :, 1))
642 !
643 IF (grad) THEN
644 DO i = 1, 3
645 gradient(i, :, 2) = matmul(qdr(i, :, :), dedq(:))
646 stress(i, :, 2) = matmul(qdl(i, :, :), dedq(:))
647 END DO
648 END IF
649 !
650 IF (doabc) THEN
651 ALLOCATE (q(natom), qq(natom))
652 q(:) = 0.0_dp; qq(:) = 0.0_dp
653 ALLOCATE (enerd3(natom))
654 enerd3(:) = 0.0_dp
655 CALL disp%weight_references(mol, cnd, q, gwvec, gwdcn, gwdq)
656 CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
657 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, lattr)
658 CALL param%get_dispersion3(mol, lattr, cutoff%disp3, cutoff%width3, disp%r4r2, c6, &
659 dc6dcn, dc6dq, enerd3, dedcn, qq, gradient(:, :, 3), &
660 stress(:, :, 3))
661 END IF
662 IF (grad) THEN
663 DO i = 1, 3
664 gradient(i, :, 4) = matmul(cndr(i, :, :), dedcn(:))
665 stress(i, :, 4) = matmul(cndl(i, :, :), dedcn(:))
666 END DO
667 END IF
668
669 END SUBROUTINE refd4_debug
670
671#else
672
673! **************************************************************************************************
674!> \brief ...
675!> \param qs_env ...
676!> \param dispersion_env ...
677!> \param evdw ...
678!> \param calculate_forces ...
679!> \param iw ...
680!> \param atomic_energy ...
681! **************************************************************************************************
682 SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
683 iw, atomic_energy)
684 TYPE(qs_environment_type), POINTER :: qs_env
685 TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
686 REAL(kind=dp), INTENT(INOUT) :: evdw
687 LOGICAL, INTENT(IN) :: calculate_forces
688 INTEGER, INTENT(IN) :: iw
689 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atomic_energy
690
691 mark_used(qs_env)
692 mark_used(dispersion_env)
693 mark_used(evdw)
694 mark_used(calculate_forces)
695 mark_used(iw)
696 mark_used(atomic_energy)
697
698 cpabort("CP2K build without DFTD4")
699
700 END SUBROUTINE calculate_dispersion_d4_pairpot
701
702#endif
703
704! **************************************************************************************************
705!> \brief ...
706!> \param dispersion_env ...
707!> \param cutoff ...
708!> \param r4r2 ...
709!> \param gwvec ...
710!> \param gwdcn ...
711!> \param gwdq ...
712!> \param c6ref ...
713!> \param mrefs ...
714!> \param energies ...
715!> \param dEdcn ...
716!> \param dEdq ...
717!> \param calculate_forces ...
718!> \param gradient ...
719!> \param stress ...
720!> \param atom_map ...
721!> \param species ...
722! **************************************************************************************************
723 SUBROUTINE dispersion_2b(dispersion_env, cutoff, r4r2, &
724 gwvec, gwdcn, gwdq, c6ref, mrefs, &
725 energies, dEdcn, dEdq, &
726 calculate_forces, gradient, stress, &
727 atom_map, species)
728 TYPE(qs_dispersion_type), POINTER :: dispersion_env
729 REAL(kind=dp), INTENT(IN) :: cutoff
730 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r4r2
731 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
732 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
733 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
734 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: energies, dedcn, dedq
735 LOGICAL, INTENT(IN) :: calculate_forces
736 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient, stress
737 INTEGER, DIMENSION(:), INTENT(IN) :: atom_map, species
738
739 INTEGER :: ia, iatom, ik, ikind, ja, jatom, jk, &
740 jkind, mepos, num_pe
741 REAL(kind=dp) :: a1, a2, c6ij, cutoff2, d6, d8, de, dr2, &
742 edisp, fac, gdisp, r0ij, rrij, s6, s8, &
743 t6, t8
744 REAL(kind=dp), DIMENSION(2) :: dcdcn, dcdq
745 REAL(kind=dp), DIMENSION(3) :: dg, rij
746 REAL(kind=dp), DIMENSION(3, 3) :: ds
747 TYPE(neighbor_list_iterator_p_type), &
748 DIMENSION(:), POINTER :: nl_iterator
749 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
750 POINTER :: sab_vdw
751
752 a1 = dispersion_env%a1
753 a2 = dispersion_env%a2
754 s6 = dispersion_env%s6
755 s8 = dispersion_env%s8
756 cutoff2 = cutoff*cutoff
757
758 sab_vdw => dispersion_env%sab_vdw
759
760 num_pe = 1
761 CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
762
763 mepos = 0
764 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
765 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
766 iatom=iatom, jatom=jatom, r=rij)
767 ! Skip ghost/floating atoms
768 ia = atom_map(iatom)
769 ja = atom_map(jatom)
770 IF (ia == 0 .OR. ja == 0) cycle
771 ! D4 species indices
772 ik = species(iatom)
773 jk = species(jatom)
774 ! vdW potential
775 dr2 = sum(rij(:)**2)
776 IF (dr2 <= cutoff2 .AND. dr2 > 0.0000001_dp) THEN
777 rrij = 3._dp*r4r2(ik)*r4r2(jk)
778 r0ij = a1*sqrt(rrij) + a2
779 IF (calculate_forces) THEN
780 CALL get_c6derivs(c6ij, dcdcn, dcdq, ia, ja, ik, jk, &
781 gwvec, gwdcn, gwdq, c6ref, mrefs)
782 ELSE
783 CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
784 END IF
785 fac = 1._dp
786 IF (iatom == jatom) fac = 0.5_dp
787 t6 = 1.0_dp/(dr2**3 + r0ij**6)
788 t8 = 1.0_dp/(dr2**4 + r0ij**8)
789
790 edisp = (s6*t6 + s8*rrij*t8)*fac
791 de = -c6ij*edisp
792 energies(iatom) = energies(iatom) + de*0.5_dp
793 energies(jatom) = energies(jatom) + de*0.5_dp
794
795 IF (calculate_forces) THEN
796 d6 = -6.0_dp*dr2**2*t6**2
797 d8 = -8.0_dp*dr2**3*t8**2
798 gdisp = (s6*d6 + s8*rrij*d8)*fac
799 dg(:) = -c6ij*gdisp*rij(:)
800 gradient(:, iatom) = gradient(:, iatom) - dg
801 gradient(:, jatom) = gradient(:, jatom) + dg
802 ds(:, :) = spread(dg, 1, 3)*spread(rij, 2, 3)
803 stress(:, :) = stress(:, :) + ds(:, :)
804 dedcn(iatom) = dedcn(iatom) - dcdcn(1)*edisp
805 dedq(iatom) = dedq(iatom) - dcdq(1)*edisp
806 dedcn(jatom) = dedcn(jatom) - dcdcn(2)*edisp
807 dedq(jatom) = dedq(jatom) - dcdq(2)*edisp
808 END IF
809 END IF
810 END DO
811
812 CALL neighbor_list_iterator_release(nl_iterator)
813
814 END SUBROUTINE dispersion_2b
815
816! **************************************************************************************************
817!> \brief ...
818!> \param qs_env ...
819!> \param dispersion_env ...
820!> \param tvec ...
821!> \param cutoff ...
822!> \param r4r2 ...
823!> \param gwvec ...
824!> \param gwdcn ...
825!> \param gwdq ...
826!> \param c6ref ...
827!> \param mrefs ...
828!> \param energies ...
829!> \param dEdcn ...
830!> \param dEdq ...
831!> \param calculate_forces ...
832!> \param gradient ...
833!> \param stress ...
834!> \param atom_map ...
835!> \param species ...
836! **************************************************************************************************
837 SUBROUTINE dispersion_3b(qs_env, dispersion_env, tvec, cutoff, r4r2, &
838 gwvec, gwdcn, gwdq, c6ref, mrefs, &
839 energies, dEdcn, dEdq, &
840 calculate_forces, gradient, stress, &
841 atom_map, species)
842 TYPE(qs_environment_type), POINTER :: qs_env
843 TYPE(qs_dispersion_type), POINTER :: dispersion_env
844 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: tvec
845 REAL(kind=dp), INTENT(IN) :: cutoff
846 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r4r2
847 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
848 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
849 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
850 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: energies, dedcn, dedq
851 LOGICAL, INTENT(IN) :: calculate_forces
852 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient, stress
853 INTEGER, DIMENSION(:), INTENT(IN) :: atom_map, species
854
855 INTEGER :: ia, iatom, ik, ikind, ja, jatom, jk, &
856 jkind, ka, katom, kk, ktr, mepos, &
857 natom, num_pe
858 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
859 INTEGER, DIMENSION(3) :: cell_b
860 REAL(kind=dp) :: a1, a2, alp, ang, c6ij, c6ik, c6jk, c9, &
861 cutoff2, dang, de, dfdmp, fac, fdmp, &
862 r0, r0ij, r0ik, r0jk, r1, r2, r2ij, &
863 r2ik, r2jk, r3, r5, rr, s6, s8, s9
864 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
865 REAL(kind=dp), DIMENSION(2) :: dc6dcnij, dc6dcnik, dc6dcnjk, dc6dqij, &
866 dc6dqik, dc6dqjk
867 REAL(kind=dp), DIMENSION(3) :: dgij, dgik, dgjk, ra, rb, rb0, rij, vij, &
868 vik, vjk
869 REAL(kind=dp), DIMENSION(3, 3) :: ds
870 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
871 TYPE(cell_type), POINTER :: cell
872 TYPE(neighbor_list_iterator_p_type), &
873 DIMENSION(:), POINTER :: nl_iterator
874 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
875 POINTER :: sab_vdw
876 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
877
878 CALL get_qs_env(qs_env=qs_env, natom=natom, cell=cell, &
879 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
880
881 ALLOCATE (rcpbc(3, natom))
882 DO iatom = 1, natom
883 rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
884 END DO
885 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
886
887 a1 = dispersion_env%a1
888 a2 = dispersion_env%a2
889 s6 = dispersion_env%s6
890 s8 = dispersion_env%s8
891 s9 = dispersion_env%s9
892 alp = dispersion_env%alp
893
894 cutoff2 = cutoff**2
895
896 sab_vdw => dispersion_env%sab_vdw
897
898 num_pe = 1
899 CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
900
901 mepos = 0
902 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
903 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
904
905 ! Skip ghost/floating atoms
906 ia = atom_map(iatom)
907 ja = atom_map(jatom)
908 IF (ia == 0 .OR. ja == 0) cycle
909 ik = species(iatom)
910 jk = species(jatom)
911
912 r2ij = sum(rij(:)**2)
913 IF (calculate_forces) THEN
914 CALL get_c6derivs(c6ij, dc6dcnij, dc6dqij, ia, ja, ik, jk, &
915 gwvec, gwdcn, gwdq, c6ref, mrefs)
916 ELSE
917 CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
918 END IF
919 r0ij = a1*sqrt(3._dp*r4r2(jk)*r4r2(ik)) + a2
920 IF (r2ij <= cutoff2 .AND. r2ij > epsilon(1._dp)) THEN
921 CALL get_iterator_info(nl_iterator, cell=cell_b)
922 rb0(:) = matmul(cell%hmat, cell_b)
923 ra(:) = rcpbc(:, iatom)
924 rb(:) = rcpbc(:, jatom) + rb0
925 vij(:) = rb(:) - ra(:)
926
927 DO katom = 1, min(iatom, jatom)
928 ka = atom_map(katom)
929 IF (ka == 0) cycle
930 kk = species(katom)
931 IF (calculate_forces) THEN
932 CALL get_c6derivs(c6ik, dc6dcnik, dc6dqik, ka, ia, kk, ik, &
933 gwvec, gwdcn, gwdq, c6ref, mrefs)
934 CALL get_c6derivs(c6jk, dc6dcnjk, dc6dqjk, ka, ja, kk, jk, &
935 gwvec, gwdcn, gwdq, c6ref, mrefs)
936 ELSE
937 CALL get_c6value(c6ik, ka, ia, kk, ik, gwvec, c6ref, mrefs)
938 CALL get_c6value(c6jk, ka, ja, kk, jk, gwvec, c6ref, mrefs)
939 END IF
940 c9 = -s9*sqrt(abs(c6ij*c6ik*c6jk))
941 r0ik = a1*sqrt(3._dp*r4r2(kk)*r4r2(ik)) + a2
942 r0jk = a1*sqrt(3._dp*r4r2(kk)*r4r2(jk)) + a2
943 r0 = r0ij*r0ik*r0jk
944 fac = triple_scale(iatom, jatom, katom)
945 DO ktr = 1, SIZE(tvec, 2)
946 vik(:) = rcpbc(:, katom) + tvec(:, ktr) - rcpbc(:, iatom)
947 r2ik = vik(1)*vik(1) + vik(2)*vik(2) + vik(3)*vik(3)
948 IF (r2ik > cutoff2 .OR. r2ik < epsilon(1.0_dp)) cycle
949 vjk(:) = rcpbc(:, katom) + tvec(:, ktr) - rb(:)
950 r2jk = vjk(1)*vjk(1) + vjk(2)*vjk(2) + vjk(3)*vjk(3)
951 IF (r2jk > cutoff2 .OR. r2jk < epsilon(1.0_dp)) cycle
952 r2 = r2ij*r2ik*r2jk
953 r1 = sqrt(r2)
954 r3 = r2*r1
955 r5 = r3*r2
956
957 fdmp = 1.0_dp/(1.0_dp + 6.0_dp*(r0/r1)**(alp/3.0_dp))
958 ang = 0.375_dp*(r2ij + r2jk - r2ik)*(r2ij - r2jk + r2ik)* &
959 (-r2ij + r2jk + r2ik)/r5 + 1.0_dp/r3
960
961 rr = ang*fdmp
962 de = rr*c9*fac
963 energies(iatom) = energies(iatom) - de/3._dp
964 energies(jatom) = energies(jatom) - de/3._dp
965 energies(katom) = energies(katom) - de/3._dp
966
967 IF (calculate_forces) THEN
968
969 dfdmp = -2.0_dp*alp*(r0/r1)**(alp/3.0_dp)*fdmp**2
970
971 ! d/drij
972 dang = -0.375_dp*(r2ij**3 + r2ij**2*(r2jk + r2ik) &
973 + r2ij*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ik &
974 + 3.0_dp*r2ik**2) &
975 - 5.0_dp*(r2jk - r2ik)**2*(r2jk + r2ik))/r5
976 dgij(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ij*vij
977
978 ! d/drik
979 dang = -0.375_dp*(r2ik**3 + r2ik**2*(r2jk + r2ij) &
980 + r2ik*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ij &
981 + 3.0_dp*r2ij**2) &
982 - 5.0_dp*(r2jk - r2ij)**2*(r2jk + r2ij))/r5
983 dgik(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ik*vik
984
985 ! d/drjk
986 dang = -0.375_dp*(r2jk**3 + r2jk**2*(r2ik + r2ij) &
987 + r2jk*(3.0_dp*r2ik**2 + 2.0_dp*r2ik*r2ij &
988 + 3.0_dp*r2ij**2) &
989 - 5.0_dp*(r2ik - r2ij)**2*(r2ik + r2ij))/r5
990 dgjk(:) = c9*(-dang*fdmp + ang*dfdmp)/r2jk*vjk
991
992 gradient(:, iatom) = gradient(:, iatom) - dgij - dgik
993 gradient(:, jatom) = gradient(:, jatom) + dgij - dgjk
994 gradient(:, katom) = gradient(:, katom) + dgik + dgjk
995
996 ds(:, :) = spread(dgij, 1, 3)*spread(vij, 2, 3) &
997 + spread(dgik, 1, 3)*spread(vik, 2, 3) &
998 + spread(dgjk, 1, 3)*spread(vjk, 2, 3)
999
1000 stress(:, :) = stress + ds*fac
1001
1002 dedcn(iatom) = dedcn(iatom) - de*0.5_dp &
1003 *(dc6dcnij(1)/c6ij + dc6dcnik(2)/c6ik)
1004 dedcn(jatom) = dedcn(jatom) - de*0.5_dp &
1005 *(dc6dcnij(2)/c6ij + dc6dcnjk(2)/c6jk)
1006 dedcn(katom) = dedcn(katom) - de*0.5_dp &
1007 *(dc6dcnik(1)/c6ik + dc6dcnjk(1)/c6jk)
1008
1009 dedq(iatom) = dedq(iatom) - de*0.5_dp &
1010 *(dc6dqij(1)/c6ij + dc6dqik(2)/c6ik)
1011 dedq(jatom) = dedq(jatom) - de*0.5_dp &
1012 *(dc6dqij(2)/c6ij + dc6dqjk(2)/c6jk)
1013 dedq(katom) = dedq(katom) - de*0.5_dp &
1014 *(dc6dqik(1)/c6ik + dc6dqjk(1)/c6jk)
1015
1016 END IF
1017
1018 END DO
1019 END DO
1020 END IF
1021 END DO
1022
1023 CALL neighbor_list_iterator_release(nl_iterator)
1024
1025 DEALLOCATE (rcpbc)
1026
1027 END SUBROUTINE dispersion_3b
1028
1029! **************************************************************************************************
1030!> \brief ...
1031!> \param ii ...
1032!> \param jj ...
1033!> \param kk ...
1034!> \return ...
1035! **************************************************************************************************
1036 FUNCTION triple_scale(ii, jj, kk) RESULT(triple)
1037 INTEGER, INTENT(IN) :: ii, jj, kk
1038 REAL(kind=dp) :: triple
1039
1040 IF (ii == jj) THEN
1041 IF (ii == kk) THEN
1042 ! ii'i" -> 1/6
1043 triple = 1.0_dp/6.0_dp
1044 ELSE
1045 ! ii'j -> 1/2
1046 triple = 0.5_dp
1047 END IF
1048 ELSE
1049 IF (ii /= kk .AND. jj /= kk) THEN
1050 ! ijk -> 1 (full)
1051 triple = 1.0_dp
1052 ELSE
1053 ! ijj' and iji' -> 1/2
1054 triple = 0.5_dp
1055 END IF
1056 END IF
1057
1058 END FUNCTION triple_scale
1059
1060! **************************************************************************************************
1061!> \brief ...
1062!> \param qs_env ...
1063!> \param dEdcn ...
1064!> \param dcnum ...
1065!> \param gradient ...
1066!> \param stress ...
1067! **************************************************************************************************
1068 SUBROUTINE dedcn_force(qs_env, dEdcn, dcnum, gradient, stress)
1069 TYPE(qs_environment_type), POINTER :: qs_env
1070 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: dedcn
1071 TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
1072 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient
1073 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: stress
1074
1075 CHARACTER(len=*), PARAMETER :: routinen = 'dEdcn_force'
1076
1077 INTEGER :: handle, i, ia, iatom, ikind, katom, &
1078 natom, nkind
1079 LOGICAL :: use_virial
1080 REAL(kind=dp) :: drk
1081 REAL(kind=dp), DIMENSION(3) :: fdik, rik
1082 TYPE(distribution_1d_type), POINTER :: local_particles
1083 TYPE(virial_type), POINTER :: virial
1084
1085 CALL timeset(routinen, handle)
1086
1087 CALL get_qs_env(qs_env, nkind=nkind, natom=natom, &
1088 local_particles=local_particles, &
1089 virial=virial)
1090 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1091
1092 DO ikind = 1, nkind
1093 DO ia = 1, local_particles%n_el(ikind)
1094 iatom = local_particles%list(ikind)%array(ia)
1095 DO i = 1, dcnum(iatom)%neighbors
1096 katom = dcnum(iatom)%nlist(i)
1097 rik = dcnum(iatom)%rik(:, i)
1098 drk = sqrt(sum(rik(:)**2))
1099 fdik(:) = -(dedcn(iatom) + dedcn(katom))*dcnum(iatom)%dvals(i)*rik(:)/drk
1100 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
1101 IF (use_virial) THEN
1102 CALL virial_pair_force(stress, -0.5_dp, fdik, rik)
1103 END IF
1104 END DO
1105 END DO
1106 END DO
1107
1108 CALL timestop(handle)
1109
1110 END SUBROUTINE dedcn_force
1111
1112! **************************************************************************************************
1113!> \brief ...
1114!> \param c6ij ...
1115!> \param ia ...
1116!> \param ja ...
1117!> \param ik ...
1118!> \param jk ...
1119!> \param gwvec ...
1120!> \param c6ref ...
1121!> \param mrefs ...
1122! **************************************************************************************************
1123 SUBROUTINE get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
1124 REAL(kind=dp), INTENT(OUT) :: c6ij
1125 INTEGER, INTENT(IN) :: ia, ja, ik, jk
1126 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec
1127 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
1128 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
1129
1130 INTEGER :: iref, jref
1131 REAL(kind=dp) :: refc6
1132
1133 c6ij = 0.0_dp
1134 DO jref = 1, mrefs(jk)
1135 DO iref = 1, mrefs(ik)
1136 refc6 = c6ref(iref, jref, ik, jk)
1137 c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1138 END DO
1139 END DO
1140
1141 END SUBROUTINE get_c6value
1142
1143! **************************************************************************************************
1144!> \brief ...
1145!> \param c6ij ...
1146!> \param dc6dcn ...
1147!> \param dc6dq ...
1148!> \param ia ...
1149!> \param ja ...
1150!> \param ik ...
1151!> \param jk ...
1152!> \param gwvec ...
1153!> \param gwdcn ...
1154!> \param gwdq ...
1155!> \param c6ref ...
1156!> \param mrefs ...
1157! **************************************************************************************************
1158 SUBROUTINE get_c6derivs(c6ij, dc6dcn, dc6dq, ia, ja, ik, jk, &
1159 gwvec, gwdcn, gwdq, c6ref, mrefs)
1160 REAL(kind=dp), INTENT(OUT) :: c6ij
1161 REAL(kind=dp), DIMENSION(2), INTENT(OUT) :: dc6dcn, dc6dq
1162 INTEGER, INTENT(IN) :: ia, ja, ik, jk
1163 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
1164 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
1165 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
1166
1167 INTEGER :: iref, jref
1168 REAL(kind=dp) :: refc6
1169
1170 c6ij = 0.0_dp
1171 dc6dcn = 0.0_dp
1172 dc6dq = 0.0_dp
1173 DO jref = 1, mrefs(jk)
1174 DO iref = 1, mrefs(ik)
1175 refc6 = c6ref(iref, jref, ik, jk)
1176 c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1177 dc6dcn(1) = dc6dcn(1) + gwdcn(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1178 dc6dcn(2) = dc6dcn(2) + gwvec(iref, ia, 1)*gwdcn(jref, ja, 1)*refc6
1179 dc6dq(1) = dc6dq(1) + gwdq(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1180 dc6dq(2) = dc6dq(2) + gwvec(iref, ia, 1)*gwdq(jref, ja, 1)*refc6
1181 END DO
1182 END DO
1183
1184 END SUBROUTINE get_c6derivs
1185
1186! **************************************************************************************************
1187!> \brief ...
1188!> \param ga ...
1189!> \param gd ...
1190!> \param ev1 ...
1191!> \param ev2 ...
1192!> \param ev3 ...
1193!> \param ev4 ...
1194! **************************************************************************************************
1195 SUBROUTINE gerror(ga, gd, ev1, ev2, ev3, ev4)
1196 REAL(kind=dp), DIMENSION(:, :) :: ga, gd
1197 REAL(kind=dp), INTENT(OUT) :: ev1, ev2, ev3, ev4
1198
1199 INTEGER :: na, np(2)
1200
1201 na = SIZE(ga, 2)
1202 ev1 = sqrt(sum((gd - ga)**2)/na)
1203 ev2 = ev1/sqrt(sum(gd**2)/na)*100._dp
1204 np = maxloc(abs(gd - ga))
1205 ev3 = abs(gd(np(1), np(2)) - ga(np(1), np(2)))
1206 ev4 = abs(gd(np(1), np(2)))
1207 IF (ev4 > 1.e-6_dp) THEN
1208 ev4 = ev3/ev4*100._dp
1209 ELSE
1210 ev4 = 100.0_dp
1211 END IF
1212
1213 END SUBROUTINE gerror
1214
1215! **************************************************************************************************
1216!> \brief ...
1217!> \param sa ...
1218!> \param sd ...
1219!> \param ev1 ...
1220!> \param ev2 ...
1221! **************************************************************************************************
1222 SUBROUTINE serror(sa, sd, ev1, ev2)
1223 REAL(kind=dp), DIMENSION(3, 3) :: sa, sd
1224 REAL(kind=dp), INTENT(OUT) :: ev1, ev2
1225
1226 INTEGER :: i, j
1227 REAL(kind=dp) :: rel
1228
1229 ev1 = maxval(abs(sd - sa))
1230 ev2 = 0.0_dp
1231 DO i = 1, 3
1232 DO j = 1, 3
1233 IF (abs(sd(i, j)) > 1.e-6_dp) THEN
1234 rel = abs(sd(i, j) - sa(i, j))/abs(sd(i, j))*100._dp
1235 ev2 = max(ev2, rel)
1236 END IF
1237 END DO
1238 END DO
1239
1240 END SUBROUTINE serror
1241
1242! **************************************************************************************************
1243!> \brief ...
1244!> \param va ...
1245!> \param vd ...
1246!> \param ev1 ...
1247!> \param ev2 ...
1248! **************************************************************************************************
1249 SUBROUTINE verror(va, vd, ev1, ev2)
1250 REAL(kind=dp), DIMENSION(:) :: va, vd
1251 REAL(kind=dp), INTENT(OUT) :: ev1, ev2
1252
1253 INTEGER :: i, na
1254 REAL(kind=dp) :: rel
1255
1256 na = SIZE(va)
1257 ev1 = maxval(abs(vd - va))
1258 ev2 = 0.0_dp
1259 DO i = 1, na
1260 IF (abs(vd(i)) > 1.e-8_dp) THEN
1261 rel = abs(vd(i) - va(i))/abs(vd(i))*100._dp
1262 ev2 = max(ev2, rel)
1263 END IF
1264 END DO
1265
1266 END SUBROUTINE verror
1267
1268END MODULE qs_dispersion_d4
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:56
cell_types::pbc
Definition cell_types.F:104
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:324
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
eeq_method
Calculation of charge equilibration method.
Definition eeq_method.F:12
eeq_method::eeq_charges
subroutine, public eeq_charges(qs_env, charges, eeq_sparam, eeq_model, enshift_type, exclude, cn_max)
...
Definition eeq_method.F:203
eeq_method::eeq_forces
subroutine, public eeq_forces(qs_env, charges, dcharges, gradient, stress, eeq_sparam, eeq_model, enshift_type, response_only, exclude, cn_max)
...
Definition eeq_method.F:388
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
qs_dispersion_cnum::cnumber_release
subroutine, public cnumber_release(cnumbers, dcnum, derivatives)
...
Definition qs_dispersion_cnum.F:1340
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1278
qs_dispersion_d4
Calculation of dispersion using pair potentials.
Definition qs_dispersion_d4.F:12
qs_dispersion_d4::calculate_dispersion_d4_pairpot
subroutine, public calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, atomic_energy)
...
Definition qs_dispersion_d4.F:684
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_dispersion_utils
Set disperson types for DFT calculations.
Definition qs_dispersion_utils.F:12
qs_dispersion_utils::cellhash
integer function, public cellhash(cell, ncell)
...
Definition qs_dispersion_utils.F:1409
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::set_qs_kind
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, hund_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3323
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26