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