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 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 1702 : 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 1702 : 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 1702 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chia, efr, gam
95 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
96 1702 : 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 1702 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 1702 : 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 1702 : CALL timeset(routineN, handle)
110 :
111 1702 : 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 1702 : dft_control=dft_control)
120 1702 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
121 :
122 1702 : xtb_control => dft_control%qs_control%xtb_control
123 :
124 1702 : totalcharge = dft_control%charge
125 :
126 1702 : IF (atprop%energy) THEN
127 0 : CALL atprop_array_init(atprop%atecoul, natom)
128 : END IF
129 :
130 : ! gamma[a,b]
131 10212 : ALLOCATE (gab(nkind, nkind), gam(nkind))
132 17090 : gab = 0.0_dp
133 5954 : gam = 0.0_dp
134 5954 : DO ikind = 1, nkind
135 4252 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
136 4252 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
137 4252 : IF (.NOT. defined) CYCLE
138 4252 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
139 4252 : gam(ikind) = gama
140 21342 : DO jkind = 1, nkind
141 11136 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
142 11136 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
143 11136 : IF (.NOT. defined) CYCLE
144 11136 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
145 : !
146 26524 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
147 : !
148 : END DO
149 : END DO
150 :
151 : ! Chi[a,a]
152 1702 : enshift_type = xtb_control%enshift_type
153 1702 : IF (enshift_type == 0) THEN
154 0 : enshift_type = 2
155 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
156 : END IF
157 1702 : sgamma = 8.0_dp ! see D4 for periodic systems paper
158 1702 : esg = 1.0_dp + EXP(sgamma)
159 5106 : ALLOCATE (chia(natom))
160 1702 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
161 9676 : DO iatom = 1, natom
162 7974 : ikind = kind_of(iatom)
163 7974 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
164 7974 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
165 : !
166 7974 : IF (enshift_type == 1) THEN
167 7974 : 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 17650 : chia(iatom) = xi - kappa*scn
175 : !
176 : END DO
177 :
178 1702 : ef_energy = 0.0_dp
179 1702 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
180 : dft_control%apply_efield_field) THEN
181 984 : ALLOCATE (efr(natom))
182 1640 : efr(1:natom) = 0.0_dp
183 328 : CALL eeq_efield_pot(qs_env, efr)
184 3014 : chia(1:natom) = chia(1:natom) + efr(1:natom)
185 : END IF
186 :
187 1702 : do_ewald = xtb_control%do_ewald
188 :
189 1702 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
190 1702 : IF (do_ewald) THEN
191 : CALL get_qs_env(qs_env=qs_env, &
192 756 : 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 756 : 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 946 : totalcharge=totalcharge, iounit=iunit)
203 : END IF
204 :
205 1702 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
206 : dft_control%apply_efield_field) THEN
207 328 : CALL eeq_efield_energy(qs_env, charges, ef_energy)
208 1640 : eeq_energy = eeq_energy - SUM(charges*efr)
209 328 : DEALLOCATE (efr)
210 : END IF
211 :
212 1702 : DEALLOCATE (gab, gam, chia)
213 :
214 1702 : CALL timestop(handle)
215 :
216 3404 : 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 42 : 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 42 : 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 42 : 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 42 : REAL(KIND=dp), DIMENSION(:), POINTER :: chrgx, dchia, qlag
254 42 : 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 42 : DIMENSION(:), POINTER :: nl_iterator
265 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
266 42 : POINTER :: sab_tbe
267 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
268 42 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
269 42 : 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 42 : CALL timeset(routineN, handle)
275 :
276 42 : 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 42 : dft_control=dft_control)
287 42 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
288 42 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
289 :
290 42 : xtb_control => dft_control%qs_control%xtb_control
291 :
292 42 : totalcharge = dft_control%charge
293 :
294 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
295 42 : atom_of_kind=atom_of_kind, kind_of=kind_of)
296 :
297 : ! gamma[a,b]
298 252 : ALLOCATE (gab(nkind, nkind), gam(nkind))
299 498 : gab = 0.0_dp
300 160 : DO ikind = 1, nkind
301 118 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
302 118 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
303 118 : IF (.NOT. defined) CYCLE
304 118 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
305 118 : gam(ikind) = gama
306 616 : DO jkind = 1, nkind
307 338 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
308 338 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
309 338 : IF (.NOT. defined) CYCLE
310 338 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
311 : !
312 794 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
313 : !
314 : END DO
315 : END DO
316 :
317 126 : ALLOCATE (qlag(natom))
318 :
319 42 : do_ewald = xtb_control%do_ewald
320 :
321 42 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
322 42 : IF (do_ewald) THEN
323 : CALL get_qs_env(qs_env=qs_env, &
324 28 : 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 172 : 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 70 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
333 : END IF
334 :
335 42 : enshift_type = xtb_control%enshift_type
336 42 : IF (enshift_type == 0) THEN
337 0 : enshift_type = 2
338 0 : IF (ALL(cell%perd == 0)) enshift_type = 1
339 : END IF
340 42 : sgamma = 8.0_dp ! see D4 for periodic systems paper
341 42 : esg = 1.0_dp + EXP(sgamma)
342 168 : ALLOCATE (chrgx(natom), dchia(natom))
343 242 : DO iatom = 1, natom
344 200 : ikind = kind_of(iatom)
345 200 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
346 200 : CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
347 : !
348 200 : ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
349 242 : IF (enshift_type == 1) THEN
350 200 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
351 200 : 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 42 : IF (dft_control%apply_period_efield) THEN
363 8 : CALL eeq_efield_force_periodic(qs_env, charges, qlag)
364 34 : ELSE IF (dft_control%apply_efield) THEN
365 8 : CALL eeq_efield_force_loc(qs_env, charges, qlag)
366 26 : 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 42 : local_particles=local_particles)
373 160 : DO ikind = 1, nkind
374 260 : DO ia = 1, local_particles%n_el(ikind)
375 100 : iatom = local_particles%list(ikind)%array(ia)
376 100 : atom_a = atom_of_kind(iatom)
377 688 : DO i = 1, dcnum(iatom)%neighbors
378 470 : katom = dcnum(iatom)%nlist(i)
379 470 : kkind = kind_of(katom)
380 470 : atom_c = atom_of_kind(katom)
381 1880 : rik = dcnum(iatom)%rik(:, i)
382 1880 : drk = SQRT(SUM(rik(:)**2))
383 570 : IF (drk > 1.e-3_dp) THEN
384 1880 : fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
385 1880 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
386 1880 : force(kkind)%rho_elec(:, atom_c) = force(kkind)%rho_elec(:, atom_c) + fdik(:)
387 470 : IF (use_virial) THEN
388 374 : 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 42 : IF (do_ewald) THEN
397 28 : CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
398 28 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
399 28 : rcut = 1.0_dp*rcut
400 28 : CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
401 19968 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
402 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
403 19940 : iatom=iatom, jatom=jatom, r=rij)
404 19940 : atom_a = atom_of_kind(iatom)
405 19940 : atom_b = atom_of_kind(jatom)
406 : !
407 79760 : dr2 = SUM(rij**2)
408 19940 : dr = SQRT(dr2)
409 19940 : IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
410 2309 : fe = 1.0_dp
411 2309 : IF (iatom == jatom) fe = 0.5_dp
412 2309 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
413 2309 : gama = gab(ikind, jkind)
414 2309 : gam2 = gama*gama
415 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
416 2309 : - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
417 2309 : qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
418 2309 : qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
419 9236 : fdik(:) = -qq1*grc*rij(:)/dr
420 9236 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
421 9236 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
422 2309 : IF (use_virial) THEN
423 1554 : CALL virial_pair_force(virial%pv_virial, fe, fdik, rij)
424 : END IF
425 9236 : fdik(:) = qq2*grc*rij(:)/dr
426 9236 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
427 9236 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
428 2337 : IF (use_virial) THEN
429 1554 : CALL virial_pair_force(virial%pv_virial, -fe, fdik, rij)
430 : END IF
431 : END DO
432 28 : CALL neighbor_list_iterator_release(nl_iterator)
433 : ELSE
434 56 : DO ikind = 1, nkind
435 84 : DO ia = 1, local_particles%n_el(ikind)
436 28 : iatom = local_particles%list(ikind)%array(ia)
437 28 : atom_a = atom_of_kind(iatom)
438 112 : ri(1:3) = particle_set(iatom)%r(1:3)
439 182 : DO jatom = 1, natom
440 112 : IF (iatom == jatom) CYCLE
441 84 : jkind = kind_of(jatom)
442 84 : atom_b = atom_of_kind(jatom)
443 84 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
444 336 : rj(1:3) = particle_set(jatom)%r(1:3)
445 336 : rij(1:3) = ri(1:3) - rj(1:3)
446 336 : rij = pbc(rij, cell)
447 336 : dr2 = SUM(rij**2)
448 84 : dr = SQRT(dr2)
449 84 : gama = gab(ikind, jkind)
450 84 : gam2 = gama*gama
451 84 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
452 336 : fdik(:) = qq*grc*rij(:)/dr
453 336 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
454 364 : force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
455 : END DO
456 : END DO
457 : END DO
458 : END IF
459 :
460 : ! Forces from Ewald potential: (q+l)*A*q
461 42 : IF (do_ewald) THEN
462 84 : ALLOCATE (epforce(3, natom))
463 604 : epforce = 0.0_dp
464 316 : dchia = -charges + qlag
465 172 : chrgx = charges
466 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
467 28 : particle_set, dchia, epforce)
468 172 : dchia = charges
469 316 : chrgx = qlag
470 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
471 28 : particle_set, dchia, epforce)
472 172 : DO iatom = 1, natom
473 144 : ikind = kind_of(iatom)
474 144 : i = atom_of_kind(iatom)
475 604 : force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + epforce(:, iatom)
476 : END DO
477 28 : DEALLOCATE (epforce)
478 :
479 : ! virial
480 28 : IF (use_virial) THEN
481 154 : chrgx = charges - qlag
482 10 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
483 130 : virial%pv_virial = virial%pv_virial + pvir
484 154 : chrgx = qlag
485 10 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
486 130 : virial%pv_virial = virial%pv_virial - pvir
487 : END IF
488 : END IF
489 :
490 : ! return Lagrange multipliers
491 242 : qlagrange(1:natom) = qlag(1:natom)
492 :
493 42 : DEALLOCATE (gab, chrgx, dchia, qlag)
494 :
495 42 : CALL timestop(handle)
496 :
497 84 : END SUBROUTINE xtb_eeq_forces
498 :
499 : ! **************************************************************************************************
500 : !> \brief ...
501 : !> \param qs_env ...
502 : !> \param dcharges ...
503 : !> \param qlagrange ...
504 : !> \param eeq_sparam ...
505 : ! **************************************************************************************************
506 56 : SUBROUTINE xtb_eeq_lagrange(qs_env, dcharges, qlagrange, eeq_sparam)
507 :
508 : TYPE(qs_environment_type), POINTER :: qs_env
509 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dcharges
510 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: qlagrange
511 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
512 :
513 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_lagrange'
514 :
515 : INTEGER :: handle, ikind, iunit, jkind, natom, nkind
516 56 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
517 : LOGICAL :: defined, do_ewald
518 : REAL(KIND=dp) :: ala, alb, elag, gama, qlam, totalcharge
519 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
520 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
521 56 : REAL(KIND=dp), DIMENSION(:), POINTER :: qlag
522 56 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
523 : TYPE(cell_type), POINTER :: cell
524 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
525 : TYPE(dft_control_type), POINTER :: dft_control
526 : TYPE(ewald_environment_type), POINTER :: ewald_env
527 : TYPE(ewald_pw_type), POINTER :: ewald_pw
528 : TYPE(mp_para_env_type), POINTER :: para_env
529 56 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
530 56 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
531 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
532 : TYPE(xtb_control_type), POINTER :: xtb_control
533 :
534 56 : CALL timeset(routineN, handle)
535 :
536 56 : iunit = cp_logger_get_default_unit_nr()
537 :
538 : CALL get_qs_env(qs_env, &
539 : qs_kind_set=qs_kind_set, &
540 : atomic_kind_set=atomic_kind_set, &
541 : particle_set=particle_set, &
542 : cell=cell, &
543 56 : dft_control=dft_control)
544 56 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
545 :
546 56 : xtb_control => dft_control%qs_control%xtb_control
547 :
548 56 : totalcharge = dft_control%charge
549 :
550 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
551 56 : atom_of_kind=atom_of_kind, kind_of=kind_of)
552 :
553 : ! gamma[a,b]
554 336 : ALLOCATE (gab(nkind, nkind), gam(nkind))
555 728 : gab = 0.0_dp
556 224 : DO ikind = 1, nkind
557 168 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
558 168 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
559 168 : IF (.NOT. defined) CYCLE
560 168 : CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
561 168 : gam(ikind) = gama
562 896 : DO jkind = 1, nkind
563 504 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
564 504 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
565 504 : IF (.NOT. defined) CYCLE
566 504 : CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
567 : !
568 1176 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
569 : !
570 : END DO
571 : END DO
572 :
573 168 : ALLOCATE (qlag(natom))
574 :
575 56 : do_ewald = xtb_control%do_ewald
576 :
577 56 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
578 56 : IF (do_ewald) THEN
579 : CALL get_qs_env(qs_env=qs_env, &
580 28 : ewald_env=ewald_env, ewald_pw=ewald_pw)
581 : CALL eeq_solver(qlag, qlam, elag, &
582 : particle_set, kind_of, cell, -dcharges, gam, gab, &
583 : para_env, blacs_env, dft_control, eeq_sparam, &
584 140 : ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
585 : ELSE
586 : CALL eeq_solver(qlag, qlam, elag, &
587 : particle_set, kind_of, cell, -dcharges, gam, gab, &
588 140 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
589 : END IF
590 :
591 : ! return Lagrange multipliers
592 280 : qlagrange(1:natom) = qlag(1:natom)
593 :
594 56 : DEALLOCATE (gab, qlag)
595 :
596 56 : CALL timestop(handle)
597 :
598 112 : END SUBROUTINE xtb_eeq_lagrange
599 :
600 : END MODULE xtb_eeq
|
