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 dispersion in DFTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dftb_dispersion
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE atprop_types, ONLY: atprop_array_init,&
16 : atprop_type
17 : USE cp_control_types, ONLY: dft_control_type,&
18 : dftb_control_type
19 : USE input_constants, ONLY: dispersion_d2,&
20 : dispersion_d3,&
21 : dispersion_d3bj,&
22 : dispersion_uff
23 : USE kinds, ONLY: dp
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE particle_types, ONLY: particle_type
26 : USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
27 : qs_dftb_pairpot_type
28 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
29 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
30 : USE qs_dispersion_types, ONLY: qs_dispersion_type
31 : USE qs_energy_types, ONLY: qs_energy_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_force_types, ONLY: qs_force_type
35 : USE qs_kind_types, ONLY: get_qs_kind,&
36 : qs_kind_type
37 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
38 : neighbor_list_iterate,&
39 : neighbor_list_iterator_create,&
40 : neighbor_list_iterator_p_type,&
41 : neighbor_list_iterator_release,&
42 : neighbor_list_set_p_type
43 : USE virial_methods, ONLY: virial_pair_force
44 : USE virial_types, ONLY: virial_type
45 : #include "./base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_dispersion'
52 :
53 : PUBLIC :: calculate_dftb_dispersion
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param qs_env ...
60 : !> \param para_env ...
61 : !> \param calculate_forces ...
62 : ! **************************************************************************************************
63 2822 : SUBROUTINE calculate_dftb_dispersion(qs_env, para_env, calculate_forces)
64 :
65 : TYPE(qs_environment_type), POINTER :: qs_env
66 : TYPE(mp_para_env_type), POINTER :: para_env
67 : LOGICAL, INTENT(IN) :: calculate_forces
68 :
69 : TYPE(dft_control_type), POINTER :: dft_control
70 : TYPE(dftb_control_type), POINTER :: dftb_control
71 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
72 : TYPE(qs_energy_type), POINTER :: energy
73 :
74 : CALL get_qs_env(qs_env=qs_env, &
75 : energy=energy, &
76 2822 : dft_control=dft_control)
77 :
78 2822 : energy%dispersion = 0._dp
79 :
80 2822 : dftb_control => dft_control%qs_control%dftb_control
81 2822 : IF (dftb_control%dispersion) THEN
82 2926 : SELECT CASE (dftb_control%dispersion_type)
83 : CASE (dispersion_uff)
84 1376 : CALL calculate_dispersion_uff(qs_env, para_env, calculate_forces)
85 : CASE (dispersion_d3, dispersion_d3bj, dispersion_d2)
86 174 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
87 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
88 174 : energy%dispersion, calculate_forces)
89 : CASE DEFAULT
90 1550 : CPABORT("")
91 : END SELECT
92 : END IF
93 :
94 2822 : END SUBROUTINE calculate_dftb_dispersion
95 :
96 : ! **************************************************************************************************
97 : !> \brief ...
98 : !> \param qs_env ...
99 : !> \param para_env ...
100 : !> \param calculate_forces ...
101 : ! **************************************************************************************************
102 1376 : SUBROUTINE calculate_dispersion_uff(qs_env, para_env, calculate_forces)
103 :
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 : TYPE(mp_para_env_type), POINTER :: para_env
106 : LOGICAL, INTENT(IN) :: calculate_forces
107 :
108 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_uff'
109 :
110 : INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
111 : jatom, jkind, nkind
112 1376 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
113 : LOGICAL :: use_virial
114 1376 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: define_kind
115 1376 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: rc_kind
116 : REAL(KIND=dp) :: a, b, c, devdw, dij, dr, eij, evdw, fac, &
117 : rc, x0ij, xij, xp
118 : REAL(KIND=dp), DIMENSION(3) :: fdij, rij
119 1376 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120 : TYPE(atprop_type), POINTER :: atprop
121 : TYPE(dft_control_type), POINTER :: dft_control
122 : TYPE(dftb_control_type), POINTER :: dftb_control
123 : TYPE(neighbor_list_iterator_p_type), &
124 1376 : DIMENSION(:), POINTER :: nl_iterator
125 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
126 1376 : POINTER :: sab_vdw
127 1376 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
128 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a
129 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
130 1376 : POINTER :: dftb_potential
131 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij
132 : TYPE(qs_energy_type), POINTER :: energy
133 1376 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
134 1376 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135 : TYPE(virial_type), POINTER :: virial
136 :
137 1376 : CALL timeset(routineN, handle)
138 :
139 1376 : NULLIFY (atomic_kind_set, sab_vdw, atprop)
140 :
141 : CALL get_qs_env(qs_env=qs_env, &
142 : energy=energy, &
143 : atomic_kind_set=atomic_kind_set, &
144 : qs_kind_set=qs_kind_set, &
145 : virial=virial, atprop=atprop, &
146 1376 : dft_control=dft_control)
147 :
148 1376 : energy%dispersion = 0._dp
149 :
150 1376 : dftb_control => dft_control%qs_control%dftb_control
151 :
152 1376 : IF (dftb_control%dispersion) THEN
153 :
154 1376 : NULLIFY (dftb_potential)
155 1376 : CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
156 1376 : IF (calculate_forces) THEN
157 450 : NULLIFY (force, particle_set)
158 : CALL get_qs_env(qs_env=qs_env, &
159 : particle_set=particle_set, &
160 450 : force=force)
161 450 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
162 450 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
163 : ELSE
164 : use_virial = .FALSE.
165 : END IF
166 1376 : nkind = SIZE(atomic_kind_set)
167 6880 : ALLOCATE (define_kind(nkind), rc_kind(nkind))
168 4156 : DO ikind = 1, nkind
169 2780 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
170 4156 : CALL get_dftb_atom_param(dftb_kind_a, defined=define_kind(ikind), rcdisp=rc_kind(ikind))
171 : END DO
172 :
173 1376 : evdw = 0._dp
174 :
175 1376 : IF (atprop%energy) THEN
176 58 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
177 58 : CALL atprop_array_init(atprop%atevdw, natom=SIZE(particle_set))
178 : END IF
179 :
180 1376 : CALL get_qs_env(qs_env=qs_env, sab_vdw=sab_vdw)
181 1376 : CALL neighbor_list_iterator_create(nl_iterator, sab_vdw)
182 8663828 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
183 8662452 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
184 8662452 : IF ((.NOT. define_kind(ikind)) .OR. (.NOT. define_kind(jkind))) CYCLE
185 8662452 : rc = rc_kind(ikind) + rc_kind(jkind)
186 : ! vdW potential
187 34649808 : dr = SQRT(SUM(rij(:)**2))
188 8663828 : IF (dr <= rc .AND. dr > 0.001_dp) THEN
189 8655844 : fac = 1._dp
190 8655844 : IF (iatom == jatom) fac = 0.5_dp
191 : ! retrieve information on potential
192 8655844 : dftb_param_ij => dftb_potential(ikind, jkind)
193 : ! vdW parameter
194 8655844 : xij = dftb_param_ij%xij
195 8655844 : dij = dftb_param_ij%dij
196 8655844 : x0ij = dftb_param_ij%x0ij
197 8655844 : a = dftb_param_ij%a
198 8655844 : b = dftb_param_ij%b
199 8655844 : c = dftb_param_ij%c
200 8655844 : IF (dr > x0ij) THEN
201 : ! This is the standard London contribution.
202 : ! UFF1 - Eq. 20 (long-range)
203 8640421 : xp = xij/dr
204 8640421 : eij = dij*(-2._dp*xp**6 + xp**12)*fac
205 8640421 : evdw = evdw + eij
206 8640421 : IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
207 2058704 : devdw = dij*12._dp*(xp**6 - xp**12)/dr*fac
208 2058704 : atom_a = atom_of_kind(iatom)
209 2058704 : atom_b = atom_of_kind(jatom)
210 8234816 : fdij(:) = devdw*rij(:)/dr
211 : force(ikind)%dispersion(:, atom_a) = &
212 8234816 : force(ikind)%dispersion(:, atom_a) - fdij(:)
213 : force(jkind)%dispersion(:, atom_b) = &
214 8234816 : force(jkind)%dispersion(:, atom_b) + fdij(:)
215 : END IF
216 : ELSE
217 : ! Shorter distance.
218 : ! London contribution should converge to a finite value.
219 : ! Using a parabola of the form (y = A - Bx**5 -Cx**10).
220 : ! Analytic parameters by forcing energy, first and second
221 : ! derivatives to be continuous.
222 15423 : eij = (A - B*dr**5 - C*dr**10)*fac
223 15423 : evdw = evdw + eij
224 15423 : IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
225 6139 : atom_a = atom_of_kind(iatom)
226 6139 : atom_b = atom_of_kind(jatom)
227 6139 : devdw = (-5*B*dr**4 - 10*C*dr**9)*fac
228 24556 : fdij(:) = devdw*rij(:)/dr
229 : force(ikind)%dispersion(:, atom_a) = &
230 24556 : force(ikind)%dispersion(:, atom_a) - fdij(:)
231 : force(jkind)%dispersion(:, atom_b) = &
232 24556 : force(jkind)%dispersion(:, atom_b) + fdij(:)
233 : END IF
234 : END IF
235 8655844 : IF (atprop%energy) THEN
236 1368540 : atprop%atevdw(iatom) = atprop%atevdw(iatom) + 0.5_dp*eij
237 1368540 : atprop%atevdw(jatom) = atprop%atevdw(jatom) + 0.5_dp*eij
238 : END IF
239 8655844 : IF (calculate_forces .AND. (dr > 0.001_dp) .AND. use_virial) THEN
240 1094831 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
241 1094831 : IF (atprop%stress) THEN
242 684224 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdij, rij)
243 684224 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fdij, rij)
244 : END IF
245 : END IF
246 : END IF
247 : END DO
248 1376 : CALL neighbor_list_iterator_release(nl_iterator)
249 :
250 : ! set dispersion energy
251 1376 : CALL para_env%sum(evdw)
252 1376 : energy%dispersion = evdw
253 :
254 : END IF
255 :
256 1376 : CALL timestop(handle)
257 :
258 2752 : END SUBROUTINE calculate_dispersion_uff
259 :
260 : END MODULE qs_dftb_dispersion
261 :
|
