(git:3db43b4)
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#if defined(__DFTD4_V3)
59 USE dftd4_charge, ONLY: get_charges
60#else
61 USE multicharge, ONLY: get_charges
62 USE mctc_env, ONLY: error_type
63#endif
64!&>
65#endif
66#include "./base/base_uses.f90"
67
68 IMPLICIT NONE
69
70 PRIVATE
71
72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
73
74 PUBLIC :: calculate_dispersion_d4_pairpot
75
76! **************************************************************************************************
77
78CONTAINS
79
80#if defined(__DFTD4)
81! **************************************************************************************************
82!> \brief ...
83!> \param qs_env ...
84!> \param dispersion_env ...
85!> \param evdw ...
86!> \param calculate_forces ...
87!> \param iw ...
88!> \param atomic_energy ...
89! **************************************************************************************************
90 SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, &
91 atomic_energy)
92 TYPE(qs_environment_type), POINTER :: qs_env
93 TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
94 REAL(KIND=dp), INTENT(INOUT) :: evdw
95 LOGICAL, INTENT(IN) :: calculate_forces
96 INTEGER, INTENT(IN) :: iw
97 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atomic_energy
98
99 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
100
101 INTEGER :: atoma, cnfun, enshift, handle, i, iatom
102#if defined(__DFTD4_V3)
103 INTEGER :: ifull, ikind, mref, natom, &
104 natom_full, nghost
105#else
106 INTEGER :: ifull, ikind, mref, natom, ncoup, &
107 natom_full, nghost
108#endif
109 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_map, atom_map_back, atom_of_kind, &
110 atomtype, kind_of, species, t_atomtype
111 INTEGER, DIMENSION(3) :: periodic
112 LOGICAL :: debug, grad, ifloating, ighost, &
113 use_virial
114 LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_ghost, exclude_ghost
115 LOGICAL, DIMENSION(3) :: lperiod
116 REAL(KIND=dp) :: ed2, ed3, ev1, ev2, ev3, ev4, pd2, pd3, &
117 ta, tb, tc, td, te, ts
118 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cn, cn_red, cnd, dedcn, dedq, &
119 edcn, edq, enerd2, &
120 enerd3, energies, energies3, q_red, qd
121#if defined(__DFTD4_V3)
122 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ga, gradient, gwdcn, gwdq, &
123 gwvec, t_xyz, tvec, xyz
124 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gdeb
125#else
126 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ga, gradient, t_xyz, tvec, xyz
127 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gdeb, gwdcn, gwdq, gwvec
128#endif
129 REAL(KIND=dp), DIMENSION(3, 3) :: sigma, stress
130 REAL(KIND=dp), DIMENSION(3, 3, 4) :: sdeb
131 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
132 TYPE(cell_type), POINTER :: cell
133 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
134 TYPE(mp_para_env_type), POINTER :: para_env
135 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
136 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
137 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
138 TYPE(virial_type), POINTER :: virial
139
140 CLASS(damping_param), ALLOCATABLE :: param
141 TYPE(d4_model) :: disp
142 TYPE(structure_type) :: mol
143 TYPE(realspace_cutoff) :: cutoff
144
145#if !defined(__DFTD4_V3)
146 TYPE(error_type), ALLOCATABLE :: error
147#endif
148
149 CALL timeset(routinen, handle)
150
151 debug = dispersion_env%d4_debug
152
153 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
154 cell=cell, force=force, virial=virial, para_env=para_env)
155 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
156
157 !get information about particles
158 natom_full = SIZE(particle_set)
159 nghost = 0
160 ALLOCATE (t_xyz(3, natom_full), t_atomtype(natom_full), a_ghost(natom_full))
161 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
162 DO iatom = 1, natom_full
163 t_xyz(:, iatom) = particle_set(iatom)%r(:)
164 ikind = kind_of(iatom)
165 CALL get_qs_kind(qs_kind_set(ikind), zatom=t_atomtype(iatom), ghost=ighost, floating=ifloating)
166 a_ghost(iatom) = ighost .OR. ifloating
167 IF (a_ghost(iatom)) nghost = nghost + 1
168 END DO
169
170 natom = natom_full - nghost
171 ! Build atom mapping: full index -> reduced index (0 for ghost)
172 ALLOCATE (atom_map(natom_full), atom_map_back(natom))
173 atom_map = 0
174 iatom = 0
175 ALLOCATE (xyz(3, natom), atomtype(natom))
176 DO i = 1, natom_full
177 IF (.NOT. a_ghost(i)) THEN
178 iatom = iatom + 1
179 atom_map(i) = iatom
180 atom_map_back(iatom) = i
181 xyz(:, iatom) = t_xyz(:, i)
182 atomtype(iatom) = t_atomtype(i)
183 END IF
184 END DO
185 DEALLOCATE (a_ghost, t_xyz, t_atomtype)
186
187 !get information about cell / lattice
188 CALL get_cell(cell=cell, periodic=periodic)
189 lperiod(1) = periodic(1) == 1
190 lperiod(2) = periodic(2) == 1
191 lperiod(3) = periodic(3) == 1
192 ! enforce en shift method 1 (original/molecular)
193 ! method 2 from paper on PBC seems not to work
194 enshift = 1
195 !IF (ALL(periodic == 0)) enshift = 1
196
197 !prepare for the call to the dispersion function
198 CALL new(mol, atomtype, xyz, lattice=cell%hmat, periodic=lperiod)
199#if defined(__DFTD4_V3)
200 CALL new_d4_model(disp, mol)
201#else
202 CALL new_d4_model(error, disp, mol)
203 IF (ALLOCATED(error)) THEN
204 cpabort(error%message)
205 END IF
206
207 ! Number of coupling
208 ncoup = disp%ncoup
209#endif
210
211 ! Build species mapping: full atom index -> D4 species ID (0 for ghost)
212 ALLOCATE (species(natom_full))
213 species = 0
214 DO i = 1, natom
215 species(atom_map_back(i)) = mol%id(i)
216 END DO
217
218 ! Build per-kind exclusion mask for EEQ (ghost/floating kinds excluded)
219 ALLOCATE (exclude_ghost(SIZE(qs_kind_set)))
220 exclude_ghost = .false.
221 DO i = 1, SIZE(qs_kind_set)
222 CALL get_qs_kind(qs_kind_set(i), ghost=ighost, floating=ifloating)
223 exclude_ghost(i) = ighost .OR. ifloating
224 END DO
225
226 IF (dispersion_env%ref_functional == "none") THEN
227 CALL get_rational_damping("pbe", param, s9=0.0_dp)
228 SELECT TYPE (param)
229 TYPE is (rational_damping_param)
230 param%s6 = dispersion_env%s6
231 param%s8 = dispersion_env%s8
232 param%a1 = dispersion_env%a1
233 param%a2 = dispersion_env%a2
234 param%alp = dispersion_env%alp
235 END SELECT
236 ELSE
237 CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
238 SELECT TYPE (param)
239 TYPE is (rational_damping_param)
240 dispersion_env%s6 = param%s6
241 dispersion_env%s8 = param%s8
242 dispersion_env%a1 = param%a1
243 dispersion_env%a2 = param%a2
244 dispersion_env%alp = param%alp
245 END SELECT
246 END IF
247
248 ! Coordination number cutoff
249 cutoff%cn = dispersion_env%rc_cn
250 ! Two-body interaction cutoff
251 cutoff%disp2 = dispersion_env%rc_d4*2._dp
252 ! Three-body interaction cutoff
253 cutoff%disp3 = dispersion_env%rc_disp*2._dp
254 IF (cutoff%disp3 > cutoff%disp2) THEN
255 cpabort("D4: Three-body cutoff should be smaller than two-body cutoff")
256 END IF
257
258 IF (calculate_forces) THEN
259 grad = .true.
260 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
261 ELSE
262 grad = .false.
263 use_virial = .false.
264 END IF
265
266 IF (dispersion_env%d4_reference_code) THEN
267
268 !> Wrapper to handle the evaluation of dispersion energy and derivatives
269 IF (.NOT. dispersion_env%doabc) THEN
270 cpwarn("Using D4_REFERENCE_CODE enforces calculation of C9 term.")
271 END IF
272 IF (grad) THEN
273 ALLOCATE (gradient(3, natom))
274 CALL get_dispersion(mol, disp, param, cutoff, evdw, gradient, stress)
275 IF (calculate_forces) THEN
276 IF (use_virial) THEN
277 virial%pv_virial = virial%pv_virial - stress/para_env%num_pe
278 END IF
279 DO iatom = 1, natom
280 ifull = atom_map_back(iatom)
281 ikind = kind_of(ifull)
282 atoma = atom_of_kind(ifull)
283 force(ikind)%dispersion(:, atoma) = &
284 force(ikind)%dispersion(:, atoma) + gradient(:, iatom)/para_env%num_pe
285 END DO
286 END IF
287 DEALLOCATE (gradient)
288 ELSE
289 CALL get_dispersion(mol, disp, param, cutoff, evdw)
290 END IF
291 !dispersion energy is computed by every MPI process
292 evdw = evdw/para_env%num_pe
293 IF (dispersion_env%ext_charges) dispersion_env%dcharges = 0.0_dp
294 IF (PRESENT(atomic_energy)) THEN
295 cpwarn("Atomic energies not available for D4 reference code")
296 atomic_energy = 0.0_dp
297 END IF
298
299 ELSE
300
301 IF (iw > 0) THEN
302 WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
303 WRITE (iw, fmt="(T32,A)") "DEBUG D4 DISPERSION"
304 WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
305 WRITE (iw, '(A,T71,A10)') " DEBUG D4| Reference functional ", trim(dispersion_env%ref_functional)
306 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s6) ", dispersion_env%s6
307 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s8) ", dispersion_env%s8
308 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a1) ", dispersion_env%a1
309 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a2) ", dispersion_env%a2
310 WRITE (iw, '(A,T71,E10.4)') " DEBUG D4| Cutoff value coordination numbers ", dispersion_env%eps_cn
311 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius coordination numbers ", dispersion_env%rc_cn
312 WRITE (iw, '(A,T71,I10)') " DEBUG D4| Coordination number function type ", dispersion_env%cnfun
313 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 2-body terms [bohr]", 2._dp*dispersion_env%rc_d4
314 WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 3-body terms [bohr]", 2._dp*dispersion_env%rc_disp
315 END IF
316
317 td = 0.0_dp
318 IF (debug .AND. iw > 0) THEN
319 ts = m_walltime()
320 CALL refd4_debug(param, disp, mol, cutoff, grad, dispersion_env%doabc, &
321 enerd2, enerd3, cnd, qd, edcn, edq, gdeb, sdeb)
322 te = m_walltime()
323 td = te - ts
324 END IF
325
326 tc = 0.0_dp
327 ts = m_walltime()
328
329 mref = maxval(disp%ref)
330 ! Coordination numbers (full-size from qs_env; ghosts excluded internally)
331 cnfun = dispersion_env%cnfun
332 CALL cnumber_init(qs_env, cn, dcnum, cnfun, grad)
333 ! cn has size natom_full; ghost entries are 0
334
335 ! Filter CN to reduced space for D4 model
336 ALLOCATE (cn_red(natom))
337 DO i = 1, natom
338 cn_red(i) = cn(atom_map_back(i))
339 END DO
340 IF (debug .AND. iw > 0) THEN
341 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (max)", maxval(abs(cn_red - cnd))
342 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (ave)", sum(abs(cn_red - cnd))/natom
343 END IF
344
345 ! EEQ charges
346 ! Use CP2K's MPI-parallel EEQ solver with ghost exclusion.
347 ! Ghost kinds get huge hardness + zero coupling -> q_ghost = 0 exactly,
348 ! without influencing real atom charges. Preserves full MPI parallelism.
349 IF (dispersion_env%ext_charges) THEN
350 ALLOCATE (q_red(natom))
351 q_red(1:natom) = dispersion_env%charges(1:natom)
352 ELSE
353 ! eeq_charges writes natom_full entries (from qs_env), so allocate full-size
354 ALLOCATE (q_red(natom_full))
355 CALL eeq_charges(qs_env, q_red, dispersion_env%eeq_sparam, 2, enshift, &
356 exclude=exclude_ghost, cn_max=8.0_dp)
357 ! Filter to reduced size for D4 model
358 block
359 REAL(KIND=dp), ALLOCATABLE :: q_tmp(:)
360 ALLOCATE (q_tmp(natom))
361 DO i = 1, natom
362 q_tmp(i) = q_red(atom_map_back(i))
363 END DO
364 DEALLOCATE (q_red)
365 CALL move_alloc(q_tmp, q_red)
366 END block
367 END IF
368 IF (debug .AND. iw > 0) THEN
369 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (max)", maxval(abs(q_red - qd))
370 WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (ave)", sum(abs(q_red - qd))/natom
371 END IF
372 ! Weights for C6 calculation (reduced space)
373#if defined(__DFTD4_V3)
374 ALLOCATE (gwvec(mref, natom))
375 IF (grad) ALLOCATE (gwdcn(mref, natom), gwdq(mref, natom))
376#else
377 ALLOCATE (gwvec(mref, natom, ncoup))
378 IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
379#endif
380 CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
381
382 ! Energies and derivatives (full-size for CP2K infrastructure compatibility)
383 ALLOCATE (energies(natom_full))
384 energies(:) = 0.0_dp
385 IF (grad) THEN
386 ALLOCATE (gradient(3, natom_full), ga(3, natom_full))
387 ALLOCATE (dedcn(natom_full), dedq(natom_full))
388 dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
389 ga(:, :) = 0.0_dp
390 sigma(:, :) = 0.0_dp
391 END IF
392 CALL dispersion_2b(dispersion_env, cutoff%disp2, disp%r4r2, &
393 gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
394 energies, dedcn, dedq, grad, ga, sigma, &
395 atom_map, species)
396 IF (grad) THEN
397 gradient(1:3, 1:natom_full) = ga(1:3, 1:natom_full)
398 stress = sigma
399 IF (debug) THEN
400 CALL para_env%sum(ga)
401 CALL para_env%sum(sigma)
402 IF (iw > 0) THEN
403 CALL gerror(ga, gdeb(:, :, 1), ev1, ev2, ev3, ev4)
404 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [2B]", ev1, ev2, " %"
405 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [2B]", ev3, ev4, " %"
406 IF (use_virial) THEN
407 CALL serror(sigma, sdeb(:, :, 1), ev1, ev2)
408 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [2B]", ev1, ev2, " %"
409 END IF
410 END IF
411 END IF
412 END IF
413 ! no contribution from dispersion_3b as q=0 (but q is changed!)
414 ! so we calculate this here
415 IF (grad) THEN
416 IF (dispersion_env%ext_charges) THEN
417 dispersion_env%dcharges = dedq
418 ELSE
419 CALL para_env%sum(dedq)
420 ! Reconstruct full-sized charges for eeq_forces
421 block
422 REAL(KIND=dp), ALLOCATABLE :: q_full(:)
423 ALLOCATE (q_full(natom_full))
424 q_full = 0.0_dp
425 DO i = 1, natom
426 q_full(atom_map_back(i)) = q_red(i)
427 END DO
428 ga(:, :) = 0.0_dp
429 sigma = 0.0_dp
430 CALL eeq_forces(qs_env, q_full, dedq, ga, sigma, dispersion_env%eeq_sparam, &
431 2, enshift, response_only=.true., exclude=exclude_ghost, &
432 cn_max=8.0_dp)
433 DEALLOCATE (q_full)
434 END block
435 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
436 stress = stress + sigma
437 IF (debug) THEN
438 CALL para_env%sum(ga)
439 CALL para_env%sum(sigma)
440 IF (iw > 0) THEN
441 CALL verror(dedq, edq, ev1, ev2)
442 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdq", ev1, ev2, " %"
443 CALL gerror(ga, gdeb(:, :, 2), ev1, ev2, ev3, ev4)
444 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdq]", ev1, ev2, " %"
445 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdq]", ev3, ev4, " %"
446 IF (use_virial) THEN
447 CALL serror(sigma, sdeb(:, :, 2), ev1, ev2)
448 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdq]", ev1, ev2, " %"
449 END IF
450 END IF
451 END IF
452 END IF
453 END IF
454
455 IF (dispersion_env%doabc) THEN
456 ALLOCATE (energies3(natom_full))
457 energies3(:) = 0.0_dp
458 q_red(:) = 0.0_dp
459 ! i.e. dc6dq = dEdq = 0
460 CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
461 !
462 IF (grad) THEN
463 gwdq = 0.0_dp
464 ga(:, :) = 0.0_dp
465 sigma = 0.0_dp
466 END IF
467 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, tvec)
468 CALL dispersion_3b(qs_env, dispersion_env, tvec, cutoff%disp3, disp%r4r2, &
469 gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
470 energies3, dedcn, dedq, grad, ga, sigma, &
471 atom_map, species)
472 IF (grad) THEN
473 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
474 stress = stress + sigma
475 IF (debug) THEN
476 CALL para_env%sum(ga)
477 CALL para_env%sum(sigma)
478 IF (iw > 0) THEN
479 CALL gerror(ga, gdeb(:, :, 3), ev1, ev2, ev3, ev4)
480 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [3B]", ev1, ev2, " %"
481 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [3B]", ev3, ev4, " %"
482 IF (use_virial) THEN
483 CALL serror(sigma, sdeb(:, :, 3), ev1, ev2)
484 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [3B]", ev1, ev2, " %"
485 END IF
486 END IF
487 END IF
488 END IF
489 END IF
490
491 IF (grad) THEN
492 CALL para_env%sum(dedcn)
493 ga(:, :) = 0.0_dp
494 sigma = 0.0_dp
495 CALL dedcn_force(qs_env, dedcn, dcnum, ga, sigma)
496 gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
497 stress = stress + sigma
498 IF (debug) THEN
499 CALL para_env%sum(ga)
500 CALL para_env%sum(sigma)
501 IF (iw > 0) THEN
502 CALL verror(dedcn, edcn, ev1, ev2)
503 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdcn", ev1, ev2, " %"
504 CALL gerror(ga, gdeb(:, :, 4), ev1, ev2, ev3, ev4)
505 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdcn]", ev1, ev2, " %"
506 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdcn]", ev3, ev4, " %"
507 IF (use_virial) THEN
508 CALL serror(sigma, sdeb(:, :, 4), ev1, ev2)
509 WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdcn]", ev1, ev2, " %"
510 END IF
511 END IF
512 END IF
513 END IF
514 DEALLOCATE (q_red, cn_red)
515 CALL cnumber_release(cn, dcnum, grad)
516 te = m_walltime()
517 tc = tc + te - ts
518
519 IF (debug) THEN
520 ta = sum(energies)
521 CALL para_env%sum(ta)
522 IF (iw > 0) THEN
523 tb = sum(enerd2)
524 ed2 = ta - tb
525 pd2 = abs(ed2)/abs(tb)*100.
526 WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 2-body", ed2, pd2, " %"
527 END IF
528 IF (dispersion_env%doabc) THEN
529 ta = sum(energies3)
530 CALL para_env%sum(ta)
531 IF (iw > 0) THEN
532 tb = sum(enerd3)
533 ed3 = ta - tb
534 pd3 = abs(ed3)/abs(tb)*100.
535 WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 3-body", ed3, pd3, " %"
536 END IF
537 END IF
538 IF (iw > 0) THEN
539 WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for reference code [s]", td
540 WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for production code [s]", tc
541 END IF
542 END IF
543
544 IF (dispersion_env%doabc) THEN
545 energies(:) = energies(:) + energies3(:)
546 END IF
547 evdw = sum(energies)
548 IF (PRESENT(atomic_energy)) THEN
549 atomic_energy(1:natom_full) = energies(1:natom_full)
550 END IF
551
552 IF (use_virial .AND. calculate_forces) THEN
553 virial%pv_virial = virial%pv_virial - stress
554 END IF
555 IF (calculate_forces) THEN
556 DO iatom = 1, natom_full
557 IF (atom_map(iatom) == 0) cycle
558 ikind = kind_of(iatom)
559 atoma = atom_of_kind(iatom)
560 force(ikind)%dispersion(:, atoma) = &
561 force(ikind)%dispersion(:, atoma) + gradient(:, iatom)
562 END DO
563 END IF
564
565 DEALLOCATE (energies)
566 IF (dispersion_env%doabc) DEALLOCATE (energies3)
567 IF (grad) THEN
568 DEALLOCATE (gradient, ga)
569 END IF
570
571 END IF
572
573 DEALLOCATE (xyz, atomtype, atom_map, atom_map_back, species, exclude_ghost)
574
575 CALL timestop(handle)
576
577 END SUBROUTINE calculate_dispersion_d4_pairpot
578
579! **************************************************************************************************
580!> \brief ...
581!> \param param ...
582!> \param disp ...
583!> \param mol ...
584!> \param cutoff ...
585!> \param grad ...
586!> \param doabc ...
587!> \param enerd2 ...
588!> \param enerd3 ...
589!> \param cnd ...
590!> \param qd ...
591!> \param dEdcn ...
592!> \param dEdq ...
593!> \param gradient ...
594!> \param stress ...
595! **************************************************************************************************
596 SUBROUTINE refd4_debug(param, disp, mol, cutoff, grad, doabc, &
597 enerd2, enerd3, cnd, qd, dEdcn, dEdq, gradient, stress)
598 CLASS(damping_param) :: param
599 TYPE(d4_model) :: disp
600 TYPE(structure_type) :: mol
601 TYPE(realspace_cutoff) :: cutoff
602#if !defined(__DFTD4_V3)
603 TYPE(error_type), ALLOCATABLE :: error
604#endif
605 LOGICAL, INTENT(IN) :: grad, doabc
606 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: enerd2, enerd3, cnd, qd, dedcn, dedq
607 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gradient
608 REAL(KIND=dp), DIMENSION(3, 3, 4) :: stress
609
610#if defined(__DFTD4_V3)
611 INTEGER :: mref, natom, i
612 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q, qq
613 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: lattr, gwdcn, gwdq, gwvec, &
614 c6, dc6dcn, dc6dq
615 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cndr, cndl, qdr, qdl
616#else
617 INTEGER :: mref, natom, i, ncoup
618 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q, qq
619 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: lattr, c6, dc6dcn, dc6dq
620 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cndr, cndl, qdr, qdl, gwdcn, gwdq, gwvec
621#endif
622
623 mref = maxval(disp%ref)
624 natom = mol%nat
625#if !defined(__DFTD4_V3)
626 ncoup = disp%ncoup
627#endif
628
629 ! Coordination numbers
630 ALLOCATE (cnd(natom))
631 IF (grad) ALLOCATE (cndr(3, natom, natom), cndl(3, 3, natom))
632 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
633 CALL get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, disp%en, &
634 cnd, cndr, cndl)
635 ! EEQ charges
636 ALLOCATE (qd(natom))
637 IF (grad) ALLOCATE (qdr(3, natom, natom), qdl(3, 3, natom))
638#if defined(__DFTD4_V3)
639 CALL get_charges(mol, qd, qdr, qdl)
640#else
641 CALL get_charges(disp%mchrg, mol, error, qd, qdr, qdl)
642 IF (ALLOCATED(error)) THEN
643 cpabort(error%message)
644 END IF
645#endif
646 ! C6 interpolation
647#if defined(__DFTD4_V3)
648 ALLOCATE (gwvec(mref, natom))
649 IF (grad) ALLOCATE (gwdcn(mref, natom), gwdq(mref, natom))
650#else
651 ALLOCATE (gwvec(mref, natom, ncoup))
652 IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
653#endif
654 CALL disp%weight_references(mol, cnd, qd, gwvec, gwdcn, gwdq)
655 ALLOCATE (c6(natom, natom))
656 IF (grad) ALLOCATE (dc6dcn(natom, natom), dc6dq(natom, natom))
657 CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
658 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp2, lattr)
659 !
660 IF (grad) THEN
661 ALLOCATE (gradient(3, natom, 4))
662 gradient = 0.0_dp
663 stress = 0.0_dp
664 END IF
665 !
666 ALLOCATE (enerd2(natom))
667 enerd2(:) = 0.0_dp
668 IF (grad) THEN
669 ALLOCATE (dedcn(natom), dedq(natom))
670 dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
671 END IF
672 CALL param%get_dispersion2(mol, lattr, cutoff%disp2, disp%r4r2, c6, dc6dcn, dc6dq, &
673 enerd2, dedcn, dedq, gradient(:, :, 1), stress(:, :, 1))
674 !
675 IF (grad) THEN
676 DO i = 1, 3
677 gradient(i, :, 2) = matmul(qdr(i, :, :), dedq(:))
678 stress(i, :, 2) = matmul(qdl(i, :, :), dedq(:))
679 END DO
680 END IF
681 !
682 IF (doabc) THEN
683 ALLOCATE (q(natom), qq(natom))
684 q(:) = 0.0_dp; qq(:) = 0.0_dp
685 ALLOCATE (enerd3(natom))
686 enerd3(:) = 0.0_dp
687 CALL disp%weight_references(mol, cnd, q, gwvec, gwdcn, gwdq)
688 CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
689 CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, lattr)
690 CALL param%get_dispersion3(mol, lattr, cutoff%disp3, disp%r4r2, c6, dc6dcn, dc6dq, &
691 enerd3, dedcn, qq, gradient(:, :, 3), stress(:, :, 3))
692 END IF
693 IF (grad) THEN
694 DO i = 1, 3
695 gradient(i, :, 4) = matmul(cndr(i, :, :), dedcn(:))
696 stress(i, :, 4) = matmul(cndl(i, :, :), dedcn(:))
697 END DO
698 END IF
699
700 END SUBROUTINE refd4_debug
701
702#else
703
704! **************************************************************************************************
705!> \brief ...
706!> \param qs_env ...
707!> \param dispersion_env ...
708!> \param evdw ...
709!> \param calculate_forces ...
710!> \param iw ...
711!> \param atomic_energy ...
712! **************************************************************************************************
713 SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
714 iw, atomic_energy)
715 TYPE(qs_environment_type), POINTER :: qs_env
716 TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
717 REAL(kind=dp), INTENT(INOUT) :: evdw
718 LOGICAL, INTENT(IN) :: calculate_forces
719 INTEGER, INTENT(IN) :: iw
720 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atomic_energy
721
722 mark_used(qs_env)
723 mark_used(dispersion_env)
724 mark_used(evdw)
725 mark_used(calculate_forces)
726 mark_used(iw)
727 mark_used(atomic_energy)
728
729 cpabort("CP2K build without DFTD4")
730
731 END SUBROUTINE calculate_dispersion_d4_pairpot
732
733#endif
734
735! **************************************************************************************************
736!> \brief ...
737!> \param dispersion_env ...
738!> \param cutoff ...
739!> \param r4r2 ...
740!> \param gwvec ...
741!> \param gwdcn ...
742!> \param gwdq ...
743!> \param c6ref ...
744!> \param mrefs ...
745!> \param energies ...
746!> \param dEdcn ...
747!> \param dEdq ...
748!> \param calculate_forces ...
749!> \param gradient ...
750!> \param stress ...
751!> \param atom_map ...
752!> \param species ...
753! **************************************************************************************************
754 SUBROUTINE dispersion_2b(dispersion_env, cutoff, r4r2, &
755 gwvec, gwdcn, gwdq, c6ref, mrefs, &
756 energies, dEdcn, dEdq, &
757 calculate_forces, gradient, stress, &
758 atom_map, species)
759 TYPE(qs_dispersion_type), POINTER :: dispersion_env
760 REAL(kind=dp), INTENT(IN) :: cutoff
761 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r4r2
762#if defined(__DFTD4_V3)
763 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gwvec, gwdcn, gwdq
764#else
765 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
766#endif
767 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
768 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
769 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: energies, dedcn, dedq
770 LOGICAL, INTENT(IN) :: calculate_forces
771 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient, stress
772 INTEGER, DIMENSION(:), INTENT(IN) :: atom_map, species
773
774 INTEGER :: ia, iatom, ik, ikind, ja, jatom, jk, &
775 jkind, mepos, num_pe
776 REAL(kind=dp) :: a1, a2, c6ij, cutoff2, d6, d8, de, dr2, &
777 edisp, fac, gdisp, r0ij, rrij, s6, s8, &
778 t6, t8
779 REAL(kind=dp), DIMENSION(2) :: dcdcn, dcdq
780 REAL(kind=dp), DIMENSION(3) :: dg, rij
781 REAL(kind=dp), DIMENSION(3, 3) :: ds
782 TYPE(neighbor_list_iterator_p_type), &
783 DIMENSION(:), POINTER :: nl_iterator
784 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
785 POINTER :: sab_vdw
786
787 a1 = dispersion_env%a1
788 a2 = dispersion_env%a2
789 s6 = dispersion_env%s6
790 s8 = dispersion_env%s8
791 cutoff2 = cutoff*cutoff
792
793 sab_vdw => dispersion_env%sab_vdw
794
795 num_pe = 1
796 CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
797
798 mepos = 0
799 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
800 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
801 iatom=iatom, jatom=jatom, r=rij)
802 ! Skip ghost/floating atoms
803 ia = atom_map(iatom)
804 ja = atom_map(jatom)
805 IF (ia == 0 .OR. ja == 0) cycle
806 ! D4 species indices
807 ik = species(iatom)
808 jk = species(jatom)
809 ! vdW potential
810 dr2 = sum(rij(:)**2)
811 IF (dr2 <= cutoff2 .AND. dr2 > 0.0000001_dp) THEN
812 rrij = 3._dp*r4r2(ik)*r4r2(jk)
813 r0ij = a1*sqrt(rrij) + a2
814 IF (calculate_forces) THEN
815 CALL get_c6derivs(c6ij, dcdcn, dcdq, ia, ja, ik, jk, &
816 gwvec, gwdcn, gwdq, c6ref, mrefs)
817 ELSE
818 CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
819 END IF
820 fac = 1._dp
821 IF (iatom == jatom) fac = 0.5_dp
822 t6 = 1.0_dp/(dr2**3 + r0ij**6)
823 t8 = 1.0_dp/(dr2**4 + r0ij**8)
824
825 edisp = (s6*t6 + s8*rrij*t8)*fac
826 de = -c6ij*edisp
827 energies(iatom) = energies(iatom) + de*0.5_dp
828 energies(jatom) = energies(jatom) + de*0.5_dp
829
830 IF (calculate_forces) THEN
831 d6 = -6.0_dp*dr2**2*t6**2
832 d8 = -8.0_dp*dr2**3*t8**2
833 gdisp = (s6*d6 + s8*rrij*d8)*fac
834 dg(:) = -c6ij*gdisp*rij(:)
835 gradient(:, iatom) = gradient(:, iatom) - dg
836 gradient(:, jatom) = gradient(:, jatom) + dg
837 ds(:, :) = spread(dg, 1, 3)*spread(rij, 2, 3)
838 stress(:, :) = stress(:, :) + ds(:, :)
839 dedcn(iatom) = dedcn(iatom) - dcdcn(1)*edisp
840 dedq(iatom) = dedq(iatom) - dcdq(1)*edisp
841 dedcn(jatom) = dedcn(jatom) - dcdcn(2)*edisp
842 dedq(jatom) = dedq(jatom) - dcdq(2)*edisp
843 END IF
844 END IF
845 END DO
846
847 CALL neighbor_list_iterator_release(nl_iterator)
848
849 END SUBROUTINE dispersion_2b
850
851! **************************************************************************************************
852!> \brief ...
853!> \param qs_env ...
854!> \param dispersion_env ...
855!> \param tvec ...
856!> \param cutoff ...
857!> \param r4r2 ...
858!> \param gwvec ...
859!> \param gwdcn ...
860!> \param gwdq ...
861!> \param c6ref ...
862!> \param mrefs ...
863!> \param energies ...
864!> \param dEdcn ...
865!> \param dEdq ...
866!> \param calculate_forces ...
867!> \param gradient ...
868!> \param stress ...
869!> \param atom_map ...
870!> \param species ...
871! **************************************************************************************************
872 SUBROUTINE dispersion_3b(qs_env, dispersion_env, tvec, cutoff, r4r2, &
873 gwvec, gwdcn, gwdq, c6ref, mrefs, &
874 energies, dEdcn, dEdq, &
875 calculate_forces, gradient, stress, &
876 atom_map, species)
877 TYPE(qs_environment_type), POINTER :: qs_env
878 TYPE(qs_dispersion_type), POINTER :: dispersion_env
879 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: tvec
880 REAL(kind=dp), INTENT(IN) :: cutoff
881 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r4r2
882#if defined(__DFTD4_V3)
883 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gwvec, gwdcn, gwdq
884#else
885 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
886#endif
887 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
888 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
889 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: energies, dedcn, dedq
890 LOGICAL, INTENT(IN) :: calculate_forces
891 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient, stress
892 INTEGER, DIMENSION(:), INTENT(IN) :: atom_map, species
893
894 INTEGER :: ia, iatom, ik, ikind, ja, jatom, jk, &
895 jkind, ka, katom, kk, ktr, mepos, &
896 natom, num_pe
897 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
898 INTEGER, DIMENSION(3) :: cell_b
899 REAL(kind=dp) :: a1, a2, alp, ang, c6ij, c6ik, c6jk, c9, &
900 cutoff2, dang, de, dfdmp, fac, fdmp, &
901 r0, r0ij, r0ik, r0jk, r1, r2, r2ij, &
902 r2ik, r2jk, r3, r5, rr, s6, s8, s9
903 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
904 REAL(kind=dp), DIMENSION(2) :: dc6dcnij, dc6dcnik, dc6dcnjk, dc6dqij, &
905 dc6dqik, dc6dqjk
906 REAL(kind=dp), DIMENSION(3) :: dgij, dgik, dgjk, ra, rb, rb0, rij, vij, &
907 vik, vjk
908 REAL(kind=dp), DIMENSION(3, 3) :: ds
909 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
910 TYPE(cell_type), POINTER :: cell
911 TYPE(neighbor_list_iterator_p_type), &
912 DIMENSION(:), POINTER :: nl_iterator
913 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
914 POINTER :: sab_vdw
915 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
916
917 CALL get_qs_env(qs_env=qs_env, natom=natom, cell=cell, &
918 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
919
920 ALLOCATE (rcpbc(3, natom))
921 DO iatom = 1, natom
922 rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
923 END DO
924 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
925
926 a1 = dispersion_env%a1
927 a2 = dispersion_env%a2
928 s6 = dispersion_env%s6
929 s8 = dispersion_env%s8
930 s9 = dispersion_env%s9
931 alp = dispersion_env%alp
932
933 cutoff2 = cutoff**2
934
935 sab_vdw => dispersion_env%sab_vdw
936
937 num_pe = 1
938 CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
939
940 mepos = 0
941 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
942 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
943
944 ! Skip ghost/floating atoms
945 ia = atom_map(iatom)
946 ja = atom_map(jatom)
947 IF (ia == 0 .OR. ja == 0) cycle
948 ik = species(iatom)
949 jk = species(jatom)
950
951 r2ij = sum(rij(:)**2)
952 IF (calculate_forces) THEN
953 CALL get_c6derivs(c6ij, dc6dcnij, dc6dqij, ia, ja, ik, jk, &
954 gwvec, gwdcn, gwdq, c6ref, mrefs)
955 ELSE
956 CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
957 END IF
958 r0ij = a1*sqrt(3._dp*r4r2(jk)*r4r2(ik)) + a2
959 IF (r2ij <= cutoff2 .AND. r2ij > epsilon(1._dp)) THEN
960 CALL get_iterator_info(nl_iterator, cell=cell_b)
961 rb0(:) = matmul(cell%hmat, cell_b)
962 ra(:) = rcpbc(:, iatom)
963 rb(:) = rcpbc(:, jatom) + rb0
964 vij(:) = rb(:) - ra(:)
965
966 DO katom = 1, min(iatom, jatom)
967 ka = atom_map(katom)
968 IF (ka == 0) cycle
969 kk = species(katom)
970 IF (calculate_forces) THEN
971 CALL get_c6derivs(c6ik, dc6dcnik, dc6dqik, ka, ia, kk, ik, &
972 gwvec, gwdcn, gwdq, c6ref, mrefs)
973 CALL get_c6derivs(c6jk, dc6dcnjk, dc6dqjk, ka, ja, kk, jk, &
974 gwvec, gwdcn, gwdq, c6ref, mrefs)
975 ELSE
976 CALL get_c6value(c6ik, ka, ia, kk, ik, gwvec, c6ref, mrefs)
977 CALL get_c6value(c6jk, ka, ja, kk, jk, gwvec, c6ref, mrefs)
978 END IF
979 c9 = -s9*sqrt(abs(c6ij*c6ik*c6jk))
980 r0ik = a1*sqrt(3._dp*r4r2(kk)*r4r2(ik)) + a2
981 r0jk = a1*sqrt(3._dp*r4r2(kk)*r4r2(jk)) + a2
982 r0 = r0ij*r0ik*r0jk
983 fac = triple_scale(iatom, jatom, katom)
984 DO ktr = 1, SIZE(tvec, 2)
985 vik(:) = rcpbc(:, katom) + tvec(:, ktr) - rcpbc(:, iatom)
986 r2ik = vik(1)*vik(1) + vik(2)*vik(2) + vik(3)*vik(3)
987 IF (r2ik > cutoff2 .OR. r2ik < epsilon(1.0_dp)) cycle
988 vjk(:) = rcpbc(:, katom) + tvec(:, ktr) - rb(:)
989 r2jk = vjk(1)*vjk(1) + vjk(2)*vjk(2) + vjk(3)*vjk(3)
990 IF (r2jk > cutoff2 .OR. r2jk < epsilon(1.0_dp)) cycle
991 r2 = r2ij*r2ik*r2jk
992 r1 = sqrt(r2)
993 r3 = r2*r1
994 r5 = r3*r2
995
996 fdmp = 1.0_dp/(1.0_dp + 6.0_dp*(r0/r1)**(alp/3.0_dp))
997 ang = 0.375_dp*(r2ij + r2jk - r2ik)*(r2ij - r2jk + r2ik)* &
998 (-r2ij + r2jk + r2ik)/r5 + 1.0_dp/r3
999
1000 rr = ang*fdmp
1001 de = rr*c9*fac
1002 energies(iatom) = energies(iatom) - de/3._dp
1003 energies(jatom) = energies(jatom) - de/3._dp
1004 energies(katom) = energies(katom) - de/3._dp
1005
1006 IF (calculate_forces) THEN
1007
1008 dfdmp = -2.0_dp*alp*(r0/r1)**(alp/3.0_dp)*fdmp**2
1009
1010 ! d/drij
1011 dang = -0.375_dp*(r2ij**3 + r2ij**2*(r2jk + r2ik) &
1012 + r2ij*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ik &
1013 + 3.0_dp*r2ik**2) &
1014 - 5.0_dp*(r2jk - r2ik)**2*(r2jk + r2ik))/r5
1015 dgij(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ij*vij
1016
1017 ! d/drik
1018 dang = -0.375_dp*(r2ik**3 + r2ik**2*(r2jk + r2ij) &
1019 + r2ik*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ij &
1020 + 3.0_dp*r2ij**2) &
1021 - 5.0_dp*(r2jk - r2ij)**2*(r2jk + r2ij))/r5
1022 dgik(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ik*vik
1023
1024 ! d/drjk
1025 dang = -0.375_dp*(r2jk**3 + r2jk**2*(r2ik + r2ij) &
1026 + r2jk*(3.0_dp*r2ik**2 + 2.0_dp*r2ik*r2ij &
1027 + 3.0_dp*r2ij**2) &
1028 - 5.0_dp*(r2ik - r2ij)**2*(r2ik + r2ij))/r5
1029 dgjk(:) = c9*(-dang*fdmp + ang*dfdmp)/r2jk*vjk
1030
1031 gradient(:, iatom) = gradient(:, iatom) - dgij - dgik
1032 gradient(:, jatom) = gradient(:, jatom) + dgij - dgjk
1033 gradient(:, katom) = gradient(:, katom) + dgik + dgjk
1034
1035 ds(:, :) = spread(dgij, 1, 3)*spread(vij, 2, 3) &
1036 + spread(dgik, 1, 3)*spread(vik, 2, 3) &
1037 + spread(dgjk, 1, 3)*spread(vjk, 2, 3)
1038
1039 stress(:, :) = stress + ds*fac
1040
1041 dedcn(iatom) = dedcn(iatom) - de*0.5_dp &
1042 *(dc6dcnij(1)/c6ij + dc6dcnik(2)/c6ik)
1043 dedcn(jatom) = dedcn(jatom) - de*0.5_dp &
1044 *(dc6dcnij(2)/c6ij + dc6dcnjk(2)/c6jk)
1045 dedcn(katom) = dedcn(katom) - de*0.5_dp &
1046 *(dc6dcnik(1)/c6ik + dc6dcnjk(1)/c6jk)
1047
1048 dedq(iatom) = dedq(iatom) - de*0.5_dp &
1049 *(dc6dqij(1)/c6ij + dc6dqik(2)/c6ik)
1050 dedq(jatom) = dedq(jatom) - de*0.5_dp &
1051 *(dc6dqij(2)/c6ij + dc6dqjk(2)/c6jk)
1052 dedq(katom) = dedq(katom) - de*0.5_dp &
1053 *(dc6dqik(1)/c6ik + dc6dqjk(1)/c6jk)
1054
1055 END IF
1056
1057 END DO
1058 END DO
1059 END IF
1060 END DO
1061
1062 CALL neighbor_list_iterator_release(nl_iterator)
1063
1064 DEALLOCATE (rcpbc)
1065
1066 END SUBROUTINE dispersion_3b
1067
1068! **************************************************************************************************
1069!> \brief ...
1070!> \param ii ...
1071!> \param jj ...
1072!> \param kk ...
1073!> \return ...
1074! **************************************************************************************************
1075 FUNCTION triple_scale(ii, jj, kk) RESULT(triple)
1076 INTEGER, INTENT(IN) :: ii, jj, kk
1077 REAL(kind=dp) :: triple
1078
1079 IF (ii == jj) THEN
1080 IF (ii == kk) THEN
1081 ! ii'i" -> 1/6
1082 triple = 1.0_dp/6.0_dp
1083 ELSE
1084 ! ii'j -> 1/2
1085 triple = 0.5_dp
1086 END IF
1087 ELSE
1088 IF (ii /= kk .AND. jj /= kk) THEN
1089 ! ijk -> 1 (full)
1090 triple = 1.0_dp
1091 ELSE
1092 ! ijj' and iji' -> 1/2
1093 triple = 0.5_dp
1094 END IF
1095 END IF
1096
1097 END FUNCTION triple_scale
1098
1099! **************************************************************************************************
1100!> \brief ...
1101!> \param qs_env ...
1102!> \param dEdcn ...
1103!> \param dcnum ...
1104!> \param gradient ...
1105!> \param stress ...
1106! **************************************************************************************************
1107 SUBROUTINE dedcn_force(qs_env, dEdcn, dcnum, gradient, stress)
1108 TYPE(qs_environment_type), POINTER :: qs_env
1109 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: dedcn
1110 TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
1111 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient
1112 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: stress
1113
1114 CHARACTER(len=*), PARAMETER :: routinen = 'dEdcn_force'
1115
1116 INTEGER :: handle, i, ia, iatom, ikind, katom, &
1117 natom, nkind
1118 LOGICAL :: use_virial
1119 REAL(kind=dp) :: drk
1120 REAL(kind=dp), DIMENSION(3) :: fdik, rik
1121 TYPE(distribution_1d_type), POINTER :: local_particles
1122 TYPE(virial_type), POINTER :: virial
1123
1124 CALL timeset(routinen, handle)
1125
1126 CALL get_qs_env(qs_env, nkind=nkind, natom=natom, &
1127 local_particles=local_particles, &
1128 virial=virial)
1129 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1130
1131 DO ikind = 1, nkind
1132 DO ia = 1, local_particles%n_el(ikind)
1133 iatom = local_particles%list(ikind)%array(ia)
1134 DO i = 1, dcnum(iatom)%neighbors
1135 katom = dcnum(iatom)%nlist(i)
1136 rik = dcnum(iatom)%rik(:, i)
1137 drk = sqrt(sum(rik(:)**2))
1138 fdik(:) = -(dedcn(iatom) + dedcn(katom))*dcnum(iatom)%dvals(i)*rik(:)/drk
1139 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
1140 IF (use_virial) THEN
1141 CALL virial_pair_force(stress, -0.5_dp, fdik, rik)
1142 END IF
1143 END DO
1144 END DO
1145 END DO
1146
1147 CALL timestop(handle)
1148
1149 END SUBROUTINE dedcn_force
1150
1151! **************************************************************************************************
1152!> \brief ...
1153!> \param c6ij ...
1154!> \param ia ...
1155!> \param ja ...
1156!> \param ik ...
1157!> \param jk ...
1158!> \param gwvec ...
1159!> \param c6ref ...
1160!> \param mrefs ...
1161! **************************************************************************************************
1162 SUBROUTINE get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
1163 REAL(kind=dp), INTENT(OUT) :: c6ij
1164 INTEGER, INTENT(IN) :: ia, ja, ik, jk
1165#if defined(__DFTD4_V3)
1166 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gwvec
1167#else
1168 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec
1169#endif
1170 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
1171 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
1172
1173 INTEGER :: iref, jref
1174 REAL(kind=dp) :: refc6
1175
1176 c6ij = 0.0_dp
1177 DO jref = 1, mrefs(jk)
1178 DO iref = 1, mrefs(ik)
1179 refc6 = c6ref(iref, jref, ik, jk)
1180#if defined(__DFTD4_V3)
1181 c6ij = c6ij + gwvec(iref, ia)*gwvec(jref, ja)*refc6
1182#else
1183 c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1184#endif
1185 END DO
1186 END DO
1187
1188 END SUBROUTINE get_c6value
1189
1190! **************************************************************************************************
1191!> \brief ...
1192!> \param c6ij ...
1193!> \param dc6dcn ...
1194!> \param dc6dq ...
1195!> \param ia ...
1196!> \param ja ...
1197!> \param ik ...
1198!> \param jk ...
1199!> \param gwvec ...
1200!> \param gwdcn ...
1201!> \param gwdq ...
1202!> \param c6ref ...
1203!> \param mrefs ...
1204! **************************************************************************************************
1205 SUBROUTINE get_c6derivs(c6ij, dc6dcn, dc6dq, ia, ja, ik, jk, &
1206 gwvec, gwdcn, gwdq, c6ref, mrefs)
1207 REAL(kind=dp), INTENT(OUT) :: c6ij
1208 REAL(kind=dp), DIMENSION(2), INTENT(OUT) :: dc6dcn, dc6dq
1209 INTEGER, INTENT(IN) :: ia, ja, ik, jk
1210#if defined(__DFTD4_V3)
1211 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gwvec, gwdcn, gwdq
1212#else
1213 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gwvec, gwdcn, gwdq
1214#endif
1215 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: c6ref
1216 INTEGER, DIMENSION(:), INTENT(IN) :: mrefs
1217
1218 INTEGER :: iref, jref
1219 REAL(kind=dp) :: refc6
1220
1221 c6ij = 0.0_dp
1222 dc6dcn = 0.0_dp
1223 dc6dq = 0.0_dp
1224 DO jref = 1, mrefs(jk)
1225 DO iref = 1, mrefs(ik)
1226 refc6 = c6ref(iref, jref, ik, jk)
1227#if defined(__DFTD4_V3)
1228 c6ij = c6ij + gwvec(iref, ia)*gwvec(jref, ja)*refc6
1229 dc6dcn(1) = dc6dcn(1) + gwdcn(iref, ia)*gwvec(jref, ja)*refc6
1230 dc6dcn(2) = dc6dcn(2) + gwvec(iref, ia)*gwdcn(jref, ja)*refc6
1231 dc6dq(1) = dc6dq(1) + gwdq(iref, ia)*gwvec(jref, ja)*refc6
1232 dc6dq(2) = dc6dq(2) + gwvec(iref, ia)*gwdq(jref, ja)*refc6
1233#else
1234 c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1235 dc6dcn(1) = dc6dcn(1) + gwdcn(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1236 dc6dcn(2) = dc6dcn(2) + gwvec(iref, ia, 1)*gwdcn(jref, ja, 1)*refc6
1237 dc6dq(1) = dc6dq(1) + gwdq(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
1238 dc6dq(2) = dc6dq(2) + gwvec(iref, ia, 1)*gwdq(jref, ja, 1)*refc6
1239#endif
1240 END DO
1241 END DO
1242
1243 END SUBROUTINE get_c6derivs
1244
1245! **************************************************************************************************
1246!> \brief ...
1247!> \param ga ...
1248!> \param gd ...
1249!> \param ev1 ...
1250!> \param ev2 ...
1251!> \param ev3 ...
1252!> \param ev4 ...
1253! **************************************************************************************************
1254 SUBROUTINE gerror(ga, gd, ev1, ev2, ev3, ev4)
1255 REAL(kind=dp), DIMENSION(:, :) :: ga, gd
1256 REAL(kind=dp), INTENT(OUT) :: ev1, ev2, ev3, ev4
1257
1258 INTEGER :: na, np(2)
1259
1260 na = SIZE(ga, 2)
1261 ev1 = sqrt(sum((gd - ga)**2)/na)
1262 ev2 = ev1/sqrt(sum(gd**2)/na)*100._dp
1263 np = maxloc(abs(gd - ga))
1264 ev3 = abs(gd(np(1), np(2)) - ga(np(1), np(2)))
1265 ev4 = abs(gd(np(1), np(2)))
1266 IF (ev4 > 1.e-6) THEN
1267 ev4 = ev3/ev4*100._dp
1268 ELSE
1269 ev4 = 100.0_dp
1270 END IF
1271
1272 END SUBROUTINE gerror
1273
1274! **************************************************************************************************
1275!> \brief ...
1276!> \param sa ...
1277!> \param sd ...
1278!> \param ev1 ...
1279!> \param ev2 ...
1280! **************************************************************************************************
1281 SUBROUTINE serror(sa, sd, ev1, ev2)
1282 REAL(kind=dp), DIMENSION(3, 3) :: sa, sd
1283 REAL(kind=dp), INTENT(OUT) :: ev1, ev2
1284
1285 INTEGER :: i, j
1286 REAL(kind=dp) :: rel
1287
1288 ev1 = maxval(abs(sd - sa))
1289 ev2 = 0.0_dp
1290 DO i = 1, 3
1291 DO j = 1, 3
1292 IF (abs(sd(i, j)) > 1.e-6_dp) THEN
1293 rel = abs(sd(i, j) - sa(i, j))/abs(sd(i, j))*100._dp
1294 ev2 = max(ev2, rel)
1295 END IF
1296 END DO
1297 END DO
1298
1299 END SUBROUTINE serror
1300
1301! **************************************************************************************************
1302!> \brief ...
1303!> \param va ...
1304!> \param vd ...
1305!> \param ev1 ...
1306!> \param ev2 ...
1307! **************************************************************************************************
1308 SUBROUTINE verror(va, vd, ev1, ev2)
1309 REAL(kind=dp), DIMENSION(:) :: va, vd
1310 REAL(kind=dp), INTENT(OUT) :: ev1, ev2
1311
1312 INTEGER :: i, na
1313 REAL(kind=dp) :: rel
1314
1315 na = SIZE(va)
1316 ev1 = maxval(abs(vd - va))
1317 ev2 = 0.0_dp
1318 DO i = 1, na
1319 IF (abs(vd(i)) > 1.e-8_dp) THEN
1320 rel = abs(vd(i) - va(i))/abs(vd(i))*100._dp
1321 ev2 = max(ev2, rel)
1322 END IF
1323 END DO
1324
1325 END SUBROUTINE verror
1326
1327END 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:51
cell_types::pbc
Definition cell_types.F:97
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:200
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:257
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:195
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:380
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:136
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:153
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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:1426
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1364
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:715
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:1202
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, 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:530
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, 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, 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:468
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, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3026
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:549
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:743
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:220
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:177
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