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