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