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 charge equilibration in xTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE xtb_eeq
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 cell_types, ONLY: cell_type,&
18 : pbc
19 : USE cp_blacs_env, ONLY: cp_blacs_env_type
20 : USE cp_control_types, ONLY: dft_control_type,&
21 : xtb_control_type
22 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
23 : USE distribution_1d_types, ONLY: distribution_1d_type
24 : USE eeq_input, ONLY: eeq_solver_type
25 : USE eeq_method, ONLY: eeq_efield_energy,&
26 : eeq_efield_force_loc,&
27 : eeq_efield_force_periodic,&
28 : eeq_efield_pot,&
29 : eeq_solver
30 : USE ewald_environment_types, ONLY: ewald_env_get,&
31 : ewald_environment_type
32 : USE ewald_pw_types, ONLY: ewald_pw_type
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: oorootpi
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE particle_types, ONLY: particle_type
37 : USE qs_dispersion_cnum, ONLY: dcnum_type
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_force_types, ONLY: qs_force_type
41 : USE qs_kind_types, ONLY: get_qs_kind,&
42 : qs_kind_type
43 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
44 : neighbor_list_iterate,&
45 : neighbor_list_iterator_create,&
46 : neighbor_list_iterator_p_type,&
47 : neighbor_list_iterator_release,&
48 : neighbor_list_set_p_type
49 : USE spme, ONLY: spme_forces,&
50 : spme_virial
51 : USE virial_methods, ONLY: virial_pair_force
52 : USE virial_types, ONLY: virial_type
53 : USE xtb_types, ONLY: get_xtb_atom_param,&
54 : xtb_atom_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_eeq'
62 :
63 : PUBLIC :: xtb_eeq_calculation, xtb_eeq_forces
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param qs_env ...
70 : !> \param charges ...
71 : !> \param cnumbers ...
72 : !> \param eeq_sparam ...
73 : !> \param eeq_energy ...
74 : !> \param ef_energy ...
75 : !> \param lambda ...
76 : ! **************************************************************************************************
77 1436 : SUBROUTINE xtb_eeq_calculation(qs_env, charges, cnumbers, &
78 : eeq_sparam, eeq_energy, ef_energy, lambda)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
82 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cnumbers
83 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
84 : REAL(KIND=dp), INTENT(INOUT) :: eeq_energy, ef_energy, lambda
85 :
86 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_calculation'
87 :
88 : INTEGER :: enshift_type, handle, iatom, ikind, &
89 : iunit, jkind, natom, nkind
90 1436 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
91 : LOGICAL :: defined, do_ewald
92 : REAL(KIND=dp) :: ala, alb, cn, esg, gama, kappa, scn, &
93 : sgamma, totalcharge, xi
94 1436 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chia, efr, gam
95 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
96 1436 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
97 : TYPE(atprop_type), POINTER :: atprop
98 : TYPE(cell_type), POINTER :: cell
99 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
100 : TYPE(dft_control_type), POINTER :: dft_control
101 : TYPE(ewald_environment_type), POINTER :: ewald_env
102 : TYPE(ewald_pw_type), POINTER :: ewald_pw
103 : TYPE(mp_para_env_type), POINTER :: para_env
104 1436 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 1436 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
106 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
107 : TYPE(xtb_control_type), POINTER :: xtb_control
108 :
109 1436 : CALL timeset(routineN, handle)
110 :
111 1436 : iunit = cp_logger_get_default_unit_nr()
112 :
113 : CALL get_qs_env(qs_env, &
114 : qs_kind_set=qs_kind_set, &
115 : atomic_kind_set=atomic_kind_set, &
116 : particle_set=particle_set, &
117 : cell=cell, &
118 : atprop=atprop, &
119 1436 : dft_control=dft_control)
120 1436 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
121 :
122 1436 : xtb_control => dft_control%qs_control%xtb_control
123 :
124 1436 : totalcharge = dft_control%charge
125 :
126 1436 : IF (atprop%energy) THEN
127 0 : CALL atprop_array_init(atprop%atecoul, natom)
128 : END IF
129 :
130 : ! gamma[a,b]
131 8616 : ALLOCATE (gab(nkind, nkind), gam(nkind))
132 13632 : gab = 0.0_dp
133 4890 : gam = 0.0_dp
134 4890 : DO ikind = 1, nkind
135 3454 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
136 3454 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
137 3454 : IF (.NOT. defined) CYCLE
138 3454 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
139 3454 : gam(ikind) = gama
140 17086 : DO jkind = 1, nkind
141 8742 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
142 8742 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
143 8742 : IF (.NOT. defined) CYCLE
144 8742 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
145 : !
146 20938 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
147 : !
148 : END DO
149 : END DO
150 :
151 : ! Chi[a,a]
152 1436 : enshift_type = xtb_control%enshift_type
153 1436 : IF (enshift_type == 0) THEN
154 0 : enshift_type = 2
155 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
156 : END IF
157 1436 : sgamma = 8.0_dp ! see D4 for periodic systems paper
158 1436 : esg = 1.0_dp + EXP(sgamma)
159 4308 : ALLOCATE (chia(natom))
160 1436 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
161 8346 : DO iatom = 1, natom
162 6910 : ikind = kind_of(iatom)
163 6910 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
164 6910 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
165 : !
166 6910 : IF (enshift_type == 1) THEN
167 6910 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
168 0 : ELSE IF (enshift_type == 2) THEN
169 0 : cn = cnumbers(iatom)/esg
170 0 : scn = LOG(esg/(esg - cnumbers(iatom)))
171 : ELSE
172 0 : CPABORT("Unknown enshift_type")
173 : END IF
174 15256 : chia(iatom) = xi - kappa*scn
175 : !
176 : END DO
177 :
178 1436 : ef_energy = 0.0_dp
179 1436 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
180 : dft_control%apply_efield_field) THEN
181 312 : ALLOCATE (efr(natom))
182 520 : efr(1:natom) = 0.0_dp
183 104 : CALL eeq_efield_pot(qs_env, efr)
184 1852 : chia(1:natom) = chia(1:natom) + efr(1:natom)
185 : END IF
186 :
187 1436 : do_ewald = xtb_control%do_ewald
188 :
189 1436 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
190 1436 : IF (do_ewald) THEN
191 : CALL get_qs_env(qs_env=qs_env, &
192 550 : ewald_env=ewald_env, ewald_pw=ewald_pw)
193 : CALL eeq_solver(charges, lambda, eeq_energy, &
194 : particle_set, kind_of, cell, chia, gam, gab, &
195 : para_env, blacs_env, dft_control, eeq_sparam, &
196 : totalcharge=totalcharge, ewald=do_ewald, &
197 550 : ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
198 : ELSE
199 : CALL eeq_solver(charges, lambda, eeq_energy, &
200 : particle_set, kind_of, cell, chia, gam, gab, &
201 : para_env, blacs_env, dft_control, eeq_sparam, &
202 886 : totalcharge=totalcharge, iounit=iunit)
203 : END IF
204 :
205 1436 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
206 : dft_control%apply_efield_field) THEN
207 104 : CALL eeq_efield_energy(qs_env, charges, ef_energy)
208 520 : eeq_energy = eeq_energy - SUM(charges*efr)
209 104 : DEALLOCATE (efr)
210 : END IF
211 :
212 1436 : DEALLOCATE (gab, gam, chia)
213 :
214 1436 : CALL timestop(handle)
215 :
216 2872 : END SUBROUTINE xtb_eeq_calculation
217 :
218 : ! **************************************************************************************************
219 : !> \brief ...
220 : !> \param qs_env ...
221 : !> \param charges ...
222 : !> \param dcharges ...
223 : !> \param cnumbers ...
224 : !> \param dcnum ...
225 : !> \param eeq_sparam ...
226 : ! **************************************************************************************************
227 28 : SUBROUTINE xtb_eeq_forces(qs_env, charges, dcharges, cnumbers, dcnum, eeq_sparam)
228 :
229 : TYPE(qs_environment_type), POINTER :: qs_env
230 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges, cnumbers
231 : TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
232 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
233 :
234 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_forces'
235 :
236 : INTEGER :: atom_a, atom_b, atom_c, enshift_type, &
237 : handle, i, ia, iatom, ikind, iunit, &
238 : jatom, jkind, katom, kkind, natom, &
239 : nkind
240 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
241 : LOGICAL :: defined, do_ewald, use_virial
242 : REAL(KIND=dp) :: ala, alb, alpha, cn, ctot, dr, dr2, drk, &
243 : elag, esg, gam2, gama, grc, kappa, &
244 : qlam, qq, qq1, qq2, rcut, scn, sgamma, &
245 : totalcharge, xi
246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
247 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: epforce, gab
248 : REAL(KIND=dp), DIMENSION(3) :: fdik, ri, rij, rik, rj
249 : REAL(KIND=dp), DIMENSION(3, 3) :: pvir
250 28 : REAL(KIND=dp), DIMENSION(:), POINTER :: chrgx, dchia, qlag
251 28 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
252 : TYPE(atprop_type), POINTER :: atprop
253 : TYPE(cell_type), POINTER :: cell
254 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
255 : TYPE(dft_control_type), POINTER :: dft_control
256 : TYPE(distribution_1d_type), POINTER :: local_particles
257 : TYPE(ewald_environment_type), POINTER :: ewald_env
258 : TYPE(ewald_pw_type), POINTER :: ewald_pw
259 : TYPE(mp_para_env_type), POINTER :: para_env
260 : TYPE(neighbor_list_iterator_p_type), &
261 28 : DIMENSION(:), POINTER :: nl_iterator
262 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
263 28 : POINTER :: sab_tbe
264 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
265 28 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
266 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
267 : TYPE(virial_type), POINTER :: virial
268 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
269 : TYPE(xtb_control_type), POINTER :: xtb_control
270 :
271 28 : CALL timeset(routineN, handle)
272 :
273 28 : iunit = cp_logger_get_default_unit_nr()
274 :
275 : CALL get_qs_env(qs_env, &
276 : qs_kind_set=qs_kind_set, &
277 : atomic_kind_set=atomic_kind_set, &
278 : particle_set=particle_set, &
279 : atprop=atprop, &
280 : force=force, &
281 : virial=virial, &
282 : cell=cell, &
283 28 : dft_control=dft_control)
284 28 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
285 28 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
286 :
287 28 : xtb_control => dft_control%qs_control%xtb_control
288 :
289 28 : totalcharge = dft_control%charge
290 :
291 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
292 28 : atom_of_kind=atom_of_kind, kind_of=kind_of)
293 :
294 : ! gamma[a,b]
295 168 : ALLOCATE (gab(nkind, nkind), gam(nkind))
296 316 : gab = 0.0_dp
297 104 : DO ikind = 1, nkind
298 76 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
299 76 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
300 76 : IF (.NOT. defined) CYCLE
301 76 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
302 76 : gam(ikind) = gama
303 392 : DO jkind = 1, nkind
304 212 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
305 212 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
306 212 : IF (.NOT. defined) CYCLE
307 212 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
308 : !
309 500 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
310 : !
311 : END DO
312 : END DO
313 :
314 84 : ALLOCATE (qlag(natom))
315 :
316 28 : do_ewald = xtb_control%do_ewald
317 :
318 28 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
319 28 : IF (do_ewald) THEN
320 : CALL get_qs_env(qs_env=qs_env, &
321 18 : ewald_env=ewald_env, ewald_pw=ewald_pw)
322 : CALL eeq_solver(qlag, qlam, elag, &
323 : particle_set, kind_of, cell, -dcharges, gam, gab, &
324 : para_env, blacs_env, dft_control, eeq_sparam, &
325 122 : ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
326 : ELSE
327 : CALL eeq_solver(qlag, qlam, elag, &
328 : particle_set, kind_of, cell, -dcharges, gam, gab, &
329 50 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
330 : END IF
331 :
332 28 : enshift_type = xtb_control%enshift_type
333 28 : IF (enshift_type == 0) THEN
334 0 : enshift_type = 2
335 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
336 : END IF
337 28 : sgamma = 8.0_dp ! see D4 for periodic systems paper
338 28 : esg = 1.0_dp + EXP(sgamma)
339 112 : ALLOCATE (chrgx(natom), dchia(natom))
340 172 : DO iatom = 1, natom
341 144 : ikind = kind_of(iatom)
342 144 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
343 144 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
344 : !
345 144 : ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
346 172 : IF (enshift_type == 1) THEN
347 144 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
348 144 : dchia(iatom) = -ctot*kappa/scn
349 0 : ELSE IF (enshift_type == 2) THEN
350 0 : cn = cnumbers(iatom)
351 0 : scn = 1.0_dp/(esg - cn)
352 0 : dchia(iatom) = -ctot*kappa*scn
353 : ELSE
354 0 : CPABORT("Unknown enshift_type")
355 : END IF
356 : END DO
357 :
358 : ! Efield
359 28 : IF (dft_control%apply_period_efield) THEN
360 2 : CALL eeq_efield_force_periodic(qs_env, charges, qlag)
361 26 : ELSE IF (dft_control%apply_efield) THEN
362 2 : CALL eeq_efield_force_loc(qs_env, charges, qlag)
363 24 : ELSE IF (dft_control%apply_efield_field) THEN
364 0 : CPABORT("apply field")
365 : END IF
366 :
367 : ! Forces from q*X
368 : CALL get_qs_env(qs_env=qs_env, &
369 28 : local_particles=local_particles)
370 104 : DO ikind = 1, nkind
371 176 : DO ia = 1, local_particles%n_el(ikind)
372 72 : iatom = local_particles%list(ikind)%array(ia)
373 72 : atom_a = atom_of_kind(iatom)
374 576 : DO i = 1, dcnum(iatom)%neighbors
375 428 : katom = dcnum(iatom)%nlist(i)
376 428 : kkind = kind_of(katom)
377 428 : atom_c = atom_of_kind(katom)
378 1712 : rik = dcnum(iatom)%rik(:, i)
379 1712 : drk = SQRT(SUM(rik(:)**2))
380 500 : IF (drk > 1.e-3_dp) THEN
381 1712 : fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
382 1712 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
383 1712 : force(kkind)%rho_elec(:, atom_c) = force(kkind)%rho_elec(:, atom_c) + fdik(:)
384 428 : IF (use_virial) THEN
385 368 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
386 : END IF
387 : END IF
388 : END DO
389 : END DO
390 : END DO
391 :
392 : ! Forces from (0.5*q+l)*dA/dR*q
393 28 : IF (do_ewald) THEN
394 18 : CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
395 18 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
396 18 : rcut = 2.0_dp*rcut
397 18 : CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
398 15879 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
399 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
400 15861 : iatom=iatom, jatom=jatom, r=rij)
401 15861 : atom_a = atom_of_kind(iatom)
402 15861 : atom_b = atom_of_kind(jatom)
403 : !
404 63444 : dr2 = SUM(rij**2)
405 15861 : dr = SQRT(dr2)
406 15861 : IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
407 15809 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
408 15809 : gama = gab(ikind, jkind)
409 15809 : gam2 = gama*gama
410 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
411 15809 : - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
412 15809 : qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
413 15809 : qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
414 63236 : fdik(:) = -qq1*grc*rij(:)/dr
415 63236 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
416 63236 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
417 15809 : IF (use_virial) THEN
418 11768 : CALL virial_pair_force(virial%pv_virial, 1._dp, fdik, rij)
419 : END IF
420 63236 : fdik(:) = qq2*grc*rij(:)/dr
421 63236 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
422 63236 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
423 15827 : IF (use_virial) THEN
424 11768 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rij)
425 : END IF
426 : END DO
427 18 : CALL neighbor_list_iterator_release(nl_iterator)
428 : ELSE
429 40 : DO ikind = 1, nkind
430 60 : DO ia = 1, local_particles%n_el(ikind)
431 20 : iatom = local_particles%list(ikind)%array(ia)
432 20 : atom_a = atom_of_kind(iatom)
433 80 : ri(1:3) = particle_set(iatom)%r(1:3)
434 130 : DO jatom = 1, natom
435 80 : IF (iatom == jatom) CYCLE
436 60 : jkind = kind_of(jatom)
437 60 : atom_b = atom_of_kind(jatom)
438 60 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
439 240 : rj(1:3) = particle_set(jatom)%r(1:3)
440 240 : rij(1:3) = ri(1:3) - rj(1:3)
441 240 : rij = pbc(rij, cell)
442 240 : dr2 = SUM(rij**2)
443 60 : dr = SQRT(dr2)
444 60 : gama = gab(ikind, jkind)
445 60 : gam2 = gama*gama
446 60 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
447 240 : fdik(:) = qq*grc*rij(:)/dr
448 240 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
449 260 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
450 : END DO
451 : END DO
452 : END DO
453 : END IF
454 :
455 : ! Forces from Ewald potential: (q+l)*A*q
456 28 : IF (do_ewald) THEN
457 54 : ALLOCATE (epforce(3, natom))
458 434 : epforce = 0.0_dp
459 226 : dchia = -charges + qlag
460 122 : chrgx = charges
461 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
462 18 : particle_set, dchia, epforce)
463 122 : dchia = charges
464 226 : chrgx = qlag
465 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
466 18 : particle_set, dchia, epforce)
467 122 : DO iatom = 1, natom
468 104 : ikind = kind_of(iatom)
469 104 : i = atom_of_kind(iatom)
470 434 : force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + epforce(:, iatom)
471 : END DO
472 18 : DEALLOCATE (epforce)
473 :
474 : ! virial
475 18 : IF (use_virial) THEN
476 136 : chrgx = charges - qlag
477 8 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
478 104 : virial%pv_virial = virial%pv_virial + pvir
479 136 : chrgx = qlag
480 8 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
481 104 : virial%pv_virial = virial%pv_virial - pvir
482 : END IF
483 : END IF
484 :
485 28 : DEALLOCATE (gab, chrgx, dchia, qlag)
486 :
487 28 : CALL timestop(handle)
488 :
489 56 : END SUBROUTINE xtb_eeq_forces
490 :
491 : END MODULE xtb_eeq
|
