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