Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of D2 dispersion
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 kinds, ONLY: dp
21 : USE physcon, ONLY: bohr,&
22 : kjmol
23 : USE qs_dispersion_types, ONLY: qs_atom_dispersion_type,&
24 : qs_dispersion_type
25 : USE qs_environment_types, ONLY: get_qs_env,&
26 : qs_environment_type
27 : USE qs_force_types, ONLY: qs_force_type
28 : USE qs_kind_types, ONLY: get_qs_kind,&
29 : qs_kind_type
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 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d2'
45 :
46 : PUBLIC :: calculate_dispersion_d2_pairpot, dftd2_param
47 :
48 : ! **************************************************************************************************
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief ...
54 : !> \param qs_env ...
55 : !> \param dispersion_env ...
56 : !> \param evdw ...
57 : !> \param calculate_forces ...
58 : !> \param atevdw ...
59 : ! **************************************************************************************************
60 58 : SUBROUTINE calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
61 :
62 : TYPE(qs_environment_type), POINTER :: qs_env
63 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
64 : REAL(KIND=dp), INTENT(OUT) :: evdw
65 : LOGICAL, INTENT(IN) :: calculate_forces
66 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atevdw
67 :
68 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d2_pairpot'
69 :
70 : INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
71 : jatom, jkind, mepos, natom, nkind, &
72 : num_pe, za, zb
73 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
74 : LOGICAL :: atenergy, atex, floating_a, ghost_a, &
75 : use_virial
76 58 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: dodisp, floating, ghost
77 : REAL(KIND=dp) :: c6, dd, devdw, dfdmp, dr, er, fac, fdmp, &
78 : rcc, rcut, s6, xp
79 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c6d2, radd2
80 : REAL(KIND=dp), DIMENSION(3) :: fdij, rij
81 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread
82 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
83 58 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
84 : TYPE(atprop_type), POINTER :: atprop
85 : TYPE(cell_type), POINTER :: cell
86 : TYPE(neighbor_list_iterator_p_type), &
87 58 : DIMENSION(:), POINTER :: nl_iterator
88 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
89 58 : POINTER :: sab_vdw
90 : TYPE(qs_atom_dispersion_type), POINTER :: disp_a
91 58 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
92 58 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
93 : TYPE(virial_type), POINTER :: virial
94 :
95 58 : CALL timeset(routineN, handle)
96 :
97 58 : evdw = 0._dp
98 :
99 58 : NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
100 :
101 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
102 58 : qs_kind_set=qs_kind_set, cell=cell, virial=virial, atprop=atprop)
103 :
104 : ! atomic energy and stress arrays
105 58 : atenergy = atprop%energy
106 58 : IF (atenergy) THEN
107 0 : CALL atprop_array_init(atprop%atevdw, natom)
108 0 : atener => atprop%atevdw
109 : END IF
110 : ! external atomic energy
111 58 : atex = .FALSE.
112 58 : IF (PRESENT(atevdw)) THEN
113 0 : atex = .TRUE.
114 : END IF
115 :
116 58 : NULLIFY (force)
117 58 : CALL get_qs_env(qs_env=qs_env, force=force)
118 58 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
119 58 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
120 58 : pv_virial_thread(:, :) = 0._dp
121 :
122 522 : ALLOCATE (dodisp(nkind), ghost(nkind), floating(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
123 188 : DO ikind = 1, nkind
124 130 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
125 130 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
126 130 : dodisp(ikind) = disp_a%defined
127 130 : ghost(ikind) = ghost_a
128 130 : floating(ikind) = floating_a
129 130 : atomnumber(ikind) = za
130 130 : c6d2(ikind) = disp_a%c6
131 318 : radd2(ikind) = disp_a%vdw_radii
132 : END DO
133 :
134 58 : rcut = 2._dp*dispersion_env%rc_disp
135 58 : s6 = dispersion_env%scaling
136 58 : dd = dispersion_env%exp_pre
137 :
138 58 : sab_vdw => dispersion_env%sab_vdw
139 58 : num_pe = 1
140 58 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
141 :
142 58 : mepos = 0
143 32703 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
144 32645 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
145 :
146 32645 : IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) CYCLE
147 :
148 32147 : IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
149 :
150 32147 : za = atomnumber(ikind)
151 32147 : zb = atomnumber(jkind)
152 : ! vdW potential
153 128588 : dr = SQRT(SUM(rij(:)**2))
154 32205 : IF (dr <= rcut) THEN
155 32147 : fac = 1._dp
156 32147 : IF (iatom == jatom) fac = 0.5_dp
157 32147 : IF (dr > 0.001_dp) THEN
158 31997 : c6 = SQRT(c6d2(ikind)*c6d2(jkind))
159 31997 : rcc = radd2(ikind) + radd2(jkind)
160 31997 : er = EXP(-dd*(dr/rcc - 1._dp))
161 31997 : fdmp = 1._dp/(1._dp + er)
162 31997 : xp = s6*c6/dr**6
163 31997 : evdw = evdw - xp*fdmp*fac
164 31997 : IF (calculate_forces) THEN
165 31432 : dfdmp = dd/rcc*er*fdmp*fdmp
166 31432 : devdw = -xp*(-6._dp*fdmp/dr + dfdmp)
167 125728 : fdij(:) = devdw*rij(:)/dr*fac
168 31432 : atom_a = atom_of_kind(iatom)
169 31432 : atom_b = atom_of_kind(jatom)
170 125728 : force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
171 125728 : force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
172 31432 : IF (use_virial) THEN
173 9842 : CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
174 : END IF
175 : END IF
176 31997 : IF (atenergy) THEN
177 0 : atener(iatom) = atener(iatom) - 0.5_dp*xp*fdmp*fac
178 0 : atener(jatom) = atener(jatom) - 0.5_dp*xp*fdmp*fac
179 : END IF
180 31997 : IF (atex) THEN
181 0 : atevdw(iatom) = atevdw(iatom) - 0.5_dp*xp*fdmp*fac
182 0 : atevdw(jatom) = atevdw(jatom) - 0.5_dp*xp*fdmp*fac
183 : END IF
184 : END IF
185 : END IF
186 :
187 : END DO
188 :
189 754 : virial%pv_virial = virial%pv_virial + pv_virial_thread
190 :
191 58 : CALL neighbor_list_iterator_release(nl_iterator)
192 :
193 58 : DEALLOCATE (dodisp, ghost, floating, atomnumber, radd2, c6d2)
194 :
195 58 : CALL timestop(handle)
196 :
197 116 : END SUBROUTINE calculate_dispersion_d2_pairpot
198 :
199 : ! **************************************************************************************************
200 : !> \brief ...
201 : !> \param z ...
202 : !> \param c6 ...
203 : !> \param r ...
204 : !> \param found ...
205 : ! **************************************************************************************************
206 38 : SUBROUTINE dftd2_param(z, c6, r, found)
207 :
208 : INTEGER, INTENT(in) :: z
209 : REAL(KIND=dp), INTENT(inout) :: c6, r
210 : LOGICAL, INTENT(inout) :: found
211 :
212 : REAL(KIND=dp), DIMENSION(54), PARAMETER :: c6val = (/0.14_dp, 0.08_dp, 1.61_dp, 1.61_dp, &
213 : 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,&
214 : 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, &
215 : 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, &
216 : 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, &
217 : 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, &
218 : 38.44_dp, 31.74_dp, 31.50_dp, 29.99_dp/)
219 : REAL(KIND=dp), DIMENSION(54), PARAMETER :: rval = (/1.001_dp, 1.012_dp, 0.825_dp, 1.408_dp, &
220 : 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, &
221 : 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, &
222 : 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, &
223 : 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, &
224 : 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, &
225 : 1.804_dp, 1.881_dp, 1.892_dp, 1.892_dp, 1.881_dp/)
226 :
227 : !
228 : ! GRIMME DISPERSION PARAMETERS
229 : ! Stefan Grimme, Semiempirical GGA-Type Density Functional Constructed
230 : ! with a Long-Range Dispersion Correction, J. Comp. Chem. 27: 1787-1799 (2006)
231 : ! doi:10.1002/jcc.20495
232 : !
233 : ! Conversion factor [Jnm^6mol^-1] -> [a.u.] : 17.34527758021901
234 : ! Conversion factor [A] -> [a.u.] : 1.889726132885643
235 : !
236 : ! C6 values in [Jnm^6/mol]
237 : ! vdW radii [A]
238 :
239 38 : IF (z > 0 .AND. z <= 54) THEN
240 38 : found = .TRUE.
241 38 : c6 = c6val(z)*1000._dp*bohr**6/kjmol
242 38 : r = rval(z)*bohr
243 : ELSE
244 0 : found = .FALSE.
245 : END IF
246 :
247 38 : END SUBROUTINE dftd2_param
248 :
249 : ! **************************************************************************************************
250 :
251 : END MODULE qs_dispersion_d2
|
