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