(git:6d276e9)
Loading...
Searching...
No Matches
eeq_method.F
Go to the documentation of this file.
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 method
10!> \author JGH
11! **************************************************************************************************
12MODULE eeq_method
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind,&
15 get_atomic_kind_set
16 USE atprop_types, ONLY: atprop_type
17 USE cell_types, ONLY: cell_type,&
18 get_cell,&
19 pbc,&
20 plane_distance
21 USE cp_blacs_env, ONLY: cp_blacs_env_type
22 USE cp_control_types, ONLY: dft_control_type
23 USE cp_fm_basic_linalg, ONLY: cp_fm_invert,&
24 cp_fm_matvec,&
25 cp_fm_solve
26 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 cp_fm_struct_release,&
28 cp_fm_struct_type
29 USE cp_fm_types, ONLY: cp_fm_create,&
30 cp_fm_get_info,&
31 cp_fm_release,&
32 cp_fm_set_all,&
33 cp_fm_type
34 USE cp_log_handling, ONLY: cp_get_default_logger,&
35 cp_logger_get_default_unit_nr,&
36 cp_logger_type
37 USE cp_output_handling, ONLY: medium_print_level
38 USE distribution_1d_types, ONLY: distribution_1d_type
39 USE distribution_2d_types, ONLY: distribution_2d_type
40 USE eeq_data, ONLY: get_eeq_data
41 USE eeq_input, ONLY: eeq_solver_type
42 USE ewald_environment_types, ONLY: ewald_env_create,&
43 ewald_env_get,&
44 ewald_env_release,&
45 ewald_env_set,&
46 ewald_environment_type,&
47 read_ewald_section_tb
48 USE ewald_pw_types, ONLY: ewald_pw_create,&
49 ewald_pw_release,&
50 ewald_pw_type
51 USE input_section_types, ONLY: section_vals_get_subs_vals,&
52 section_vals_type
53 USE kinds, ONLY: dp,&
54 int_8
55 USE machine, ONLY: m_walltime
56 USE mathconstants, ONLY: oorootpi,&
57 twopi
58 USE mathlib, ONLY: invmat
59 USE message_passing, ONLY: mp_para_env_type
60 USE molecule_types, ONLY: molecule_type
61 USE particle_types, ONLY: particle_type
62 USE physcon, ONLY: bohr
63 USE pw_poisson_types, ONLY: do_ewald_spme
64 USE qs_dispersion_cnum, ONLY: cnumber_init,&
65 cnumber_release,&
66 dcnum_type
67 USE qs_dispersion_types, ONLY: qs_dispersion_release,&
68 qs_dispersion_type
69 USE qs_environment_types, ONLY: get_qs_env,&
70 qs_environment_type
71 USE qs_force_types, ONLY: qs_force_type
72 USE qs_kind_types, ONLY: get_qs_kind,&
73 qs_kind_type
74 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
75 neighbor_list_iterate,&
76 neighbor_list_iterator_create,&
77 neighbor_list_iterator_p_type,&
78 neighbor_list_iterator_release,&
79 neighbor_list_set_p_type,&
80 release_neighbor_list_sets
81 USE qs_neighbor_lists, ONLY: atom2d_build,&
82 atom2d_cleanup,&
83 build_neighbor_lists,&
84 local_atoms_type,&
85 pair_radius_setup
86 USE spme, ONLY: spme_forces,&
87 spme_potential,&
88 spme_virial
89 USE util, ONLY: sort
90 USE virial_methods, ONLY: virial_pair_force
91 USE virial_types, ONLY: virial_type
92#include "./base/base_uses.f90"
93
94 IMPLICIT NONE
95
96 PRIVATE
97
98 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eeq_method'
99
100 INTEGER, PARAMETER :: maxElem = 86
101
102 TYPE eeq_sparse_matrix_type
103 INTEGER, ALLOCATABLE, DIMENSION(:) :: col, row
104 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag, value
105 INTEGER :: nactive = 0
106 END TYPE eeq_sparse_matrix_type
107 ! covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 188-197)
108 ! values for metals decreased by 10 %
109 REAL(kind=dp), PARAMETER :: rcov(1:maxelem) = [&
110 & 0.32_dp, 0.46_dp, 1.20_dp, 0.94_dp, 0.77_dp, 0.75_dp, 0.71_dp, 0.63_dp, &
111 & 0.64_dp, 0.67_dp, 1.40_dp, 1.25_dp, 1.13_dp, 1.04_dp, 1.10_dp, 1.02_dp, &
112 & 0.99_dp, 0.96_dp, 1.76_dp, 1.54_dp, 1.33_dp, 1.22_dp, 1.21_dp, 1.10_dp, &
113 & 1.07_dp, 1.04_dp, 1.00_dp, 0.99_dp, 1.01_dp, 1.09_dp, 1.12_dp, 1.09_dp, &
114 & 1.15_dp, 1.10_dp, 1.14_dp, 1.17_dp, 1.89_dp, 1.67_dp, 1.47_dp, 1.39_dp, &
115 & 1.32_dp, 1.24_dp, 1.15_dp, 1.13_dp, 1.13_dp, 1.08_dp, 1.15_dp, 1.23_dp, &
116 & 1.28_dp, 1.26_dp, 1.26_dp, 1.23_dp, 1.32_dp, 1.31_dp, 2.09_dp, 1.76_dp, &
117 & 1.62_dp, 1.47_dp, 1.58_dp, 1.57_dp, 1.56_dp, 1.55_dp, 1.51_dp, 1.52_dp, &
118 & 1.51_dp, 1.50_dp, 1.49_dp, 1.49_dp, 1.48_dp, 1.53_dp, 1.46_dp, 1.37_dp, &
119 & 1.31_dp, 1.23_dp, 1.18_dp, 1.16_dp, 1.11_dp, 1.12_dp, 1.13_dp, 1.32_dp, &
120 & 1.30_dp, 1.30_dp, 1.36_dp, 1.31_dp, 1.38_dp, 1.42_dp]
121
122 PUBLIC :: eeq_solver, eeq_print, eeq_charges, eeq_forces, &
123 eeq_efield_energy, eeq_efield_pot, eeq_efield_force_loc, eeq_efield_force_periodic
124
125CONTAINS
126
127! **************************************************************************************************
128!> \brief ...
129!> \param qs_env ...
130!> \param iounit ...
131!> \param print_level ...
132!> \param ext ...
133! **************************************************************************************************
134 SUBROUTINE eeq_print(qs_env, iounit, print_level, ext)
135
136 TYPE(qs_environment_type), POINTER :: qs_env
137 INTEGER, INTENT(IN) :: iounit, print_level
138 LOGICAL, INTENT(IN) :: ext
139
140 CHARACTER(LEN=2) :: element_symbol
141 INTEGER :: enshift_type, iatom, ikind, natom
142 REAL(kind=dp), DIMENSION(:), POINTER :: charges
143 TYPE(cell_type), POINTER :: cell
144 TYPE(eeq_solver_type) :: eeq_sparam
145 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
146
147 mark_used(print_level)
148
149 CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, cell=cell)
150 IF (ext) THEN
151 NULLIFY (charges)
152 CALL get_qs_env(qs_env, eeq=charges)
153 cpassert(ASSOCIATED(charges))
154 enshift_type = 0
155 ELSE
156 ALLOCATE (charges(natom))
157 ! enforce en shift method 1 (original/molecular)
158 ! method 2 from paper on PBC seems not to work
159 enshift_type = 1
160 !IF (ALL(cell%perd == 0)) enshift_type = 1
161 CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
162 END IF
163
164 IF (iounit > 0) THEN
165
166 IF (enshift_type == 0) THEN
167 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (External)"
168 ELSE IF (enshift_type == 1) THEN
169 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Molecules))"
170 ELSE IF (enshift_type == 2) THEN
171 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Crystals))"
172 ELSE
173 cpabort("Unknown enshift_type")
174 END IF
175 WRITE (unit=iounit, fmt="(/,T2,A)") &
176 "# Atom Element Kind Atomic Charge"
177
178 DO iatom = 1, natom
179 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
180 element_symbol=element_symbol, &
181 kind_number=ikind)
182 WRITE (unit=iounit, fmt="(T4,I8,T18,A2,I10,T43,F12.4)") &
183 iatom, element_symbol, ikind, charges(iatom)
184 END DO
185
186 END IF
187
188 IF (.NOT. ext) DEALLOCATE (charges)
189
190 END SUBROUTINE eeq_print
191
192! **************************************************************************************************
193!> \brief ...
194!> \param qs_env ...
195!> \param charges ...
196!> \param eeq_sparam ...
197!> \param eeq_model ...
198!> \param enshift_type ...
199!> \param exclude ...
200!> \param cn_max ...
201! **************************************************************************************************
202 SUBROUTINE eeq_charges(qs_env, charges, eeq_sparam, eeq_model, enshift_type, exclude, cn_max)
203
204 TYPE(qs_environment_type), POINTER :: qs_env
205 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
206 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
207 INTEGER, INTENT(IN) :: eeq_model, enshift_type
208 LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: exclude
209 REAL(kind=dp), INTENT(IN), OPTIONAL :: cn_max
210
211 CHARACTER(len=*), PARAMETER :: routinen = 'eeq_charges'
212
213 INTEGER :: handle, iatom, ikind, iunit, jkind, &
214 natom, nkind, za, zb
215 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
216 INTEGER, DIMENSION(3) :: periodic
217 LOGICAL :: do_ewald
218 REAL(kind=dp) :: ala, alb, eeq_energy, esg, kappa, &
219 lambda, scn, sgamma, totalcharge, xi
220 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chia, cnumbers, efr, gam
221 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
222 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
223 TYPE(cell_type), POINTER :: cell, cell_ref
224 TYPE(cp_blacs_env_type), POINTER :: blacs_env
225 TYPE(cp_logger_type), POINTER :: logger
226 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
227 TYPE(dft_control_type), POINTER :: dft_control
228 TYPE(ewald_environment_type), POINTER :: ewald_env
229 TYPE(ewald_pw_type), POINTER :: ewald_pw
230 TYPE(mp_para_env_type), POINTER :: para_env
231 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
232 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
233 TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
234 print_section
235
236 CALL timeset(routinen, handle)
237
238 CALL get_qs_env(qs_env, &
239 qs_kind_set=qs_kind_set, &
240 atomic_kind_set=atomic_kind_set, &
241 particle_set=particle_set, &
242 para_env=para_env, &
243 blacs_env=blacs_env, &
244 cell=cell, &
245 dft_control=dft_control)
246 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
247
248 logger => cp_get_default_logger()
249 IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
250 iunit = cp_logger_get_default_unit_nr()
251 ELSE
252 iunit = -1
253 END IF
254
255 totalcharge = dft_control%charge
256
257 CALL get_cnumbers(qs_env, cnumbers, dcnum, .false.)
258
259 ! Apply smooth logistic CN cutoff to match multicharge/D4 behavior
260 IF (PRESENT(cn_max)) THEN
261 DO iatom = 1, natom
262 cnumbers(iatom) = log(1.0_dp + exp(cn_max)) - log(1.0_dp + exp(cn_max - cnumbers(iatom)))
263 END DO
264 END IF
265
266 ! gamma[a,b]
267 ALLOCATE (gab(nkind, nkind), gam(nkind))
268 gab = 0.0_dp
269 gam = 0.0_dp
270 DO ikind = 1, nkind
271 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
272 CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
273 DO jkind = 1, nkind
274 CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
275 CALL get_eeq_data(zb, eeq_model, rad=alb)
276 !
277 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
278 !
279 END DO
280 END DO
281
282 ! Override parameters for excluded kinds (ghost/floating atoms in BSSE):
283 ! huge hardness + zero coupling -> q = 0, no influence on other atoms.
284 IF (PRESENT(exclude)) THEN
285 DO ikind = 1, nkind
286 IF (exclude(ikind)) THEN
287 gam(ikind) = 1.0e30_dp
288 gab(ikind, :) = 0.0_dp
289 gab(:, ikind) = 0.0_dp
290 END IF
291 END DO
292 END IF
293
294 ! Chi[a,a]
295 sgamma = 8.0_dp ! see D4 for periodic systems paper
296 esg = 1.0_dp + exp(sgamma)
297 ALLOCATE (chia(natom))
298 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
299 DO iatom = 1, natom
300 ikind = kind_of(iatom)
301 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
302 CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
303 !
304 IF (enshift_type == 1) THEN
305 scn = cnumbers(iatom)/sqrt(cnumbers(iatom) + 1.0e-14_dp)
306 ELSE IF (enshift_type == 2) THEN
307 scn = log(esg/(esg - cnumbers(iatom)))
308 ELSE
309 cpabort("Unknown enshift_type")
310 END IF
311 chia(iatom) = xi - kappa*scn
312 !
313 END DO
314
315 ! Zero electronegativity for excluded atoms (ghost/floating in BSSE)
316 IF (PRESENT(exclude)) THEN
317 DO iatom = 1, natom
318 ikind = kind_of(iatom)
319 IF (exclude(ikind)) chia(iatom) = 0.0_dp
320 END DO
321 END IF
322
323 ! efield
324 IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
325 dft_control%apply_efield_field) THEN
326 ALLOCATE (efr(natom))
327 efr(1:natom) = 0.0_dp
328 CALL eeq_efield_pot(qs_env, efr)
329 chia(1:natom) = chia(1:natom) + efr(1:natom)
330 DEALLOCATE (efr)
331 END IF
332
333 CALL cnumber_release(cnumbers, dcnum, .false.)
334
335 CALL get_cell(cell, periodic=periodic)
336 do_ewald = .NOT. all(periodic == 0)
337 IF (do_ewald) THEN
338 ALLOCATE (ewald_env)
339 CALL ewald_env_create(ewald_env, para_env)
340 poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
341 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
342 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
343 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
344 CALL get_qs_env(qs_env, cell_ref=cell_ref)
345 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
346 silent=.true., pset="EEQ", cell_periodic=cell%perd)
347 ALLOCATE (ewald_pw)
348 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
349 !
350 CALL eeq_solver(charges, lambda, eeq_energy, &
351 particle_set, kind_of, cell, chia, gam, gab, &
352 para_env, blacs_env, dft_control, eeq_sparam, &
353 totalcharge=totalcharge, ewald=do_ewald, &
354 ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit, qs_env=qs_env)
355 !
356 CALL ewald_env_release(ewald_env)
357 CALL ewald_pw_release(ewald_pw)
358 DEALLOCATE (ewald_env, ewald_pw)
359 ELSE
360 CALL eeq_solver(charges, lambda, eeq_energy, &
361 particle_set, kind_of, cell, chia, gam, gab, &
362 para_env, blacs_env, dft_control, eeq_sparam, &
363 totalcharge=totalcharge, iounit=iunit, qs_env=qs_env)
364 END IF
365
366 DEALLOCATE (gab, gam, chia)
367
368 CALL timestop(handle)
369
370 END SUBROUTINE eeq_charges
371
372! **************************************************************************************************
373!> \brief ...
374!> \param qs_env ...
375!> \param charges ...
376!> \param dcharges ...
377!> \param gradient ...
378!> \param stress ...
379!> \param eeq_sparam ...
380!> \param eeq_model ...
381!> \param enshift_type ...
382!> \param response_only ...
383!> \param exclude ...
384!> \param cn_max ...
385! **************************************************************************************************
386 SUBROUTINE eeq_forces(qs_env, charges, dcharges, gradient, stress, &
387 eeq_sparam, eeq_model, enshift_type, response_only, exclude, cn_max)
388
389 TYPE(qs_environment_type), POINTER :: qs_env
390 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges
391 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient
392 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: stress
393 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
394 INTEGER, INTENT(IN) :: eeq_model, enshift_type
395 LOGICAL, INTENT(IN) :: response_only
396 LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: exclude
397 REAL(kind=dp), INTENT(IN), OPTIONAL :: cn_max
398
399 CHARACTER(len=*), PARAMETER :: routinen = 'eeq_forces'
400
401 INTEGER :: handle, i, ia, iatom, ikind, iunit, &
402 jatom, jkind, katom, natom, nkind, za, &
403 zb
404 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
405 INTEGER, DIMENSION(3) :: periodic
406 LOGICAL :: do_ewald, use_virial
407 LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
408 REAL(kind=dp) :: ala, alb, alpha, cn, ctot, dcnpdcn, dr, dr2, drk, elag, esg, fe, gam2, &
409 gama, grc, kappa, qlam, qq, qq1, qq2, rcut, scn, sgamma, subcells, totalcharge, xi
410 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: c_radius, cnumbers, gam, qlag
411 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: epforce, gab, pair_radius
412 REAL(kind=dp), DIMENSION(3) :: fdik, ri, rij, rik, rj
413 REAL(kind=dp), DIMENSION(3, 3) :: pvir
414 REAL(kind=dp), DIMENSION(:), POINTER :: chrgx, dchia
415 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
416 TYPE(atprop_type), POINTER :: atprop
417 TYPE(cell_type), POINTER :: cell, cell_ref
418 TYPE(cp_blacs_env_type), POINTER :: blacs_env
419 TYPE(cp_logger_type), POINTER :: logger
420 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
421 TYPE(dft_control_type), POINTER :: dft_control
422 TYPE(distribution_1d_type), POINTER :: distribution_1d, local_particles
423 TYPE(distribution_2d_type), POINTER :: distribution_2d
424 TYPE(ewald_environment_type), POINTER :: ewald_env
425 TYPE(ewald_pw_type), POINTER :: ewald_pw
426 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
427 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
428 TYPE(mp_para_env_type), POINTER :: para_env
429 TYPE(neighbor_list_iterator_p_type), &
430 DIMENSION(:), POINTER :: nl_iterator
431 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
432 POINTER :: sab_ew
433 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
434 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
435 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
436 TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
437 print_section
438 TYPE(virial_type), POINTER :: virial
439
440 CALL timeset(routinen, handle)
441
442 CALL get_qs_env(qs_env, &
443 qs_kind_set=qs_kind_set, &
444 atomic_kind_set=atomic_kind_set, &
445 particle_set=particle_set, &
446 para_env=para_env, &
447 blacs_env=blacs_env, &
448 cell=cell, &
449 force=force, &
450 virial=virial, &
451 atprop=atprop, &
452 dft_control=dft_control)
453
454 logger => cp_get_default_logger()
455 IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
456 iunit = cp_logger_get_default_unit_nr()
457 ELSE
458 iunit = -1
459 END IF
460
461 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
462 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
463
464 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
465
466 totalcharge = dft_control%charge
467
468 CALL get_cnumbers(qs_env, cnumbers, dcnum, .true.)
469
470 ! Apply smooth logistic CN cutoff to match multicharge/D4 behavior
471 ! Chain rule: dcn_cut/dR = (dcn_cut/dcn) * (dcn/dR)
472 IF (PRESENT(cn_max)) THEN
473 DO iatom = 1, natom
474 dcnpdcn = exp(cn_max)/(exp(cn_max) + exp(cnumbers(iatom)))
475 cnumbers(iatom) = log(1.0_dp + exp(cn_max)) - log(1.0_dp + exp(cn_max - cnumbers(iatom)))
476 DO i = 1, dcnum(iatom)%neighbors
477 dcnum(iatom)%dvals(i) = dcnum(iatom)%dvals(i)*dcnpdcn
478 END DO
479 END DO
480 END IF
481
482 ! gamma[a,b]
483 ALLOCATE (gab(nkind, nkind), gam(nkind))
484 gab = 0.0_dp
485 gam = 0.0_dp
486 DO ikind = 1, nkind
487 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
488 CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
489 DO jkind = 1, nkind
490 CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
491 CALL get_eeq_data(zb, eeq_model, rad=alb)
492 !
493 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
494 !
495 END DO
496 END DO
497
498 ! Override parameters for excluded kinds (ghost/floating atoms in BSSE)
499 IF (PRESENT(exclude)) THEN
500 DO ikind = 1, nkind
501 IF (exclude(ikind)) THEN
502 gam(ikind) = 1.0e30_dp
503 gab(ikind, :) = 0.0_dp
504 gab(:, ikind) = 0.0_dp
505 END IF
506 END DO
507 END IF
508
509 ALLOCATE (qlag(natom))
510
511 CALL get_cell(cell, periodic=periodic)
512 do_ewald = .NOT. all(periodic == 0)
513 IF (do_ewald) THEN
514 ALLOCATE (ewald_env)
515 CALL ewald_env_create(ewald_env, para_env)
516 poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
517 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
518 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
519 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
520 CALL get_qs_env(qs_env, cell_ref=cell_ref)
521 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
522 silent=.true., pset="EEQ", cell_periodic=cell%perd)
523 ALLOCATE (ewald_pw)
524 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
525 !
526 CALL eeq_solver(qlag, qlam, elag, &
527 particle_set, kind_of, cell, -dcharges, gam, gab, &
528 para_env, blacs_env, dft_control, eeq_sparam, &
529 ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit, qs_env=qs_env)
530 ELSE
531 CALL eeq_solver(qlag, qlam, elag, &
532 particle_set, kind_of, cell, -dcharges, gam, gab, &
533 para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit, qs_env=qs_env)
534 END IF
535
536 sgamma = 8.0_dp ! see D4 for periodic systems paper
537 esg = 1.0_dp + exp(sgamma)
538 ALLOCATE (chrgx(natom), dchia(natom))
539 DO iatom = 1, natom
540 ikind = kind_of(iatom)
541 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
542 CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
543 !
544 IF (response_only) THEN
545 ctot = -0.5_dp*qlag(iatom)
546 ELSE
547 ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
548 END IF
549 IF (enshift_type == 1) THEN
550 scn = sqrt(cnumbers(iatom)) + 1.0e-14_dp
551 dchia(iatom) = -ctot*kappa/scn
552 ELSE IF (enshift_type == 2) THEN
553 cn = cnumbers(iatom)
554 scn = 1.0_dp/(esg - cn)
555 dchia(iatom) = -ctot*kappa*scn
556 ELSE
557 cpabort("Unknown enshift_type")
558 END IF
559 END DO
560
561 ! Efield
562 IF (dft_control%apply_period_efield) THEN
563 CALL eeq_efield_force_periodic(qs_env, charges, qlag)
564 ELSE IF (dft_control%apply_efield) THEN
565 CALL eeq_efield_force_loc(qs_env, charges, qlag)
566 ELSE IF (dft_control%apply_efield_field) THEN
567 cpabort("apply field")
568 END IF
569
570 ! Forces from q*X
571 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
572 DO ikind = 1, nkind
573 DO ia = 1, local_particles%n_el(ikind)
574 iatom = local_particles%list(ikind)%array(ia)
575 DO i = 1, dcnum(iatom)%neighbors
576 katom = dcnum(iatom)%nlist(i)
577 rik = dcnum(iatom)%rik(:, i)
578 drk = norm2(rik)
579 IF (drk > 1.e-3_dp) THEN
580 fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
581 gradient(:, iatom) = gradient(:, iatom) - fdik(:)
582 gradient(:, katom) = gradient(:, katom) + fdik(:)
583 IF (use_virial) THEN
584 CALL virial_pair_force(stress, 1._dp, fdik, rik)
585 END IF
586 END IF
587 END DO
588 END DO
589 END DO
590
591 ! Forces from (0.5*q+l)*dA/dR*q
592 IF (do_ewald) THEN
593
594 ! Build the neighbor lists for the CN
595 CALL get_qs_env(qs_env, &
596 distribution_2d=distribution_2d, &
597 local_particles=distribution_1d, &
598 molecule_set=molecule_set)
599 subcells = 2.0_dp
600 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
601 rcut = 2.0_dp*rcut
602 NULLIFY (sab_ew)
603 ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
604 c_radius(:) = rcut
605 default_present = .true.
606 ALLOCATE (atom2d(nkind))
607 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
608 molecule_set, .false., particle_set=particle_set)
609 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
610 CALL build_neighbor_lists(sab_ew, particle_set, atom2d, cell, pair_radius, &
611 subcells=subcells, operator_type="PP", nlname="sab_ew")
612 DEALLOCATE (c_radius, pair_radius, default_present)
613 CALL atom2d_cleanup(atom2d)
614 !
615 CALL neighbor_list_iterator_create(nl_iterator, sab_ew)
616 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
617 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
618 iatom=iatom, jatom=jatom, r=rij)
619 !
620 dr2 = sum(rij**2)
621 dr = sqrt(dr2)
622 IF (dr > rcut .OR. dr < 1.e-6_dp) cycle
623 fe = 1.0_dp
624 IF (iatom == jatom) fe = 0.5_dp
625 IF (response_only) THEN
626 qq = -qlag(iatom)*charges(jatom)
627 ELSE
628 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
629 END IF
630 gama = gab(ikind, jkind)
631 gam2 = gama*gama
632 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
633 - 2._dp*alpha*exp(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
634 IF (response_only) THEN
635 qq1 = -qlag(iatom)*charges(jatom)
636 qq2 = -qlag(jatom)*charges(iatom)
637 ELSE
638 qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
639 qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
640 END IF
641 fdik(:) = -qq1*grc*rij(:)/dr
642 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
643 gradient(:, jatom) = gradient(:, jatom) - fdik(:)
644 IF (use_virial) THEN
645 CALL virial_pair_force(stress, -fe, fdik, rij)
646 END IF
647 fdik(:) = qq2*grc*rij(:)/dr
648 gradient(:, iatom) = gradient(:, iatom) - fdik(:)
649 gradient(:, jatom) = gradient(:, jatom) + fdik(:)
650 IF (use_virial) THEN
651 CALL virial_pair_force(stress, fe, fdik, rij)
652 END IF
653 END DO
654 CALL neighbor_list_iterator_release(nl_iterator)
655 !
656 CALL release_neighbor_list_sets(sab_ew)
657 ELSE
658 DO ikind = 1, nkind
659 DO ia = 1, local_particles%n_el(ikind)
660 iatom = local_particles%list(ikind)%array(ia)
661 ri(1:3) = particle_set(iatom)%r(1:3)
662 DO jatom = 1, natom
663 IF (iatom == jatom) cycle
664 jkind = kind_of(jatom)
665 IF (response_only) THEN
666 qq = -qlag(iatom)*charges(jatom)
667 ELSE
668 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
669 END IF
670 rj(1:3) = particle_set(jatom)%r(1:3)
671 rij(1:3) = ri(1:3) - rj(1:3)
672 rij = pbc(rij, cell)
673 dr2 = sum(rij**2)
674 dr = sqrt(dr2)
675 gama = gab(ikind, jkind)
676 gam2 = gama*gama
677 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
678 fdik(:) = qq*grc*rij(:)/dr
679 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
680 gradient(:, jatom) = gradient(:, jatom) - fdik(:)
681 END DO
682 END DO
683 END DO
684 END IF
685
686 ! Forces from Ewald potential: (q+l)*A*q
687 IF (do_ewald) THEN
688 ALLOCATE (epforce(3, natom))
689 epforce = 0.0_dp
690 IF (response_only) THEN
691 dchia(1:natom) = qlag(1:natom)
692 ELSE
693 dchia(1:natom) = -charges(1:natom) + qlag(1:natom)
694 END IF
695 chrgx(1:natom) = charges(1:natom)
696 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
697 particle_set, dchia, epforce)
698 dchia(1:natom) = charges(1:natom)
699 chrgx(1:natom) = qlag(1:natom)
700 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
701 particle_set, dchia, epforce)
702 gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + epforce(1:3, 1:natom)
703 DEALLOCATE (epforce)
704
705 ! virial
706 IF (use_virial) THEN
707 chrgx(1:natom) = charges(1:natom) - qlag(1:natom)
708 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
709 stress = stress - pvir
710 chrgx(1:natom) = qlag(1:natom)
711 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
712 stress = stress + pvir
713 IF (response_only) THEN
714 chrgx(1:natom) = charges(1:natom)
715 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
716 stress = stress + pvir
717 END IF
718 END IF
719 !
720 CALL ewald_env_release(ewald_env)
721 CALL ewald_pw_release(ewald_pw)
722 DEALLOCATE (ewald_env, ewald_pw)
723 END IF
724
725 CALL cnumber_release(cnumbers, dcnum, .true.)
726
727 DEALLOCATE (gab, gam, qlag, chrgx, dchia)
728
729 CALL timestop(handle)
730
731 END SUBROUTINE eeq_forces
732
733! **************************************************************************************************
734!> \brief ...
735!> \param qs_env ...
736!> \param cnumbers ...
737!> \param dcnum ...
738!> \param calculate_forces ...
739! **************************************************************************************************
740 SUBROUTINE get_cnumbers(qs_env, cnumbers, dcnum, calculate_forces)
741
742 TYPE(qs_environment_type), POINTER :: qs_env
743 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
744 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
745 LOGICAL, INTENT(IN) :: calculate_forces
746
747 INTEGER :: ikind, natom, nkind, za
748 LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
749 REAL(kind=dp) :: subcells
750 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: c_radius
751 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
752 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
753 TYPE(cell_type), POINTER :: cell
754 TYPE(distribution_1d_type), POINTER :: distribution_1d
755 TYPE(distribution_2d_type), POINTER :: distribution_2d
756 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
757 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
758 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
759 POINTER :: sab_cn
760 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
761 TYPE(qs_dispersion_type), POINTER :: disp
762 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
763
764 CALL get_qs_env(qs_env, &
765 qs_kind_set=qs_kind_set, &
766 atomic_kind_set=atomic_kind_set, &
767 particle_set=particle_set, &
768 cell=cell, &
769 distribution_2d=distribution_2d, &
770 local_particles=distribution_1d, &
771 molecule_set=molecule_set)
772 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
773
774 ! Check for dispersion_env and sab_cn needed for cnumbers
775 ALLOCATE (disp)
776 disp%k1 = 16.0_dp
777 disp%k2 = 4._dp/3._dp
778 disp%eps_cn = 1.e-6_dp
779 disp%max_elem = maxelem
780 ALLOCATE (disp%rcov(maxelem))
781 disp%rcov(1:maxelem) = bohr*disp%k2*rcov(1:maxelem)
782 subcells = 2.0_dp
783 ! Build the neighbor lists for the CN
784 NULLIFY (sab_cn)
785 ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
786 c_radius(:) = 0.0_dp
787 default_present = .true.
788 DO ikind = 1, nkind
789 CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
790 c_radius(ikind) = 4._dp*rcov(za)*bohr
791 END DO
792 ALLOCATE (atom2d(nkind))
793 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
794 molecule_set, .false., particle_set=particle_set)
795 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
796 CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
797 subcells=subcells, operator_type="PP", nlname="sab_cn")
798 disp%sab_cn => sab_cn
799 DEALLOCATE (c_radius, pair_radius, default_present)
800 CALL atom2d_cleanup(atom2d)
801
802 ! Calculate coordination numbers
803 CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces, disp_env=disp)
804
805 CALL qs_dispersion_release(disp)
806
807 END SUBROUTINE get_cnumbers
808
809! **************************************************************************************************
810!> \brief ...
811!> \param charges ...
812!> \param lambda ...
813!> \param eeq_energy ...
814!> \param particle_set ...
815!> \param kind_of ...
816!> \param cell ...
817!> \param chia ...
818!> \param gam ...
819!> \param gab ...
820!> \param para_env ...
821!> \param blacs_env ...
822!> \param dft_control ...
823!> \param eeq_sparam ...
824!> \param totalcharge ...
825!> \param ewald ...
826!> \param ewald_env ...
827!> \param ewald_pw ...
828!> \param iounit ...
829!> \param qs_env environment used to build the sparse real-space neighbor list
830! **************************************************************************************************
831 SUBROUTINE eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, &
832 chia, gam, gab, para_env, blacs_env, dft_control, eeq_sparam, &
833 totalcharge, ewald, ewald_env, ewald_pw, iounit, qs_env)
834
835 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
836 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
837 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
838 POINTER :: particle_set
839 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
840 TYPE(cell_type), POINTER :: cell
841 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
842 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
843 TYPE(mp_para_env_type), POINTER :: para_env
844 TYPE(cp_blacs_env_type), POINTER :: blacs_env
845 TYPE(dft_control_type), POINTER :: dft_control
846 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
847 REAL(kind=dp), INTENT(IN), OPTIONAL :: totalcharge
848 LOGICAL, INTENT(IN), OPTIONAL :: ewald
849 TYPE(ewald_environment_type), OPTIONAL, POINTER :: ewald_env
850 TYPE(ewald_pw_type), OPTIONAL, POINTER :: ewald_pw
851 INTEGER, INTENT(IN), OPTIONAL :: iounit
852 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
853
854 CHARACTER(len=*), PARAMETER :: routinen = 'eeq_solver'
855
856 INTEGER :: handle, ierror, iunit, natom, nkind, ns
857 LOGICAL :: do_direct, do_displ, do_ewald, &
858 do_sparse, do_sparse_auto, fm_created
859 REAL(kind=dp) :: alpha, deth, ftime, qtot
860 TYPE(cp_fm_struct_type), POINTER :: mat_struct
861 TYPE(cp_fm_type) :: eeq_mat
862
863 CALL timeset(routinen, handle)
864
865 do_ewald = .false.
866 IF (PRESENT(ewald)) do_ewald = ewald
867 !
868 qtot = 0.0_dp
869 IF (PRESENT(totalcharge)) qtot = totalcharge
870 !
871 iunit = -1
872 IF (PRESENT(iounit)) iunit = iounit
873
874 ! EEQ solver parameters
875 do_direct = eeq_sparam%direct
876 natom = SIZE(particle_set)
877 do_displ = .false.
878 IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
879 do_displ = dft_control%period_efield%displacement_field
880 END IF
881 do_sparse_auto = eeq_sparam%sparse_threshold > 0 .AND. natom >= eeq_sparam%sparse_threshold
882 ! An explicit request diagnoses unsupported displacement-field use below.
883 ! Automatic selection must retain the historical dense path instead of
884 ! turning a system-size threshold into a new input-dependent abort.
885 do_sparse = eeq_sparam%sparse .OR. (do_sparse_auto .AND. .NOT. do_displ)
886 fm_created = .false.
887
888 nkind = SIZE(gam)
889 ns = natom + 1
890 IF (.NOT. (do_ewald .AND. do_sparse .AND. .NOT. do_direct)) THEN
891 CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
892 nrow_global=ns, ncol_global=ns)
893 CALL cp_fm_create(eeq_mat, mat_struct)
894 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
895 fm_created = .true.
896 END IF
897 !
898 IF (do_ewald) THEN
899 cpassert(PRESENT(ewald_env))
900 cpassert(PRESENT(ewald_pw))
901 IF (do_direct) THEN
902 IF (do_displ) THEN
903 cpabort("NYA")
904 ELSE
905 CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
906 kind_of, cell, chia, gam, gab, qtot, &
907 ewald_env, ewald_pw, iounit)
908 END IF
909 ELSE IF (do_sparse) THEN
910 IF (do_displ) THEN
911 cpabort("Sparse periodic EEQ is not available with a displacement field")
912 ELSE
913 cpassert(PRESENT(qs_env))
914 ierror = 0
915 CALL pbc_sparse_solver(charges, lambda, eeq_energy, particle_set, &
916 kind_of, cell, chia, gam, gab, qtot, &
917 ewald_env, ewald_pw, eeq_sparam, qs_env, ierror, iounit)
918 IF (ierror /= 0) THEN
919 CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
920 nrow_global=ns, ncol_global=ns)
921 CALL cp_fm_create(eeq_mat, mat_struct)
922 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
923 fm_created = .true.
924 CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
925 kind_of, cell, chia, gam, gab, qtot, &
926 ewald_env, ewald_pw, iounit)
927 END IF
928 END IF
929 ELSE
930 IF (do_displ) THEN
931 cpabort("NYA")
932 ELSE
933 ierror = 0
934 CALL pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
935 kind_of, cell, chia, gam, gab, qtot, &
936 ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
937 IF (ierror /= 0) THEN
938 ! backup to non-iterative method
939 CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
940 kind_of, cell, chia, gam, gab, qtot, &
941 ewald_env, ewald_pw, iounit)
942 END IF
943 END IF
944 END IF
945 IF (qtot /= 0._dp) THEN
946 CALL get_cell(cell=cell, deth=deth)
947 CALL ewald_env_get(ewald_env, alpha=alpha)
948 eeq_energy = eeq_energy - 0.5_dp*qtot**2/alpha**2/deth
949 END IF
950 ELSE
951 CALL mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, &
952 cell, chia, gam, gab, qtot, ftime)
953 IF (iounit > 0) THEN
954 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
955 END IF
956 END IF
957 IF (fm_created) THEN
958 CALL cp_fm_struct_release(mat_struct)
959 CALL cp_fm_release(eeq_mat)
960 END IF
961
962 CALL timestop(handle)
963
964 END SUBROUTINE eeq_solver
965
966! **************************************************************************************************
967!> \brief ...
968!> \param charges ...
969!> \param lambda ...
970!> \param eeq_energy ...
971!> \param eeq_mat ...
972!> \param particle_set ...
973!> \param kind_of ...
974!> \param cell ...
975!> \param chia ...
976!> \param gam ...
977!> \param gab ...
978!> \param qtot ...
979!> \param ftime ...
980! **************************************************************************************************
981 SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
982 chia, gam, gab, qtot, ftime)
983
984 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
985 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
986 TYPE(cp_fm_type) :: eeq_mat
987 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
988 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
989 TYPE(cell_type), POINTER :: cell
990 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
991 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
992 REAL(kind=dp), INTENT(IN) :: qtot
993 REAL(kind=dp), INTENT(OUT) :: ftime
994
995 CHARACTER(len=*), PARAMETER :: routinen = 'mi_solver'
996
997 INTEGER :: handle, ia, iac, iar, ic, ikind, ir, &
998 jkind, natom, ncloc, ncvloc, nkind, &
999 nrloc, nrvloc, ns
1000 INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
1001 REAL(kind=dp) :: dr, grc, te, ti, xr
1002 REAL(kind=dp), DIMENSION(3) :: ri, rij, rj
1003 TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
1004 TYPE(cp_fm_type) :: rhs_vec
1005 TYPE(mp_para_env_type), POINTER :: para_env
1006
1007 CALL timeset(routinen, handle)
1008 ti = m_walltime()
1009
1010 natom = SIZE(particle_set)
1011 nkind = SIZE(gam)
1012 !
1013 ns = natom + 1
1014 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1015 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1016 row_indices=rind, col_indices=cind)
1017 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1018 nrow_global=ns, ncol_global=1)
1019 CALL cp_fm_create(rhs_vec, vec_struct)
1020 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1021 row_indices=rvind, col_indices=cvind)
1022 !
1023 ! set up matrix
1024 CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
1025 CALL cp_fm_set_all(rhs_vec, 0.0_dp)
1026 DO ir = 1, nrloc
1027 iar = rind(ir)
1028 IF (iar > natom) cycle
1029 ikind = kind_of(iar)
1030 ri(1:3) = particle_set(iar)%r(1:3)
1031 DO ic = 1, ncloc
1032 iac = cind(ic)
1033 IF (iac > natom) cycle
1034 jkind = kind_of(iac)
1035 rj(1:3) = particle_set(iac)%r(1:3)
1036 IF (iar == iac) THEN
1037 grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1038 ELSE
1039 rij(1:3) = ri(1:3) - rj(1:3)
1040 rij = pbc(rij, cell)
1041 dr = norm2(rij)
1042 grc = erf(gab(ikind, jkind)*dr)/dr
1043 END IF
1044 eeq_mat%local_data(ir, ic) = grc
1045 END DO
1046 END DO
1047 ! set up rhs vector
1048 DO ir = 1, nrvloc
1049 iar = rvind(ir)
1050 DO ic = 1, ncvloc
1051 iac = cvind(ic)
1052 ia = max(iar, iac)
1053 IF (ia > natom) THEN
1054 xr = qtot
1055 ELSE
1056 xr = -chia(ia)
1057 END IF
1058 rhs_vec%local_data(ir, ic) = xr
1059 END DO
1060 END DO
1061 !
1062 CALL cp_fm_solve(eeq_mat, rhs_vec)
1063 !
1064 charges = 0.0_dp
1065 lambda = 0.0_dp
1066 DO ir = 1, nrvloc
1067 iar = rvind(ir)
1068 DO ic = 1, ncvloc
1069 iac = cvind(ic)
1070 ia = max(iar, iac)
1071 IF (ia <= natom) THEN
1072 xr = rhs_vec%local_data(ir, ic)
1073 charges(ia) = xr
1074 ELSE
1075 lambda = rhs_vec%local_data(ir, ic)
1076 END IF
1077 END DO
1078 END DO
1079 CALL para_env%sum(lambda)
1080 CALL para_env%sum(charges)
1081 !
1082 ! energy: 0.5*(q^T.X - lambda*totalcharge)
1083 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1084
1085 CALL cp_fm_struct_release(vec_struct)
1086 CALL cp_fm_release(rhs_vec)
1087
1088 te = m_walltime()
1089 ftime = te - ti
1090 CALL timestop(handle)
1091
1092 END SUBROUTINE mi_solver
1093
1094! **************************************************************************************************
1095!> \brief Solve the periodic EEQ equations in the fixed-charge subspace with projected PCG.
1096!>
1097!> The short-range Gaussian/Ewald correction is stored as a distributed sparse pair list.
1098!> The long-range Coulomb response is applied matrix-free with SPME. A constrained diagonal
1099!> hardness preconditioner preserves the total charge exactly.
1100!> \param charges converged atomic charges
1101!> \param lambda charge-constraint multiplier
1102!> \param eeq_energy EEQ energy
1103!> \param particle_set particles
1104!> \param kind_of atomic kind for every particle
1105!> \param cell simulation cell
1106!> \param chia electronegativities
1107!> \param gam atomic hardness parameters
1108!> \param gab pair Gaussian exponents
1109!> \param qtot requested total charge
1110!> \param ewald_env Ewald environment
1111!> \param ewald_pw Ewald plane-wave data
1112!> \param eeq_sparam solver settings
1113!> \param qs_env environment used for the cutoff neighbor list
1114!> \param ierror zero on convergence
1115!> \param iounit output unit
1116! **************************************************************************************************
1117 SUBROUTINE pbc_sparse_solver(charges, lambda, eeq_energy, particle_set, &
1118 kind_of, cell, chia, gam, gab, qtot, &
1119 ewald_env, ewald_pw, eeq_sparam, qs_env, ierror, iounit)
1120
1121 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
1122 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
1123 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1124 POINTER :: particle_set
1125 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1126 TYPE(cell_type), POINTER :: cell
1127 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1128 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1129 REAL(kind=dp), INTENT(IN) :: qtot
1130 TYPE(ewald_environment_type), POINTER :: ewald_env
1131 TYPE(ewald_pw_type), POINTER :: ewald_pw
1132 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
1133 TYPE(qs_environment_type), POINTER :: qs_env
1134 INTEGER, INTENT(OUT) :: ierror
1135 INTEGER, OPTIONAL :: iounit
1136
1137 CHARACTER(len=*), PARAMETER :: routinen = 'pbc_sparse_solver'
1138
1139 INTEGER :: handle, i, iter, iunit, max_iter, natom
1140 REAL(kind=dp) :: alpha_cg, beta_cg, denom, eps_solver, &
1141 ftime, res, res_initial, rz, rz_new, &
1142 te, ti
1143 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ap, aq, gradient, p, precond_diag, &
1144 residual, work, z
1145 TYPE(eeq_sparse_matrix_type) :: response
1146
1147 CALL timeset(routinen, handle)
1148 ti = m_walltime()
1149
1150 iunit = -1
1151 IF (PRESENT(iounit)) iunit = iounit
1152 natom = SIZE(charges)
1153 max_iter = eeq_sparam%max_diis
1154 eps_solver = eeq_sparam%eps_diis
1155 ierror = 0
1156
1157 CALL build_sparse_pbc_response(response, qs_env, particle_set, kind_of, cell, gam, gab, &
1158 ewald_env, eps_solver, iunit)
1159
1160 ALLOCATE (ap(natom), aq(natom), gradient(natom), p(natom), precond_diag(natom), &
1161 residual(natom), work(natom), z(natom))
1162 DO i = 1, natom
1163 precond_diag(i) = gam(kind_of(i)) + 2.0_dp*gab(kind_of(i), kind_of(i))*oorootpi
1164 precond_diag(i) = max(precond_diag(i), sqrt(epsilon(1.0_dp)))
1165 END DO
1166
1167 ! Diagonal constrained minimizer as a deterministic charge-conserving initial guess.
1168 lambda = (-qtot - sum(chia/precond_diag))/sum(1.0_dp/precond_diag)
1169 charges = -(chia + lambda)/precond_diag
1170 charges = charges + (qtot - sum(charges))/real(natom, kind=dp)
1171
1172 CALL apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1173 charges, aq, work)
1174 gradient(:) = aq + chia
1175 residual(:) = -gradient
1176 CALL project_eeq_charge_tangent(residual)
1177 res = norm2(residual)
1178 res_initial = max(res, sqrt(epsilon(1.0_dp)))
1179 iter = 0
1180
1181 IF (res >= eps_solver) THEN
1182 CALL apply_eeq_constrained_diagonal_preconditioner(residual, precond_diag, z)
1183 p(:) = z
1184 rz = sum(residual*z)
1185
1186 DO iter = 1, max_iter
1187 CALL apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1188 p, ap, work)
1189 CALL project_eeq_charge_tangent(ap)
1190 denom = sum(p*ap)
1191 IF (.NOT. (denom > sqrt(tiny(1.0_dp)) .AND. abs(denom) < huge(denom))) THEN
1192 ierror = 1
1193 EXIT
1194 END IF
1195
1196 alpha_cg = rz/denom
1197 IF (.NOT. (abs(alpha_cg) < huge(alpha_cg))) THEN
1198 ierror = 1
1199 EXIT
1200 END IF
1201 charges = charges + alpha_cg*p
1202 residual(:) = residual - alpha_cg*ap
1203 ! Remove accumulated roundoff from both constrained spaces.
1204 charges = charges + (qtot - sum(charges))/real(natom, kind=dp)
1205 CALL project_eeq_charge_tangent(residual)
1206 res = norm2(residual)
1207 IF (res < eps_solver) EXIT
1208 IF (res > 100.0_dp*res_initial .OR. .NOT. (res < huge(res))) THEN
1209 ierror = 1
1210 EXIT
1211 END IF
1212
1213 CALL apply_eeq_constrained_diagonal_preconditioner(residual, precond_diag, z)
1214 rz_new = sum(residual*z)
1215 IF (.NOT. (rz_new > 0.0_dp .AND. abs(rz_new) < huge(rz_new))) THEN
1216 ierror = 1
1217 EXIT
1218 END IF
1219 beta_cg = rz_new/rz
1220 p(:) = z + beta_cg*p
1221 rz = rz_new
1222 END DO
1223 IF (iter > max_iter) ierror = 1
1224 END IF
1225
1226 IF (ierror == 0) THEN
1227 CALL apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1228 charges, aq, work)
1229 gradient(:) = aq + chia
1230 lambda = -sum(gradient)/real(natom, kind=dp)
1231 eeq_energy = 0.5_dp*sum(charges*aq) + sum(charges*chia)
1232 END IF
1233
1234 te = m_walltime()
1235 ftime = te - ti
1236 IF (iunit > 0) THEN
1237 IF (ierror /= 0) THEN
1238 WRITE (iunit, '(A,T50,I6,T61,E20.5)') &
1239 " EEQ| Sparse PBC solver failed: iterations/accuracy ", iter, res
1240 ELSE
1241 WRITE (iunit, '(A,T50,I6,T61,E20.5)') &
1242 " EEQ| Sparse PBC solver: iterations/accuracy ", iter, res
1243 END IF
1244 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Sparse PBC solver: time[s]", ftime
1245 END IF
1246
1247 DEALLOCATE (ap, aq, gradient, p, precond_diag, residual, work, z)
1248 CALL release_sparse_pbc_response(response)
1249 CALL timestop(handle)
1250
1251 END SUBROUTINE pbc_sparse_solver
1252
1253! **************************************************************************************************
1254!> \brief Build the distributed short-range part of the periodic EEQ response.
1255!> \param response ...
1256!> \param qs_env ...
1257!> \param particle_set ...
1258!> \param kind_of ...
1259!> \param cell ...
1260!> \param gam ...
1261!> \param gab ...
1262!> \param ewald_env ...
1263!> \param eps_solver ...
1264!> \param iunit ...
1265! **************************************************************************************************
1266 SUBROUTINE build_sparse_pbc_response(response, qs_env, particle_set, kind_of, cell, gam, gab, &
1267 ewald_env, eps_solver, iunit)
1268
1269 TYPE(eeq_sparse_matrix_type), INTENT(OUT) :: response
1270 TYPE(qs_environment_type), POINTER :: qs_env
1271 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1272 POINTER :: particle_set
1273 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1274 TYPE(cell_type), POINTER :: cell
1275 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: gam
1276 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1277 TYPE(ewald_environment_type), POINTER :: ewald_env
1278 REAL(kind=dp), INTENT(IN) :: eps_solver
1279 INTEGER, INTENT(IN) :: iunit
1280
1281 INTEGER :: iatom, ix, iy, iz, jatom, &
1282 n_candidates_global, n_pairs, &
1283 n_pairs_alloc, n_pairs_global, &
1284 n_unique, natom, owner
1285 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: pair_key
1286 INTEGER, ALLOCATABLE, DIMENSION(:) :: pair_order
1287 INTEGER, DIMENSION(3) :: periodic
1288 LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
1289 REAL(kind=dp) :: drop_tolerance, neighbor_cutoff, rmax, &
1290 subcells, value
1291 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: c_radius
1292 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
1293 REAL(kind=dp), DIMENSION(3) :: cell_extent, vertex
1294 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1295 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1296 TYPE(distribution_1d_type), POINTER :: distribution_1d
1297 TYPE(distribution_2d_type), POINTER :: distribution_2d
1298 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
1299 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1300 TYPE(mp_para_env_type), POINTER :: para_env
1301 TYPE(neighbor_list_iterator_p_type), &
1302 DIMENSION(:), POINTER :: nl_iterator
1303 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1304 POINTER :: sab_eeq
1305
1306 NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, nl_iterator, sab_eeq)
1307 CALL ewald_env_get(ewald_env, rcut=rmax, para_env=para_env)
1308 rmax = 2.0_dp*rmax
1309 natom = SIZE(particle_set)
1310 drop_tolerance = max(100.0_dp*epsilon(1.0_dp), &
1311 eps_solver/(1000.0_dp*real(max(1, natom), kind=dp)))
1312
1313 ! A minimum-image atom pair can never be farther away than the most distant vertex of the
1314 ! centered cell. Limiting the neighbor-list radius to that bound avoids materializing many
1315 ! redundant periodic images when the Ewald real-space cutoff exceeds a small simulation cell.
1316 CALL get_cell(cell, h=hmat, periodic=periodic)
1317 WHERE (periodic /= 0)
1318 cell_extent = 0.5_dp
1319 ELSE WHERE
1320 cell_extent = 1.0_dp
1321 END WHERE
1322 neighbor_cutoff = 0.0_dp
1323 DO ix = -1, 1, 2
1324 DO iy = -1, 1, 2
1325 DO iz = -1, 1, 2
1326 vertex = matmul(hmat, cell_extent*real([ix, iy, iz], kind=dp))
1327 neighbor_cutoff = max(neighbor_cutoff, norm2(vertex))
1328 END DO
1329 END DO
1330 END DO
1331 neighbor_cutoff = min(rmax, neighbor_cutoff + &
1332 1.0e-10_dp*max(1.0_dp, neighbor_cutoff))
1333
1334 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
1335 distribution_2d=distribution_2d, local_particles=distribution_1d, &
1336 molecule_set=molecule_set)
1337 ALLOCATE (c_radius(SIZE(atomic_kind_set)), default_present(SIZE(atomic_kind_set)), &
1338 pair_radius(SIZE(atomic_kind_set), SIZE(atomic_kind_set)), &
1339 atom2d(SIZE(atomic_kind_set)))
1340 c_radius = 0.5_dp*neighbor_cutoff
1341 default_present = .true.
1342 subcells = 2.0_dp
1343 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
1344 molecule_set, .false., particle_set=particle_set)
1345 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
1346 CALL build_neighbor_lists(sab_eeq, particle_set, atom2d, cell, pair_radius, &
1347 subcells=subcells, mic=.false., symmetric=.true., &
1348 operator_type="PP", nlname="sab_eeq_sparse")
1349 DEALLOCATE (c_radius, pair_radius, default_present)
1350 CALL atom2d_cleanup(atom2d)
1351
1352 ALLOCATE (response%diag(natom), source=0.0_dp)
1353 DO iatom = 1, natom
1354 owner = mod(iatom - 1, para_env%num_pe)
1355 IF (owner == para_env%mepos) THEN
1356 response%diag(iatom) = eeq_short_range_element(iatom, iatom, particle_set, kind_of, &
1357 cell, gam, gab, ewald_env)
1358 END IF
1359 END DO
1360
1361 ! Collect the canonical pair keys from the bounded image list. Sorting removes image-cell
1362 ! duplicates without an O(N_atom**2) lookup table. The symmetric list assigns every pair to
1363 ! one MPI owner; the full Ewald-cutoff image sum is evaluated once for each unique key.
1364 n_pairs_alloc = 0
1365 CALL neighbor_list_iterator_create(nl_iterator, sab_eeq)
1366 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1367 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
1368 IF (iatom /= jatom) n_pairs_alloc = n_pairs_alloc + 1
1369 END DO
1370 CALL neighbor_list_iterator_release(nl_iterator)
1371
1372 ALLOCATE (pair_key(n_pairs_alloc), pair_order(n_pairs_alloc))
1373 n_pairs_alloc = 0
1374 CALL neighbor_list_iterator_create(nl_iterator, sab_eeq)
1375 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1376 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
1377 IF (iatom == jatom) cycle
1378 n_pairs_alloc = n_pairs_alloc + 1
1379 pair_key(n_pairs_alloc) = int(min(iatom, jatom) - 1, int_8)*int(natom, int_8) + &
1380 int(max(iatom, jatom), int_8)
1381 END DO
1382 CALL neighbor_list_iterator_release(nl_iterator)
1383 CALL release_neighbor_list_sets(sab_eeq)
1384
1385 IF (n_pairs_alloc > 0) CALL sort(pair_key, n_pairs_alloc, pair_order)
1386 n_unique = 0
1387 DO iatom = 1, n_pairs_alloc
1388 IF (iatom == 1) THEN
1389 n_unique = n_unique + 1
1390 ELSE IF (pair_key(iatom) /= pair_key(iatom - 1)) THEN
1391 n_unique = n_unique + 1
1392 END IF
1393 END DO
1394 ALLOCATE (response%row(n_unique), response%col(n_unique), response%value(n_unique))
1395 n_pairs = 0
1396 DO iatom = 1, n_pairs_alloc
1397 IF (iatom > 1) THEN
1398 IF (pair_key(iatom) == pair_key(iatom - 1)) cycle
1399 END IF
1400 jatom = int(modulo(pair_key(iatom) - 1_int_8, int(natom, int_8))) + 1
1401 owner = int((pair_key(iatom) - 1_int_8)/int(natom, int_8)) + 1
1402 value = eeq_short_range_element(owner, jatom, particle_set, kind_of, cell, &
1403 gam, gab, ewald_env)
1404 IF (abs(value) <= drop_tolerance) cycle
1405 n_pairs = n_pairs + 1
1406 response%row(n_pairs) = owner
1407 response%col(n_pairs) = jatom
1408 response%value(n_pairs) = value
1409 END DO
1410 DEALLOCATE (pair_key, pair_order)
1411 response%nactive = n_pairs
1412
1413 n_pairs_global = n_pairs
1414 CALL para_env%sum(n_pairs_global)
1415 n_candidates_global = n_unique
1416 CALL para_env%sum(n_candidates_global)
1417 IF (iunit > 0) THEN
1418 WRITE (iunit, '(A,T54,I12)') " EEQ| Sparse short-range candidate atom pairs", n_candidates_global
1419 WRITE (iunit, '(A,T54,I12)') " EEQ| Sparse short-range retained atom pairs", n_pairs_global
1420 WRITE (iunit, '(A,T61,E20.5)') " EEQ| Sparse short-range drop tolerance", drop_tolerance
1421 WRITE (iunit, '(A,T61,E20.5)') " EEQ| Sparse short-range MIC cutoff", neighbor_cutoff
1422 END IF
1423
1424 END SUBROUTINE build_sparse_pbc_response
1425
1426! **************************************************************************************************
1427!> \brief Evaluate one aggregated short-range periodic EEQ response element.
1428!> \param iatom ...
1429!> \param jatom ...
1430!> \param particle_set ...
1431!> \param kind_of ...
1432!> \param cell ...
1433!> \param gam ...
1434!> \param gab ...
1435!> \param ewald_env ...
1436!> \return ...
1437! **************************************************************************************************
1438 FUNCTION eeq_short_range_element(iatom, jatom, particle_set, kind_of, cell, gam, gab, &
1439 ewald_env) RESULT(element)
1440
1441 INTEGER, INTENT(IN) :: iatom, jatom
1442 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1443 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1444 TYPE(cell_type), POINTER :: cell
1445 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: gam
1446 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1447 TYPE(ewald_environment_type), POINTER :: ewald_env
1448 REAL(kind=dp) :: element
1449
1450 INTEGER :: ix, iy, iz
1451 INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1452 REAL(kind=dp) :: alpha, dr, rcut, rmax
1453 REAL(kind=dp), DIMENSION(3) :: rij, rijl
1454 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1455
1456 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
1457 rmax = 2.0_dp*rcut
1458 CALL get_cell(cell, h=hmat, periodic=periodic)
1459 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1460 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1461 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1462 WHERE (periodic == 0) ncell = 0
1463
1464 rij = pbc(particle_set(iatom)%r - particle_set(jatom)%r, cell)
1465 element = 0.0_dp
1466 DO ix = -ncell(1), ncell(1)
1467 DO iy = -ncell(2), ncell(2)
1468 DO iz = -ncell(3), ncell(3)
1469 cvec = [ix, iy, iz]
1470 rijl = rij + matmul(hmat, cvec)
1471 dr = norm2(rijl)
1472 IF (dr > rmax) cycle
1473 IF (iatom == jatom .AND. dr < 1.0e-5_dp) THEN
1474 element = element + gam(kind_of(iatom)) + &
1475 2.0_dp*gab(kind_of(iatom), kind_of(iatom))*oorootpi - &
1476 2.0_dp*alpha*oorootpi
1477 ELSE IF (dr > 1.0e-12_dp) THEN
1478 element = element + &
1479 (erf(gab(kind_of(iatom), kind_of(jatom))*dr) - erf(alpha*dr))/dr
1480 END IF
1481 END DO
1482 END DO
1483 END DO
1484
1485 END FUNCTION eeq_short_range_element
1486
1487! **************************************************************************************************
1488!> \brief Apply sparse short range plus the matrix-free SPME long-range response.
1489!> \param response ...
1490!> \param ewald_env ...
1491!> \param ewald_pw ...
1492!> \param cell ...
1493!> \param particle_set ...
1494!> \param charges ...
1495!> \param potential ...
1496!> \param work ...
1497! **************************************************************************************************
1498 SUBROUTINE apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1499 charges, potential, work)
1500
1501 TYPE(eeq_sparse_matrix_type), INTENT(IN) :: response
1502 TYPE(ewald_environment_type), POINTER :: ewald_env
1503 TYPE(ewald_pw_type), POINTER :: ewald_pw
1504 TYPE(cell_type), POINTER :: cell
1505 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1506 REAL(kind=dp), DIMENSION(:), INTENT(IN), TARGET :: charges
1507 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: potential, work
1508
1509 INTEGER :: i, ia, ja
1510 TYPE(mp_para_env_type), POINTER :: para_env
1511
1512 CALL ewald_env_get(ewald_env, para_env=para_env)
1513 potential = response%diag*charges
1514 DO i = 1, response%nactive
1515 ia = response%row(i)
1516 ja = response%col(i)
1517 potential(ia) = potential(ia) + response%value(i)*charges(ja)
1518 potential(ja) = potential(ja) + response%value(i)*charges(ia)
1519 END DO
1520 CALL para_env%sum(potential)
1521
1522 CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, work)
1523 potential = potential + work
1524
1525 END SUBROUTINE apply_sparse_pbc_response
1526
1527! **************************************************************************************************
1528!> \brief Apply a positive diagonal EEQ preconditioner in the zero-total-charge subspace.
1529!> \param residual ...
1530!> \param diagonal ...
1531!> \param RESULT ...
1532! **************************************************************************************************
1533 SUBROUTINE apply_eeq_constrained_diagonal_preconditioner(residual, diagonal, RESULT)
1534
1535 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: residual, diagonal
1536 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: result
1537
1538 REAL(kind=dp) :: constraint_shift
1539
1540 result = residual/diagonal
1541 constraint_shift = sum(result)/sum(1.0_dp/diagonal)
1542 result = result - constraint_shift/diagonal
1543
1544 END SUBROUTINE apply_eeq_constrained_diagonal_preconditioner
1545
1546! **************************************************************************************************
1547!> \brief Project a vector onto the tangent space of the EEQ total-charge constraint.
1548!> \param vector ...
1549! **************************************************************************************************
1550 SUBROUTINE project_eeq_charge_tangent(vector)
1551
1552 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: vector
1553
1554 vector = vector - sum(vector)/real(SIZE(vector), kind=dp)
1555
1556 END SUBROUTINE project_eeq_charge_tangent
1557
1558! **************************************************************************************************
1559!> \brief Release a sparse periodic EEQ response.
1560!> \param response ...
1561! **************************************************************************************************
1562 SUBROUTINE release_sparse_pbc_response(response)
1563
1564 TYPE(eeq_sparse_matrix_type), INTENT(INOUT) :: response
1565
1566 IF (ALLOCATED(response%col)) DEALLOCATE (response%col)
1567 IF (ALLOCATED(response%row)) DEALLOCATE (response%row)
1568 IF (ALLOCATED(response%diag)) DEALLOCATE (response%diag)
1569 IF (ALLOCATED(response%value)) DEALLOCATE (response%value)
1570 response%nactive = 0
1571
1572 END SUBROUTINE release_sparse_pbc_response
1573
1574! **************************************************************************************************
1575!> \brief ...
1576!> \param charges ...
1577!> \param lambda ...
1578!> \param eeq_energy ...
1579!> \param eeq_mat ...
1580!> \param particle_set ...
1581!> \param kind_of ...
1582!> \param cell ...
1583!> \param chia ...
1584!> \param gam ...
1585!> \param gab ...
1586!> \param qtot ...
1587!> \param ewald_env ...
1588!> \param ewald_pw ...
1589!> \param eeq_sparam ...
1590!> \param ierror ...
1591!> \param iounit ...
1592! **************************************************************************************************
1593 SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1594 kind_of, cell, chia, gam, gab, qtot, &
1595 ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
1596
1597 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
1598 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
1599 TYPE(cp_fm_type) :: eeq_mat
1600 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1601 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1602 TYPE(cell_type), POINTER :: cell
1603 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1604 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1605 REAL(kind=dp), INTENT(IN) :: qtot
1606 TYPE(ewald_environment_type), POINTER :: ewald_env
1607 TYPE(ewald_pw_type), POINTER :: ewald_pw
1608 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
1609 INTEGER, INTENT(OUT) :: ierror
1610 INTEGER, OPTIONAL :: iounit
1611
1612 CHARACTER(len=*), PARAMETER :: routinen = 'pbc_solver'
1613
1614 INTEGER :: ewald_type, handle, i, iac, iar, ic, ikind, info, ir, iunit, iv, ix, iy, iz, &
1615 jkind, max_diis, mdiis, natom, ncloc, ndiis, nkind, now, nrloc, ns, sdiis
1616 INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1617 INTEGER, DIMENSION(:), POINTER :: cind, rind
1618 REAL(kind=dp) :: ad, alpha, astep, deth, dr, eeqn, &
1619 eps_diis, ftime, grc1, grc2, rcut, &
1620 res, resin, rmax, te, ti
1621 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1622 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat, fvec, vmat, xvec
1623 REAL(kind=dp), DIMENSION(3) :: ri, rij, rijl, rj
1624 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1625 REAL(kind=dp), DIMENSION(:), POINTER :: rhs, rv0, xv0
1626 TYPE(cp_fm_struct_type), POINTER :: mat_struct
1627 TYPE(cp_fm_type) :: mmat, pmat
1628 TYPE(mp_para_env_type), POINTER :: para_env
1629
1630 CALL timeset(routinen, handle)
1631 ti = m_walltime()
1632
1633 ierror = 0
1634
1635 iunit = -1
1636 IF (PRESENT(iounit)) iunit = iounit
1637
1638 natom = SIZE(particle_set)
1639 nkind = SIZE(gam)
1640 !
1641 CALL get_cell(cell=cell, deth=deth)
1642 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1643 ad = 2.0_dp*alpha*oorootpi
1644 IF (ewald_type /= do_ewald_spme) THEN
1645 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1646 END IF
1647 !
1648 rmax = 2.0_dp*rcut
1649 ! max cells used
1650 CALL get_cell(cell, h=hmat, periodic=periodic)
1651 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1652 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1653 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1654 IF (periodic(1) == 0) ncell(1) = 0
1655 IF (periodic(2) == 0) ncell(2) = 0
1656 IF (periodic(3) == 0) ncell(3) = 0
1657 !
1658 CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1659 chia, gam, gab, qtot, ftime)
1660 IF (iunit > 0) THEN
1661 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1662 END IF
1663 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1664 CALL cp_fm_create(pmat, mat_struct)
1665 CALL cp_fm_create(mmat, mat_struct)
1666 !
1667 ! response matrix
1668 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1669 row_indices=rind, col_indices=cind)
1670 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1671 DO ir = 1, nrloc
1672 iar = rind(ir)
1673 ri = 0.0_dp
1674 IF (iar <= natom) THEN
1675 ikind = kind_of(iar)
1676 ri(1:3) = particle_set(iar)%r(1:3)
1677 END IF
1678 DO ic = 1, ncloc
1679 iac = cind(ic)
1680 IF (iac > natom .AND. iar > natom) THEN
1681 eeq_mat%local_data(ir, ic) = 0.0_dp
1682 cycle
1683 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1684 eeq_mat%local_data(ir, ic) = 1.0_dp
1685 cycle
1686 END IF
1687 jkind = kind_of(iac)
1688 rj(1:3) = particle_set(iac)%r(1:3)
1689 rij(1:3) = ri(1:3) - rj(1:3)
1690 rij = pbc(rij, cell)
1691 DO ix = -ncell(1), ncell(1)
1692 DO iy = -ncell(2), ncell(2)
1693 DO iz = -ncell(3), ncell(3)
1694 cvec = [ix, iy, iz]
1695 rijl = rij + matmul(hmat, cvec)
1696 dr = norm2(rijl)
1697 IF (dr > rmax) cycle
1698 IF (iar == iac .AND. dr < 0.00001_dp) THEN
1699 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1700 ELSE
1701 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1702 END IF
1703 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1704 END DO
1705 END DO
1706 END DO
1707 END DO
1708 END DO
1709 !
1710 ! preconditioner
1711 CALL cp_fm_get_info(pmat, nrow_local=nrloc, ncol_local=ncloc, &
1712 row_indices=rind, col_indices=cind)
1713 CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1714 DO ir = 1, nrloc
1715 iar = rind(ir)
1716 ri = 0.0_dp
1717 IF (iar <= natom) THEN
1718 ikind = kind_of(iar)
1719 ri(1:3) = particle_set(iar)%r(1:3)
1720 END IF
1721 DO ic = 1, ncloc
1722 iac = cind(ic)
1723 IF (iac > natom .AND. iar > natom) THEN
1724 pmat%local_data(ir, ic) = 0.0_dp
1725 cycle
1726 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1727 pmat%local_data(ir, ic) = 1.0_dp
1728 cycle
1729 END IF
1730 jkind = kind_of(iac)
1731 rj(1:3) = particle_set(iac)%r(1:3)
1732 rij(1:3) = ri(1:3) - rj(1:3)
1733 rij = pbc(rij, cell)
1734 IF (iar == iac) THEN
1735 grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1736 ELSE
1737 grc2 = erf(gab(ikind, jkind)*dr)/dr
1738 END IF
1739 pmat%local_data(ir, ic) = grc2
1740 END DO
1741 END DO
1742 CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1743 ! preconditioner invers
1744 CALL cp_fm_invert(pmat, mmat)
1745 !
1746 ! rhs
1747 ns = natom + 1
1748 ALLOCATE (rhs(ns))
1749 rhs(1:natom) = chia(1:natom)
1750 rhs(ns) = -qtot
1751 !
1752 ALLOCATE (xv0(ns), rv0(ns))
1753 ! initial guess
1754 xv0(1:natom) = charges(1:natom)
1755 xv0(ns) = 0.0_dp
1756 ! DIIS optimizer
1757 max_diis = eeq_sparam%max_diis
1758 mdiis = eeq_sparam%mdiis
1759 sdiis = eeq_sparam%sdiis
1760 eps_diis = eeq_sparam%eps_diis
1761 astep = eeq_sparam%alpha
1762 ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1763 xvec = 0.0_dp; fvec = 0.0_dp
1764 ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1765 dmat = 0.0_dp; dvec = 0.0_dp
1766 ndiis = 1
1767 now = 1
1768 CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1769 cell, particle_set, xv0, rhs, rv0)
1770 resin = norm2(rv0)
1771 !
1772 DO iv = 1, max_diis
1773 res = norm2(rv0)
1774 IF (res > 10._dp*resin) EXIT
1775 IF (res < eps_diis) EXIT
1776 !
1777 now = mod(iv - 1, mdiis) + 1
1778 ndiis = min(iv, mdiis)
1779 xvec(1:ns, now) = xv0(1:ns)
1780 fvec(1:ns, now) = rv0(1:ns)
1781 DO i = 1, ndiis
1782 vmat(now, i) = sum(fvec(:, now)*fvec(:, i))
1783 vmat(i, now) = vmat(now, i)
1784 END DO
1785 IF (ndiis < sdiis) THEN
1786 xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1787 ELSE
1788 dvec = 0.0_dp
1789 dvec(ndiis + 1) = 1.0_dp
1790 dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1791 dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1792 dmat(1:ndiis, ndiis + 1) = 1.0_dp
1793 dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1794 CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1795 dvec(1:ndiis + 1) = matmul(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1796 xv0(1:ns) = matmul(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1797 xv0(1:ns) = xv0(1:ns) + matmul(fvec(1:ns, 1:ndiis), dvec(1:ndiis))
1798 END IF
1799 !
1800 CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1801 cell, particle_set, xv0, rhs, rv0)
1802 END DO
1803 charges(1:natom) = xv0(1:natom)
1804 lambda = xv0(ns)
1805 eeq_energy = eeqn
1806 IF (res > eps_diis) ierror = 1
1807 !
1808 DEALLOCATE (xvec, fvec, bvec)
1809 DEALLOCATE (vmat, dmat, dvec)
1810 DEALLOCATE (xv0, rv0)
1811 DEALLOCATE (rhs)
1812 CALL cp_fm_release(pmat)
1813 CALL cp_fm_release(mmat)
1814
1815 te = m_walltime()
1816 IF (iunit > 0) THEN
1817 IF (ierror == 1) THEN
1818 WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1819 ELSE
1820 WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1821 END IF
1822 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1823 END IF
1824 CALL timestop(handle)
1825
1826 END SUBROUTINE pbc_solver
1827
1828! **************************************************************************************************
1829!> \brief ...
1830!> \param charges ...
1831!> \param lambda ...
1832!> \param eeq_energy ...
1833!> \param eeq_mat ...
1834!> \param particle_set ...
1835!> \param kind_of ...
1836!> \param cell ...
1837!> \param chia ...
1838!> \param gam ...
1839!> \param gab ...
1840!> \param qtot ...
1841!> \param ewald_env ...
1842!> \param ewald_pw ...
1843!> \param iounit ...
1844! **************************************************************************************************
1845 SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1846 kind_of, cell, chia, gam, gab, qtot, ewald_env, ewald_pw, iounit)
1847
1848 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
1849 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
1850 TYPE(cp_fm_type) :: eeq_mat
1851 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1852 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1853 TYPE(cell_type), POINTER :: cell
1854 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1855 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1856 REAL(kind=dp), INTENT(IN) :: qtot
1857 TYPE(ewald_environment_type), POINTER :: ewald_env
1858 TYPE(ewald_pw_type), POINTER :: ewald_pw
1859 INTEGER, INTENT(IN), OPTIONAL :: iounit
1860
1861 CHARACTER(len=*), PARAMETER :: routinen = 'fpbc_solver'
1862
1863 INTEGER :: ewald_type, handle, ia, iac, iar, ic, &
1864 ikind, ir, iunit, ix, iy, iz, jkind, &
1865 natom, ncloc, ncvloc, nkind, nrloc, &
1866 nrvloc, ns
1867 INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1868 INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
1869 REAL(kind=dp) :: ad, alpha, deth, dr, grc1, rcut, rmax, &
1870 te, ti, xr
1871 REAL(kind=dp), DIMENSION(3) :: ri, rij, rijl, rj
1872 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1873 REAL(kind=dp), DIMENSION(:), POINTER :: pval, xval
1874 TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
1875 TYPE(cp_fm_type) :: rhs_vec
1876 TYPE(mp_para_env_type), POINTER :: para_env
1877
1878 CALL timeset(routinen, handle)
1879 ti = m_walltime()
1880
1881 iunit = -1
1882 IF (PRESENT(iounit)) iunit = iounit
1883
1884 natom = SIZE(particle_set)
1885 nkind = SIZE(gam)
1886 ns = natom + 1
1887 !
1888 CALL get_cell(cell=cell, deth=deth)
1889 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1890 ad = 2.0_dp*alpha*oorootpi
1891 IF (ewald_type /= do_ewald_spme) THEN
1892 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1893 END IF
1894 !
1895 rmax = 2.0_dp*rcut
1896 ! max cells used
1897 CALL get_cell(cell, h=hmat, periodic=periodic)
1898 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1899 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1900 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1901 IF (periodic(1) == 0) ncell(1) = 0
1902 IF (periodic(2) == 0) ncell(2) = 0
1903 IF (periodic(3) == 0) ncell(3) = 0
1904 !
1905 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1906 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1907 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1908 row_indices=rind, col_indices=cind)
1909 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1910 nrow_global=ns, ncol_global=1)
1911 CALL cp_fm_create(rhs_vec, vec_struct)
1912 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1913 row_indices=rvind, col_indices=cvind)
1914 ! response matrix
1915 DO ir = 1, nrloc
1916 iar = rind(ir)
1917 ri = 0.0_dp
1918 IF (iar <= natom) THEN
1919 ikind = kind_of(iar)
1920 ri(1:3) = particle_set(iar)%r(1:3)
1921 END IF
1922 DO ic = 1, ncloc
1923 iac = cind(ic)
1924 IF (iac > natom .AND. iar > natom) THEN
1925 eeq_mat%local_data(ir, ic) = 0.0_dp
1926 cycle
1927 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1928 eeq_mat%local_data(ir, ic) = 1.0_dp
1929 cycle
1930 END IF
1931 jkind = kind_of(iac)
1932 rj(1:3) = particle_set(iac)%r(1:3)
1933 rij(1:3) = ri(1:3) - rj(1:3)
1934 rij = pbc(rij, cell)
1935 DO ix = -ncell(1), ncell(1)
1936 DO iy = -ncell(2), ncell(2)
1937 DO iz = -ncell(3), ncell(3)
1938 cvec = [ix, iy, iz]
1939 rijl = rij + matmul(hmat, cvec)
1940 dr = norm2(rijl)
1941 IF (dr > rmax) cycle
1942 IF (iar == iac .AND. dr < 0.0001_dp) THEN
1943 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1944 ELSE
1945 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1946 END IF
1947 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1948 END DO
1949 END DO
1950 END DO
1951 END DO
1952 END DO
1953 !
1954 ALLOCATE (xval(natom), pval(natom))
1955 DO ia = 1, natom
1956 xval = 0.0_dp
1957 xval(ia) = 1.0_dp
1958 CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1959 !
1960 DO ir = 1, nrloc
1961 iar = rind(ir)
1962 IF (iar /= ia) cycle
1963 DO ic = 1, ncloc
1964 iac = cind(ic)
1965 IF (iac > natom) cycle
1966 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1967 END DO
1968 END DO
1969 END DO
1970 DEALLOCATE (xval, pval)
1971 !
1972 ! set up rhs vector
1973 DO ir = 1, nrvloc
1974 iar = rvind(ir)
1975 DO ic = 1, ncvloc
1976 iac = cvind(ic)
1977 ia = max(iar, iac)
1978 IF (ia > natom) THEN
1979 xr = qtot
1980 ELSE
1981 xr = -chia(ia)
1982 END IF
1983 rhs_vec%local_data(ir, ic) = xr
1984 END DO
1985 END DO
1986 !
1987 CALL cp_fm_solve(eeq_mat, rhs_vec)
1988 !
1989 charges = 0.0_dp
1990 lambda = 0.0_dp
1991 DO ir = 1, nrvloc
1992 iar = rvind(ir)
1993 DO ic = 1, ncvloc
1994 iac = cvind(ic)
1995 ia = max(iar, iac)
1996 IF (ia <= natom) THEN
1997 xr = rhs_vec%local_data(ir, ic)
1998 charges(ia) = xr
1999 ELSE
2000 lambda = rhs_vec%local_data(ir, ic)
2001 END IF
2002 END DO
2003 END DO
2004 CALL para_env%sum(lambda)
2005 CALL para_env%sum(charges)
2006 !
2007 ! energy: 0.5*(q^T.X - lambda*totalcharge)
2008 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
2009
2010 CALL cp_fm_struct_release(vec_struct)
2011 CALL cp_fm_release(rhs_vec)
2012
2013 te = m_walltime()
2014 IF (iunit > 0) THEN
2015 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
2016 END IF
2017 CALL timestop(handle)
2018
2019 END SUBROUTINE fpbc_solver
2020
2021! **************************************************************************************************
2022!> \brief ...
2023!> \param ewald_env ...
2024!> \param ewald_pw ...
2025!> \param cell ...
2026!> \param particle_set ...
2027!> \param charges ...
2028!> \param potential ...
2029! **************************************************************************************************
2030 SUBROUTINE apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, potential)
2031 TYPE(ewald_environment_type), POINTER :: ewald_env
2032 TYPE(ewald_pw_type), POINTER :: ewald_pw
2033 TYPE(cell_type), POINTER :: cell
2034 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
2035 REAL(kind=dp), DIMENSION(:), INTENT(IN), TARGET :: charges
2036 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: potential
2037
2038 TYPE(mp_para_env_type), POINTER :: para_env
2039
2040 CALL ewald_env_get(ewald_env, para_env=para_env)
2041 potential = 0.0_dp
2042 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
2043 particle_set, potential)
2044 CALL para_env%sum(potential)
2045
2046 END SUBROUTINE apply_potential
2047
2048! **************************************************************************************************
2049!> \brief ...
2050!> \param eeqn ...
2051!> \param fm_mat ...
2052!> \param mmat ...
2053!> \param ewald_env ...
2054!> \param ewald_pw ...
2055!> \param cell ...
2056!> \param particle_set ...
2057!> \param charges ...
2058!> \param rhs ...
2059!> \param potential ...
2060! **************************************************************************************************
2061 SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
2062 cell, particle_set, charges, rhs, potential)
2063 REAL(kind=dp), INTENT(INOUT) :: eeqn
2064 TYPE(cp_fm_type), INTENT(IN) :: fm_mat, mmat
2065 TYPE(ewald_environment_type), POINTER :: ewald_env
2066 TYPE(ewald_pw_type), POINTER :: ewald_pw
2067 TYPE(cell_type), POINTER :: cell
2068 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
2069 REAL(kind=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
2070 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rhs
2071 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: potential
2072
2073 INTEGER :: na, ns
2074 REAL(kind=dp) :: lambda
2075 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mvec
2076 TYPE(mp_para_env_type), POINTER :: para_env
2077
2078 ns = SIZE(charges)
2079 na = ns - 1
2080 CALL ewald_env_get(ewald_env, para_env=para_env)
2081 potential = 0.0_dp
2082 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
2083 particle_set, potential(1:na))
2084 CALL para_env%sum(potential(1:na))
2085 CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
2086 eeqn = 0.5_dp*sum(charges(1:na)*potential(1:na)) + sum(charges(1:na)*rhs(1:na))
2087 potential(1:ns) = potential(1:ns) + rhs(1:ns)
2088 ALLOCATE (mvec(ns))
2089 CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
2090 lambda = -sum(mvec(1:na))/real(na, kind=dp)
2091 potential(1:na) = mvec(1:na) + lambda
2092 DEALLOCATE (mvec)
2093
2094 END SUBROUTINE get_energy_gradient
2095
2096! **************************************************************************************************
2097!> \brief ...
2098!> \param qs_env ...
2099!> \param charges ...
2100!> \param ef_energy ...
2101! **************************************************************************************************
2102 SUBROUTINE eeq_efield_energy(qs_env, charges, ef_energy)
2103 TYPE(qs_environment_type), POINTER :: qs_env
2104 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges
2105 REAL(kind=dp), INTENT(OUT) :: ef_energy
2106
2107 COMPLEX(KIND=dp) :: zdeta
2108 COMPLEX(KIND=dp), DIMENSION(3) :: zi
2109 INTEGER :: ia, idir, natom
2110 LOGICAL :: dfield
2111 REAL(kind=dp) :: kr, omega, q
2112 REAL(kind=dp), DIMENSION(3) :: ci, dfilter, fieldpol, fpolvec, kvec, &
2113 qi, ria
2114 REAL(kind=dp), DIMENSION(3, 3) :: hmat
2115 TYPE(cell_type), POINTER :: cell
2116 TYPE(dft_control_type), POINTER :: dft_control
2117 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2118
2119 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
2120 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2121
2122 IF (dft_control%apply_period_efield) THEN
2123 dfield = dft_control%period_efield%displacement_field
2124
2125 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2126 cpabort("use of strength_list not implemented for eeq_efield_energy")
2127 END IF
2128
2129 fieldpol = dft_control%period_efield%polarisation
2130 fieldpol = fieldpol/norm2(fieldpol)
2131 fieldpol = -fieldpol*dft_control%period_efield%strength
2132 hmat = cell%hmat(:, :)/twopi
2133 DO idir = 1, 3
2134 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
2135 + fieldpol(3)*hmat(3, idir)
2136 END DO
2137
2138 zi(:) = cmplx(1._dp, 0._dp, dp)
2139 DO ia = 1, natom
2140 q = charges(ia)
2141 ria = particle_set(ia)%r
2142 ria = pbc(ria, cell)
2143 DO idir = 1, 3
2144 kvec(:) = twopi*cell%h_inv(idir, :)
2145 kr = sum(kvec(:)*ria(:))
2146 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
2147 zi(idir) = zi(idir)*zdeta
2148 END DO
2149 END DO
2150 qi = aimag(log(zi))
2151 IF (dfield) THEN
2152 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
2153 omega = cell%deth
2154 ci = matmul(hmat, qi)/omega
2155 ef_energy = 0.0_dp
2156 DO idir = 1, 3
2157 ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
2158 END DO
2159 ef_energy = -0.25_dp*omega/twopi*ef_energy
2160 ELSE
2161 ef_energy = sum(fpolvec(:)*qi(:))
2162 END IF
2163
2164 ELSE IF (dft_control%apply_efield) THEN
2165
2166 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2167 dft_control%efield_fields(1)%efield%strength
2168
2169 ef_energy = 0.0_dp
2170 DO ia = 1, natom
2171 ria = particle_set(ia)%r
2172 ria = pbc(ria, cell)
2173 q = charges(ia)
2174 ef_energy = ef_energy - q*sum(fieldpol*ria)
2175 END DO
2176
2177 ELSE
2178 cpabort("apply field")
2179 END IF
2180
2181 END SUBROUTINE eeq_efield_energy
2182
2183! **************************************************************************************************
2184!> \brief ...
2185!> \param qs_env ...
2186!> \param efr ...
2187! **************************************************************************************************
2188 SUBROUTINE eeq_efield_pot(qs_env, efr)
2189 TYPE(qs_environment_type), POINTER :: qs_env
2190 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: efr
2191
2192 INTEGER :: ia, idir, natom
2193 LOGICAL :: dfield
2194 REAL(kind=dp) :: kr
2195 REAL(kind=dp), DIMENSION(3) :: fieldpol, fpolvec, kvec, ria
2196 REAL(kind=dp), DIMENSION(3, 3) :: hmat
2197 TYPE(cell_type), POINTER :: cell
2198 TYPE(dft_control_type), POINTER :: dft_control
2199 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2200
2201 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
2202 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2203
2204 IF (dft_control%apply_period_efield) THEN
2205 dfield = dft_control%period_efield%displacement_field
2206
2207 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2208 cpabort("use of strength_list not implemented for eeq_efield_pot")
2209 END IF
2210
2211 fieldpol = dft_control%period_efield%polarisation
2212 fieldpol = fieldpol/norm2(fieldpol)
2213 fieldpol = -fieldpol*dft_control%period_efield%strength
2214 hmat = cell%hmat(:, :)/twopi
2215 DO idir = 1, 3
2216 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
2217 + fieldpol(3)*hmat(3, idir)
2218 END DO
2219
2220 IF (dfield) THEN
2221 ! dE/dq depends on q, postpone calculation
2222 efr = 0.0_dp
2223 ELSE
2224 efr = 0.0_dp
2225 DO ia = 1, natom
2226 ria = particle_set(ia)%r
2227 ria = pbc(ria, cell)
2228 DO idir = 1, 3
2229 kvec(:) = twopi*cell%h_inv(idir, :)
2230 kr = sum(kvec(:)*ria(:))
2231 efr(ia) = efr(ia) + kr*fpolvec(idir)
2232 END DO
2233 END DO
2234 END IF
2235
2236 ELSE IF (dft_control%apply_efield) THEN
2237
2238 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2239 dft_control%efield_fields(1)%efield%strength
2240
2241 DO ia = 1, natom
2242 ria = particle_set(ia)%r
2243 ria = pbc(ria, cell)
2244 efr(ia) = -sum(fieldpol*ria)
2245 END DO
2246
2247 ELSE
2248 cpabort("apply field")
2249 END IF
2250
2251 END SUBROUTINE eeq_efield_pot
2252
2253! **************************************************************************************************
2254!> \brief ...
2255!> \param charges ...
2256!> \param dft_control ...
2257!> \param particle_set ...
2258!> \param cell ...
2259!> \param efr ...
2260! **************************************************************************************************
2261 SUBROUTINE eeq_dfield_pot(charges, dft_control, particle_set, cell, efr)
2262 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges
2263 TYPE(dft_control_type), POINTER :: dft_control
2264 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
2265 TYPE(cell_type), POINTER :: cell
2266 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: efr
2267
2268 COMPLEX(KIND=dp) :: zdeta
2269 COMPLEX(KIND=dp), DIMENSION(3) :: zi
2270 INTEGER :: ia, idir, natom
2271 REAL(kind=dp) :: kr, omega, q
2272 REAL(kind=dp), DIMENSION(3) :: ci, dfilter, fieldpol, kvec, qi, ria
2273 REAL(kind=dp), DIMENSION(3, 3) :: hmat
2274
2275 natom = SIZE(particle_set)
2276
2277 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2278 cpabort("use of strength_list not implemented for eeq_dfield_pot")
2279 END IF
2280
2281 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
2282 fieldpol = dft_control%period_efield%polarisation
2283 fieldpol = fieldpol/norm2(fieldpol)
2284 fieldpol = -fieldpol*dft_control%period_efield%strength
2285 hmat = cell%hmat(:, :)/twopi
2286 omega = cell%deth
2287 !
2288 zi(:) = cmplx(1._dp, 0._dp, dp)
2289 DO ia = 1, natom
2290 q = charges(ia)
2291 ria = particle_set(ia)%r
2292 ria = pbc(ria, cell)
2293 DO idir = 1, 3
2294 kvec(:) = twopi*cell%h_inv(idir, :)
2295 kr = sum(kvec(:)*ria(:))
2296 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
2297 zi(idir) = zi(idir)*zdeta
2298 END DO
2299 END DO
2300 qi = aimag(log(zi))
2301 ci = matmul(hmat, qi)/omega
2302 ci = dfilter*(fieldpol - 2._dp*twopi*ci)
2303 DO ia = 1, natom
2304 ria = particle_set(ia)%r
2305 ria = pbc(ria, cell)
2306 efr(ia) = efr(ia) - sum(ci*ria)
2307 END DO
2308
2309 END SUBROUTINE eeq_dfield_pot
2310
2311! **************************************************************************************************
2312!> \brief ...
2313!> \param qs_env ...
2314!> \param charges ...
2315!> \param qlag ...
2316! **************************************************************************************************
2317 SUBROUTINE eeq_efield_force_loc(qs_env, charges, qlag)
2318 TYPE(qs_environment_type), POINTER :: qs_env
2319 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
2320
2321 INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
2322 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
2323 REAL(kind=dp) :: q
2324 REAL(kind=dp), DIMENSION(3) :: fieldpol
2325 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2326 TYPE(cell_type), POINTER :: cell
2327 TYPE(dft_control_type), POINTER :: dft_control
2328 TYPE(distribution_1d_type), POINTER :: local_particles
2329 TYPE(mp_para_env_type), POINTER :: para_env
2330 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2331 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2332
2333 CALL get_qs_env(qs_env=qs_env, &
2334 dft_control=dft_control, &
2335 cell=cell, particle_set=particle_set, &
2336 nkind=nkind, natom=natom, &
2337 para_env=para_env, &
2338 local_particles=local_particles)
2339
2340 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2341 dft_control%efield_fields(1)%efield%strength
2342
2343 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
2344 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
2345 CALL get_qs_env(qs_env=qs_env, force=force)
2346
2347 DO ikind = 1, nkind
2348 force(ikind)%efield = 0.0_dp
2349 DO ia = 1, local_particles%n_el(ikind)
2350 iatom = local_particles%list(ikind)%array(ia)
2351 q = charges(iatom) - qlag(iatom)
2352 atom_a = atom_of_kind(iatom)
2353 force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
2354 END DO
2355 CALL para_env%sum(force(ikind)%efield)
2356 END DO
2357
2358 END SUBROUTINE eeq_efield_force_loc
2359
2360! **************************************************************************************************
2361!> \brief ...
2362!> \param qs_env ...
2363!> \param charges ...
2364!> \param qlag ...
2365! **************************************************************************************************
2366 SUBROUTINE eeq_efield_force_periodic(qs_env, charges, qlag)
2367 TYPE(qs_environment_type), POINTER :: qs_env
2368 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
2369
2370 INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
2371 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
2372 LOGICAL :: dfield, use_virial
2373 REAL(kind=dp) :: q
2374 REAL(kind=dp), DIMENSION(3) :: fa, fieldpol, ria
2375 REAL(kind=dp), DIMENSION(3, 3) :: pve
2376 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2377 TYPE(cell_type), POINTER :: cell
2378 TYPE(dft_control_type), POINTER :: dft_control
2379 TYPE(distribution_1d_type), POINTER :: local_particles
2380 TYPE(mp_para_env_type), POINTER :: para_env
2381 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2382 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2383 TYPE(virial_type), POINTER :: virial
2384
2385 CALL get_qs_env(qs_env=qs_env, &
2386 dft_control=dft_control, &
2387 cell=cell, particle_set=particle_set, &
2388 virial=virial, &
2389 nkind=nkind, natom=natom, &
2390 para_env=para_env, &
2391 local_particles=local_particles)
2392
2393 dfield = dft_control%period_efield%displacement_field
2394 cpassert(.NOT. dfield)
2395
2396 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2397 cpabort("use of strength_list not implemented for eeq_efield_force_periodic")
2398 END IF
2399
2400 fieldpol = dft_control%period_efield%polarisation
2401 fieldpol = fieldpol/norm2(fieldpol)
2402 fieldpol = -fieldpol*dft_control%period_efield%strength
2403
2404 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
2405
2406 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
2407 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
2408 CALL get_qs_env(qs_env=qs_env, force=force)
2409
2410 pve = 0.0_dp
2411 DO ikind = 1, nkind
2412 force(ikind)%efield = 0.0_dp
2413 DO ia = 1, local_particles%n_el(ikind)
2414 iatom = local_particles%list(ikind)%array(ia)
2415 q = charges(iatom) - qlag(iatom)
2416 fa(1:3) = q*fieldpol(1:3)
2417 atom_a = atom_of_kind(iatom)
2418 force(ikind)%efield(1:3, atom_a) = fa
2419 IF (use_virial) THEN
2420 ria = particle_set(ia)%r
2421 ria = pbc(ria, cell)
2422 CALL virial_pair_force(pve, -0.5_dp, fa, ria)
2423 CALL virial_pair_force(pve, -0.5_dp, ria, fa)
2424 END IF
2425 END DO
2426 CALL para_env%sum(force(ikind)%efield)
2427 END DO
2428 virial%pv_virial = virial%pv_virial + pve
2429
2430 END SUBROUTINE eeq_efield_force_periodic
2431
2432END MODULE eeq_method
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cell_types::pbc
Definition cell_types.F:104
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
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
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:324
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_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_solve
subroutine, public cp_fm_solve(matrix_a, general_a)
computes the the solution to A*b=A_general using lu decomposition
Definition cp_fm_basic_linalg.F:2109
cp_fm_basic_linalg::cp_fm_invert
subroutine, public cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix.
Definition cp_fm_basic_linalg.F:1790
cp_fm_basic_linalg::cp_fm_matvec
subroutine, public cp_fm_matvec(amat, xv, yv, alpha, beta)
Calculates yv = alpha*amat*xv + beta*yv where amat: fm matrix xv : vector replicated yv : vector repl...
Definition cp_fm_basic_linalg.F:2878
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
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:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:73
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
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
eeq_data
EEQ data from different sources.
Definition eeq_data.F:12
eeq_data::get_eeq_data
subroutine, public get_eeq_data(za, model, chi, eta, kcn, rad)
...
Definition eeq_data.F:240
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:2103
eeq_method::eeq_charges
subroutine, public eeq_charges(qs_env, charges, eeq_sparam, eeq_model, enshift_type, exclude, cn_max)
...
Definition eeq_method.F:203
eeq_method::eeq_efield_force_periodic
subroutine, public eeq_efield_force_periodic(qs_env, charges, qlag)
...
Definition eeq_method.F:2367
eeq_method::rcov
real(kind=dp), dimension(1:maxelem), parameter rcov
Definition eeq_method.F:109
eeq_method::eeq_efield_pot
subroutine, public eeq_efield_pot(qs_env, efr)
...
Definition eeq_method.F:2189
eeq_method::eeq_print
subroutine, public eeq_print(qs_env, iounit, print_level, ext)
...
Definition eeq_method.F:135
eeq_method::eeq_forces
subroutine, public eeq_forces(qs_env, charges, dcharges, gradient, stress, eeq_sparam, eeq_model, enshift_type, response_only, exclude, cn_max)
...
Definition eeq_method.F:388
eeq_method::eeq_efield_force_loc
subroutine, public eeq_efield_force_loc(qs_env, charges, qlag)
...
Definition eeq_method.F:2318
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, qs_env)
...
Definition eeq_method.F:834
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_set
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
Definition ewald_environment_types.F:189
ewald_environment_types::ewald_env_create
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
Definition ewald_environment_types.F:230
ewald_environment_types::read_ewald_section_tb
subroutine, public read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset, cell_periodic)
Purpose: read the EWALD section for TB methods.
Definition ewald_environment_types.F:398
ewald_environment_types::ewald_env_release
subroutine, public ewald_env_release(ewald_env)
releases the given ewald_env (see doc/ReferenceCounting.html)
Definition ewald_environment_types.F:247
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:132
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
ewald_pw_types::ewald_pw_release
subroutine, public ewald_pw_release(ewald_pw)
releases the memory used by the ewald_pw
Definition ewald_pw_types.F:102
ewald_pw_types::ewald_pw_create
subroutine, public ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
creates the structure ewald_pw_type
Definition ewald_pw_types.F:83
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::invmat
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition mathlib.F:551
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1633
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:76
qs_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
qs_dispersion_cnum::cnumber_release
subroutine, public cnumber_release(cnumbers, dcnum, derivatives)
...
Definition qs_dispersion_cnum.F:1598
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1536
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_dispersion_types::qs_dispersion_release
subroutine, public qs_dispersion_release(dispersion_env)
...
Definition qs_dispersion_types.F:143
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, mimic, 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, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
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, cneo_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, 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:491
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::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1096
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:550
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition qs_neighbor_lists.F:19
qs_neighbor_lists::atom2d_cleanup
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
Definition qs_neighbor_lists.F:137
qs_neighbor_lists::pair_radius_setup
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
Definition qs_neighbor_lists.F:1689
qs_neighbor_lists::build_neighbor_lists
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep, stable_images)
Build simple pair neighbor lists.
Definition qs_neighbor_lists.F:1083
qs_neighbor_lists::atom2d_build
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Definition qs_neighbor_lists.F:179
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_potential
subroutine, public spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, potential)
Calculate the electrostatic potential from particles A (charge A) at positions of particles B.
Definition spme.F:427
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
util
All kind of helpful little routines.
Definition util.F:14
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
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:60
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:758
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_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:64
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
molecule_types::molecule_type
Definition molecule_types.F:159
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
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:190
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
qs_neighbor_lists::local_atoms_type
Definition qs_neighbor_lists.F:106
virial_types::virial_type
Definition virial_types.F:26