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 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, xtb_eeq_lagrange
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 2104 : 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 2104 : 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 2104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chia, efr, gam
95 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
96 2104 : 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 2104 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 2104 : 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 2104 : CALL timeset(routineN, handle)
110 :
111 2104 : 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 2104 : dft_control=dft_control)
120 2104 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
121 :
122 2104 : xtb_control => dft_control%qs_control%xtb_control
123 :
124 2104 : totalcharge = dft_control%charge
125 :
126 2104 : IF (atprop%energy) THEN
127 0 : CALL atprop_array_init(atprop%atecoul, natom)
128 : END IF
129 :
130 : ! gamma[a,b]
131 12624 : ALLOCATE (gab(nkind, nkind), gam(nkind))
132 2104 : gab = 0.0_dp
133 2104 : gam = 0.0_dp
134 7012 : DO ikind = 1, nkind
135 4908 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
136 4908 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
137 4908 : IF (.NOT. defined) CYCLE
138 4908 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
139 4908 : gam(ikind) = gama
140 24228 : DO jkind = 1, nkind
141 12308 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
142 12308 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
143 12308 : IF (.NOT. defined) CYCLE
144 12308 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
145 : !
146 29524 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
147 : !
148 : END DO
149 : END DO
150 :
151 : ! Chi[a,a]
152 2104 : enshift_type = xtb_control%enshift_type
153 2104 : IF (enshift_type == 0) THEN
154 0 : enshift_type = 2
155 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
156 : END IF
157 2104 : sgamma = 8.0_dp ! see D4 for periodic systems paper
158 2104 : esg = 1.0_dp + EXP(sgamma)
159 6312 : ALLOCATE (chia(natom))
160 2104 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
161 12048 : DO iatom = 1, natom
162 9944 : ikind = kind_of(iatom)
163 9944 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
164 9944 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
165 : !
166 9944 : IF (enshift_type == 1) THEN
167 9944 : 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 21992 : chia(iatom) = xi - kappa*scn
175 : !
176 : END DO
177 :
178 2104 : ef_energy = 0.0_dp
179 2104 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
180 : dft_control%apply_efield_field) THEN
181 996 : ALLOCATE (efr(natom))
182 1660 : efr(1:natom) = 0.0_dp
183 332 : CALL eeq_efield_pot(qs_env, efr)
184 3432 : chia(1:natom) = chia(1:natom) + efr(1:natom)
185 : END IF
186 :
187 2104 : do_ewald = xtb_control%do_ewald
188 :
189 2104 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
190 2104 : IF (do_ewald) THEN
191 : CALL get_qs_env(qs_env=qs_env, &
192 990 : 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 990 : 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 1114 : totalcharge=totalcharge, iounit=iunit)
203 : END IF
204 :
205 2104 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
206 : dft_control%apply_efield_field) THEN
207 332 : CALL eeq_efield_energy(qs_env, charges, ef_energy)
208 1660 : eeq_energy = eeq_energy - SUM(charges*efr)
209 332 : DEALLOCATE (efr)
210 : END IF
211 :
212 2104 : DEALLOCATE (gab, gam, chia)
213 :
214 2104 : CALL timestop(handle)
215 :
216 4208 : END SUBROUTINE xtb_eeq_calculation
217 :
218 : ! **************************************************************************************************
219 : !> \brief ...
220 : !> \param qs_env ...
221 : !> \param charges ...
222 : !> \param dcharges ...
223 : !> \param qlagrange ...
224 : !> \param cnumbers ...
225 : !> \param dcnum ...
226 : !> \param eeq_sparam ...
227 : ! **************************************************************************************************
228 68 : SUBROUTINE xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
229 :
230 : TYPE(qs_environment_type), POINTER :: qs_env
231 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges
232 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: qlagrange
233 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cnumbers
234 : TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
235 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
236 :
237 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_forces'
238 :
239 : INTEGER :: atom_a, atom_b, atom_c, enshift_type, &
240 : handle, i, ia, iatom, ikind, iunit, &
241 : jatom, jkind, katom, kkind, natom, &
242 : nkind
243 68 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
244 : LOGICAL :: defined, do_ewald, use_virial
245 : REAL(KIND=dp) :: ala, alb, alpha, cn, ctot, dr, dr2, drk, &
246 : elag, esg, fe, gam2, gama, grc, kappa, &
247 : qlam, qq, qq1, qq2, rcut, scn, sgamma, &
248 : totalcharge, xi
249 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
250 68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: epforce, gab
251 : REAL(KIND=dp), DIMENSION(3) :: fdik, ri, rij, rik, rj
252 : REAL(KIND=dp), DIMENSION(3, 3) :: pvir
253 68 : REAL(KIND=dp), DIMENSION(:), POINTER :: chrgx, dchia, qlag
254 68 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
255 : TYPE(atprop_type), POINTER :: atprop
256 : TYPE(cell_type), POINTER :: cell
257 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
258 : TYPE(dft_control_type), POINTER :: dft_control
259 : TYPE(distribution_1d_type), POINTER :: local_particles
260 : TYPE(ewald_environment_type), POINTER :: ewald_env
261 : TYPE(ewald_pw_type), POINTER :: ewald_pw
262 : TYPE(mp_para_env_type), POINTER :: para_env
263 : TYPE(neighbor_list_iterator_p_type), &
264 68 : DIMENSION(:), POINTER :: nl_iterator
265 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
266 68 : POINTER :: sab_tbe
267 68 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
268 68 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
269 68 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
270 : TYPE(virial_type), POINTER :: virial
271 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
272 : TYPE(xtb_control_type), POINTER :: xtb_control
273 :
274 68 : CALL timeset(routineN, handle)
275 :
276 68 : iunit = cp_logger_get_default_unit_nr()
277 :
278 : CALL get_qs_env(qs_env, &
279 : qs_kind_set=qs_kind_set, &
280 : atomic_kind_set=atomic_kind_set, &
281 : particle_set=particle_set, &
282 : atprop=atprop, &
283 : force=force, &
284 : virial=virial, &
285 : cell=cell, &
286 68 : dft_control=dft_control)
287 68 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
288 68 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
289 :
290 68 : xtb_control => dft_control%qs_control%xtb_control
291 :
292 68 : totalcharge = dft_control%charge
293 :
294 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
295 68 : atom_of_kind=atom_of_kind, kind_of=kind_of)
296 :
297 : ! gamma[a,b]
298 408 : ALLOCATE (gab(nkind, nkind), gam(nkind))
299 68 : gab = 0.0_dp
300 224 : DO ikind = 1, nkind
301 156 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
302 156 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
303 156 : IF (.NOT. defined) CYCLE
304 156 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
305 156 : gam(ikind) = gama
306 780 : DO jkind = 1, nkind
307 400 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
308 400 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
309 400 : IF (.NOT. defined) CYCLE
310 400 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
311 : !
312 956 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
313 : !
314 : END DO
315 : END DO
316 :
317 204 : ALLOCATE (qlag(natom))
318 :
319 68 : do_ewald = xtb_control%do_ewald
320 :
321 68 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
322 68 : IF (do_ewald) THEN
323 : CALL get_qs_env(qs_env=qs_env, &
324 46 : ewald_env=ewald_env, ewald_pw=ewald_pw)
325 : CALL eeq_solver(qlag, qlam, elag, &
326 : particle_set, kind_of, cell, -dcharges, gam, gab, &
327 : para_env, blacs_env, dft_control, eeq_sparam, &
328 314 : ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
329 : ELSE
330 : CALL eeq_solver(qlag, qlam, elag, &
331 : particle_set, kind_of, cell, -dcharges, gam, gab, &
332 102 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
333 : END IF
334 :
335 68 : enshift_type = xtb_control%enshift_type
336 68 : IF (enshift_type == 0) THEN
337 0 : enshift_type = 2
338 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
339 : END IF
340 68 : sgamma = 8.0_dp ! see D4 for periodic systems paper
341 68 : esg = 1.0_dp + EXP(sgamma)
342 272 : ALLOCATE (chrgx(natom), dchia(natom))
343 416 : DO iatom = 1, natom
344 348 : ikind = kind_of(iatom)
345 348 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
346 348 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
347 : !
348 348 : ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
349 416 : IF (enshift_type == 1) THEN
350 348 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
351 348 : dchia(iatom) = -ctot*kappa/scn
352 0 : ELSE IF (enshift_type == 2) THEN
353 0 : cn = cnumbers(iatom)
354 0 : scn = 1.0_dp/(esg - cn)
355 0 : dchia(iatom) = -ctot*kappa*scn
356 : ELSE
357 0 : CPABORT("Unknown enshift_type")
358 : END IF
359 : END DO
360 :
361 : ! Efield
362 68 : IF (dft_control%apply_period_efield) THEN
363 8 : CALL eeq_efield_force_periodic(qs_env, charges, qlag)
364 60 : ELSE IF (dft_control%apply_efield) THEN
365 8 : CALL eeq_efield_force_loc(qs_env, charges, qlag)
366 52 : ELSE IF (dft_control%apply_efield_field) THEN
367 0 : CPABORT("apply field")
368 : END IF
369 :
370 : ! Forces from q*X
371 : CALL get_qs_env(qs_env=qs_env, &
372 68 : local_particles=local_particles)
373 224 : DO ikind = 1, nkind
374 398 : DO ia = 1, local_particles%n_el(ikind)
375 174 : iatom = local_particles%list(ikind)%array(ia)
376 174 : atom_a = atom_of_kind(iatom)
377 1580 : DO i = 1, dcnum(iatom)%neighbors
378 1250 : katom = dcnum(iatom)%nlist(i)
379 1250 : kkind = kind_of(katom)
380 1250 : atom_c = atom_of_kind(katom)
381 5000 : rik = dcnum(iatom)%rik(:, i)
382 5000 : drk = SQRT(SUM(rik(:)**2))
383 1424 : IF (drk > 1.e-3_dp) THEN
384 5000 : fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
385 5000 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
386 5000 : force(kkind)%rho_elec(:, atom_c) = force(kkind)%rho_elec(:, atom_c) + fdik(:)
387 1250 : IF (use_virial) THEN
388 1084 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
389 : END IF
390 : END IF
391 : END DO
392 : END DO
393 : END DO
394 :
395 : ! Forces from (0.5*q+l)*dA/dR*q
396 68 : IF (do_ewald) THEN
397 46 : CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
398 46 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
399 46 : rcut = 1.0_dp*rcut
400 46 : CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
401 54750 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
402 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
403 54704 : iatom=iatom, jatom=jatom, r=rij)
404 54704 : atom_a = atom_of_kind(iatom)
405 54704 : atom_b = atom_of_kind(jatom)
406 : !
407 218816 : dr2 = SUM(rij**2)
408 54704 : dr = SQRT(dr2)
409 54704 : IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
410 6533 : fe = 1.0_dp
411 6533 : IF (iatom == jatom) fe = 0.5_dp
412 6533 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
413 6533 : gama = gab(ikind, jkind)
414 6533 : gam2 = gama*gama
415 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
416 6533 : - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
417 6533 : qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
418 6533 : qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
419 26132 : fdik(:) = -qq1*grc*rij(:)/dr
420 26132 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
421 26132 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
422 6533 : IF (use_virial) THEN
423 5586 : CALL virial_pair_force(virial%pv_virial, fe, fdik, rij)
424 : END IF
425 26132 : fdik(:) = qq2*grc*rij(:)/dr
426 26132 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
427 26132 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
428 6579 : IF (use_virial) THEN
429 5586 : CALL virial_pair_force(virial%pv_virial, -fe, fdik, rij)
430 : END IF
431 : END DO
432 46 : CALL neighbor_list_iterator_release(nl_iterator)
433 : ELSE
434 80 : DO ikind = 1, nkind
435 120 : DO ia = 1, local_particles%n_el(ikind)
436 40 : iatom = local_particles%list(ikind)%array(ia)
437 40 : atom_a = atom_of_kind(iatom)
438 160 : ri(1:3) = particle_set(iatom)%r(1:3)
439 246 : DO jatom = 1, natom
440 148 : IF (iatom == jatom) CYCLE
441 108 : jkind = kind_of(jatom)
442 108 : atom_b = atom_of_kind(jatom)
443 108 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
444 432 : rj(1:3) = particle_set(jatom)%r(1:3)
445 432 : rij(1:3) = ri(1:3) - rj(1:3)
446 432 : rij = pbc(rij, cell)
447 432 : dr2 = SUM(rij**2)
448 108 : dr = SQRT(dr2)
449 108 : gama = gab(ikind, jkind)
450 108 : gam2 = gama*gama
451 108 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
452 432 : fdik(:) = qq*grc*rij(:)/dr
453 432 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
454 432 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
455 148 : IF (use_virial) THEN
456 12 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rij)
457 : END IF
458 : END DO
459 : END DO
460 : END DO
461 : END IF
462 :
463 : ! Forces from Ewald potential: (q+l)*A*q
464 68 : IF (do_ewald) THEN
465 138 : ALLOCATE (epforce(3, natom))
466 46 : epforce = 0.0_dp
467 582 : dchia = -charges + qlag
468 314 : chrgx = charges
469 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
470 46 : particle_set, dchia, epforce)
471 314 : dchia = charges
472 582 : chrgx = qlag
473 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
474 46 : particle_set, dchia, epforce)
475 314 : DO iatom = 1, natom
476 268 : ikind = kind_of(iatom)
477 268 : i = atom_of_kind(iatom)
478 1118 : force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + epforce(:, iatom)
479 : END DO
480 46 : DEALLOCATE (epforce)
481 :
482 : ! virial
483 46 : IF (use_virial) THEN
484 372 : chrgx = charges - qlag
485 24 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
486 312 : virial%pv_virial = virial%pv_virial + pvir
487 372 : chrgx = qlag
488 24 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
489 312 : virial%pv_virial = virial%pv_virial - pvir
490 : END IF
491 : END IF
492 :
493 : ! return Lagrange multipliers
494 416 : qlagrange(1:natom) = qlag(1:natom)
495 :
496 68 : DEALLOCATE (gab, chrgx, dchia, qlag)
497 :
498 68 : CALL timestop(handle)
499 :
500 136 : END SUBROUTINE xtb_eeq_forces
501 :
502 : ! **************************************************************************************************
503 : !> \brief ...
504 : !> \param qs_env ...
505 : !> \param dcharges ...
506 : !> \param qlagrange ...
507 : !> \param eeq_sparam ...
508 : ! **************************************************************************************************
509 58 : SUBROUTINE xtb_eeq_lagrange(qs_env, dcharges, qlagrange, eeq_sparam)
510 :
511 : TYPE(qs_environment_type), POINTER :: qs_env
512 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dcharges
513 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: qlagrange
514 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
515 :
516 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_lagrange'
517 :
518 : INTEGER :: handle, ikind, iunit, jkind, natom, nkind
519 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
520 : LOGICAL :: defined, do_ewald
521 : REAL(KIND=dp) :: ala, alb, elag, gama, qlam, totalcharge
522 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
523 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
524 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: qlag
525 58 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
526 : TYPE(cell_type), POINTER :: cell
527 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
528 : TYPE(dft_control_type), POINTER :: dft_control
529 : TYPE(ewald_environment_type), POINTER :: ewald_env
530 : TYPE(ewald_pw_type), POINTER :: ewald_pw
531 : TYPE(mp_para_env_type), POINTER :: para_env
532 58 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
533 58 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
534 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
535 : TYPE(xtb_control_type), POINTER :: xtb_control
536 :
537 58 : CALL timeset(routineN, handle)
538 :
539 58 : iunit = cp_logger_get_default_unit_nr()
540 :
541 : CALL get_qs_env(qs_env, &
542 : qs_kind_set=qs_kind_set, &
543 : atomic_kind_set=atomic_kind_set, &
544 : particle_set=particle_set, &
545 : cell=cell, &
546 58 : dft_control=dft_control)
547 58 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
548 :
549 58 : xtb_control => dft_control%qs_control%xtb_control
550 :
551 58 : totalcharge = dft_control%charge
552 :
553 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
554 58 : atom_of_kind=atom_of_kind, kind_of=kind_of)
555 :
556 : ! gamma[a,b]
557 348 : ALLOCATE (gab(nkind, nkind), gam(nkind))
558 58 : gab = 0.0_dp
559 232 : DO ikind = 1, nkind
560 174 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
561 174 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
562 174 : IF (.NOT. defined) CYCLE
563 174 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
564 174 : gam(ikind) = gama
565 928 : DO jkind = 1, nkind
566 522 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
567 522 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
568 522 : IF (.NOT. defined) CYCLE
569 522 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
570 : !
571 1218 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
572 : !
573 : END DO
574 : END DO
575 :
576 174 : ALLOCATE (qlag(natom))
577 :
578 58 : do_ewald = xtb_control%do_ewald
579 :
580 58 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
581 58 : IF (do_ewald) THEN
582 : CALL get_qs_env(qs_env=qs_env, &
583 28 : ewald_env=ewald_env, ewald_pw=ewald_pw)
584 : CALL eeq_solver(qlag, qlam, elag, &
585 : particle_set, kind_of, cell, -dcharges, gam, gab, &
586 : para_env, blacs_env, dft_control, eeq_sparam, &
587 140 : ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
588 : ELSE
589 : CALL eeq_solver(qlag, qlam, elag, &
590 : particle_set, kind_of, cell, -dcharges, gam, gab, &
591 150 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
592 : END IF
593 :
594 : ! return Lagrange multipliers
595 290 : qlagrange(1:natom) = qlag(1:natom)
596 :
597 58 : DEALLOCATE (gab, qlag)
598 :
599 58 : CALL timestop(handle)
600 :
601 116 : END SUBROUTINE xtb_eeq_lagrange
602 :
603 : END MODULE xtb_eeq
|
