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 D2 dispersion with Grimme or Chai-Head-Gordon damping
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dispersion_d2
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind,&
16 : get_atomic_kind_set
17 : USE atprop_types, ONLY: atprop_array_init,&
18 : atprop_type
19 : USE cell_types, ONLY: cell_type
20 : USE input_constants, ONLY: vdw_pairpot_dftd2chg
21 : USE kinds, ONLY: dp
22 : USE physcon, ONLY: bohr,&
23 : kjmol
24 : USE qs_dispersion_types, ONLY: qs_atom_dispersion_type,&
25 : qs_dispersion_type
26 : USE qs_environment_types, ONLY: get_qs_env,&
27 : qs_environment_type
28 : USE qs_force_types, ONLY: qs_force_type
29 : USE qs_kind_types, ONLY: get_qs_kind,&
30 : qs_kind_type
31 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
32 : neighbor_list_iterate,&
33 : neighbor_list_iterator_create,&
34 : neighbor_list_iterator_p_type,&
35 : neighbor_list_iterator_release,&
36 : neighbor_list_set_p_type
37 : USE virial_methods, ONLY: virial_pair_force
38 : USE virial_types, ONLY: virial_type
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d2'
46 :
47 : PUBLIC :: calculate_dispersion_d2_pairpot, dftd2_param
48 :
49 : ! **************************************************************************************************
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief ...
55 : !> \param qs_env ...
56 : !> \param dispersion_env ...
57 : !> \param evdw ...
58 : !> \param calculate_forces ...
59 : !> \param atevdw ...
60 : ! **************************************************************************************************
61 74 : SUBROUTINE calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
62 :
63 : TYPE(qs_environment_type), POINTER :: qs_env
64 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
65 : REAL(KIND=dp), INTENT(OUT) :: evdw
66 : LOGICAL, INTENT(IN) :: calculate_forces
67 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atevdw
68 :
69 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d2_pairpot'
70 :
71 : INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
72 : jatom, jkind, mepos, natom, nkind, &
73 : num_pe, za, zb
74 74 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
75 : LOGICAL :: atenergy, atex, chg_damping, floating_a, &
76 : ghost_a, use_virial
77 74 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: dodisp, floating, ghost
78 : REAL(KIND=dp) :: c6, dd, devdw, dfdmp, dr, er, fac, fdmp, &
79 : rcc, rcut, s6, xp
80 74 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c6d2, radd2
81 : REAL(KIND=dp), DIMENSION(3) :: fdij, rij
82 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread
83 74 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
84 74 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
85 : TYPE(atprop_type), POINTER :: atprop
86 : TYPE(cell_type), POINTER :: cell
87 : TYPE(neighbor_list_iterator_p_type), &
88 74 : DIMENSION(:), POINTER :: nl_iterator
89 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
90 74 : POINTER :: sab_vdw
91 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a
92 74 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
93 74 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
94 : TYPE(virial_type), POINTER :: virial
95 :
96 74 : CALL timeset(routineN, handle)
97 :
98 74 : evdw = 0._dp
99 :
100 74 : NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
101 :
102 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
103 74 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, atprop=atprop)
104 :
105 : ! atomic energy and stress arrays
106 74 : atenergy = atprop%energy
107 74 : IF (atenergy) THEN
108 0 : CALL atprop_array_init(atprop%atevdw, natom)
109 0 : atener => atprop%atevdw
110 : END IF
111 : ! external atomic energy
112 74 : atex = .FALSE.
113 74 : IF (PRESENT(atevdw)) THEN
114 0 : atex = .TRUE.
115 : END IF
116 :
117 74 : NULLIFY (force)
118 74 : CALL get_qs_env(qs_env=qs_env, force=force)
119 74 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
120 74 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
121 74 : pv_virial_thread(:, :) = 0._dp
122 :
123 666 : ALLOCATE (dodisp(nkind), ghost(nkind), floating(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
124 236 : DO ikind = 1, nkind
125 162 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
126 162 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
127 162 : dodisp(ikind) = disp_a%defined
128 162 : ghost(ikind) = ghost_a
129 162 : floating(ikind) = floating_a
130 162 : atomnumber(ikind) = za
131 162 : c6d2(ikind) = disp_a%c6
132 398 : radd2(ikind) = disp_a%vdw_radii
133 : END DO
134 :
135 74 : rcut = 2._dp*dispersion_env%rc_disp
136 74 : s6 = dispersion_env%scaling
137 74 : dd = dispersion_env%exp_pre
138 74 : chg_damping = dispersion_env%pp_type == vdw_pairpot_dftd2chg
139 :
140 74 : sab_vdw => dispersion_env%sab_vdw
141 74 : num_pe = 1
142 74 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
143 :
144 74 : mepos = 0
145 33110 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
146 33036 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
147 :
148 33036 : IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) CYCLE
149 :
150 32538 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
151 :
152 : za = atomnumber(ikind)
153 : zb = atomnumber(jkind)
154 : ! vdW potential
155 130152 : dr = SQRT(SUM(rij(:)**2))
156 32612 : IF (dr <= rcut) THEN
157 32538 : fac = 1._dp
158 32538 : IF (iatom == jatom) fac = 0.5_dp
159 32538 : IF (dr > 0.001_dp) THEN
160 32315 : c6 = SQRT(c6d2(ikind)*c6d2(jkind))
161 32315 : rcc = radd2(ikind) + radd2(jkind)
162 32315 : IF (chg_damping) THEN
163 : ! Chai-Head-Gordon damping, using the same radii as Grimme D2.
164 3 : er = 6._dp*(rcc/dr)**12
165 : ELSE
166 32312 : er = EXP(-dd*(dr/rcc - 1._dp))
167 : END IF
168 32315 : fdmp = 1._dp/(1._dp + er)
169 32315 : xp = s6*c6/dr**6
170 32315 : evdw = evdw - xp*fdmp*fac
171 32315 : IF (calculate_forces) THEN
172 31435 : IF (chg_damping) THEN
173 3 : dfdmp = 12._dp/dr*er*fdmp*fdmp
174 : ELSE
175 31432 : dfdmp = dd/rcc*er*fdmp*fdmp
176 : END IF
177 31435 : devdw = -xp*(-6._dp*fdmp/dr + dfdmp)
178 125740 : fdij(:) = devdw*rij(:)/dr*fac
179 31435 : atom_a = atom_of_kind(iatom)
180 31435 : atom_b = atom_of_kind(jatom)
181 125740 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
182 125740 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
183 31435 : IF (use_virial) THEN
184 9842 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
185 : END IF
186 : END IF
187 32315 : IF (atenergy) THEN
188 0 : atener(iatom) = atener(iatom) - 0.5_dp*xp*fdmp*fac
189 0 : atener(jatom) = atener(jatom) - 0.5_dp*xp*fdmp*fac
190 : END IF
191 32315 : IF (atex) THEN
192 0 : atevdw(iatom) = atevdw(iatom) - 0.5_dp*xp*fdmp*fac
193 0 : atevdw(jatom) = atevdw(jatom) - 0.5_dp*xp*fdmp*fac
194 : END IF
195 : END IF
196 : END IF
197 :
198 : END DO
199 :
200 962 : virial%pv_virial = virial%pv_virial + pv_virial_thread
201 :
202 74 : CALL neighbor_list_iterator_release(nl_iterator)
203 :
204 74 : DEALLOCATE (dodisp, ghost, floating, atomnumber, radd2, c6d2)
205 :
206 74 : CALL timestop(handle)
207 :
208 148 : END SUBROUTINE calculate_dispersion_d2_pairpot
209 :
210 : ! **************************************************************************************************
211 : !> \brief ...
212 : !> \param z ...
213 : !> \param c6 ...
214 : !> \param r ...
215 : !> \param found ...
216 : ! **************************************************************************************************
217 70 : SUBROUTINE dftd2_param(z, c6, r, found)
218 :
219 : INTEGER, INTENT(in) :: z
220 : REAL(KIND=dp), INTENT(inout) :: c6, r
221 : LOGICAL, INTENT(inout) :: found
222 :
223 : REAL(KIND=dp), DIMENSION(54), PARAMETER :: c6val = [0.14_dp, 0.08_dp, 1.61_dp, 1.61_dp, &
224 : 3.13_dp, 1.75_dp, 1.23_dp, 0.70_dp, 0.75_dp, 0.63_dp, 5.71_dp, 5.71_dp, 10.79_dp, 9.23_dp,&
225 : 7.84_dp, 5.57_dp, 5.07_dp, 4.61_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, &
226 : 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 16.99_dp, 17.10_dp, &
227 : 16.37_dp, 12.64_dp, 12.47_dp, 12.01_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, &
228 : 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 37.32_dp, 38.71_dp, &
229 : 38.44_dp, 31.74_dp, 31.50_dp, 29.99_dp]
230 : REAL(KIND=dp), DIMENSION(54), PARAMETER :: rval = [1.001_dp, 1.012_dp, 0.825_dp, 1.408_dp, &
231 : 1.485_dp, 1.452_dp, 1.397_dp, 1.342_dp, 1.287_dp, 1.243_dp, 1.144_dp, 1.364_dp, 1.639_dp, &
232 : 1.716_dp, 1.705_dp, 1.683_dp, 1.639_dp, 1.595_dp, 1.485_dp, 1.474_dp, 1.562_dp, 1.562_dp, &
233 : 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.650_dp, &
234 : 1.727_dp, 1.760_dp, 1.771_dp, 1.749_dp, 1.727_dp, 1.628_dp, 1.606_dp, 1.639_dp, 1.639_dp, &
235 : 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.672_dp, &
236 : 1.804_dp, 1.881_dp, 1.892_dp, 1.892_dp, 1.881_dp]
237 :
238 : !
239 : ! GRIMME DISPERSION PARAMETERS
240 : ! Stefan Grimme, Semiempirical GGA-Type Density Functional Constructed
241 : ! with a Long-Range Dispersion Correction, J. Comp. Chem. 27: 1787-1799 (2006)
242 : ! doi:10.1002/jcc.20495
243 : !
244 : ! Conversion factor [Jnm^6mol^-1] -> [a.u.] : 17.34527758021901
245 : ! Conversion factor [A] -> [a.u.] : 1.889726132885643
246 : !
247 : ! C6 values in [Jnm^6/mol]
248 : ! vdW radii [A]
249 :
250 70 : IF (z > 0 .AND. z <= 54) THEN
251 70 : found = .TRUE.
252 70 : c6 = c6val(z)*1000._dp*bohr**6/kjmol
253 70 : r = rval(z)*bohr
254 : ELSE
255 0 : found = .FALSE.
256 : END IF
257 :
258 70 : END SUBROUTINE dftd2_param
259 :
260 : ! **************************************************************************************************
261 :
262 : END MODULE qs_dispersion_d2
|
