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