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 : USE machine, ONLY: m_walltime
55 : USE mathconstants, ONLY: oorootpi,&
56 : twopi
57 : USE mathlib, ONLY: invmat
58 : USE message_passing, ONLY: mp_para_env_type
59 : USE molecule_types, ONLY: molecule_type
60 : USE particle_types, ONLY: particle_type
61 : USE physcon, ONLY: bohr
62 : USE pw_poisson_types, ONLY: do_ewald_spme
63 : USE qs_dispersion_cnum, ONLY: cnumber_init,&
64 : cnumber_release,&
65 : dcnum_type
66 : USE qs_dispersion_types, ONLY: qs_dispersion_release,&
67 : qs_dispersion_type
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
70 : USE qs_force_types, ONLY: qs_force_type
71 : USE qs_kind_types, ONLY: get_qs_kind,&
72 : qs_kind_type
73 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
74 : neighbor_list_iterate,&
75 : neighbor_list_iterator_create,&
76 : neighbor_list_iterator_p_type,&
77 : neighbor_list_iterator_release,&
78 : neighbor_list_set_p_type,&
79 : release_neighbor_list_sets
80 : USE qs_neighbor_lists, ONLY: atom2d_build,&
81 : atom2d_cleanup,&
82 : build_neighbor_lists,&
83 : local_atoms_type,&
84 : pair_radius_setup
85 : USE spme, ONLY: spme_forces,&
86 : spme_potential,&
87 : spme_virial
88 : USE virial_methods, ONLY: virial_pair_force
89 : USE virial_types, ONLY: virial_type
90 : #include "./base/base_uses.f90"
91 :
92 : IMPLICIT NONE
93 :
94 : PRIVATE
95 :
96 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eeq_method'
97 :
98 : INTEGER, PARAMETER :: maxElem = 86
99 : ! covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 188-197)
100 : ! values for metals decreased by 10 %
101 : REAL(KIND=dp), PARAMETER :: rcov(1:maxElem) = [&
102 : & 0.32_dp, 0.46_dp, 1.20_dp, 0.94_dp, 0.77_dp, 0.75_dp, 0.71_dp, 0.63_dp, &
103 : & 0.64_dp, 0.67_dp, 1.40_dp, 1.25_dp, 1.13_dp, 1.04_dp, 1.10_dp, 1.02_dp, &
104 : & 0.99_dp, 0.96_dp, 1.76_dp, 1.54_dp, 1.33_dp, 1.22_dp, 1.21_dp, 1.10_dp, &
105 : & 1.07_dp, 1.04_dp, 1.00_dp, 0.99_dp, 1.01_dp, 1.09_dp, 1.12_dp, 1.09_dp, &
106 : & 1.15_dp, 1.10_dp, 1.14_dp, 1.17_dp, 1.89_dp, 1.67_dp, 1.47_dp, 1.39_dp, &
107 : & 1.32_dp, 1.24_dp, 1.15_dp, 1.13_dp, 1.13_dp, 1.08_dp, 1.15_dp, 1.23_dp, &
108 : & 1.28_dp, 1.26_dp, 1.26_dp, 1.23_dp, 1.32_dp, 1.31_dp, 2.09_dp, 1.76_dp, &
109 : & 1.62_dp, 1.47_dp, 1.58_dp, 1.57_dp, 1.56_dp, 1.55_dp, 1.51_dp, 1.52_dp, &
110 : & 1.51_dp, 1.50_dp, 1.49_dp, 1.49_dp, 1.48_dp, 1.53_dp, 1.46_dp, 1.37_dp, &
111 : & 1.31_dp, 1.23_dp, 1.18_dp, 1.16_dp, 1.11_dp, 1.12_dp, 1.13_dp, 1.32_dp, &
112 : & 1.30_dp, 1.30_dp, 1.36_dp, 1.31_dp, 1.38_dp, 1.42_dp]
113 :
114 : PUBLIC :: eeq_solver, eeq_print, eeq_charges, eeq_forces, &
115 : eeq_efield_energy, eeq_efield_pot, eeq_efield_force_loc, eeq_efield_force_periodic
116 :
117 : CONTAINS
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param qs_env ...
122 : !> \param iounit ...
123 : !> \param print_level ...
124 : !> \param ext ...
125 : ! **************************************************************************************************
126 38 : SUBROUTINE eeq_print(qs_env, iounit, print_level, ext)
127 :
128 : TYPE(qs_environment_type), POINTER :: qs_env
129 : INTEGER, INTENT(IN) :: iounit, print_level
130 : LOGICAL, INTENT(IN) :: ext
131 :
132 : CHARACTER(LEN=2) :: element_symbol
133 : INTEGER :: enshift_type, iatom, ikind, natom
134 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
135 : TYPE(cell_type), POINTER :: cell
136 : TYPE(eeq_solver_type) :: eeq_sparam
137 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
138 :
139 : MARK_USED(print_level)
140 :
141 38 : CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, cell=cell)
142 38 : IF (ext) THEN
143 0 : NULLIFY (charges)
144 0 : CALL get_qs_env(qs_env, eeq=charges)
145 0 : CPASSERT(ASSOCIATED(charges))
146 0 : enshift_type = 0
147 : ELSE
148 114 : ALLOCATE (charges(natom))
149 : ! enforce en shift method 1 (original/molecular)
150 : ! method 2 from paper on PBC seems not to work
151 38 : enshift_type = 1
152 : !IF (ALL(cell%perd == 0)) enshift_type = 1
153 38 : CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
154 : END IF
155 :
156 38 : IF (iounit > 0) THEN
157 :
158 19 : IF (enshift_type == 0) THEN
159 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (External)"
160 19 : ELSE IF (enshift_type == 1) THEN
161 19 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Molecules))"
162 0 : ELSE IF (enshift_type == 2) THEN
163 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Crystals))"
164 : ELSE
165 0 : CPABORT("Unknown enshift_type")
166 : END IF
167 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
168 19 : "# Atom Element Kind Atomic Charge"
169 :
170 140 : DO iatom = 1, natom
171 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
172 : element_symbol=element_symbol, &
173 121 : kind_number=ikind)
174 : WRITE (UNIT=iounit, FMT="(T4,I8,T18,A2,I10,T43,F12.4)") &
175 140 : iatom, element_symbol, ikind, charges(iatom)
176 : END DO
177 :
178 : END IF
179 :
180 38 : IF (.NOT. ext) DEALLOCATE (charges)
181 :
182 38 : END SUBROUTINE eeq_print
183 :
184 : ! **************************************************************************************************
185 : !> \brief ...
186 : !> \param qs_env ...
187 : !> \param charges ...
188 : !> \param eeq_sparam ...
189 : !> \param eeq_model ...
190 : !> \param enshift_type ...
191 : !> \param exclude ...
192 : !> \param cn_max ...
193 : ! **************************************************************************************************
194 68 : SUBROUTINE eeq_charges(qs_env, charges, eeq_sparam, eeq_model, enshift_type, exclude, cn_max)
195 :
196 : TYPE(qs_environment_type), POINTER :: qs_env
197 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
198 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
199 : INTEGER, INTENT(IN) :: eeq_model, enshift_type
200 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: exclude
201 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cn_max
202 :
203 : CHARACTER(len=*), PARAMETER :: routineN = 'eeq_charges'
204 :
205 : INTEGER :: handle, iatom, ikind, iunit, jkind, &
206 : natom, nkind, za, zb
207 68 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
208 : INTEGER, DIMENSION(3) :: periodic
209 : LOGICAL :: do_ewald
210 : REAL(KIND=dp) :: ala, alb, eeq_energy, esg, kappa, &
211 : lambda, scn, sgamma, totalcharge, xi
212 68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chia, cnumbers, efr, gam
213 68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
214 68 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
215 : TYPE(cell_type), POINTER :: cell, cell_ref
216 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
217 : TYPE(cp_logger_type), POINTER :: logger
218 68 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
219 : TYPE(dft_control_type), POINTER :: dft_control
220 : TYPE(ewald_environment_type), POINTER :: ewald_env
221 : TYPE(ewald_pw_type), POINTER :: ewald_pw
222 : TYPE(mp_para_env_type), POINTER :: para_env
223 68 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
224 68 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
225 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
226 : print_section
227 :
228 68 : CALL timeset(routineN, handle)
229 :
230 : CALL get_qs_env(qs_env, &
231 : qs_kind_set=qs_kind_set, &
232 : atomic_kind_set=atomic_kind_set, &
233 : particle_set=particle_set, &
234 : para_env=para_env, &
235 : blacs_env=blacs_env, &
236 : cell=cell, &
237 68 : dft_control=dft_control)
238 68 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
239 :
240 68 : logger => cp_get_default_logger()
241 68 : IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
242 26 : iunit = cp_logger_get_default_unit_nr()
243 : ELSE
244 42 : iunit = -1
245 : END IF
246 :
247 68 : totalcharge = dft_control%charge
248 :
249 68 : CALL get_cnumbers(qs_env, cnumbers, dcnum, .FALSE.)
250 :
251 : ! Apply smooth logistic CN cutoff to match multicharge/D4 behavior
252 68 : IF (PRESENT(cn_max)) THEN
253 266 : DO iatom = 1, natom
254 266 : cnumbers(iatom) = LOG(1.0_dp + EXP(cn_max)) - LOG(1.0_dp + EXP(cn_max - cnumbers(iatom)))
255 : END DO
256 : END IF
257 :
258 : ! gamma[a,b]
259 408 : ALLOCATE (gab(nkind, nkind), gam(nkind))
260 560 : gab = 0.0_dp
261 214 : gam = 0.0_dp
262 214 : DO ikind = 1, nkind
263 146 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
264 146 : CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
265 560 : DO jkind = 1, nkind
266 346 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
267 346 : CALL get_eeq_data(zb, eeq_model, rad=alb)
268 : !
269 492 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
270 : !
271 : END DO
272 : END DO
273 :
274 : ! Override parameters for excluded kinds (ghost/floating atoms in BSSE):
275 : ! huge hardness + zero coupling -> q = 0, no influence on other atoms.
276 68 : IF (PRESENT(exclude)) THEN
277 96 : DO ikind = 1, nkind
278 96 : IF (exclude(ikind)) THEN
279 8 : gam(ikind) = 1.0e30_dp
280 40 : gab(ikind, :) = 0.0_dp
281 40 : gab(:, ikind) = 0.0_dp
282 : END IF
283 : END DO
284 : END IF
285 :
286 : ! Chi[a,a]
287 68 : sgamma = 8.0_dp ! see D4 for periodic systems paper
288 68 : esg = 1.0_dp + EXP(sgamma)
289 204 : ALLOCATE (chia(natom))
290 68 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
291 546 : DO iatom = 1, natom
292 478 : ikind = kind_of(iatom)
293 478 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
294 478 : CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
295 : !
296 478 : IF (enshift_type == 1) THEN
297 478 : scn = cnumbers(iatom)/SQRT(cnumbers(iatom) + 1.0e-14_dp)
298 0 : ELSE IF (enshift_type == 2) THEN
299 0 : scn = LOG(esg/(esg - cnumbers(iatom)))
300 : ELSE
301 0 : CPABORT("Unknown enshift_type")
302 : END IF
303 1024 : chia(iatom) = xi - kappa*scn
304 : !
305 : END DO
306 :
307 : ! Zero electronegativity for excluded atoms (ghost/floating in BSSE)
308 68 : IF (PRESENT(exclude)) THEN
309 266 : DO iatom = 1, natom
310 236 : ikind = kind_of(iatom)
311 266 : IF (exclude(ikind)) chia(iatom) = 0.0_dp
312 : END DO
313 : END IF
314 :
315 : ! efield
316 68 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
317 : dft_control%apply_efield_field) THEN
318 0 : ALLOCATE (efr(natom))
319 0 : efr(1:natom) = 0.0_dp
320 0 : CALL eeq_efield_pot(qs_env, efr)
321 0 : chia(1:natom) = chia(1:natom) + efr(1:natom)
322 0 : DEALLOCATE (efr)
323 : END IF
324 :
325 68 : CALL cnumber_release(cnumbers, dcnum, .FALSE.)
326 :
327 68 : CALL get_cell(cell, periodic=periodic)
328 134 : do_ewald = .NOT. ALL(periodic == 0)
329 68 : IF (do_ewald) THEN
330 736 : ALLOCATE (ewald_env)
331 46 : CALL ewald_env_create(ewald_env, para_env)
332 46 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
333 46 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
334 46 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
335 46 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
336 46 : CALL get_qs_env(qs_env, cell_ref=cell_ref)
337 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
338 46 : silent=.TRUE., pset="EEQ")
339 46 : ALLOCATE (ewald_pw)
340 46 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
341 : !
342 : CALL eeq_solver(charges, lambda, eeq_energy, &
343 : particle_set, kind_of, cell, chia, gam, gab, &
344 : para_env, blacs_env, dft_control, eeq_sparam, &
345 : totalcharge=totalcharge, ewald=do_ewald, &
346 46 : ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
347 : !
348 46 : CALL ewald_env_release(ewald_env)
349 46 : CALL ewald_pw_release(ewald_pw)
350 46 : DEALLOCATE (ewald_env, ewald_pw)
351 : ELSE
352 : CALL eeq_solver(charges, lambda, eeq_energy, &
353 : particle_set, kind_of, cell, chia, gam, gab, &
354 : para_env, blacs_env, dft_control, eeq_sparam, &
355 22 : totalcharge=totalcharge, iounit=iunit)
356 : END IF
357 :
358 68 : DEALLOCATE (gab, gam, chia)
359 :
360 68 : CALL timestop(handle)
361 :
362 136 : END SUBROUTINE eeq_charges
363 :
364 : ! **************************************************************************************************
365 : !> \brief ...
366 : !> \param qs_env ...
367 : !> \param charges ...
368 : !> \param dcharges ...
369 : !> \param gradient ...
370 : !> \param stress ...
371 : !> \param eeq_sparam ...
372 : !> \param eeq_model ...
373 : !> \param enshift_type ...
374 : !> \param response_only ...
375 : !> \param exclude ...
376 : !> \param cn_max ...
377 : ! **************************************************************************************************
378 10 : SUBROUTINE eeq_forces(qs_env, charges, dcharges, gradient, stress, &
379 10 : eeq_sparam, eeq_model, enshift_type, response_only, exclude, cn_max)
380 :
381 : TYPE(qs_environment_type), POINTER :: qs_env
382 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges
383 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: gradient
384 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: stress
385 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
386 : INTEGER, INTENT(IN) :: eeq_model, enshift_type
387 : LOGICAL, INTENT(IN) :: response_only
388 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: exclude
389 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cn_max
390 :
391 : CHARACTER(len=*), PARAMETER :: routineN = 'eeq_forces'
392 :
393 : INTEGER :: handle, i, ia, iatom, ikind, iunit, &
394 : jatom, jkind, katom, natom, nkind, za, &
395 : zb
396 10 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
397 : INTEGER, DIMENSION(3) :: periodic
398 : LOGICAL :: do_ewald, use_virial
399 10 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
400 : REAL(KIND=dp) :: ala, alb, alpha, cn, ctot, dcnpdcn, dr, dr2, drk, elag, esg, fe, gam2, &
401 : gama, grc, kappa, qlam, qq, qq1, qq2, rcut, scn, sgamma, subcells, totalcharge, xi
402 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c_radius, cnumbers, gam, qlag
403 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: epforce, gab, pair_radius
404 : REAL(KIND=dp), DIMENSION(3) :: fdik, ri, rij, rik, rj
405 : REAL(KIND=dp), DIMENSION(3, 3) :: pvir
406 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: chrgx, dchia
407 10 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
408 : TYPE(atprop_type), POINTER :: atprop
409 : TYPE(cell_type), POINTER :: cell, cell_ref
410 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
411 : TYPE(cp_logger_type), POINTER :: logger
412 10 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
413 : TYPE(dft_control_type), POINTER :: dft_control
414 : TYPE(distribution_1d_type), POINTER :: distribution_1d, local_particles
415 : TYPE(distribution_2d_type), POINTER :: distribution_2d
416 : TYPE(ewald_environment_type), POINTER :: ewald_env
417 : TYPE(ewald_pw_type), POINTER :: ewald_pw
418 10 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
419 10 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
420 : TYPE(mp_para_env_type), POINTER :: para_env
421 : TYPE(neighbor_list_iterator_p_type), &
422 10 : DIMENSION(:), POINTER :: nl_iterator
423 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
424 10 : POINTER :: sab_ew
425 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426 10 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
427 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
428 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
429 : print_section
430 : TYPE(virial_type), POINTER :: virial
431 :
432 10 : CALL timeset(routineN, handle)
433 :
434 : CALL get_qs_env(qs_env, &
435 : qs_kind_set=qs_kind_set, &
436 : atomic_kind_set=atomic_kind_set, &
437 : particle_set=particle_set, &
438 : para_env=para_env, &
439 : blacs_env=blacs_env, &
440 : cell=cell, &
441 : force=force, &
442 : virial=virial, &
443 : atprop=atprop, &
444 10 : dft_control=dft_control)
445 :
446 10 : logger => cp_get_default_logger()
447 10 : IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
448 2 : iunit = cp_logger_get_default_unit_nr()
449 : ELSE
450 8 : iunit = -1
451 : END IF
452 :
453 10 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
454 10 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
455 :
456 10 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
457 :
458 10 : totalcharge = dft_control%charge
459 :
460 10 : CALL get_cnumbers(qs_env, cnumbers, dcnum, .TRUE.)
461 :
462 : ! Apply smooth logistic CN cutoff to match multicharge/D4 behavior
463 : ! Chain rule: dcn_cut/dR = (dcn_cut/dcn) * (dcn/dR)
464 10 : IF (PRESENT(cn_max)) THEN
465 98 : DO iatom = 1, natom
466 88 : dcnpdcn = EXP(cn_max)/(EXP(cn_max) + EXP(cnumbers(iatom)))
467 88 : cnumbers(iatom) = LOG(1.0_dp + EXP(cn_max)) - LOG(1.0_dp + EXP(cn_max - cnumbers(iatom)))
468 226 : DO i = 1, dcnum(iatom)%neighbors
469 216 : dcnum(iatom)%dvals(i) = dcnum(iatom)%dvals(i)*dcnpdcn
470 : END DO
471 : END DO
472 : END IF
473 :
474 : ! gamma[a,b]
475 60 : ALLOCATE (gab(nkind, nkind), gam(nkind))
476 62 : gab = 0.0_dp
477 28 : gam = 0.0_dp
478 28 : DO ikind = 1, nkind
479 18 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
480 18 : CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
481 62 : DO jkind = 1, nkind
482 34 : CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
483 34 : CALL get_eeq_data(zb, eeq_model, rad=alb)
484 : !
485 52 : gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
486 : !
487 : END DO
488 : END DO
489 :
490 : ! Override parameters for excluded kinds (ghost/floating atoms in BSSE)
491 10 : IF (PRESENT(exclude)) THEN
492 28 : DO ikind = 1, nkind
493 28 : IF (exclude(ikind)) THEN
494 0 : gam(ikind) = 1.0e30_dp
495 0 : gab(ikind, :) = 0.0_dp
496 0 : gab(:, ikind) = 0.0_dp
497 : END IF
498 : END DO
499 : END IF
500 :
501 30 : ALLOCATE (qlag(natom))
502 :
503 10 : CALL get_cell(cell, periodic=periodic)
504 22 : do_ewald = .NOT. ALL(periodic == 0)
505 10 : IF (do_ewald) THEN
506 96 : ALLOCATE (ewald_env)
507 6 : CALL ewald_env_create(ewald_env, para_env)
508 6 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
509 6 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
510 6 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
511 6 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
512 6 : CALL get_qs_env(qs_env, cell_ref=cell_ref)
513 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
514 6 : silent=.TRUE., pset="EEQ")
515 6 : ALLOCATE (ewald_pw)
516 6 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
517 : !
518 : CALL eeq_solver(qlag, qlam, elag, &
519 : particle_set, kind_of, cell, -dcharges, gam, gab, &
520 : para_env, blacs_env, dft_control, eeq_sparam, &
521 54 : ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
522 : ELSE
523 : CALL eeq_solver(qlag, qlam, elag, &
524 : particle_set, kind_of, cell, -dcharges, gam, gab, &
525 44 : para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
526 : END IF
527 :
528 10 : sgamma = 8.0_dp ! see D4 for periodic systems paper
529 10 : esg = 1.0_dp + EXP(sgamma)
530 40 : ALLOCATE (chrgx(natom), dchia(natom))
531 98 : DO iatom = 1, natom
532 88 : ikind = kind_of(iatom)
533 88 : CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
534 88 : CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
535 : !
536 88 : IF (response_only) THEN
537 88 : ctot = -0.5_dp*qlag(iatom)
538 : ELSE
539 0 : ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
540 : END IF
541 186 : IF (enshift_type == 1) THEN
542 88 : scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
543 88 : dchia(iatom) = -ctot*kappa/scn
544 0 : ELSE IF (enshift_type == 2) THEN
545 0 : cn = cnumbers(iatom)
546 0 : scn = 1.0_dp/(esg - cn)
547 0 : dchia(iatom) = -ctot*kappa*scn
548 : ELSE
549 0 : CPABORT("Unknown enshift_type")
550 : END IF
551 : END DO
552 :
553 : ! Efield
554 10 : IF (dft_control%apply_period_efield) THEN
555 0 : CALL eeq_efield_force_periodic(qs_env, charges, qlag)
556 10 : ELSE IF (dft_control%apply_efield) THEN
557 0 : CALL eeq_efield_force_loc(qs_env, charges, qlag)
558 10 : ELSE IF (dft_control%apply_efield_field) THEN
559 0 : CPABORT("apply field")
560 : END IF
561 :
562 : ! Forces from q*X
563 10 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
564 28 : DO ikind = 1, nkind
565 72 : DO ia = 1, local_particles%n_el(ikind)
566 44 : iatom = local_particles%list(ikind)%array(ia)
567 126 : DO i = 1, dcnum(iatom)%neighbors
568 64 : katom = dcnum(iatom)%nlist(i)
569 256 : rik = dcnum(iatom)%rik(:, i)
570 256 : drk = SQRT(SUM(rik(:)**2))
571 108 : IF (drk > 1.e-3_dp) THEN
572 256 : fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
573 256 : gradient(:, iatom) = gradient(:, iatom) - fdik(:)
574 256 : gradient(:, katom) = gradient(:, katom) + fdik(:)
575 64 : IF (use_virial) THEN
576 16 : CALL virial_pair_force(stress, 1._dp, fdik, rik)
577 : END IF
578 : END IF
579 : END DO
580 : END DO
581 : END DO
582 :
583 : ! Forces from (0.5*q+l)*dA/dR*q
584 10 : IF (do_ewald) THEN
585 :
586 : ! Build the neighbor lists for the CN
587 : CALL get_qs_env(qs_env, &
588 : distribution_2d=distribution_2d, &
589 : local_particles=distribution_1d, &
590 6 : molecule_set=molecule_set)
591 6 : subcells = 2.0_dp
592 6 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
593 6 : rcut = 2.0_dp*rcut
594 6 : NULLIFY (sab_ew)
595 48 : ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
596 16 : c_radius(:) = rcut
597 16 : default_present = .TRUE.
598 28 : ALLOCATE (atom2d(nkind))
599 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
600 6 : molecule_set, .FALSE., particle_set=particle_set)
601 6 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
602 : CALL build_neighbor_lists(sab_ew, particle_set, atom2d, cell, pair_radius, &
603 6 : subcells=subcells, operator_type="PP", nlname="sab_ew")
604 6 : DEALLOCATE (c_radius, pair_radius, default_present)
605 6 : CALL atom2d_cleanup(atom2d)
606 : !
607 6 : CALL neighbor_list_iterator_create(nl_iterator, sab_ew)
608 77994 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
609 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
610 77988 : iatom=iatom, jatom=jatom, r=rij)
611 : !
612 311952 : dr2 = SUM(rij**2)
613 77988 : dr = SQRT(dr2)
614 77988 : IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
615 9758 : fe = 1.0_dp
616 9758 : IF (iatom == jatom) fe = 0.5_dp
617 9758 : IF (response_only) THEN
618 9758 : qq = -qlag(iatom)*charges(jatom)
619 : ELSE
620 9758 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
621 : END IF
622 9758 : gama = gab(ikind, jkind)
623 9758 : gam2 = gama*gama
624 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
625 9758 : - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
626 9758 : IF (response_only) THEN
627 9758 : qq1 = -qlag(iatom)*charges(jatom)
628 9758 : qq2 = -qlag(jatom)*charges(iatom)
629 : ELSE
630 0 : qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
631 0 : qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
632 : END IF
633 39032 : fdik(:) = -qq1*grc*rij(:)/dr
634 39032 : gradient(:, iatom) = gradient(:, iatom) + fdik(:)
635 39032 : gradient(:, jatom) = gradient(:, jatom) - fdik(:)
636 9758 : IF (use_virial) THEN
637 5279 : CALL virial_pair_force(stress, -fe, fdik, rij)
638 : END IF
639 39032 : fdik(:) = qq2*grc*rij(:)/dr
640 39032 : gradient(:, iatom) = gradient(:, iatom) - fdik(:)
641 39032 : gradient(:, jatom) = gradient(:, jatom) + fdik(:)
642 9764 : IF (use_virial) THEN
643 5279 : CALL virial_pair_force(stress, fe, fdik, rij)
644 : END IF
645 : END DO
646 6 : CALL neighbor_list_iterator_release(nl_iterator)
647 : !
648 6 : CALL release_neighbor_list_sets(sab_ew)
649 : ELSE
650 12 : DO ikind = 1, nkind
651 32 : DO ia = 1, local_particles%n_el(ikind)
652 20 : iatom = local_particles%list(ikind)%array(ia)
653 80 : ri(1:3) = particle_set(iatom)%r(1:3)
654 228 : DO jatom = 1, natom
655 200 : IF (iatom == jatom) CYCLE
656 180 : jkind = kind_of(jatom)
657 180 : IF (response_only) THEN
658 180 : qq = -qlag(iatom)*charges(jatom)
659 : ELSE
660 0 : qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
661 : END IF
662 720 : rj(1:3) = particle_set(jatom)%r(1:3)
663 720 : rij(1:3) = ri(1:3) - rj(1:3)
664 720 : rij = pbc(rij, cell)
665 720 : dr2 = SUM(rij**2)
666 180 : dr = SQRT(dr2)
667 180 : gama = gab(ikind, jkind)
668 180 : gam2 = gama*gama
669 180 : grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
670 720 : fdik(:) = qq*grc*rij(:)/dr
671 720 : gradient(:, iatom) = gradient(:, iatom) + fdik(:)
672 740 : gradient(:, jatom) = gradient(:, jatom) - fdik(:)
673 : END DO
674 : END DO
675 : END DO
676 : END IF
677 :
678 : ! Forces from Ewald potential: (q+l)*A*q
679 10 : IF (do_ewald) THEN
680 18 : ALLOCATE (epforce(3, natom))
681 198 : epforce = 0.0_dp
682 6 : IF (response_only) THEN
683 54 : dchia(1:natom) = qlag(1:natom)
684 : ELSE
685 0 : dchia(1:natom) = -charges(1:natom) + qlag(1:natom)
686 : END IF
687 54 : chrgx(1:natom) = charges(1:natom)
688 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
689 6 : particle_set, dchia, epforce)
690 54 : dchia(1:natom) = charges(1:natom)
691 54 : chrgx(1:natom) = qlag(1:natom)
692 : CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
693 6 : particle_set, dchia, epforce)
694 198 : gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + epforce(1:3, 1:natom)
695 6 : DEALLOCATE (epforce)
696 :
697 : ! virial
698 6 : IF (use_virial) THEN
699 32 : chrgx(1:natom) = charges(1:natom) - qlag(1:natom)
700 4 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
701 52 : stress = stress - pvir
702 32 : chrgx(1:natom) = qlag(1:natom)
703 4 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
704 52 : stress = stress + pvir
705 4 : IF (response_only) THEN
706 32 : chrgx(1:natom) = charges(1:natom)
707 4 : CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
708 52 : stress = stress + pvir
709 : END IF
710 : END IF
711 : !
712 6 : CALL ewald_env_release(ewald_env)
713 6 : CALL ewald_pw_release(ewald_pw)
714 6 : DEALLOCATE (ewald_env, ewald_pw)
715 : END IF
716 :
717 10 : CALL cnumber_release(cnumbers, dcnum, .TRUE.)
718 :
719 10 : DEALLOCATE (gab, gam, qlag, chrgx, dchia)
720 :
721 10 : CALL timestop(handle)
722 :
723 20 : END SUBROUTINE eeq_forces
724 :
725 : ! **************************************************************************************************
726 : !> \brief ...
727 : !> \param qs_env ...
728 : !> \param cnumbers ...
729 : !> \param dcnum ...
730 : !> \param calculate_forces ...
731 : ! **************************************************************************************************
732 78 : SUBROUTINE get_cnumbers(qs_env, cnumbers, dcnum, calculate_forces)
733 :
734 : TYPE(qs_environment_type), POINTER :: qs_env
735 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
736 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
737 : LOGICAL, INTENT(IN) :: calculate_forces
738 :
739 : INTEGER :: ikind, natom, nkind, za
740 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
741 : REAL(KIND=dp) :: subcells
742 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c_radius
743 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
744 78 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
745 : TYPE(cell_type), POINTER :: cell
746 : TYPE(distribution_1d_type), POINTER :: distribution_1d
747 : TYPE(distribution_2d_type), POINTER :: distribution_2d
748 78 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
749 78 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
750 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
751 78 : POINTER :: sab_cn
752 78 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
753 : TYPE(qs_dispersion_type), POINTER :: disp
754 78 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
755 :
756 : CALL get_qs_env(qs_env, &
757 : qs_kind_set=qs_kind_set, &
758 : atomic_kind_set=atomic_kind_set, &
759 : particle_set=particle_set, &
760 : cell=cell, &
761 : distribution_2d=distribution_2d, &
762 : local_particles=distribution_1d, &
763 78 : molecule_set=molecule_set)
764 78 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
765 :
766 : ! Check for dispersion_env and sab_cn needed for cnumbers
767 390 : ALLOCATE (disp)
768 78 : disp%k1 = 16.0_dp
769 78 : disp%k2 = 4._dp/3._dp
770 78 : disp%eps_cn = 1.E-6_dp
771 78 : disp%max_elem = maxElem
772 78 : ALLOCATE (disp%rcov(maxElem))
773 6786 : disp%rcov(1:maxElem) = bohr*disp%k2*rcov(1:maxElem)
774 78 : subcells = 2.0_dp
775 : ! Build the neighbor lists for the CN
776 78 : NULLIFY (sab_cn)
777 624 : ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
778 242 : c_radius(:) = 0.0_dp
779 242 : default_present = .TRUE.
780 242 : DO ikind = 1, nkind
781 164 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
782 242 : c_radius(ikind) = 4._dp*rcov(za)*bohr
783 : END DO
784 398 : ALLOCATE (atom2d(nkind))
785 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
786 78 : molecule_set, .FALSE., particle_set=particle_set)
787 78 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
788 : CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
789 78 : subcells=subcells, operator_type="PP", nlname="sab_cn")
790 78 : disp%sab_cn => sab_cn
791 78 : DEALLOCATE (c_radius, pair_radius, default_present)
792 78 : CALL atom2d_cleanup(atom2d)
793 :
794 : ! Calculate coordination numbers
795 78 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces, disp_env=disp)
796 :
797 78 : CALL qs_dispersion_release(disp)
798 :
799 156 : END SUBROUTINE get_cnumbers
800 :
801 : ! **************************************************************************************************
802 : !> \brief ...
803 : !> \param charges ...
804 : !> \param lambda ...
805 : !> \param eeq_energy ...
806 : !> \param particle_set ...
807 : !> \param kind_of ...
808 : !> \param cell ...
809 : !> \param chia ...
810 : !> \param gam ...
811 : !> \param gab ...
812 : !> \param para_env ...
813 : !> \param blacs_env ...
814 : !> \param dft_control ...
815 : !> \param eeq_sparam ...
816 : !> \param totalcharge ...
817 : !> \param ewald ...
818 : !> \param ewald_env ...
819 : !> \param ewald_pw ...
820 : !> \param iounit ...
821 : ! **************************************************************************************************
822 3756 : SUBROUTINE eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, &
823 3756 : chia, gam, gab, para_env, blacs_env, dft_control, eeq_sparam, &
824 : totalcharge, ewald, ewald_env, ewald_pw, iounit)
825 :
826 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
827 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
828 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
829 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
830 : TYPE(cell_type), POINTER :: cell
831 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
832 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
833 : TYPE(mp_para_env_type), POINTER :: para_env
834 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
835 : TYPE(dft_control_type), POINTER :: dft_control
836 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
837 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: totalcharge
838 : LOGICAL, INTENT(IN), OPTIONAL :: ewald
839 : TYPE(ewald_environment_type), OPTIONAL, POINTER :: ewald_env
840 : TYPE(ewald_pw_type), OPTIONAL, POINTER :: ewald_pw
841 : INTEGER, INTENT(IN), OPTIONAL :: iounit
842 :
843 : CHARACTER(len=*), PARAMETER :: routineN = 'eeq_solver'
844 :
845 : INTEGER :: handle, ierror, iunit, natom, nkind, ns
846 : LOGICAL :: do_direct, do_displ, do_ewald, do_sparse
847 : REAL(KIND=dp) :: alpha, deth, ftime, qtot
848 : TYPE(cp_fm_struct_type), POINTER :: mat_struct
849 : TYPE(cp_fm_type) :: eeq_mat
850 :
851 1878 : CALL timeset(routineN, handle)
852 :
853 1878 : do_ewald = .FALSE.
854 1878 : IF (PRESENT(ewald)) do_ewald = ewald
855 : !
856 1878 : qtot = 0.0_dp
857 1878 : IF (PRESENT(totalcharge)) qtot = totalcharge
858 : !
859 1878 : iunit = -1
860 1878 : IF (PRESENT(iounit)) iunit = iounit
861 :
862 : ! EEQ solver parameters
863 1878 : do_direct = eeq_sparam%direct
864 1878 : do_sparse = eeq_sparam%sparse
865 :
866 1878 : do_displ = .FALSE.
867 1878 : IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
868 200 : do_displ = dft_control%period_efield%displacement_field
869 : END IF
870 :
871 : ! sparse option NYA
872 1878 : IF (do_sparse) THEN
873 0 : CALL cp_abort(__LOCATION__, "EEQ: Sparse option not yet available")
874 : END IF
875 : !
876 1878 : natom = SIZE(particle_set)
877 1878 : nkind = SIZE(gam)
878 1878 : ns = natom + 1
879 : CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
880 1878 : nrow_global=ns, ncol_global=ns)
881 1878 : CALL cp_fm_create(eeq_mat, mat_struct)
882 1878 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
883 : !
884 1878 : IF (do_ewald) THEN
885 864 : CPASSERT(PRESENT(ewald_env))
886 864 : CPASSERT(PRESENT(ewald_pw))
887 864 : IF (do_direct) THEN
888 0 : IF (do_displ) THEN
889 0 : CPABORT("NYA")
890 : ELSE
891 : CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
892 : kind_of, cell, chia, gam, gab, qtot, &
893 0 : ewald_env, ewald_pw, iounit)
894 : END IF
895 : ELSE
896 864 : IF (do_displ) THEN
897 0 : CPABORT("NYA")
898 : ELSE
899 : ierror = 0
900 : CALL pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
901 : kind_of, cell, chia, gam, gab, qtot, &
902 864 : ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
903 864 : IF (ierror /= 0) THEN
904 : ! backup to non-iterative method
905 : CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
906 : kind_of, cell, chia, gam, gab, qtot, &
907 730 : ewald_env, ewald_pw, iounit)
908 : END IF
909 : END IF
910 : END IF
911 864 : IF (qtot /= 0._dp) THEN
912 104 : CALL get_cell(cell=cell, deth=deth)
913 104 : CALL ewald_env_get(ewald_env, alpha=alpha)
914 104 : eeq_energy = eeq_energy - 0.5_dp*qtot**2/alpha**2/deth
915 : END IF
916 : ELSE
917 : CALL mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, &
918 1014 : cell, chia, gam, gab, qtot, ftime)
919 1014 : IF (iounit > 0) THEN
920 997 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
921 : END IF
922 : END IF
923 1878 : CALL cp_fm_struct_release(mat_struct)
924 1878 : CALL cp_fm_release(eeq_mat)
925 :
926 1878 : CALL timestop(handle)
927 :
928 1878 : END SUBROUTINE eeq_solver
929 :
930 : ! **************************************************************************************************
931 : !> \brief ...
932 : !> \param charges ...
933 : !> \param lambda ...
934 : !> \param eeq_energy ...
935 : !> \param eeq_mat ...
936 : !> \param particle_set ...
937 : !> \param kind_of ...
938 : !> \param cell ...
939 : !> \param chia ...
940 : !> \param gam ...
941 : !> \param gab ...
942 : !> \param qtot ...
943 : !> \param ftime ...
944 : ! **************************************************************************************************
945 1878 : SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
946 1878 : chia, gam, gab, qtot, ftime)
947 :
948 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
949 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
950 : TYPE(cp_fm_type) :: eeq_mat
951 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
952 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
953 : TYPE(cell_type), POINTER :: cell
954 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
955 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
956 : REAL(KIND=dp), INTENT(IN) :: qtot
957 : REAL(KIND=dp), INTENT(OUT) :: ftime
958 :
959 : CHARACTER(len=*), PARAMETER :: routineN = 'mi_solver'
960 :
961 : INTEGER :: handle, ia, iac, iar, ic, ikind, ir, &
962 : jkind, natom, ncloc, ncvloc, nkind, &
963 : nrloc, nrvloc, ns
964 1878 : INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
965 : REAL(KIND=dp) :: dr, grc, te, ti, xr
966 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
967 : TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
968 : TYPE(cp_fm_type) :: rhs_vec
969 : TYPE(mp_para_env_type), POINTER :: para_env
970 :
971 1878 : CALL timeset(routineN, handle)
972 1878 : ti = m_walltime()
973 :
974 1878 : natom = SIZE(particle_set)
975 1878 : nkind = SIZE(gam)
976 : !
977 1878 : ns = natom + 1
978 1878 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
979 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
980 1878 : row_indices=rind, col_indices=cind)
981 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
982 1878 : nrow_global=ns, ncol_global=1)
983 1878 : CALL cp_fm_create(rhs_vec, vec_struct)
984 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
985 1878 : row_indices=rvind, col_indices=cvind)
986 : !
987 : ! set up matrix
988 1878 : CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
989 1878 : CALL cp_fm_set_all(rhs_vec, 0.0_dp)
990 7299 : DO ir = 1, nrloc
991 5421 : iar = rind(ir)
992 5421 : IF (iar > natom) CYCLE
993 4482 : ikind = kind_of(iar)
994 17928 : ri(1:3) = particle_set(iar)%r(1:3)
995 40974 : DO ic = 1, ncloc
996 34614 : iac = cind(ic)
997 34614 : IF (iac > natom) CYCLE
998 30132 : jkind = kind_of(iac)
999 120528 : rj(1:3) = particle_set(iac)%r(1:3)
1000 30132 : IF (iar == iac) THEN
1001 4482 : grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1002 : ELSE
1003 102600 : rij(1:3) = ri(1:3) - rj(1:3)
1004 102600 : rij = pbc(rij, cell)
1005 102600 : dr = SQRT(SUM(rij**2))
1006 25650 : grc = erf(gab(ikind, jkind)*dr)/dr
1007 : END IF
1008 40035 : eeq_mat%local_data(ir, ic) = grc
1009 : END DO
1010 : END DO
1011 : ! set up rhs vector
1012 7299 : DO ir = 1, nrvloc
1013 5421 : iar = rvind(ir)
1014 12720 : DO ic = 1, ncvloc
1015 5421 : iac = cvind(ic)
1016 5421 : ia = MAX(iar, iac)
1017 5421 : IF (ia > natom) THEN
1018 939 : xr = qtot
1019 : ELSE
1020 4482 : xr = -chia(ia)
1021 : END IF
1022 10842 : rhs_vec%local_data(ir, ic) = xr
1023 : END DO
1024 : END DO
1025 : !
1026 1878 : CALL cp_fm_solve(eeq_mat, rhs_vec)
1027 : !
1028 10842 : charges = 0.0_dp
1029 1878 : lambda = 0.0_dp
1030 7299 : DO ir = 1, nrvloc
1031 5421 : iar = rvind(ir)
1032 12720 : DO ic = 1, ncvloc
1033 5421 : iac = cvind(ic)
1034 5421 : ia = MAX(iar, iac)
1035 10842 : IF (ia <= natom) THEN
1036 4482 : xr = rhs_vec%local_data(ir, ic)
1037 4482 : charges(ia) = xr
1038 : ELSE
1039 939 : lambda = rhs_vec%local_data(ir, ic)
1040 : END IF
1041 : END DO
1042 : END DO
1043 1878 : CALL para_env%sum(lambda)
1044 19806 : CALL para_env%sum(charges)
1045 : !
1046 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
1047 10842 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1048 :
1049 1878 : CALL cp_fm_struct_release(vec_struct)
1050 1878 : CALL cp_fm_release(rhs_vec)
1051 :
1052 1878 : te = m_walltime()
1053 1878 : ftime = te - ti
1054 1878 : CALL timestop(handle)
1055 :
1056 1878 : END SUBROUTINE mi_solver
1057 :
1058 : ! **************************************************************************************************
1059 : !> \brief ...
1060 : !> \param charges ...
1061 : !> \param lambda ...
1062 : !> \param eeq_energy ...
1063 : !> \param eeq_mat ...
1064 : !> \param particle_set ...
1065 : !> \param kind_of ...
1066 : !> \param cell ...
1067 : !> \param chia ...
1068 : !> \param gam ...
1069 : !> \param gab ...
1070 : !> \param qtot ...
1071 : !> \param ewald_env ...
1072 : !> \param ewald_pw ...
1073 : !> \param eeq_sparam ...
1074 : !> \param ierror ...
1075 : !> \param iounit ...
1076 : ! **************************************************************************************************
1077 864 : SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1078 864 : kind_of, cell, chia, gam, gab, qtot, &
1079 : ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
1080 :
1081 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
1082 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
1083 : TYPE(cp_fm_type) :: eeq_mat
1084 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1085 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1086 : TYPE(cell_type), POINTER :: cell
1087 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1088 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
1089 : REAL(KIND=dp), INTENT(IN) :: qtot
1090 : TYPE(ewald_environment_type), POINTER :: ewald_env
1091 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1092 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
1093 : INTEGER, INTENT(OUT) :: ierror
1094 : INTEGER, OPTIONAL :: iounit
1095 :
1096 : CHARACTER(len=*), PARAMETER :: routineN = 'pbc_solver'
1097 :
1098 : INTEGER :: ewald_type, handle, i, iac, iar, ic, ikind, info, ir, iunit, iv, ix, iy, iz, &
1099 : jkind, max_diis, mdiis, natom, ncloc, ndiis, nkind, now, nrloc, ns, sdiis
1100 : INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1101 864 : INTEGER, DIMENSION(:), POINTER :: cind, rind
1102 : REAL(KIND=dp) :: ad, alpha, astep, deth, dr, eeqn, &
1103 : eps_diis, ftime, grc1, grc2, rcut, &
1104 : res, resin, rmax, te, ti
1105 864 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1106 864 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat, fvec, vmat, xvec
1107 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rijl, rj
1108 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1109 864 : REAL(KIND=dp), DIMENSION(:), POINTER :: rhs, rv0, xv0
1110 : TYPE(cp_fm_struct_type), POINTER :: mat_struct
1111 : TYPE(cp_fm_type) :: mmat, pmat
1112 : TYPE(mp_para_env_type), POINTER :: para_env
1113 :
1114 864 : CALL timeset(routineN, handle)
1115 864 : ti = m_walltime()
1116 :
1117 864 : ierror = 0
1118 :
1119 864 : iunit = -1
1120 864 : IF (PRESENT(iounit)) iunit = iounit
1121 :
1122 864 : natom = SIZE(particle_set)
1123 864 : nkind = SIZE(gam)
1124 : !
1125 864 : CALL get_cell(cell=cell, deth=deth)
1126 864 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1127 864 : ad = 2.0_dp*alpha*oorootpi
1128 864 : IF (ewald_type /= do_ewald_spme) THEN
1129 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1130 : END IF
1131 : !
1132 864 : rmax = 2.0_dp*rcut
1133 : ! max cells used
1134 864 : CALL get_cell(cell, h=hmat, periodic=periodic)
1135 864 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1136 864 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1137 864 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1138 864 : IF (periodic(1) == 0) ncell(1) = 0
1139 864 : IF (periodic(2) == 0) ncell(2) = 0
1140 864 : IF (periodic(3) == 0) ncell(3) = 0
1141 : !
1142 : CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1143 864 : chia, gam, gab, qtot, ftime)
1144 864 : IF (iunit > 0) THEN
1145 831 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1146 : END IF
1147 864 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1148 864 : CALL cp_fm_create(pmat, mat_struct)
1149 864 : CALL cp_fm_create(mmat, mat_struct)
1150 : !
1151 : ! response matrix
1152 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1153 864 : row_indices=rind, col_indices=cind)
1154 864 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1155 3688 : DO ir = 1, nrloc
1156 2824 : iar = rind(ir)
1157 2824 : ri = 0.0_dp
1158 2824 : IF (iar <= natom) THEN
1159 2392 : ikind = kind_of(iar)
1160 9568 : ri(1:3) = particle_set(iar)%r(1:3)
1161 : END IF
1162 28070 : DO ic = 1, ncloc
1163 24382 : iac = cind(ic)
1164 24382 : IF (iac > natom .AND. iar > natom) THEN
1165 432 : eeq_mat%local_data(ir, ic) = 0.0_dp
1166 432 : CYCLE
1167 23950 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1168 4784 : eeq_mat%local_data(ir, ic) = 1.0_dp
1169 4784 : CYCLE
1170 : END IF
1171 19166 : jkind = kind_of(iac)
1172 76664 : rj(1:3) = particle_set(iac)%r(1:3)
1173 76664 : rij(1:3) = ri(1:3) - rj(1:3)
1174 76664 : rij = pbc(rij, cell)
1175 156152 : DO ix = -ncell(1), ncell(1)
1176 1092462 : DO iy = -ncell(2), ncell(2)
1177 7647234 : DO iz = -ncell(3), ncell(3)
1178 26295752 : cvec = [ix, iy, iz]
1179 124904822 : rijl = rij + MATMUL(hmat, cvec)
1180 26295752 : dr = SQRT(SUM(rijl**2))
1181 6573938 : IF (dr > rmax) CYCLE
1182 1570014 : IF (iar == iac .AND. dr < 0.00001_dp) THEN
1183 2392 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1184 : ELSE
1185 1567622 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1186 : END IF
1187 2509148 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1188 : END DO
1189 : END DO
1190 : END DO
1191 : END DO
1192 : END DO
1193 : !
1194 : ! preconditioner
1195 : CALL cp_fm_get_info(pmat, nrow_local=nrloc, ncol_local=ncloc, &
1196 864 : row_indices=rind, col_indices=cind)
1197 864 : CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1198 3688 : DO ir = 1, nrloc
1199 2824 : iar = rind(ir)
1200 2824 : ri = 0.0_dp
1201 2824 : IF (iar <= natom) THEN
1202 2392 : ikind = kind_of(iar)
1203 9568 : ri(1:3) = particle_set(iar)%r(1:3)
1204 : END IF
1205 28070 : DO ic = 1, ncloc
1206 24382 : iac = cind(ic)
1207 24382 : IF (iac > natom .AND. iar > natom) THEN
1208 432 : pmat%local_data(ir, ic) = 0.0_dp
1209 432 : CYCLE
1210 23950 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1211 4784 : pmat%local_data(ir, ic) = 1.0_dp
1212 4784 : CYCLE
1213 : END IF
1214 19166 : jkind = kind_of(iac)
1215 76664 : rj(1:3) = particle_set(iac)%r(1:3)
1216 76664 : rij(1:3) = ri(1:3) - rj(1:3)
1217 76664 : rij = pbc(rij, cell)
1218 19166 : IF (iar == iac) THEN
1219 2392 : grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1220 : ELSE
1221 16774 : grc2 = erf(gab(ikind, jkind)*dr)/dr
1222 : END IF
1223 21990 : pmat%local_data(ir, ic) = grc2
1224 : END DO
1225 : END DO
1226 864 : CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1227 : ! preconditioner invers
1228 864 : CALL cp_fm_invert(pmat, mmat)
1229 : !
1230 : ! rhs
1231 864 : ns = natom + 1
1232 2592 : ALLOCATE (rhs(ns))
1233 5648 : rhs(1:natom) = chia(1:natom)
1234 864 : rhs(ns) = -qtot
1235 : !
1236 2592 : ALLOCATE (xv0(ns), rv0(ns))
1237 : ! initial guess
1238 5648 : xv0(1:natom) = charges(1:natom)
1239 864 : xv0(ns) = 0.0_dp
1240 : ! DIIS optimizer
1241 864 : max_diis = eeq_sparam%max_diis
1242 864 : mdiis = eeq_sparam%mdiis
1243 864 : sdiis = eeq_sparam%sdiis
1244 864 : eps_diis = eeq_sparam%eps_diis
1245 864 : astep = eeq_sparam%alpha
1246 6048 : ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1247 157152 : xvec = 0.0_dp; fvec = 0.0_dp
1248 7776 : ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1249 169344 : dmat = 0.0_dp; dvec = 0.0_dp
1250 864 : ndiis = 1
1251 864 : now = 1
1252 : CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1253 864 : cell, particle_set, xv0, rhs, rv0)
1254 6512 : resin = SQRT(SUM(rv0*rv0))
1255 : !
1256 8578 : DO iv = 1, max_diis
1257 75160 : res = SQRT(SUM(rv0*rv0))
1258 8578 : IF (res > 10._dp*resin) EXIT
1259 7848 : IF (res < eps_diis) EXIT
1260 : !
1261 7714 : now = MOD(iv - 1, mdiis) + 1
1262 7714 : ndiis = MIN(iv, mdiis)
1263 68648 : xvec(1:ns, now) = xv0(1:ns)
1264 68648 : fvec(1:ns, now) = rv0(1:ns)
1265 56454 : DO i = 1, ndiis
1266 477304 : vmat(now, i) = SUM(fvec(:, now)*fvec(:, i))
1267 56454 : vmat(i, now) = vmat(now, i)
1268 : END DO
1269 7714 : IF (ndiis < sdiis) THEN
1270 24192 : xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1271 : ELSE
1272 84028 : dvec = 0.0_dp
1273 6002 : dvec(ndiis + 1) = 1.0_dp
1274 484946 : dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1275 52174 : dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1276 52174 : dmat(1:ndiis, ndiis + 1) = 1.0_dp
1277 6002 : dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1278 6002 : CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1279 693642 : dvec(1:ndiis + 1) = MATMUL(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1280 513572 : xv0(1:ns) = MATMUL(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1281 575270 : xv0(1:ns) = xv0(1:ns) + MATMUL(fvec(1:ns, 1:ndiis), dvec(1:ndiis))
1282 : END IF
1283 : !
1284 : CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1285 8578 : cell, particle_set, xv0, rhs, rv0)
1286 : END DO
1287 5648 : charges(1:natom) = xv0(1:natom)
1288 864 : lambda = xv0(ns)
1289 864 : eeq_energy = eeqn
1290 864 : IF (res > eps_diis) ierror = 1
1291 : !
1292 864 : DEALLOCATE (xvec, fvec, bvec)
1293 864 : DEALLOCATE (vmat, dmat, dvec)
1294 864 : DEALLOCATE (xv0, rv0)
1295 864 : DEALLOCATE (rhs)
1296 864 : CALL cp_fm_release(pmat)
1297 864 : CALL cp_fm_release(mmat)
1298 :
1299 864 : te = m_walltime()
1300 864 : IF (iunit > 0) THEN
1301 831 : IF (ierror == 1) THEN
1302 714 : WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1303 : ELSE
1304 117 : WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1305 : END IF
1306 831 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1307 : END IF
1308 864 : CALL timestop(handle)
1309 :
1310 3456 : END SUBROUTINE pbc_solver
1311 :
1312 : ! **************************************************************************************************
1313 : !> \brief ...
1314 : !> \param charges ...
1315 : !> \param lambda ...
1316 : !> \param eeq_energy ...
1317 : !> \param eeq_mat ...
1318 : !> \param particle_set ...
1319 : !> \param kind_of ...
1320 : !> \param cell ...
1321 : !> \param chia ...
1322 : !> \param gam ...
1323 : !> \param gab ...
1324 : !> \param qtot ...
1325 : !> \param ewald_env ...
1326 : !> \param ewald_pw ...
1327 : !> \param iounit ...
1328 : ! **************************************************************************************************
1329 730 : SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1330 730 : kind_of, cell, chia, gam, gab, qtot, ewald_env, ewald_pw, iounit)
1331 :
1332 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
1333 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
1334 : TYPE(cp_fm_type) :: eeq_mat
1335 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1336 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1337 : TYPE(cell_type), POINTER :: cell
1338 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1339 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
1340 : REAL(KIND=dp), INTENT(IN) :: qtot
1341 : TYPE(ewald_environment_type), POINTER :: ewald_env
1342 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1343 : INTEGER, INTENT(IN), OPTIONAL :: iounit
1344 :
1345 : CHARACTER(len=*), PARAMETER :: routineN = 'fpbc_solver'
1346 :
1347 : INTEGER :: ewald_type, handle, ia, iac, iar, ic, &
1348 : ikind, ir, iunit, ix, iy, iz, jkind, &
1349 : natom, ncloc, ncvloc, nkind, nrloc, &
1350 : nrvloc, ns
1351 : INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1352 730 : INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
1353 : REAL(KIND=dp) :: ad, alpha, deth, dr, grc1, rcut, rmax, &
1354 : te, ti, xr
1355 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rijl, rj
1356 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1357 730 : REAL(KIND=dp), DIMENSION(:), POINTER :: pval, xval
1358 : TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
1359 : TYPE(cp_fm_type) :: rhs_vec
1360 : TYPE(mp_para_env_type), POINTER :: para_env
1361 :
1362 730 : CALL timeset(routineN, handle)
1363 730 : ti = m_walltime()
1364 :
1365 730 : iunit = -1
1366 730 : IF (PRESENT(iounit)) iunit = iounit
1367 :
1368 730 : natom = SIZE(particle_set)
1369 730 : nkind = SIZE(gam)
1370 730 : ns = natom + 1
1371 : !
1372 730 : CALL get_cell(cell=cell, deth=deth)
1373 730 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1374 730 : ad = 2.0_dp*alpha*oorootpi
1375 730 : IF (ewald_type /= do_ewald_spme) THEN
1376 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1377 : END IF
1378 : !
1379 730 : rmax = 2.0_dp*rcut
1380 : ! max cells used
1381 730 : CALL get_cell(cell, h=hmat, periodic=periodic)
1382 730 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1383 730 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1384 730 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1385 730 : IF (periodic(1) == 0) ncell(1) = 0
1386 730 : IF (periodic(2) == 0) ncell(2) = 0
1387 730 : IF (periodic(3) == 0) ncell(3) = 0
1388 : !
1389 730 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1390 730 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1391 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1392 730 : row_indices=rind, col_indices=cind)
1393 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1394 730 : nrow_global=ns, ncol_global=1)
1395 730 : CALL cp_fm_create(rhs_vec, vec_struct)
1396 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1397 730 : row_indices=rvind, col_indices=cvind)
1398 : ! response matrix
1399 3051 : DO ir = 1, nrloc
1400 2321 : iar = rind(ir)
1401 2321 : ri = 0.0_dp
1402 2321 : IF (iar <= natom) THEN
1403 1956 : ikind = kind_of(iar)
1404 7824 : ri(1:3) = particle_set(iar)%r(1:3)
1405 : END IF
1406 19210 : DO ic = 1, ncloc
1407 16159 : iac = cind(ic)
1408 16159 : IF (iac > natom .AND. iar > natom) THEN
1409 365 : eeq_mat%local_data(ir, ic) = 0.0_dp
1410 365 : CYCLE
1411 15794 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1412 3912 : eeq_mat%local_data(ir, ic) = 1.0_dp
1413 3912 : CYCLE
1414 : END IF
1415 11882 : jkind = kind_of(iac)
1416 47528 : rj(1:3) = particle_set(iac)%r(1:3)
1417 47528 : rij(1:3) = ri(1:3) - rj(1:3)
1418 47528 : rij = pbc(rij, cell)
1419 97377 : DO ix = -ncell(1), ncell(1)
1420 677274 : DO iy = -ncell(2), ncell(2)
1421 4740918 : DO iz = -ncell(3), ncell(3)
1422 16302104 : cvec = [ix, iy, iz]
1423 77434994 : rijl = rij + MATMUL(hmat, cvec)
1424 16302104 : dr = SQRT(SUM(rijl**2))
1425 4075526 : IF (dr > rmax) CYCLE
1426 972848 : IF (iar == iac .AND. dr < 0.0001_dp) THEN
1427 1956 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1428 : ELSE
1429 970892 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1430 : END IF
1431 1555066 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1432 : END DO
1433 : END DO
1434 : END DO
1435 : END DO
1436 : END DO
1437 : !
1438 2920 : ALLOCATE (xval(natom), pval(natom))
1439 4642 : DO ia = 1, natom
1440 27676 : xval = 0.0_dp
1441 3912 : xval(ia) = 1.0_dp
1442 3912 : CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1443 : !
1444 18480 : DO ir = 1, nrloc
1445 13838 : iar = rind(ir)
1446 13838 : IF (iar /= ia) CYCLE
1447 19706 : DO ic = 1, ncloc
1448 13838 : iac = cind(ic)
1449 13838 : IF (iac > natom) CYCLE
1450 27676 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1451 : END DO
1452 : END DO
1453 : END DO
1454 730 : DEALLOCATE (xval, pval)
1455 : !
1456 : ! set up rhs vector
1457 3051 : DO ir = 1, nrvloc
1458 2321 : iar = rvind(ir)
1459 5372 : DO ic = 1, ncvloc
1460 2321 : iac = cvind(ic)
1461 2321 : ia = MAX(iar, iac)
1462 2321 : IF (ia > natom) THEN
1463 365 : xr = qtot
1464 : ELSE
1465 1956 : xr = -chia(ia)
1466 : END IF
1467 4642 : rhs_vec%local_data(ir, ic) = xr
1468 : END DO
1469 : END DO
1470 : !
1471 730 : CALL cp_fm_solve(eeq_mat, rhs_vec)
1472 : !
1473 4642 : charges = 0.0_dp
1474 730 : lambda = 0.0_dp
1475 3051 : DO ir = 1, nrvloc
1476 2321 : iar = rvind(ir)
1477 5372 : DO ic = 1, ncvloc
1478 2321 : iac = cvind(ic)
1479 2321 : ia = MAX(iar, iac)
1480 4642 : IF (ia <= natom) THEN
1481 1956 : xr = rhs_vec%local_data(ir, ic)
1482 1956 : charges(ia) = xr
1483 : ELSE
1484 365 : lambda = rhs_vec%local_data(ir, ic)
1485 : END IF
1486 : END DO
1487 : END DO
1488 730 : CALL para_env%sum(lambda)
1489 8554 : CALL para_env%sum(charges)
1490 : !
1491 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
1492 4642 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1493 :
1494 730 : CALL cp_fm_struct_release(vec_struct)
1495 730 : CALL cp_fm_release(rhs_vec)
1496 :
1497 730 : te = m_walltime()
1498 730 : IF (iunit > 0) THEN
1499 714 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
1500 : END IF
1501 730 : CALL timestop(handle)
1502 :
1503 2920 : END SUBROUTINE fpbc_solver
1504 :
1505 : ! **************************************************************************************************
1506 : !> \brief ...
1507 : !> \param ewald_env ...
1508 : !> \param ewald_pw ...
1509 : !> \param cell ...
1510 : !> \param particle_set ...
1511 : !> \param charges ...
1512 : !> \param potential ...
1513 : ! **************************************************************************************************
1514 3912 : SUBROUTINE apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, potential)
1515 : TYPE(ewald_environment_type), POINTER :: ewald_env
1516 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1517 : TYPE(cell_type), POINTER :: cell
1518 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1519 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1520 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential
1521 :
1522 : TYPE(mp_para_env_type), POINTER :: para_env
1523 :
1524 3912 : CALL ewald_env_get(ewald_env, para_env=para_env)
1525 27676 : potential = 0.0_dp
1526 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
1527 3912 : particle_set, potential)
1528 51440 : CALL para_env%sum(potential)
1529 :
1530 3912 : END SUBROUTINE apply_potential
1531 :
1532 : ! **************************************************************************************************
1533 : !> \brief ...
1534 : !> \param eeqn ...
1535 : !> \param fm_mat ...
1536 : !> \param mmat ...
1537 : !> \param ewald_env ...
1538 : !> \param ewald_pw ...
1539 : !> \param cell ...
1540 : !> \param particle_set ...
1541 : !> \param charges ...
1542 : !> \param rhs ...
1543 : !> \param potential ...
1544 : ! **************************************************************************************************
1545 8578 : SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
1546 8578 : cell, particle_set, charges, rhs, potential)
1547 : REAL(KIND=dp), INTENT(INOUT) :: eeqn
1548 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat, mmat
1549 : TYPE(ewald_environment_type), POINTER :: ewald_env
1550 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1551 : TYPE(cell_type), POINTER :: cell
1552 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1553 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1554 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rhs
1555 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential
1556 :
1557 : INTEGER :: na, ns
1558 : REAL(KIND=dp) :: lambda
1559 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mvec
1560 : TYPE(mp_para_env_type), POINTER :: para_env
1561 :
1562 8578 : ns = SIZE(charges)
1563 8578 : na = ns - 1
1564 8578 : CALL ewald_env_get(ewald_env, para_env=para_env)
1565 75160 : potential = 0.0_dp
1566 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
1567 8578 : particle_set, potential(1:na))
1568 124586 : CALL para_env%sum(potential(1:na))
1569 8578 : CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
1570 133164 : eeqn = 0.5_dp*SUM(charges(1:na)*potential(1:na)) + SUM(charges(1:na)*rhs(1:na))
1571 75160 : potential(1:ns) = potential(1:ns) + rhs(1:ns)
1572 25734 : ALLOCATE (mvec(ns))
1573 8578 : CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
1574 66582 : lambda = -SUM(mvec(1:na))/REAL(na, KIND=dp)
1575 66582 : potential(1:na) = mvec(1:na) + lambda
1576 8578 : DEALLOCATE (mvec)
1577 :
1578 8578 : END SUBROUTINE get_energy_gradient
1579 :
1580 : ! **************************************************************************************************
1581 : !> \brief ...
1582 : !> \param qs_env ...
1583 : !> \param charges ...
1584 : !> \param ef_energy ...
1585 : ! **************************************************************************************************
1586 328 : SUBROUTINE eeq_efield_energy(qs_env, charges, ef_energy)
1587 : TYPE(qs_environment_type), POINTER :: qs_env
1588 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges
1589 : REAL(KIND=dp), INTENT(OUT) :: ef_energy
1590 :
1591 : COMPLEX(KIND=dp) :: zdeta
1592 : COMPLEX(KIND=dp), DIMENSION(3) :: zi
1593 : INTEGER :: ia, idir, natom
1594 : LOGICAL :: dfield
1595 : REAL(KIND=dp) :: kr, omega, q
1596 : REAL(KIND=dp), DIMENSION(3) :: ci, dfilter, fieldpol, fpolvec, kvec, &
1597 : qi, ria
1598 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1599 : TYPE(cell_type), POINTER :: cell
1600 : TYPE(dft_control_type), POINTER :: dft_control
1601 328 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1602 :
1603 328 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1604 328 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1605 :
1606 328 : IF (dft_control%apply_period_efield) THEN
1607 164 : dfield = dft_control%period_efield%displacement_field
1608 :
1609 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1610 0 : CPABORT("use of strength_list not implemented for eeq_efield_energy")
1611 : END IF
1612 :
1613 656 : fieldpol = dft_control%period_efield%polarisation
1614 1148 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1615 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
1616 2132 : hmat = cell%hmat(:, :)/twopi
1617 656 : DO idir = 1, 3
1618 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1619 656 : + fieldpol(3)*hmat(3, idir)
1620 : END DO
1621 :
1622 656 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1623 820 : DO ia = 1, natom
1624 656 : q = charges(ia)
1625 2624 : ria = particle_set(ia)%r
1626 2624 : ria = pbc(ria, cell)
1627 2788 : DO idir = 1, 3
1628 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
1629 7872 : kr = SUM(kvec(:)*ria(:))
1630 1968 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
1631 2624 : zi(idir) = zi(idir)*zdeta
1632 : END DO
1633 : END DO
1634 656 : qi = AIMAG(LOG(zi))
1635 164 : IF (dfield) THEN
1636 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1637 0 : omega = cell%deth
1638 0 : ci = MATMUL(hmat, qi)/omega
1639 0 : ef_energy = 0.0_dp
1640 0 : DO idir = 1, 3
1641 0 : ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
1642 : END DO
1643 0 : ef_energy = -0.25_dp*omega/twopi*ef_energy
1644 : ELSE
1645 656 : ef_energy = SUM(fpolvec(:)*qi(:))
1646 : END IF
1647 :
1648 164 : ELSE IF (dft_control%apply_efield) THEN
1649 :
1650 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1651 656 : dft_control%efield_fields(1)%efield%strength
1652 :
1653 164 : ef_energy = 0.0_dp
1654 820 : DO ia = 1, natom
1655 2624 : ria = particle_set(ia)%r
1656 2624 : ria = pbc(ria, cell)
1657 656 : q = charges(ia)
1658 2788 : ef_energy = ef_energy - q*SUM(fieldpol*ria)
1659 : END DO
1660 :
1661 : ELSE
1662 0 : CPABORT("apply field")
1663 : END IF
1664 :
1665 328 : END SUBROUTINE eeq_efield_energy
1666 :
1667 : ! **************************************************************************************************
1668 : !> \brief ...
1669 : !> \param qs_env ...
1670 : !> \param efr ...
1671 : ! **************************************************************************************************
1672 328 : SUBROUTINE eeq_efield_pot(qs_env, efr)
1673 : TYPE(qs_environment_type), POINTER :: qs_env
1674 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: efr
1675 :
1676 : INTEGER :: ia, idir, natom
1677 : LOGICAL :: dfield
1678 : REAL(KIND=dp) :: kr
1679 : REAL(KIND=dp), DIMENSION(3) :: fieldpol, fpolvec, kvec, ria
1680 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1681 : TYPE(cell_type), POINTER :: cell
1682 : TYPE(dft_control_type), POINTER :: dft_control
1683 328 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1684 :
1685 328 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1686 328 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1687 :
1688 328 : IF (dft_control%apply_period_efield) THEN
1689 164 : dfield = dft_control%period_efield%displacement_field
1690 :
1691 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1692 0 : CPABORT("use of strength_list not implemented for eeq_efield_pot")
1693 : END IF
1694 :
1695 656 : fieldpol = dft_control%period_efield%polarisation
1696 1148 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1697 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
1698 2132 : hmat = cell%hmat(:, :)/twopi
1699 656 : DO idir = 1, 3
1700 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1701 656 : + fieldpol(3)*hmat(3, idir)
1702 : END DO
1703 :
1704 164 : IF (dfield) THEN
1705 : ! dE/dq depends on q, postpone calculation
1706 0 : efr = 0.0_dp
1707 : ELSE
1708 820 : efr = 0.0_dp
1709 820 : DO ia = 1, natom
1710 2624 : ria = particle_set(ia)%r
1711 2624 : ria = pbc(ria, cell)
1712 2788 : DO idir = 1, 3
1713 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
1714 7872 : kr = SUM(kvec(:)*ria(:))
1715 2624 : efr(ia) = efr(ia) + kr*fpolvec(idir)
1716 : END DO
1717 : END DO
1718 : END IF
1719 :
1720 164 : ELSE IF (dft_control%apply_efield) THEN
1721 :
1722 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1723 656 : dft_control%efield_fields(1)%efield%strength
1724 :
1725 820 : DO ia = 1, natom
1726 2624 : ria = particle_set(ia)%r
1727 2624 : ria = pbc(ria, cell)
1728 2788 : efr(ia) = -SUM(fieldpol*ria)
1729 : END DO
1730 :
1731 : ELSE
1732 0 : CPABORT("apply field")
1733 : END IF
1734 :
1735 328 : END SUBROUTINE eeq_efield_pot
1736 :
1737 : ! **************************************************************************************************
1738 : !> \brief ...
1739 : !> \param charges ...
1740 : !> \param dft_control ...
1741 : !> \param particle_set ...
1742 : !> \param cell ...
1743 : !> \param efr ...
1744 : ! **************************************************************************************************
1745 0 : SUBROUTINE eeq_dfield_pot(charges, dft_control, particle_set, cell, efr)
1746 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges
1747 : TYPE(dft_control_type), POINTER :: dft_control
1748 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1749 : TYPE(cell_type), POINTER :: cell
1750 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: efr
1751 :
1752 : COMPLEX(KIND=dp) :: zdeta
1753 : COMPLEX(KIND=dp), DIMENSION(3) :: zi
1754 : INTEGER :: ia, idir, natom
1755 : REAL(KIND=dp) :: kr, omega, q
1756 : REAL(KIND=dp), DIMENSION(3) :: ci, dfilter, fieldpol, kvec, qi, ria
1757 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1758 :
1759 0 : natom = SIZE(particle_set)
1760 :
1761 0 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1762 0 : CPABORT("use of strength_list not implemented for eeq_dfield_pot")
1763 : END IF
1764 :
1765 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1766 0 : fieldpol = dft_control%period_efield%polarisation
1767 0 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1768 0 : fieldpol = -fieldpol*dft_control%period_efield%strength
1769 0 : hmat = cell%hmat(:, :)/twopi
1770 0 : omega = cell%deth
1771 : !
1772 0 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1773 0 : DO ia = 1, natom
1774 0 : q = charges(ia)
1775 0 : ria = particle_set(ia)%r
1776 0 : ria = pbc(ria, cell)
1777 0 : DO idir = 1, 3
1778 0 : kvec(:) = twopi*cell%h_inv(idir, :)
1779 0 : kr = SUM(kvec(:)*ria(:))
1780 0 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
1781 0 : zi(idir) = zi(idir)*zdeta
1782 : END DO
1783 : END DO
1784 0 : qi = AIMAG(LOG(zi))
1785 0 : ci = MATMUL(hmat, qi)/omega
1786 0 : ci = dfilter*(fieldpol - 2._dp*twopi*ci)
1787 0 : DO ia = 1, natom
1788 0 : ria = particle_set(ia)%r
1789 0 : ria = pbc(ria, cell)
1790 0 : efr(ia) = efr(ia) - SUM(ci*ria)
1791 : END DO
1792 :
1793 0 : END SUBROUTINE eeq_dfield_pot
1794 :
1795 : ! **************************************************************************************************
1796 : !> \brief ...
1797 : !> \param qs_env ...
1798 : !> \param charges ...
1799 : !> \param qlag ...
1800 : ! **************************************************************************************************
1801 8 : SUBROUTINE eeq_efield_force_loc(qs_env, charges, qlag)
1802 : TYPE(qs_environment_type), POINTER :: qs_env
1803 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1804 :
1805 : INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1806 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1807 : REAL(KIND=dp) :: q
1808 : REAL(KIND=dp), DIMENSION(3) :: fieldpol
1809 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1810 : TYPE(cell_type), POINTER :: cell
1811 : TYPE(dft_control_type), POINTER :: dft_control
1812 : TYPE(distribution_1d_type), POINTER :: local_particles
1813 : TYPE(mp_para_env_type), POINTER :: para_env
1814 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1815 8 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1816 :
1817 : CALL get_qs_env(qs_env=qs_env, &
1818 : dft_control=dft_control, &
1819 : cell=cell, particle_set=particle_set, &
1820 : nkind=nkind, natom=natom, &
1821 : para_env=para_env, &
1822 8 : local_particles=local_particles)
1823 :
1824 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1825 32 : dft_control%efield_fields(1)%efield%strength
1826 :
1827 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1828 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1829 8 : CALL get_qs_env(qs_env=qs_env, force=force)
1830 :
1831 32 : DO ikind = 1, nkind
1832 152 : force(ikind)%efield = 0.0_dp
1833 40 : DO ia = 1, local_particles%n_el(ikind)
1834 16 : iatom = local_particles%list(ikind)%array(ia)
1835 16 : q = charges(iatom) - qlag(iatom)
1836 16 : atom_a = atom_of_kind(iatom)
1837 88 : force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
1838 : END DO
1839 288 : CALL para_env%sum(force(ikind)%efield)
1840 : END DO
1841 :
1842 16 : END SUBROUTINE eeq_efield_force_loc
1843 :
1844 : ! **************************************************************************************************
1845 : !> \brief ...
1846 : !> \param qs_env ...
1847 : !> \param charges ...
1848 : !> \param qlag ...
1849 : ! **************************************************************************************************
1850 8 : SUBROUTINE eeq_efield_force_periodic(qs_env, charges, qlag)
1851 : TYPE(qs_environment_type), POINTER :: qs_env
1852 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1853 :
1854 : INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1855 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1856 : LOGICAL :: dfield, use_virial
1857 : REAL(KIND=dp) :: q
1858 : REAL(KIND=dp), DIMENSION(3) :: fa, fieldpol, ria
1859 : REAL(KIND=dp), DIMENSION(3, 3) :: pve
1860 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1861 : TYPE(cell_type), POINTER :: cell
1862 : TYPE(dft_control_type), POINTER :: dft_control
1863 : TYPE(distribution_1d_type), POINTER :: local_particles
1864 : TYPE(mp_para_env_type), POINTER :: para_env
1865 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1866 8 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1867 : TYPE(virial_type), POINTER :: virial
1868 :
1869 : CALL get_qs_env(qs_env=qs_env, &
1870 : dft_control=dft_control, &
1871 : cell=cell, particle_set=particle_set, &
1872 : virial=virial, &
1873 : nkind=nkind, natom=natom, &
1874 : para_env=para_env, &
1875 8 : local_particles=local_particles)
1876 :
1877 8 : dfield = dft_control%period_efield%displacement_field
1878 8 : CPASSERT(.NOT. dfield)
1879 :
1880 8 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1881 0 : CPABORT("use of strength_list not implemented for eeq_efield_force_periodic")
1882 : END IF
1883 :
1884 32 : fieldpol = dft_control%period_efield%polarisation
1885 56 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1886 32 : fieldpol = -fieldpol*dft_control%period_efield%strength
1887 :
1888 8 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1889 :
1890 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1891 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1892 8 : CALL get_qs_env(qs_env=qs_env, force=force)
1893 :
1894 8 : pve = 0.0_dp
1895 32 : DO ikind = 1, nkind
1896 152 : force(ikind)%efield = 0.0_dp
1897 40 : DO ia = 1, local_particles%n_el(ikind)
1898 16 : iatom = local_particles%list(ikind)%array(ia)
1899 16 : q = charges(iatom) - qlag(iatom)
1900 64 : fa(1:3) = q*fieldpol(1:3)
1901 16 : atom_a = atom_of_kind(iatom)
1902 64 : force(ikind)%efield(1:3, atom_a) = fa
1903 40 : IF (use_virial) THEN
1904 0 : ria = particle_set(ia)%r
1905 0 : ria = pbc(ria, cell)
1906 0 : CALL virial_pair_force(pve, -0.5_dp, fa, ria)
1907 0 : CALL virial_pair_force(pve, -0.5_dp, ria, fa)
1908 : END IF
1909 : END DO
1910 288 : CALL para_env%sum(force(ikind)%efield)
1911 : END DO
1912 104 : virial%pv_virial = virial%pv_virial + pve
1913 :
1914 16 : END SUBROUTINE eeq_efield_force_periodic
1915 :
1916 12004 : END MODULE eeq_method
|
