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