Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of charge equilibration method
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE 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 :
125 : CONTAINS
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param qs_env ...
130 : !> \param iounit ...
131 : !> \param print_level ...
132 : !> \param ext ...
133 : ! **************************************************************************************************
134 38 : 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 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
143 : TYPE(cell_type), POINTER :: cell
144 : TYPE(eeq_solver_type) :: eeq_sparam
145 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
146 :
147 : MARK_USED(print_level)
148 :
149 38 : CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, cell=cell)
150 38 : IF (ext) THEN
151 0 : NULLIFY (charges)
152 0 : CALL get_qs_env(qs_env, eeq=charges)
153 0 : CPASSERT(ASSOCIATED(charges))
154 0 : enshift_type = 0
155 : ELSE
156 114 : ALLOCATE (charges(natom))
157 : ! enforce en shift method 1 (original/molecular)
158 : ! method 2 from paper on PBC seems not to work
159 38 : enshift_type = 1
160 : !IF (ALL(cell%perd == 0)) enshift_type = 1
161 38 : CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
162 : END IF
163 :
164 38 : IF (iounit > 0) THEN
165 :
166 19 : IF (enshift_type == 0) THEN
167 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (External)"
168 19 : ELSE IF (enshift_type == 1) THEN
169 19 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Molecules))"
170 0 : ELSE IF (enshift_type == 2) THEN
171 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Crystals))"
172 : ELSE
173 0 : CPABORT("Unknown enshift_type")
174 : END IF
175 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
176 19 : "# Atom Element Kind Atomic Charge"
177 :
178 140 : DO iatom = 1, natom
179 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
180 : element_symbol=element_symbol, &
181 121 : kind_number=ikind)
182 : WRITE (UNIT=iounit, FMT="(T4,I8,T18,A2,I10,T43,F12.4)") &
183 140 : iatom, element_symbol, ikind, charges(iatom)
184 : END DO
185 :
186 : END IF
187 :
188 38 : IF (.NOT. ext) DEALLOCATE (charges)
189 :
190 38 : 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 90 : 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 90 : 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 90 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chia, cnumbers, efr, gam
221 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
222 90 : 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 90 : 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 90 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
232 90 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
233 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
234 : print_section
235 :
236 90 : 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 90 : dft_control=dft_control)
246 90 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
247 :
248 90 : logger => cp_get_default_logger()
249 90 : IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
250 37 : iunit = cp_logger_get_default_unit_nr()
251 : ELSE
252 53 : iunit = -1
253 : END IF
254 :
255 90 : totalcharge = dft_control%charge
256 :
257 90 : CALL get_cnumbers(qs_env, cnumbers, dcnum, .FALSE.)
258 :
259 : ! Apply smooth logistic CN cutoff to match multicharge/D4 behavior
260 90 : IF (PRESENT(cn_max)) THEN
261 472 : DO iatom = 1, natom
262 472 : 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 540 : ALLOCATE (gab(nkind, nkind), gam(nkind))
268 90 : gab = 0.0_dp
269 90 : gam = 0.0_dp
270 286 : DO ikind = 1, nkind
271 196 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
272 196 : CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
273 750 : DO jkind = 1, nkind
274 464 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
275 464 : CALL get_eeq_data(zb, eeq_model, rad=alb)
276 : !
277 660 : 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 90 : IF (PRESENT(exclude)) THEN
285 168 : DO ikind = 1, nkind
286 168 : IF (exclude(ikind)) THEN
287 8 : gam(ikind) = 1.0e30_dp
288 40 : gab(ikind, :) = 0.0_dp
289 40 : gab(:, ikind) = 0.0_dp
290 : END IF
291 : END DO
292 : END IF
293 :
294 : ! Chi[a,a]
295 90 : sgamma = 8.0_dp ! see D4 for periodic systems paper
296 90 : esg = 1.0_dp + EXP(sgamma)
297 270 : ALLOCATE (chia(natom))
298 90 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
299 752 : DO iatom = 1, natom
300 662 : ikind = kind_of(iatom)
301 662 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
302 662 : CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
303 : !
304 662 : IF (enshift_type == 1) THEN
305 662 : scn = cnumbers(iatom)/SQRT(cnumbers(iatom) + 1.0e-14_dp)
306 0 : ELSE IF (enshift_type == 2) THEN
307 0 : scn = LOG(esg/(esg - cnumbers(iatom)))
308 : ELSE
309 0 : CPABORT("Unknown enshift_type")
310 : END IF
311 1414 : chia(iatom) = xi - kappa*scn
312 : !
313 : END DO
314 :
315 : ! Zero electronegativity for excluded atoms (ghost/floating in BSSE)
316 90 : IF (PRESENT(exclude)) THEN
317 472 : DO iatom = 1, natom
318 420 : ikind = kind_of(iatom)
319 472 : IF (exclude(ikind)) chia(iatom) = 0.0_dp
320 : END DO
321 : END IF
322 :
323 : ! efield
324 90 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
325 : dft_control%apply_efield_field) THEN
326 0 : ALLOCATE (efr(natom))
327 0 : efr(1:natom) = 0.0_dp
328 0 : CALL eeq_efield_pot(qs_env, efr)
329 0 : chia(1:natom) = chia(1:natom) + efr(1:natom)
330 0 : DEALLOCATE (efr)
331 : END IF
332 :
333 90 : CALL cnumber_release(cnumbers, dcnum, .FALSE.)
334 :
335 90 : CALL get_cell(cell, periodic=periodic)
336 204 : do_ewald = .NOT. ALL(periodic == 0)
337 90 : IF (do_ewald) THEN
338 832 : ALLOCATE (ewald_env)
339 52 : CALL ewald_env_create(ewald_env, para_env)
340 52 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
341 52 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
342 52 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
343 52 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
344 52 : CALL get_qs_env(qs_env, cell_ref=cell_ref)
345 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
346 52 : silent=.TRUE., pset="EEQ", cell_periodic=cell%perd)
347 52 : ALLOCATE (ewald_pw)
348 52 : 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 52 : ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit, qs_env=qs_env)
355 : !
356 52 : CALL ewald_env_release(ewald_env)
357 52 : CALL ewald_pw_release(ewald_pw)
358 52 : 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 38 : totalcharge=totalcharge, iounit=iunit, qs_env=qs_env)
364 : END IF
365 :
366 90 : DEALLOCATE (gab, gam, chia)
367 :
368 90 : CALL timestop(handle)
369 :
370 180 : 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 14 : SUBROUTINE eeq_forces(qs_env, charges, dcharges, gradient, stress, &
387 14 : 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 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
405 : INTEGER, DIMENSION(3) :: periodic
406 : LOGICAL :: do_ewald, use_virial
407 14 : 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 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c_radius, cnumbers, gam, qlag
411 14 : 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 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: chrgx, dchia
415 14 : 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 14 : 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 14 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
427 14 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
428 : TYPE(mp_para_env_type), POINTER :: para_env
429 : TYPE(neighbor_list_iterator_p_type), &
430 14 : DIMENSION(:), POINTER :: nl_iterator
431 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
432 14 : POINTER :: sab_ew
433 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
434 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
435 14 : 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 14 : 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 14 : dft_control=dft_control)
453 :
454 14 : logger => cp_get_default_logger()
455 14 : IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
456 4 : iunit = cp_logger_get_default_unit_nr()
457 : ELSE
458 10 : iunit = -1
459 : END IF
460 :
461 14 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
462 14 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
463 :
464 14 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
465 :
466 14 : totalcharge = dft_control%charge
467 :
468 14 : 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 14 : IF (PRESENT(cn_max)) THEN
473 118 : DO iatom = 1, natom
474 104 : dcnpdcn = EXP(cn_max)/(EXP(cn_max) + EXP(cnumbers(iatom)))
475 104 : cnumbers(iatom) = LOG(1.0_dp + EXP(cn_max)) - LOG(1.0_dp + EXP(cn_max - cnumbers(iatom)))
476 270 : DO i = 1, dcnum(iatom)%neighbors
477 256 : dcnum(iatom)%dvals(i) = dcnum(iatom)%dvals(i)*dcnpdcn
478 : END DO
479 : END DO
480 : END IF
481 :
482 : ! gamma[a,b]
483 84 : ALLOCATE (gab(nkind, nkind), gam(nkind))
484 14 : gab = 0.0_dp
485 14 : gam = 0.0_dp
486 44 : DO ikind = 1, nkind
487 30 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
488 30 : CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
489 114 : DO jkind = 1, nkind
490 70 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
491 70 : CALL get_eeq_data(zb, eeq_model, rad=alb)
492 : !
493 100 : 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 14 : IF (PRESENT(exclude)) THEN
500 44 : DO ikind = 1, nkind
501 44 : IF (exclude(ikind)) THEN
502 0 : gam(ikind) = 1.0e30_dp
503 0 : gab(ikind, :) = 0.0_dp
504 0 : gab(:, ikind) = 0.0_dp
505 : END IF
506 : END DO
507 : END IF
508 :
509 42 : ALLOCATE (qlag(natom))
510 :
511 14 : CALL get_cell(cell, periodic=periodic)
512 26 : do_ewald = .NOT. ALL(periodic == 0)
513 14 : IF (do_ewald) THEN
514 160 : ALLOCATE (ewald_env)
515 10 : CALL ewald_env_create(ewald_env, para_env)
516 10 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
517 10 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
518 10 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
519 10 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
520 10 : CALL get_qs_env(qs_env, cell_ref=cell_ref)
521 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
522 10 : silent=.TRUE., pset="EEQ", cell_periodic=cell%perd)
523 10 : ALLOCATE (ewald_pw)
524 10 : 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 74 : 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 44 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit, qs_env=qs_env)
534 : END IF
535 :
536 14 : sgamma = 8.0_dp ! see D4 for periodic systems paper
537 14 : esg = 1.0_dp + EXP(sgamma)
538 56 : ALLOCATE (chrgx(natom), dchia(natom))
539 118 : DO iatom = 1, natom
540 104 : ikind = kind_of(iatom)
541 104 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
542 104 : CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
543 : !
544 104 : IF (response_only) THEN
545 104 : ctot = -0.5_dp*qlag(iatom)
546 : ELSE
547 0 : ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
548 : END IF
549 222 : IF (enshift_type == 1) THEN
550 104 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
551 104 : dchia(iatom) = -ctot*kappa/scn
552 0 : ELSE IF (enshift_type == 2) THEN
553 0 : cn = cnumbers(iatom)
554 0 : scn = 1.0_dp/(esg - cn)
555 0 : dchia(iatom) = -ctot*kappa*scn
556 : ELSE
557 0 : CPABORT("Unknown enshift_type")
558 : END IF
559 : END DO
560 :
561 : ! Efield
562 14 : IF (dft_control%apply_period_efield) THEN
563 0 : CALL eeq_efield_force_periodic(qs_env, charges, qlag)
564 14 : ELSE IF (dft_control%apply_efield) THEN
565 0 : CALL eeq_efield_force_loc(qs_env, charges, qlag)
566 14 : ELSE IF (dft_control%apply_efield_field) THEN
567 0 : CPABORT("apply field")
568 : END IF
569 :
570 : ! Forces from q*X
571 14 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
572 44 : DO ikind = 1, nkind
573 96 : DO ia = 1, local_particles%n_el(ikind)
574 52 : iatom = local_particles%list(ikind)%array(ia)
575 158 : DO i = 1, dcnum(iatom)%neighbors
576 76 : katom = dcnum(iatom)%nlist(i)
577 304 : rik = dcnum(iatom)%rik(:, i)
578 304 : drk = NORM2(rik)
579 128 : IF (drk > 1.e-3_dp) THEN
580 304 : fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
581 304 : gradient(:, iatom) = gradient(:, iatom) - fdik(:)
582 304 : gradient(:, katom) = gradient(:, katom) + fdik(:)
583 76 : IF (use_virial) THEN
584 28 : 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 14 : 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 10 : molecule_set=molecule_set)
599 10 : subcells = 2.0_dp
600 10 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
601 10 : rcut = 2.0_dp*rcut
602 10 : NULLIFY (sab_ew)
603 80 : ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
604 32 : c_radius(:) = rcut
605 32 : default_present = .TRUE.
606 52 : ALLOCATE (atom2d(nkind))
607 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
608 10 : molecule_set, .FALSE., particle_set=particle_set)
609 10 : 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 10 : subcells=subcells, operator_type="PP", nlname="sab_ew")
612 10 : DEALLOCATE (c_radius, pair_radius, default_present)
613 10 : CALL atom2d_cleanup(atom2d)
614 : !
615 10 : CALL neighbor_list_iterator_create(nl_iterator, sab_ew)
616 103703 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
617 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
618 103693 : iatom=iatom, jatom=jatom, r=rij)
619 : !
620 414772 : dr2 = SUM(rij**2)
621 103693 : dr = SQRT(dr2)
622 103693 : IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
623 12994 : fe = 1.0_dp
624 12994 : IF (iatom == jatom) fe = 0.5_dp
625 12994 : IF (response_only) THEN
626 12994 : qq = -qlag(iatom)*charges(jatom)
627 : ELSE
628 12994 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
629 : END IF
630 12994 : gama = gab(ikind, jkind)
631 12994 : gam2 = gama*gama
632 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
633 12994 : - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
634 12994 : IF (response_only) THEN
635 12994 : qq1 = -qlag(iatom)*charges(jatom)
636 12994 : qq2 = -qlag(jatom)*charges(iatom)
637 : ELSE
638 0 : qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
639 0 : qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
640 : END IF
641 51976 : fdik(:) = -qq1*grc*rij(:)/dr
642 51976 : gradient(:, iatom) = gradient(:, iatom) + fdik(:)
643 51976 : gradient(:, jatom) = gradient(:, jatom) - fdik(:)
644 12994 : IF (use_virial) THEN
645 8515 : CALL virial_pair_force(stress, -fe, fdik, rij)
646 : END IF
647 51976 : fdik(:) = qq2*grc*rij(:)/dr
648 51976 : gradient(:, iatom) = gradient(:, iatom) - fdik(:)
649 51976 : gradient(:, jatom) = gradient(:, jatom) + fdik(:)
650 13004 : IF (use_virial) THEN
651 8515 : CALL virial_pair_force(stress, fe, fdik, rij)
652 : END IF
653 : END DO
654 10 : CALL neighbor_list_iterator_release(nl_iterator)
655 : !
656 10 : CALL release_neighbor_list_sets(sab_ew)
657 : ELSE
658 12 : DO ikind = 1, nkind
659 32 : DO ia = 1, local_particles%n_el(ikind)
660 20 : iatom = local_particles%list(ikind)%array(ia)
661 80 : ri(1:3) = particle_set(iatom)%r(1:3)
662 228 : DO jatom = 1, natom
663 200 : IF (iatom == jatom) CYCLE
664 180 : jkind = kind_of(jatom)
665 180 : IF (response_only) THEN
666 180 : qq = -qlag(iatom)*charges(jatom)
667 : ELSE
668 0 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
669 : END IF
670 720 : rj(1:3) = particle_set(jatom)%r(1:3)
671 720 : rij(1:3) = ri(1:3) - rj(1:3)
672 720 : rij = pbc(rij, cell)
673 720 : dr2 = SUM(rij**2)
674 180 : dr = SQRT(dr2)
675 180 : gama = gab(ikind, jkind)
676 180 : gam2 = gama*gama
677 180 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
678 720 : fdik(:) = qq*grc*rij(:)/dr
679 720 : gradient(:, iatom) = gradient(:, iatom) + fdik(:)
680 740 : 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 14 : IF (do_ewald) THEN
688 30 : ALLOCATE (epforce(3, natom))
689 10 : epforce = 0.0_dp
690 10 : IF (response_only) THEN
691 74 : dchia(1:natom) = qlag(1:natom)
692 : ELSE
693 0 : dchia(1:natom) = -charges(1:natom) + qlag(1:natom)
694 : END IF
695 74 : chrgx(1:natom) = charges(1:natom)
696 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
697 10 : particle_set, dchia, epforce)
698 74 : dchia(1:natom) = charges(1:natom)
699 74 : chrgx(1:natom) = qlag(1:natom)
700 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
701 10 : particle_set, dchia, epforce)
702 266 : gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + epforce(1:3, 1:natom)
703 10 : DEALLOCATE (epforce)
704 :
705 : ! virial
706 10 : IF (use_virial) THEN
707 52 : chrgx(1:natom) = charges(1:natom) - qlag(1:natom)
708 8 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
709 104 : stress = stress - pvir
710 52 : chrgx(1:natom) = qlag(1:natom)
711 8 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
712 104 : stress = stress + pvir
713 8 : IF (response_only) THEN
714 52 : chrgx(1:natom) = charges(1:natom)
715 8 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
716 104 : stress = stress + pvir
717 : END IF
718 : END IF
719 : !
720 10 : CALL ewald_env_release(ewald_env)
721 10 : CALL ewald_pw_release(ewald_pw)
722 10 : DEALLOCATE (ewald_env, ewald_pw)
723 : END IF
724 :
725 14 : CALL cnumber_release(cnumbers, dcnum, .TRUE.)
726 :
727 14 : DEALLOCATE (gab, gam, qlag, chrgx, dchia)
728 :
729 14 : CALL timestop(handle)
730 :
731 28 : END SUBROUTINE eeq_forces
732 :
733 : ! **************************************************************************************************
734 : !> \brief ...
735 : !> \param qs_env ...
736 : !> \param cnumbers ...
737 : !> \param dcnum ...
738 : !> \param calculate_forces ...
739 : ! **************************************************************************************************
740 104 : 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 104 : 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 104 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
757 104 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
758 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
759 104 : POINTER :: sab_cn
760 104 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
761 : TYPE(qs_dispersion_type), POINTER :: disp
762 104 : 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 104 : molecule_set=molecule_set)
772 104 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
773 :
774 : ! Check for dispersion_env and sab_cn needed for cnumbers
775 520 : ALLOCATE (disp)
776 104 : disp%k1 = 16.0_dp
777 104 : disp%k2 = 4._dp/3._dp
778 104 : disp%eps_cn = 1.E-6_dp
779 104 : disp%max_elem = maxElem
780 104 : ALLOCATE (disp%rcov(maxElem))
781 9048 : disp%rcov(1:maxElem) = bohr*disp%k2*rcov(1:maxElem)
782 104 : subcells = 2.0_dp
783 : ! Build the neighbor lists for the CN
784 104 : NULLIFY (sab_cn)
785 832 : ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
786 104 : c_radius(:) = 0.0_dp
787 330 : default_present = .TRUE.
788 330 : DO ikind = 1, nkind
789 226 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
790 330 : c_radius(ikind) = 4._dp*rcov(za)*bohr
791 : END DO
792 538 : ALLOCATE (atom2d(nkind))
793 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
794 104 : molecule_set, .FALSE., particle_set=particle_set)
795 104 : 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 104 : subcells=subcells, operator_type="PP", nlname="sab_cn")
798 104 : disp%sab_cn => sab_cn
799 104 : DEALLOCATE (c_radius, pair_radius, default_present)
800 104 : CALL atom2d_cleanup(atom2d)
801 :
802 : ! Calculate coordination numbers
803 104 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces, disp_env=disp)
804 :
805 104 : CALL qs_dispersion_release(disp)
806 :
807 208 : 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 2334 : SUBROUTINE eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, &
832 4668 : 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 2334 : CALL timeset(routineN, handle)
864 :
865 2334 : do_ewald = .FALSE.
866 2334 : IF (PRESENT(ewald)) do_ewald = ewald
867 : !
868 2334 : qtot = 0.0_dp
869 2334 : IF (PRESENT(totalcharge)) qtot = totalcharge
870 : !
871 2334 : iunit = -1
872 2334 : IF (PRESENT(iounit)) iunit = iounit
873 :
874 : ! EEQ solver parameters
875 2334 : do_direct = eeq_sparam%direct
876 2334 : natom = SIZE(particle_set)
877 2334 : do_displ = .FALSE.
878 2334 : IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
879 200 : do_displ = dft_control%period_efield%displacement_field
880 : END IF
881 2334 : 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 2334 : do_sparse = eeq_sparam%sparse .OR. (do_sparse_auto .AND. .NOT. do_displ)
886 2334 : fm_created = .FALSE.
887 :
888 2334 : nkind = SIZE(gam)
889 2334 : ns = natom + 1
890 2334 : 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 2222 : nrow_global=ns, ncol_global=ns)
893 2222 : CALL cp_fm_create(eeq_mat, mat_struct)
894 2222 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
895 2222 : fm_created = .TRUE.
896 : END IF
897 : !
898 2334 : IF (do_ewald) THEN
899 1126 : CPASSERT(PRESENT(ewald_env))
900 1126 : CPASSERT(PRESENT(ewald_pw))
901 1126 : IF (do_direct) THEN
902 0 : IF (do_displ) THEN
903 0 : 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 0 : ewald_env, ewald_pw, iounit)
908 : END IF
909 1126 : ELSE IF (do_sparse) THEN
910 112 : IF (do_displ) THEN
911 0 : CPABORT("Sparse periodic EEQ is not available with a displacement field")
912 : ELSE
913 112 : 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 112 : ewald_env, ewald_pw, eeq_sparam, qs_env, ierror, iounit)
918 112 : IF (ierror /= 0) THEN
919 : CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
920 0 : nrow_global=ns, ncol_global=ns)
921 0 : CALL cp_fm_create(eeq_mat, mat_struct)
922 0 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
923 0 : fm_created = .TRUE.
924 : CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
925 : kind_of, cell, chia, gam, gab, qtot, &
926 0 : ewald_env, ewald_pw, iounit)
927 : END IF
928 : END IF
929 : ELSE
930 1014 : IF (do_displ) THEN
931 0 : 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 1014 : ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
937 1014 : 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 892 : ewald_env, ewald_pw, iounit)
942 : END IF
943 : END IF
944 : END IF
945 1126 : IF (qtot /= 0._dp) THEN
946 104 : CALL get_cell(cell=cell, deth=deth)
947 104 : CALL ewald_env_get(ewald_env, alpha=alpha)
948 104 : 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 1208 : cell, chia, gam, gab, qtot, ftime)
953 1208 : IF (iounit > 0) THEN
954 1183 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
955 : END IF
956 : END IF
957 2334 : IF (fm_created) THEN
958 2222 : CALL cp_fm_struct_release(mat_struct)
959 2222 : CALL cp_fm_release(eeq_mat)
960 : END IF
961 :
962 2334 : CALL timestop(handle)
963 :
964 2334 : 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 2222 : SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
982 2222 : 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 2222 : 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 2222 : CALL timeset(routineN, handle)
1008 2222 : ti = m_walltime()
1009 :
1010 2222 : natom = SIZE(particle_set)
1011 2222 : nkind = SIZE(gam)
1012 : !
1013 2222 : ns = natom + 1
1014 2222 : 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 2222 : row_indices=rind, col_indices=cind)
1017 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1018 2222 : nrow_global=ns, ncol_global=1)
1019 2222 : CALL cp_fm_create(rhs_vec, vec_struct)
1020 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1021 2222 : row_indices=rvind, col_indices=cvind)
1022 : !
1023 : ! set up matrix
1024 2222 : CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
1025 2222 : CALL cp_fm_set_all(rhs_vec, 0.0_dp)
1026 8550 : DO ir = 1, nrloc
1027 6328 : iar = rind(ir)
1028 6328 : IF (iar > natom) CYCLE
1029 5217 : ikind = kind_of(iar)
1030 20868 : ri(1:3) = particle_set(iar)%r(1:3)
1031 46845 : DO ic = 1, ncloc
1032 39406 : iac = cind(ic)
1033 39406 : IF (iac > natom) CYCLE
1034 34189 : jkind = kind_of(iac)
1035 136756 : rj(1:3) = particle_set(iac)%r(1:3)
1036 34189 : IF (iar == iac) THEN
1037 5217 : grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1038 : ELSE
1039 115888 : rij(1:3) = ri(1:3) - rj(1:3)
1040 115888 : rij = pbc(rij, cell)
1041 115888 : dr = NORM2(rij)
1042 28972 : grc = erf(gab(ikind, jkind)*dr)/dr
1043 : END IF
1044 45734 : eeq_mat%local_data(ir, ic) = grc
1045 : END DO
1046 : END DO
1047 : ! set up rhs vector
1048 8550 : DO ir = 1, nrvloc
1049 6328 : iar = rvind(ir)
1050 14878 : DO ic = 1, ncvloc
1051 6328 : iac = cvind(ic)
1052 6328 : ia = MAX(iar, iac)
1053 6328 : IF (ia > natom) THEN
1054 1111 : xr = qtot
1055 : ELSE
1056 5217 : xr = -chia(ia)
1057 : END IF
1058 12656 : rhs_vec%local_data(ir, ic) = xr
1059 : END DO
1060 : END DO
1061 : !
1062 2222 : CALL cp_fm_solve(eeq_mat, rhs_vec)
1063 : !
1064 12656 : charges = 0.0_dp
1065 2222 : lambda = 0.0_dp
1066 8550 : DO ir = 1, nrvloc
1067 6328 : iar = rvind(ir)
1068 14878 : DO ic = 1, ncvloc
1069 6328 : iac = cvind(ic)
1070 6328 : ia = MAX(iar, iac)
1071 12656 : IF (ia <= natom) THEN
1072 5217 : xr = rhs_vec%local_data(ir, ic)
1073 5217 : charges(ia) = xr
1074 : ELSE
1075 1111 : lambda = rhs_vec%local_data(ir, ic)
1076 : END IF
1077 : END DO
1078 : END DO
1079 2222 : CALL para_env%sum(lambda)
1080 23090 : CALL para_env%sum(charges)
1081 : !
1082 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
1083 12656 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1084 :
1085 2222 : CALL cp_fm_struct_release(vec_struct)
1086 2222 : CALL cp_fm_release(rhs_vec)
1087 :
1088 2222 : te = m_walltime()
1089 2222 : ftime = te - ti
1090 2222 : CALL timestop(handle)
1091 :
1092 2222 : 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 112 : SUBROUTINE pbc_sparse_solver(charges, lambda, eeq_energy, particle_set, &
1118 112 : 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 112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ap, aq, gradient, p, precond_diag, &
1144 : residual, work, z
1145 112 : TYPE(eeq_sparse_matrix_type) :: response
1146 :
1147 112 : CALL timeset(routineN, handle)
1148 112 : ti = m_walltime()
1149 :
1150 112 : iunit = -1
1151 112 : IF (PRESENT(iounit)) iunit = iounit
1152 112 : natom = SIZE(charges)
1153 112 : max_iter = eeq_sparam%max_diis
1154 112 : eps_solver = eeq_sparam%eps_diis
1155 112 : ierror = 0
1156 :
1157 : CALL build_sparse_pbc_response(response, qs_env, particle_set, kind_of, cell, gam, gab, &
1158 112 : ewald_env, eps_solver, iunit)
1159 :
1160 : ALLOCATE (ap(natom), aq(natom), gradient(natom), p(natom), precond_diag(natom), &
1161 1120 : residual(natom), work(natom), z(natom))
1162 968 : DO i = 1, natom
1163 856 : precond_diag(i) = gam(kind_of(i)) + 2.0_dp*gab(kind_of(i), kind_of(i))*oorootpi
1164 968 : 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 1824 : lambda = (-qtot - SUM(chia/precond_diag))/SUM(1.0_dp/precond_diag)
1169 968 : charges = -(chia + lambda)/precond_diag
1170 1824 : 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 112 : charges, aq, work)
1174 968 : gradient(:) = aq + chia
1175 968 : residual(:) = -gradient
1176 112 : CALL project_eeq_charge_tangent(residual)
1177 968 : res = NORM2(residual)
1178 112 : res_initial = MAX(res, SQRT(EPSILON(1.0_dp)))
1179 112 : iter = 0
1180 :
1181 112 : IF (res >= eps_solver) THEN
1182 112 : CALL apply_eeq_constrained_diagonal_preconditioner(residual, precond_diag, z)
1183 968 : p(:) = z
1184 968 : rz = SUM(residual*z)
1185 :
1186 632 : DO iter = 1, max_iter
1187 : CALL apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1188 632 : p, ap, work)
1189 632 : CALL project_eeq_charge_tangent(ap)
1190 5608 : denom = SUM(p*ap)
1191 632 : IF (.NOT. (denom > SQRT(TINY(1.0_dp)) .AND. ABS(denom) < HUGE(denom))) THEN
1192 0 : ierror = 1
1193 0 : EXIT
1194 : END IF
1195 :
1196 632 : alpha_cg = rz/denom
1197 632 : IF (.NOT. (ABS(alpha_cg) < HUGE(alpha_cg))) THEN
1198 0 : ierror = 1
1199 0 : EXIT
1200 : END IF
1201 5608 : charges = charges + alpha_cg*p
1202 5608 : residual(:) = residual - alpha_cg*ap
1203 : ! Remove accumulated roundoff from both constrained spaces.
1204 10584 : charges = charges + (qtot - SUM(charges))/REAL(natom, KIND=dp)
1205 632 : CALL project_eeq_charge_tangent(residual)
1206 5608 : res = NORM2(residual)
1207 632 : IF (res < eps_solver) EXIT
1208 520 : IF (res > 100.0_dp*res_initial .OR. .NOT. (res < HUGE(res))) THEN
1209 0 : ierror = 1
1210 0 : EXIT
1211 : END IF
1212 :
1213 520 : CALL apply_eeq_constrained_diagonal_preconditioner(residual, precond_diag, z)
1214 4640 : rz_new = SUM(residual*z)
1215 520 : IF (.NOT. (rz_new > 0.0_dp .AND. ABS(rz_new) < HUGE(rz_new))) THEN
1216 0 : ierror = 1
1217 0 : EXIT
1218 : END IF
1219 520 : beta_cg = rz_new/rz
1220 4640 : p(:) = z + beta_cg*p
1221 632 : rz = rz_new
1222 : END DO
1223 112 : IF (iter > max_iter) ierror = 1
1224 : END IF
1225 :
1226 112 : IF (ierror == 0) THEN
1227 : CALL apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1228 112 : charges, aq, work)
1229 968 : gradient(:) = aq + chia
1230 968 : lambda = -SUM(gradient)/REAL(natom, KIND=dp)
1231 1824 : eeq_energy = 0.5_dp*SUM(charges*aq) + SUM(charges*chia)
1232 : END IF
1233 :
1234 112 : te = m_walltime()
1235 112 : ftime = te - ti
1236 112 : IF (iunit > 0) THEN
1237 107 : IF (ierror /= 0) THEN
1238 : WRITE (iunit, '(A,T50,I6,T61,E20.5)') &
1239 0 : " EEQ| Sparse PBC solver failed: iterations/accuracy ", iter, res
1240 : ELSE
1241 : WRITE (iunit, '(A,T50,I6,T61,E20.5)') &
1242 107 : " EEQ| Sparse PBC solver: iterations/accuracy ", iter, res
1243 : END IF
1244 107 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Sparse PBC solver: time[s]", ftime
1245 : END IF
1246 :
1247 112 : DEALLOCATE (ap, aq, gradient, p, precond_diag, residual, work, z)
1248 112 : CALL release_sparse_pbc_response(response)
1249 112 : CALL timestop(handle)
1250 :
1251 224 : 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 224 : 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 112 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: pair_key
1286 112 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pair_order
1287 : INTEGER, DIMENSION(3) :: periodic
1288 112 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
1289 : REAL(KIND=dp) :: drop_tolerance, neighbor_cutoff, rmax, &
1290 : subcells, value
1291 112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c_radius
1292 112 : 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 112 : 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 112 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
1299 112 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1300 : TYPE(mp_para_env_type), POINTER :: para_env
1301 : TYPE(neighbor_list_iterator_p_type), &
1302 112 : DIMENSION(:), POINTER :: nl_iterator
1303 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1304 112 : POINTER :: sab_eeq
1305 :
1306 112 : NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, nl_iterator, sab_eeq)
1307 112 : CALL ewald_env_get(ewald_env, rcut=rmax, para_env=para_env)
1308 112 : rmax = 2.0_dp*rmax
1309 112 : natom = SIZE(particle_set)
1310 : drop_tolerance = MAX(100.0_dp*EPSILON(1.0_dp), &
1311 112 : 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 112 : CALL get_cell(cell, h=hmat, periodic=periodic)
1317 448 : WHERE (periodic /= 0)
1318 : cell_extent = 0.5_dp
1319 : ELSE WHERE
1320 : cell_extent = 1.0_dp
1321 : END WHERE
1322 112 : neighbor_cutoff = 0.0_dp
1323 336 : DO ix = -1, 1, 2
1324 784 : DO iy = -1, 1, 2
1325 672 : DO iz = -1, 1, 2
1326 14336 : vertex = MATMUL(hmat, cell_extent*REAL([ix, iy, iz], KIND=dp))
1327 3584 : neighbor_cutoff = MAX(neighbor_cutoff, NORM2(vertex))
1328 : END DO
1329 : END DO
1330 : END DO
1331 : neighbor_cutoff = MIN(rmax, neighbor_cutoff + &
1332 112 : 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 112 : 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 1354 : atom2d(SIZE(atomic_kind_set)))
1340 346 : c_radius = 0.5_dp*neighbor_cutoff
1341 346 : default_present = .TRUE.
1342 112 : subcells = 2.0_dp
1343 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
1344 112 : molecule_set, .FALSE., particle_set=particle_set)
1345 112 : 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 112 : operator_type="PP", nlname="sab_eeq_sparse")
1349 112 : DEALLOCATE (c_radius, pair_radius, default_present)
1350 112 : CALL atom2d_cleanup(atom2d)
1351 :
1352 1192 : ALLOCATE (response%diag(natom), source=0.0_dp)
1353 968 : DO iatom = 1, natom
1354 856 : owner = MOD(iatom - 1, para_env%num_pe)
1355 968 : IF (owner == para_env%mepos) THEN
1356 : response%diag(iatom) = eeq_short_range_element(iatom, iatom, particle_set, kind_of, &
1357 428 : 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 112 : n_pairs_alloc = 0
1365 112 : CALL neighbor_list_iterator_create(nl_iterator, sab_eeq)
1366 6308 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1367 6196 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
1368 6308 : IF (iatom /= jatom) n_pairs_alloc = n_pairs_alloc + 1
1369 : END DO
1370 112 : CALL neighbor_list_iterator_release(nl_iterator)
1371 :
1372 560 : ALLOCATE (pair_key(n_pairs_alloc), pair_order(n_pairs_alloc))
1373 112 : n_pairs_alloc = 0
1374 112 : CALL neighbor_list_iterator_create(nl_iterator, sab_eeq)
1375 6308 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1376 6196 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
1377 6196 : IF (iatom == jatom) CYCLE
1378 5768 : 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 6308 : INT(MAX(iatom, jatom), int_8)
1381 : END DO
1382 112 : CALL neighbor_list_iterator_release(nl_iterator)
1383 112 : CALL release_neighbor_list_sets(sab_eeq)
1384 :
1385 112 : IF (n_pairs_alloc > 0) CALL sort(pair_key, n_pairs_alloc, pair_order)
1386 112 : n_unique = 0
1387 5880 : DO iatom = 1, n_pairs_alloc
1388 5880 : IF (iatom == 1) THEN
1389 112 : n_unique = n_unique + 1
1390 5656 : ELSE IF (pair_key(iatom) /= pair_key(iatom - 1)) THEN
1391 1346 : n_unique = n_unique + 1
1392 : END IF
1393 : END DO
1394 672 : ALLOCATE (response%row(n_unique), response%col(n_unique), response%value(n_unique))
1395 112 : n_pairs = 0
1396 5880 : DO iatom = 1, n_pairs_alloc
1397 5768 : IF (iatom > 1) THEN
1398 5656 : IF (pair_key(iatom) == pair_key(iatom - 1)) CYCLE
1399 : END IF
1400 1458 : jatom = INT(MODULO(pair_key(iatom) - 1_int_8, INT(natom, int_8))) + 1
1401 1458 : 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 1458 : gam, gab, ewald_env)
1404 1458 : IF (ABS(value) <= drop_tolerance) CYCLE
1405 1458 : n_pairs = n_pairs + 1
1406 1458 : response%row(n_pairs) = owner
1407 1458 : response%col(n_pairs) = jatom
1408 5880 : response%value(n_pairs) = value
1409 : END DO
1410 112 : DEALLOCATE (pair_key, pair_order)
1411 112 : response%nactive = n_pairs
1412 :
1413 112 : n_pairs_global = n_pairs
1414 112 : CALL para_env%sum(n_pairs_global)
1415 112 : n_candidates_global = n_unique
1416 112 : CALL para_env%sum(n_candidates_global)
1417 112 : IF (iunit > 0) THEN
1418 107 : WRITE (iunit, '(A,T54,I12)') " EEQ| Sparse short-range candidate atom pairs", n_candidates_global
1419 107 : WRITE (iunit, '(A,T54,I12)') " EEQ| Sparse short-range retained atom pairs", n_pairs_global
1420 107 : WRITE (iunit, '(A,T61,E20.5)') " EEQ| Sparse short-range drop tolerance", drop_tolerance
1421 107 : WRITE (iunit, '(A,T61,E20.5)') " EEQ| Sparse short-range MIC cutoff", neighbor_cutoff
1422 : END IF
1423 :
1424 224 : 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 1886 : 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 1886 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
1457 1886 : rmax = 2.0_dp*rcut
1458 1886 : CALL get_cell(cell, h=hmat, periodic=periodic)
1459 1886 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1460 1886 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1461 1886 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1462 7544 : WHERE (periodic == 0) ncell = 0
1463 :
1464 7544 : rij = pbc(particle_set(iatom)%r - particle_set(jatom)%r, cell)
1465 1886 : element = 0.0_dp
1466 15148 : DO ix = -ncell(1), ncell(1)
1467 108522 : DO iy = -ncell(2), ncell(2)
1468 765114 : DO iz = -ncell(3), ncell(3)
1469 2633912 : cvec = [ix, iy, iz]
1470 12511082 : rijl = rij + MATMUL(hmat, cvec)
1471 2633912 : dr = NORM2(rijl)
1472 658478 : IF (dr > rmax) CYCLE
1473 252838 : 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 428 : 2.0_dp*alpha*oorootpi
1477 159036 : ELSE IF (dr > 1.0E-12_dp) THEN
1478 : element = element + &
1479 159036 : (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 1886 : 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 2568 : SUBROUTINE apply_sparse_pbc_response(response, ewald_env, ewald_pw, cell, particle_set, &
1499 856 : 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 856 : CALL ewald_env_get(ewald_env, para_env=para_env)
1513 7544 : potential = response%diag*charges
1514 12400 : DO i = 1, response%nactive
1515 11544 : ia = response%row(i)
1516 11544 : ja = response%col(i)
1517 11544 : potential(ia) = potential(ia) + response%value(i)*charges(ja)
1518 12400 : potential(ja) = potential(ja) + response%value(i)*charges(ia)
1519 : END DO
1520 14232 : CALL para_env%sum(potential)
1521 :
1522 856 : CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, work)
1523 7544 : potential = potential + work
1524 :
1525 856 : 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 632 : 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 5608 : RESULT = residual/diagonal
1541 11216 : constraint_shift = SUM(RESULT)/SUM(1.0_dp/diagonal)
1542 5608 : RESULT = RESULT - constraint_shift/diagonal
1543 :
1544 632 : 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 1376 : SUBROUTINE project_eeq_charge_tangent(vector)
1551 :
1552 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vector
1553 :
1554 22992 : vector = vector - SUM(vector)/REAL(SIZE(vector), KIND=dp)
1555 :
1556 1376 : END SUBROUTINE project_eeq_charge_tangent
1557 :
1558 : ! **************************************************************************************************
1559 : !> \brief Release a sparse periodic EEQ response.
1560 : !> \param response ...
1561 : ! **************************************************************************************************
1562 112 : SUBROUTINE release_sparse_pbc_response(response)
1563 :
1564 : TYPE(eeq_sparse_matrix_type), INTENT(INOUT) :: response
1565 :
1566 112 : IF (ALLOCATED(response%col)) DEALLOCATE (response%col)
1567 112 : IF (ALLOCATED(response%row)) DEALLOCATE (response%row)
1568 112 : IF (ALLOCATED(response%diag)) DEALLOCATE (response%diag)
1569 112 : IF (ALLOCATED(response%value)) DEALLOCATE (response%value)
1570 112 : response%nactive = 0
1571 :
1572 112 : 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 1014 : SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1594 1014 : 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 1014 : 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 1014 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1622 1014 : 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 1014 : 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 1014 : CALL timeset(routineN, handle)
1631 1014 : ti = m_walltime()
1632 :
1633 1014 : ierror = 0
1634 :
1635 1014 : iunit = -1
1636 1014 : IF (PRESENT(iounit)) iunit = iounit
1637 :
1638 1014 : natom = SIZE(particle_set)
1639 1014 : nkind = SIZE(gam)
1640 : !
1641 1014 : CALL get_cell(cell=cell, deth=deth)
1642 1014 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1643 1014 : ad = 2.0_dp*alpha*oorootpi
1644 1014 : IF (ewald_type /= do_ewald_spme) THEN
1645 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1646 : END IF
1647 : !
1648 1014 : rmax = 2.0_dp*rcut
1649 : ! max cells used
1650 1014 : CALL get_cell(cell, h=hmat, periodic=periodic)
1651 1014 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1652 1014 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1653 1014 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1654 1014 : IF (periodic(1) == 0) ncell(1) = 0
1655 1014 : IF (periodic(2) == 0) ncell(2) = 0
1656 1014 : IF (periodic(3) == 0) ncell(3) = 0
1657 : !
1658 : CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1659 1014 : chia, gam, gab, qtot, ftime)
1660 1014 : IF (iunit > 0) THEN
1661 981 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1662 : END IF
1663 1014 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1664 1014 : CALL cp_fm_create(pmat, mat_struct)
1665 1014 : 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 1014 : row_indices=rind, col_indices=cind)
1670 1014 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1671 4298 : DO ir = 1, nrloc
1672 3284 : iar = rind(ir)
1673 3284 : ri = 0.0_dp
1674 3284 : IF (iar <= natom) THEN
1675 2777 : ikind = kind_of(iar)
1676 11108 : ri(1:3) = particle_set(iar)%r(1:3)
1677 : END IF
1678 31960 : DO ic = 1, ncloc
1679 27662 : iac = cind(ic)
1680 27662 : IF (iac > natom .AND. iar > natom) THEN
1681 507 : eeq_mat%local_data(ir, ic) = 0.0_dp
1682 507 : CYCLE
1683 27155 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1684 5554 : eeq_mat%local_data(ir, ic) = 1.0_dp
1685 5554 : CYCLE
1686 : END IF
1687 21601 : jkind = kind_of(iac)
1688 86404 : rj(1:3) = particle_set(iac)%r(1:3)
1689 86404 : rij(1:3) = ri(1:3) - rj(1:3)
1690 86404 : rij = pbc(rij, cell)
1691 192202 : DO ix = -ncell(1), ncell(1)
1692 1128591 : DO iy = -ncell(2), ncell(2)
1693 7662099 : DO iz = -ncell(3), ncell(3)
1694 26220436 : cvec = [ix, iy, iz]
1695 124547071 : rijl = rij + MATMUL(hmat, cvec)
1696 26220436 : dr = NORM2(rijl)
1697 6555109 : IF (dr > rmax) CYCLE
1698 1615133 : IF (iar == iac .AND. dr < 0.00001_dp) THEN
1699 2777 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1700 : ELSE
1701 1612356 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1702 : END IF
1703 2554806 : 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 1014 : row_indices=rind, col_indices=cind)
1713 1014 : CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1714 4298 : DO ir = 1, nrloc
1715 3284 : iar = rind(ir)
1716 3284 : ri = 0.0_dp
1717 3284 : IF (iar <= natom) THEN
1718 2777 : ikind = kind_of(iar)
1719 11108 : ri(1:3) = particle_set(iar)%r(1:3)
1720 : END IF
1721 31960 : DO ic = 1, ncloc
1722 27662 : iac = cind(ic)
1723 27662 : IF (iac > natom .AND. iar > natom) THEN
1724 507 : pmat%local_data(ir, ic) = 0.0_dp
1725 507 : CYCLE
1726 27155 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1727 5554 : pmat%local_data(ir, ic) = 1.0_dp
1728 5554 : CYCLE
1729 : END IF
1730 21601 : jkind = kind_of(iac)
1731 86404 : rj(1:3) = particle_set(iac)%r(1:3)
1732 86404 : rij(1:3) = ri(1:3) - rj(1:3)
1733 86404 : rij = pbc(rij, cell)
1734 21601 : IF (iar == iac) THEN
1735 2777 : grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1736 : ELSE
1737 18824 : grc2 = erf(gab(ikind, jkind)*dr)/dr
1738 : END IF
1739 24885 : pmat%local_data(ir, ic) = grc2
1740 : END DO
1741 : END DO
1742 1014 : CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1743 : ! preconditioner invers
1744 1014 : CALL cp_fm_invert(pmat, mmat)
1745 : !
1746 : ! rhs
1747 1014 : ns = natom + 1
1748 3042 : ALLOCATE (rhs(ns))
1749 6568 : rhs(1:natom) = chia(1:natom)
1750 1014 : rhs(ns) = -qtot
1751 : !
1752 3042 : ALLOCATE (xv0(ns), rv0(ns))
1753 : ! initial guess
1754 6568 : xv0(1:natom) = charges(1:natom)
1755 1014 : xv0(ns) = 0.0_dp
1756 : ! DIIS optimizer
1757 1014 : max_diis = eeq_sparam%max_diis
1758 1014 : mdiis = eeq_sparam%mdiis
1759 1014 : sdiis = eeq_sparam%sdiis
1760 1014 : eps_diis = eeq_sparam%eps_diis
1761 1014 : astep = eeq_sparam%alpha
1762 7098 : ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1763 1014 : xvec = 0.0_dp; fvec = 0.0_dp
1764 9126 : ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1765 1014 : dmat = 0.0_dp; dvec = 0.0_dp
1766 1014 : ndiis = 1
1767 1014 : now = 1
1768 : CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1769 1014 : cell, particle_set, xv0, rhs, rv0)
1770 7582 : resin = NORM2(rv0)
1771 : !
1772 9864 : DO iv = 1, max_diis
1773 85656 : res = NORM2(rv0)
1774 9864 : IF (res > 10._dp*resin) EXIT
1775 8972 : IF (res < eps_diis) EXIT
1776 : !
1777 8850 : now = MOD(iv - 1, mdiis) + 1
1778 8850 : ndiis = MIN(iv, mdiis)
1779 78074 : xvec(1:ns, now) = xv0(1:ns)
1780 78074 : fvec(1:ns, now) = rv0(1:ns)
1781 65160 : DO i = 1, ndiis
1782 547592 : vmat(now, i) = SUM(fvec(:, now)*fvec(:, i))
1783 65160 : vmat(i, now) = vmat(now, i)
1784 : END DO
1785 8850 : IF (ndiis < sdiis) THEN
1786 26956 : xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1787 : ELSE
1788 6902 : dvec = 0.0_dp
1789 6902 : dvec(ndiis + 1) = 1.0_dp
1790 562846 : dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1791 60290 : dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1792 60290 : dmat(1:ndiis, ndiis + 1) = 1.0_dp
1793 6902 : dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1794 6902 : CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1795 804006 : dvec(1:ndiis + 1) = MATMUL(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1796 589536 : xv0(1:ns) = MATMUL(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1797 603340 : 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 9864 : cell, particle_set, xv0, rhs, rv0)
1802 : END DO
1803 6568 : charges(1:natom) = xv0(1:natom)
1804 1014 : lambda = xv0(ns)
1805 1014 : eeq_energy = eeqn
1806 1014 : IF (res > eps_diis) ierror = 1
1807 : !
1808 1014 : DEALLOCATE (xvec, fvec, bvec)
1809 1014 : DEALLOCATE (vmat, dmat, dvec)
1810 1014 : DEALLOCATE (xv0, rv0)
1811 1014 : DEALLOCATE (rhs)
1812 1014 : CALL cp_fm_release(pmat)
1813 1014 : CALL cp_fm_release(mmat)
1814 :
1815 1014 : te = m_walltime()
1816 1014 : IF (iunit > 0) THEN
1817 981 : IF (ierror == 1) THEN
1818 876 : WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1819 : ELSE
1820 105 : WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1821 : END IF
1822 981 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1823 : END IF
1824 1014 : CALL timestop(handle)
1825 :
1826 4056 : 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 892 : SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1846 892 : 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 892 : 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 892 : CALL timeset(routineN, handle)
1879 892 : ti = m_walltime()
1880 :
1881 892 : iunit = -1
1882 892 : IF (PRESENT(iounit)) iunit = iounit
1883 :
1884 892 : natom = SIZE(particle_set)
1885 892 : nkind = SIZE(gam)
1886 892 : ns = natom + 1
1887 : !
1888 892 : CALL get_cell(cell=cell, deth=deth)
1889 892 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1890 892 : ad = 2.0_dp*alpha*oorootpi
1891 892 : IF (ewald_type /= do_ewald_spme) THEN
1892 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1893 : END IF
1894 : !
1895 892 : rmax = 2.0_dp*rcut
1896 : ! max cells used
1897 892 : CALL get_cell(cell, h=hmat, periodic=periodic)
1898 892 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1899 892 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1900 892 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1901 892 : IF (periodic(1) == 0) ncell(1) = 0
1902 892 : IF (periodic(2) == 0) ncell(2) = 0
1903 892 : IF (periodic(3) == 0) ncell(3) = 0
1904 : !
1905 892 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1906 892 : 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 892 : row_indices=rind, col_indices=cind)
1909 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1910 892 : nrow_global=ns, ncol_global=1)
1911 892 : CALL cp_fm_create(rhs_vec, vec_struct)
1912 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1913 892 : row_indices=rvind, col_indices=cvind)
1914 : ! response matrix
1915 3652 : DO ir = 1, nrloc
1916 2760 : iar = rind(ir)
1917 2760 : ri = 0.0_dp
1918 2760 : IF (iar <= natom) THEN
1919 2314 : ikind = kind_of(iar)
1920 9256 : ri(1:3) = particle_set(iar)%r(1:3)
1921 : END IF
1922 22522 : DO ic = 1, ncloc
1923 18870 : iac = cind(ic)
1924 18870 : IF (iac > natom .AND. iar > natom) THEN
1925 446 : eeq_mat%local_data(ir, ic) = 0.0_dp
1926 446 : CYCLE
1927 18424 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1928 4628 : eeq_mat%local_data(ir, ic) = 1.0_dp
1929 4628 : CYCLE
1930 : END IF
1931 13796 : jkind = kind_of(iac)
1932 55184 : rj(1:3) = particle_set(iac)%r(1:3)
1933 55184 : rij(1:3) = ri(1:3) - rj(1:3)
1934 55184 : rij = pbc(rij, cell)
1935 126476 : DO ix = -ncell(1), ncell(1)
1936 636248 : DO iy = -ncell(2), ncell(2)
1937 3959316 : DO iz = -ncell(3), ncell(3)
1938 13347456 : cvec = [ix, iy, iz]
1939 63400416 : rijl = rij + MATMUL(hmat, cvec)
1940 13347456 : dr = NORM2(rijl)
1941 3336864 : IF (dr > rmax) CYCLE
1942 826540 : IF (iar == iac .AND. dr < 0.0001_dp) THEN
1943 2314 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1944 : ELSE
1945 824226 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1946 : END IF
1947 1339072 : 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 3568 : ALLOCATE (xval(natom), pval(natom))
1955 5520 : DO ia = 1, natom
1956 32220 : xval = 0.0_dp
1957 4628 : xval(ia) = 1.0_dp
1958 4628 : CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1959 : !
1960 21630 : DO ir = 1, nrloc
1961 16110 : iar = rind(ir)
1962 16110 : IF (iar /= ia) CYCLE
1963 23052 : DO ic = 1, ncloc
1964 16110 : iac = cind(ic)
1965 16110 : IF (iac > natom) CYCLE
1966 32220 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1967 : END DO
1968 : END DO
1969 : END DO
1970 892 : DEALLOCATE (xval, pval)
1971 : !
1972 : ! set up rhs vector
1973 3652 : DO ir = 1, nrvloc
1974 2760 : iar = rvind(ir)
1975 6412 : DO ic = 1, ncvloc
1976 2760 : iac = cvind(ic)
1977 2760 : ia = MAX(iar, iac)
1978 2760 : IF (ia > natom) THEN
1979 446 : xr = qtot
1980 : ELSE
1981 2314 : xr = -chia(ia)
1982 : END IF
1983 5520 : rhs_vec%local_data(ir, ic) = xr
1984 : END DO
1985 : END DO
1986 : !
1987 892 : CALL cp_fm_solve(eeq_mat, rhs_vec)
1988 : !
1989 5520 : charges = 0.0_dp
1990 892 : lambda = 0.0_dp
1991 3652 : DO ir = 1, nrvloc
1992 2760 : iar = rvind(ir)
1993 6412 : DO ic = 1, ncvloc
1994 2760 : iac = cvind(ic)
1995 2760 : ia = MAX(iar, iac)
1996 5520 : IF (ia <= natom) THEN
1997 2314 : xr = rhs_vec%local_data(ir, ic)
1998 2314 : charges(ia) = xr
1999 : ELSE
2000 446 : lambda = rhs_vec%local_data(ir, ic)
2001 : END IF
2002 : END DO
2003 : END DO
2004 892 : CALL para_env%sum(lambda)
2005 10148 : CALL para_env%sum(charges)
2006 : !
2007 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
2008 5520 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
2009 :
2010 892 : CALL cp_fm_struct_release(vec_struct)
2011 892 : CALL cp_fm_release(rhs_vec)
2012 :
2013 892 : te = m_walltime()
2014 892 : IF (iunit > 0) THEN
2015 876 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
2016 : END IF
2017 892 : CALL timestop(handle)
2018 :
2019 3568 : 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 5484 : 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 5484 : CALL ewald_env_get(ewald_env, para_env=para_env)
2041 39764 : potential = 0.0_dp
2042 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
2043 5484 : particle_set, potential)
2044 74044 : CALL para_env%sum(potential)
2045 :
2046 5484 : 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 9864 : SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
2062 9864 : 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 9864 : ns = SIZE(charges)
2079 9864 : na = ns - 1
2080 9864 : CALL ewald_env_get(ewald_env, para_env=para_env)
2081 85656 : potential = 0.0_dp
2082 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
2083 9864 : particle_set, potential(1:na))
2084 141720 : CALL para_env%sum(potential(1:na))
2085 9864 : CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
2086 151584 : eeqn = 0.5_dp*SUM(charges(1:na)*potential(1:na)) + SUM(charges(1:na)*rhs(1:na))
2087 85656 : potential(1:ns) = potential(1:ns) + rhs(1:ns)
2088 29592 : ALLOCATE (mvec(ns))
2089 9864 : CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
2090 75792 : lambda = -SUM(mvec(1:na))/REAL(na, KIND=dp)
2091 75792 : potential(1:na) = mvec(1:na) + lambda
2092 9864 : DEALLOCATE (mvec)
2093 :
2094 9864 : END SUBROUTINE get_energy_gradient
2095 :
2096 : ! **************************************************************************************************
2097 : !> \brief ...
2098 : !> \param qs_env ...
2099 : !> \param charges ...
2100 : !> \param ef_energy ...
2101 : ! **************************************************************************************************
2102 332 : 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 332 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2118 :
2119 332 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
2120 332 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2121 :
2122 332 : IF (dft_control%apply_period_efield) THEN
2123 164 : dfield = dft_control%period_efield%displacement_field
2124 :
2125 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2126 0 : CPABORT("use of strength_list not implemented for eeq_efield_energy")
2127 : END IF
2128 :
2129 656 : fieldpol = dft_control%period_efield%polarisation
2130 1148 : fieldpol = fieldpol/NORM2(fieldpol)
2131 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
2132 2132 : hmat = cell%hmat(:, :)/twopi
2133 656 : DO idir = 1, 3
2134 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
2135 656 : + fieldpol(3)*hmat(3, idir)
2136 : END DO
2137 :
2138 656 : zi(:) = CMPLX(1._dp, 0._dp, dp)
2139 820 : DO ia = 1, natom
2140 656 : q = charges(ia)
2141 2624 : ria = particle_set(ia)%r
2142 2624 : ria = pbc(ria, cell)
2143 2788 : DO idir = 1, 3
2144 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
2145 7872 : kr = SUM(kvec(:)*ria(:))
2146 1968 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
2147 2624 : zi(idir) = zi(idir)*zdeta
2148 : END DO
2149 : END DO
2150 656 : qi = AIMAG(LOG(zi))
2151 164 : IF (dfield) THEN
2152 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
2153 0 : omega = cell%deth
2154 0 : ci = MATMUL(hmat, qi)/omega
2155 0 : ef_energy = 0.0_dp
2156 0 : DO idir = 1, 3
2157 0 : ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
2158 : END DO
2159 0 : ef_energy = -0.25_dp*omega/twopi*ef_energy
2160 : ELSE
2161 656 : ef_energy = SUM(fpolvec(:)*qi(:))
2162 : END IF
2163 :
2164 168 : ELSE IF (dft_control%apply_efield) THEN
2165 :
2166 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2167 672 : dft_control%efield_fields(1)%efield%strength
2168 :
2169 168 : ef_energy = 0.0_dp
2170 840 : DO ia = 1, natom
2171 2688 : ria = particle_set(ia)%r
2172 2688 : ria = pbc(ria, cell)
2173 672 : q = charges(ia)
2174 2856 : ef_energy = ef_energy - q*SUM(fieldpol*ria)
2175 : END DO
2176 :
2177 : ELSE
2178 0 : CPABORT("apply field")
2179 : END IF
2180 :
2181 332 : END SUBROUTINE eeq_efield_energy
2182 :
2183 : ! **************************************************************************************************
2184 : !> \brief ...
2185 : !> \param qs_env ...
2186 : !> \param efr ...
2187 : ! **************************************************************************************************
2188 332 : 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 332 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2200 :
2201 332 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
2202 332 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2203 :
2204 332 : IF (dft_control%apply_period_efield) THEN
2205 164 : dfield = dft_control%period_efield%displacement_field
2206 :
2207 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2208 0 : CPABORT("use of strength_list not implemented for eeq_efield_pot")
2209 : END IF
2210 :
2211 656 : fieldpol = dft_control%period_efield%polarisation
2212 1148 : fieldpol = fieldpol/NORM2(fieldpol)
2213 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
2214 2132 : hmat = cell%hmat(:, :)/twopi
2215 656 : DO idir = 1, 3
2216 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
2217 656 : + fieldpol(3)*hmat(3, idir)
2218 : END DO
2219 :
2220 164 : IF (dfield) THEN
2221 : ! dE/dq depends on q, postpone calculation
2222 0 : efr = 0.0_dp
2223 : ELSE
2224 820 : efr = 0.0_dp
2225 820 : DO ia = 1, natom
2226 2624 : ria = particle_set(ia)%r
2227 2624 : ria = pbc(ria, cell)
2228 2788 : DO idir = 1, 3
2229 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
2230 7872 : kr = SUM(kvec(:)*ria(:))
2231 2624 : efr(ia) = efr(ia) + kr*fpolvec(idir)
2232 : END DO
2233 : END DO
2234 : END IF
2235 :
2236 168 : ELSE IF (dft_control%apply_efield) THEN
2237 :
2238 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2239 672 : dft_control%efield_fields(1)%efield%strength
2240 :
2241 840 : DO ia = 1, natom
2242 2688 : ria = particle_set(ia)%r
2243 2688 : ria = pbc(ria, cell)
2244 2856 : efr(ia) = -SUM(fieldpol*ria)
2245 : END DO
2246 :
2247 : ELSE
2248 0 : CPABORT("apply field")
2249 : END IF
2250 :
2251 332 : 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 0 : 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 0 : natom = SIZE(particle_set)
2276 :
2277 0 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2278 0 : CPABORT("use of strength_list not implemented for eeq_dfield_pot")
2279 : END IF
2280 :
2281 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
2282 0 : fieldpol = dft_control%period_efield%polarisation
2283 0 : fieldpol = fieldpol/NORM2(fieldpol)
2284 0 : fieldpol = -fieldpol*dft_control%period_efield%strength
2285 0 : hmat = cell%hmat(:, :)/twopi
2286 0 : omega = cell%deth
2287 : !
2288 0 : zi(:) = CMPLX(1._dp, 0._dp, dp)
2289 0 : DO ia = 1, natom
2290 0 : q = charges(ia)
2291 0 : ria = particle_set(ia)%r
2292 0 : ria = pbc(ria, cell)
2293 0 : DO idir = 1, 3
2294 0 : kvec(:) = twopi*cell%h_inv(idir, :)
2295 0 : kr = SUM(kvec(:)*ria(:))
2296 0 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
2297 0 : zi(idir) = zi(idir)*zdeta
2298 : END DO
2299 : END DO
2300 0 : qi = AIMAG(LOG(zi))
2301 0 : ci = MATMUL(hmat, qi)/omega
2302 0 : ci = dfilter*(fieldpol - 2._dp*twopi*ci)
2303 0 : DO ia = 1, natom
2304 0 : ria = particle_set(ia)%r
2305 0 : ria = pbc(ria, cell)
2306 0 : efr(ia) = efr(ia) - SUM(ci*ria)
2307 : END DO
2308 :
2309 0 : END SUBROUTINE eeq_dfield_pot
2310 :
2311 : ! **************************************************************************************************
2312 : !> \brief ...
2313 : !> \param qs_env ...
2314 : !> \param charges ...
2315 : !> \param qlag ...
2316 : ! **************************************************************************************************
2317 8 : 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 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
2323 : REAL(KIND=dp) :: q
2324 : REAL(KIND=dp), DIMENSION(3) :: fieldpol
2325 8 : 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 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2331 8 : 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 8 : local_particles=local_particles)
2339 :
2340 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
2341 32 : dft_control%efield_fields(1)%efield%strength
2342 :
2343 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
2344 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
2345 8 : CALL get_qs_env(qs_env=qs_env, force=force)
2346 :
2347 32 : DO ikind = 1, nkind
2348 152 : force(ikind)%efield = 0.0_dp
2349 40 : DO ia = 1, local_particles%n_el(ikind)
2350 16 : iatom = local_particles%list(ikind)%array(ia)
2351 16 : q = charges(iatom) - qlag(iatom)
2352 16 : atom_a = atom_of_kind(iatom)
2353 88 : force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
2354 : END DO
2355 288 : CALL para_env%sum(force(ikind)%efield)
2356 : END DO
2357 :
2358 16 : END SUBROUTINE eeq_efield_force_loc
2359 :
2360 : ! **************************************************************************************************
2361 : !> \brief ...
2362 : !> \param qs_env ...
2363 : !> \param charges ...
2364 : !> \param qlag ...
2365 : ! **************************************************************************************************
2366 8 : 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 8 : 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 8 : 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 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2382 8 : 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 8 : local_particles=local_particles)
2392 :
2393 8 : dfield = dft_control%period_efield%displacement_field
2394 8 : CPASSERT(.NOT. dfield)
2395 :
2396 8 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
2397 0 : CPABORT("use of strength_list not implemented for eeq_efield_force_periodic")
2398 : END IF
2399 :
2400 32 : fieldpol = dft_control%period_efield%polarisation
2401 56 : fieldpol = fieldpol/NORM2(fieldpol)
2402 32 : fieldpol = -fieldpol*dft_control%period_efield%strength
2403 :
2404 8 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
2405 :
2406 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
2407 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
2408 8 : CALL get_qs_env(qs_env=qs_env, force=force)
2409 :
2410 8 : pve = 0.0_dp
2411 32 : DO ikind = 1, nkind
2412 152 : force(ikind)%efield = 0.0_dp
2413 40 : DO ia = 1, local_particles%n_el(ikind)
2414 16 : iatom = local_particles%list(ikind)%array(ia)
2415 16 : q = charges(iatom) - qlag(iatom)
2416 64 : fa(1:3) = q*fieldpol(1:3)
2417 16 : atom_a = atom_of_kind(iatom)
2418 64 : force(ikind)%efield(1:3, atom_a) = fa
2419 40 : IF (use_virial) THEN
2420 0 : ria = particle_set(ia)%r
2421 0 : ria = pbc(ria, cell)
2422 0 : CALL virial_pair_force(pve, -0.5_dp, fa, ria)
2423 0 : CALL virial_pair_force(pve, -0.5_dp, ria, fa)
2424 : END IF
2425 : END DO
2426 288 : CALL para_env%sum(force(ikind)%efield)
2427 : END DO
2428 104 : virial%pv_virial = virial%pv_virial + pve
2429 :
2430 16 : END SUBROUTINE eeq_efield_force_periodic
2431 :
2432 13804 : END MODULE eeq_method
|
