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