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