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