(git:ed6f26b)
Loading...
Searching...
No Matches
xtb_eeq.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
65CONTAINS
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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chia, efr, gam
95 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
96 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 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 CALL timeset(routinen, handle)
110
111 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 dft_control=dft_control)
120 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
121
122 xtb_control => dft_control%qs_control%xtb_control
123
124 totalcharge = dft_control%charge
125
126 IF (atprop%energy) THEN
127 CALL atprop_array_init(atprop%atecoul, natom)
128 END IF
129
130 ! gamma[a,b]
131 ALLOCATE (gab(nkind, nkind), gam(nkind))
132 gab = 0.0_dp
133 gam = 0.0_dp
134 DO ikind = 1, nkind
135 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
136 CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
137 IF (.NOT. defined) cycle
138 CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
139 gam(ikind) = gama
140 DO jkind = 1, nkind
141 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
142 CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
143 IF (.NOT. defined) cycle
144 CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
145 !
146 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
147 !
148 END DO
149 END DO
150
151 ! Chi[a,a]
152 enshift_type = xtb_control%enshift_type
153 IF (enshift_type == 0) THEN
154 enshift_type = 2
155 IF (all(cell%perd == 0)) enshift_type = 1
156 END IF
157 sgamma = 8.0_dp ! see D4 for periodic systems paper
158 esg = 1.0_dp + exp(sgamma)
159 ALLOCATE (chia(natom))
160 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
161 DO iatom = 1, natom
162 ikind = kind_of(iatom)
163 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
164 CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
165 !
166 IF (enshift_type == 1) THEN
167 scn = sqrt(cnumbers(iatom)) + 1.0e-14_dp
168 ELSE IF (enshift_type == 2) THEN
169 cn = cnumbers(iatom)/esg
170 scn = log(esg/(esg - cnumbers(iatom)))
171 ELSE
172 cpabort("Unknown enshift_type")
173 END IF
174 chia(iatom) = xi - kappa*scn
175 !
176 END DO
177
178 ef_energy = 0.0_dp
179 IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
180 dft_control%apply_efield_field) THEN
181 ALLOCATE (efr(natom))
182 efr(1:natom) = 0.0_dp
183 CALL eeq_efield_pot(qs_env, efr)
184 chia(1:natom) = chia(1:natom) + efr(1:natom)
185 END IF
186
187 do_ewald = xtb_control%do_ewald
188
189 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
190 IF (do_ewald) THEN
191 CALL get_qs_env(qs_env=qs_env, &
192 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 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 totalcharge=totalcharge, iounit=iunit)
203 END IF
204
205 IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
206 dft_control%apply_efield_field) THEN
207 CALL eeq_efield_energy(qs_env, charges, ef_energy)
208 eeq_energy = eeq_energy - sum(charges*efr)
209 DEALLOCATE (efr)
210 END IF
211
212 DEALLOCATE (gab, gam, chia)
213
214 CALL timestop(handle)
215
216 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 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 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: chrgx, dchia, qlag
254 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 DIMENSION(:), POINTER :: nl_iterator
265 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
266 POINTER :: sab_tbe
267 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
268 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
269 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 CALL timeset(routinen, handle)
275
276 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 dft_control=dft_control)
287 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
288 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
289
290 xtb_control => dft_control%qs_control%xtb_control
291
292 totalcharge = dft_control%charge
293
294 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
295 atom_of_kind=atom_of_kind, kind_of=kind_of)
296
297 ! gamma[a,b]
298 ALLOCATE (gab(nkind, nkind), gam(nkind))
299 gab = 0.0_dp
300 DO ikind = 1, nkind
301 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
302 CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
303 IF (.NOT. defined) cycle
304 CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
305 gam(ikind) = gama
306 DO jkind = 1, nkind
307 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
308 CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
309 IF (.NOT. defined) cycle
310 CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
311 !
312 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
313 !
314 END DO
315 END DO
316
317 ALLOCATE (qlag(natom))
318
319 do_ewald = xtb_control%do_ewald
320
321 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
322 IF (do_ewald) THEN
323 CALL get_qs_env(qs_env=qs_env, &
324 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 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 para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
333 END IF
334
335 enshift_type = xtb_control%enshift_type
336 IF (enshift_type == 0) THEN
337 enshift_type = 2
338 IF (all(cell%perd == 0)) enshift_type = 1
339 END IF
340 sgamma = 8.0_dp ! see D4 for periodic systems paper
341 esg = 1.0_dp + exp(sgamma)
342 ALLOCATE (chrgx(natom), dchia(natom))
343 DO iatom = 1, natom
344 ikind = kind_of(iatom)
345 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
346 CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
347 !
348 ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
349 IF (enshift_type == 1) THEN
350 scn = sqrt(cnumbers(iatom)) + 1.0e-14_dp
351 dchia(iatom) = -ctot*kappa/scn
352 ELSE IF (enshift_type == 2) THEN
353 cn = cnumbers(iatom)
354 scn = 1.0_dp/(esg - cn)
355 dchia(iatom) = -ctot*kappa*scn
356 ELSE
357 cpabort("Unknown enshift_type")
358 END IF
359 END DO
360
361 ! Efield
362 IF (dft_control%apply_period_efield) THEN
363 CALL eeq_efield_force_periodic(qs_env, charges, qlag)
364 ELSE IF (dft_control%apply_efield) THEN
365 CALL eeq_efield_force_loc(qs_env, charges, qlag)
366 ELSE IF (dft_control%apply_efield_field) THEN
367 cpabort("apply field")
368 END IF
369
370 ! Forces from q*X
371 CALL get_qs_env(qs_env=qs_env, &
372 local_particles=local_particles)
373 DO ikind = 1, nkind
374 DO ia = 1, local_particles%n_el(ikind)
375 iatom = local_particles%list(ikind)%array(ia)
376 atom_a = atom_of_kind(iatom)
377 DO i = 1, dcnum(iatom)%neighbors
378 katom = dcnum(iatom)%nlist(i)
379 kkind = kind_of(katom)
380 atom_c = atom_of_kind(katom)
381 rik = dcnum(iatom)%rik(:, i)
382 drk = sqrt(sum(rik(:)**2))
383 IF (drk > 1.e-3_dp) THEN
384 fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
385 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
386 force(kkind)%rho_elec(:, atom_c) = force(kkind)%rho_elec(:, atom_c) + fdik(:)
387 IF (use_virial) THEN
388 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 IF (do_ewald) THEN
397 CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
398 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
399 rcut = 1.0_dp*rcut
400 CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
401 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
402 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
403 iatom=iatom, jatom=jatom, r=rij)
404 atom_a = atom_of_kind(iatom)
405 atom_b = atom_of_kind(jatom)
406 !
407 dr2 = sum(rij**2)
408 dr = sqrt(dr2)
409 IF (dr > rcut .OR. dr < 1.e-6_dp) cycle
410 fe = 1.0_dp
411 IF (iatom == jatom) fe = 0.5_dp
412 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
413 gama = gab(ikind, jkind)
414 gam2 = gama*gama
415 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
416 - 2._dp*alpha*exp(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
417 qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
418 qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
419 fdik(:) = -qq1*grc*rij(:)/dr
420 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
421 force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
422 IF (use_virial) THEN
423 CALL virial_pair_force(virial%pv_virial, fe, fdik, rij)
424 END IF
425 fdik(:) = qq2*grc*rij(:)/dr
426 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
427 force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
428 IF (use_virial) THEN
429 CALL virial_pair_force(virial%pv_virial, -fe, fdik, rij)
430 END IF
431 END DO
432 CALL neighbor_list_iterator_release(nl_iterator)
433 ELSE
434 DO ikind = 1, nkind
435 DO ia = 1, local_particles%n_el(ikind)
436 iatom = local_particles%list(ikind)%array(ia)
437 atom_a = atom_of_kind(iatom)
438 ri(1:3) = particle_set(iatom)%r(1:3)
439 DO jatom = 1, natom
440 IF (iatom == jatom) cycle
441 jkind = kind_of(jatom)
442 atom_b = atom_of_kind(jatom)
443 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
444 rj(1:3) = particle_set(jatom)%r(1:3)
445 rij(1:3) = ri(1:3) - rj(1:3)
446 rij = pbc(rij, cell)
447 dr2 = sum(rij**2)
448 dr = sqrt(dr2)
449 gama = gab(ikind, jkind)
450 gam2 = gama*gama
451 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
452 fdik(:) = qq*grc*rij(:)/dr
453 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
454 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 IF (do_ewald) THEN
462 ALLOCATE (epforce(3, natom))
463 epforce = 0.0_dp
464 dchia = -charges + qlag
465 chrgx = charges
466 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
467 particle_set, dchia, epforce)
468 dchia = charges
469 chrgx = qlag
470 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
471 particle_set, dchia, epforce)
472 DO iatom = 1, natom
473 ikind = kind_of(iatom)
474 i = atom_of_kind(iatom)
475 force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + epforce(:, iatom)
476 END DO
477 DEALLOCATE (epforce)
478
479 ! virial
480 IF (use_virial) THEN
481 chrgx = charges - qlag
482 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
483 virial%pv_virial = virial%pv_virial + pvir
484 chrgx = qlag
485 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
486 virial%pv_virial = virial%pv_virial - pvir
487 END IF
488 END IF
489
490 ! return Lagrange multipliers
491 qlagrange(1:natom) = qlag(1:natom)
492
493 DEALLOCATE (gab, chrgx, dchia, qlag)
494
495 CALL timestop(handle)
496
497 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 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: qlag
522 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
530 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 CALL timeset(routinen, handle)
535
536 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 dft_control=dft_control)
544 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
545
546 xtb_control => dft_control%qs_control%xtb_control
547
548 totalcharge = dft_control%charge
549
550 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
551 atom_of_kind=atom_of_kind, kind_of=kind_of)
552
553 ! gamma[a,b]
554 ALLOCATE (gab(nkind, nkind), gam(nkind))
555 gab = 0.0_dp
556 DO ikind = 1, nkind
557 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
558 CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
559 IF (.NOT. defined) cycle
560 CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
561 gam(ikind) = gama
562 DO jkind = 1, nkind
563 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
564 CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
565 IF (.NOT. defined) cycle
566 CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
567 !
568 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
569 !
570 END DO
571 END DO
572
573 ALLOCATE (qlag(natom))
574
575 do_ewald = xtb_control%do_ewald
576
577 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
578 IF (do_ewald) THEN
579 CALL get_qs_env(qs_env=qs_env, &
580 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 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 para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
589 END IF
590
591 ! return Lagrange multipliers
592 qlagrange(1:natom) = qlag(1:natom)
593
594 DEALLOCATE (gab, qlag)
595
596 CALL timestop(handle)
597
598 END SUBROUTINE xtb_eeq_lagrange
599
600END MODULE xtb_eeq
cell_types::pbc
Definition cell_types.F:92
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition atprop_types.F:87
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
eeq_input
Input definition and setup for EEQ model.
Definition eeq_input.F:12
eeq_method
Calculation of charge equilibration method.
Definition eeq_method.F:12
eeq_method::eeq_efield_energy
subroutine, public eeq_efield_energy(qs_env, charges, ef_energy)
...
Definition eeq_method.F:1510
eeq_method::eeq_efield_force_periodic
subroutine, public eeq_efield_force_periodic(qs_env, charges, qlag)
...
Definition eeq_method.F:1774
eeq_method::eeq_efield_pot
subroutine, public eeq_efield_pot(qs_env, efr)
...
Definition eeq_method.F:1596
eeq_method::eeq_efield_force_loc
subroutine, public eeq_efield_force_loc(qs_env, charges, qlag)
...
Definition eeq_method.F:1725
eeq_method::eeq_solver
subroutine, public eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, chia, gam, gab, para_env, blacs_env, dft_control, eeq_sparam, totalcharge, ewald, ewald_env, ewald_pw, iounit)
...
Definition eeq_method.F:753
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:124
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
spme
Calculate the electrostatic energy by the Smooth Particle Ewald method.
Definition spme.F:14
spme::spme_forces
subroutine, public spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, charges_b, forces_b)
Calculate the forces on particles B for the electrostatic interaction betrween particles A and B.
Definition spme.F:568
spme::spme_virial
subroutine, public spme_virial(ewald_env, ewald_pw, particle_set, box, mcharge, virial)
Internal Virial for 1/2 [rho||rho] (rho=mcharge)
Definition spme.F:719
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
xtb_eeq
Calculation of charge equilibration in xTB.
Definition xtb_eeq.F:12
xtb_eeq::xtb_eeq_lagrange
subroutine, public xtb_eeq_lagrange(qs_env, dcharges, qlagrange, eeq_sparam)
...
Definition xtb_eeq.F:507
xtb_eeq::xtb_eeq_calculation
subroutine, public xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, lambda)
...
Definition xtb_eeq.F:79
xtb_eeq::xtb_eeq_forces
subroutine, public xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
...
Definition xtb_eeq.F:229
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::xtb_control_type
Definition cp_control_types.F:138
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
eeq_input::eeq_solver_type
Definition eeq_input.F:27
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42