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_displ, 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 1856 : do_displ = .FALSE.
795 1856 : IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
796 200 : do_displ = dft_control%period_efield%displacement_field
797 : END IF
798 :
799 : ! sparse option NYA
800 1856 : IF (do_sparse) THEN
801 0 : CALL cp_abort(__LOCATION__, "EEQ: Sparse option not yet available")
802 : END IF
803 : !
804 1856 : natom = SIZE(particle_set)
805 1856 : nkind = SIZE(gam)
806 1856 : ns = natom + 1
807 : CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
808 1856 : nrow_global=ns, ncol_global=ns)
809 1856 : CALL cp_fm_create(eeq_mat, mat_struct)
810 1856 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
811 : !
812 1856 : IF (do_ewald) THEN
813 858 : CPASSERT(PRESENT(ewald_env))
814 858 : CPASSERT(PRESENT(ewald_pw))
815 858 : IF (do_direct) THEN
816 0 : IF (do_displ) THEN
817 0 : CPABORT("NYA")
818 : ELSE
819 : CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
820 : kind_of, cell, chia, gam, gab, qtot, &
821 0 : ewald_env, ewald_pw, iounit)
822 : END IF
823 : ELSE
824 858 : IF (do_displ) THEN
825 0 : CPABORT("NYA")
826 : ELSE
827 : ierror = 0
828 : CALL pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
829 : kind_of, cell, chia, gam, gab, qtot, &
830 858 : ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
831 858 : IF (ierror /= 0) THEN
832 : ! backup to non-iterative method
833 : CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
834 : kind_of, cell, chia, gam, gab, qtot, &
835 730 : ewald_env, ewald_pw, iounit)
836 : END IF
837 : END IF
838 : END IF
839 858 : IF (qtot /= 0._dp) THEN
840 104 : CALL get_cell(cell=cell, deth=deth)
841 104 : CALL ewald_env_get(ewald_env, alpha=alpha)
842 104 : eeq_energy = eeq_energy - 0.5_dp*qtot**2/alpha**2/deth
843 : END IF
844 : ELSE
845 : CALL mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, &
846 998 : cell, chia, gam, gab, qtot, ftime)
847 998 : IF (iounit > 0) THEN
848 998 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
849 : END IF
850 : END IF
851 1856 : CALL cp_fm_struct_release(mat_struct)
852 1856 : CALL cp_fm_release(eeq_mat)
853 :
854 1856 : CALL timestop(handle)
855 :
856 1856 : END SUBROUTINE eeq_solver
857 :
858 : ! **************************************************************************************************
859 : !> \brief ...
860 : !> \param charges ...
861 : !> \param lambda ...
862 : !> \param eeq_energy ...
863 : !> \param eeq_mat ...
864 : !> \param particle_set ...
865 : !> \param kind_of ...
866 : !> \param cell ...
867 : !> \param chia ...
868 : !> \param gam ...
869 : !> \param gab ...
870 : !> \param qtot ...
871 : !> \param ftime ...
872 : ! **************************************************************************************************
873 1856 : SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
874 1856 : chia, gam, gab, qtot, ftime)
875 :
876 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
877 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
878 : TYPE(cp_fm_type) :: eeq_mat
879 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
880 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
881 : TYPE(cell_type), POINTER :: cell
882 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
883 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
884 : REAL(KIND=dp), INTENT(IN) :: qtot
885 : REAL(KIND=dp), INTENT(OUT) :: ftime
886 :
887 : CHARACTER(len=*), PARAMETER :: routineN = 'mi_solver'
888 :
889 : INTEGER :: handle, ia, iac, iar, ic, ikind, ir, &
890 : jkind, natom, ncloc, ncvloc, nkind, &
891 : nrloc, nrvloc, ns
892 1856 : INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
893 : REAL(KIND=dp) :: dr, grc, te, ti, xr
894 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
895 : TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
896 : TYPE(cp_fm_type) :: rhs_vec
897 : TYPE(mp_para_env_type), POINTER :: para_env
898 :
899 1856 : CALL timeset(routineN, handle)
900 1856 : ti = m_walltime()
901 :
902 1856 : natom = SIZE(particle_set)
903 1856 : nkind = SIZE(gam)
904 : !
905 1856 : ns = natom + 1
906 1856 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
907 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
908 1856 : row_indices=rind, col_indices=cind)
909 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
910 1856 : nrow_global=ns, ncol_global=1)
911 1856 : CALL cp_fm_create(rhs_vec, vec_struct)
912 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
913 1856 : row_indices=rvind, col_indices=cvind)
914 : !
915 : ! set up matrix
916 1856 : CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
917 1856 : CALL cp_fm_set_all(rhs_vec, 0.0_dp)
918 7201 : DO ir = 1, nrloc
919 5345 : iar = rind(ir)
920 5345 : IF (iar > natom) CYCLE
921 4417 : ikind = kind_of(iar)
922 17668 : ri(1:3) = particle_set(iar)%r(1:3)
923 40355 : DO ic = 1, ncloc
924 34082 : iac = cind(ic)
925 34082 : IF (iac > natom) CYCLE
926 29665 : jkind = kind_of(iac)
927 118660 : rj(1:3) = particle_set(iac)%r(1:3)
928 29665 : IF (iar == iac) THEN
929 4417 : grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
930 : ELSE
931 100992 : rij(1:3) = ri(1:3) - rj(1:3)
932 100992 : rij = pbc(rij, cell)
933 100992 : dr = SQRT(SUM(rij**2))
934 25248 : grc = erf(gab(ikind, jkind)*dr)/dr
935 : END IF
936 39427 : eeq_mat%local_data(ir, ic) = grc
937 : END DO
938 : END DO
939 : ! set up rhs vector
940 7201 : DO ir = 1, nrvloc
941 5345 : iar = rvind(ir)
942 12546 : DO ic = 1, ncvloc
943 5345 : iac = cvind(ic)
944 5345 : ia = MAX(iar, iac)
945 5345 : IF (ia > natom) THEN
946 928 : xr = qtot
947 : ELSE
948 4417 : xr = -chia(ia)
949 : END IF
950 10690 : rhs_vec%local_data(ir, ic) = xr
951 : END DO
952 : END DO
953 : !
954 1856 : CALL cp_fm_solve(eeq_mat, rhs_vec)
955 : !
956 10690 : charges = 0.0_dp
957 1856 : lambda = 0.0_dp
958 7201 : DO ir = 1, nrvloc
959 5345 : iar = rvind(ir)
960 12546 : DO ic = 1, ncvloc
961 5345 : iac = cvind(ic)
962 5345 : ia = MAX(iar, iac)
963 10690 : IF (ia <= natom) THEN
964 4417 : xr = rhs_vec%local_data(ir, ic)
965 4417 : charges(ia) = xr
966 : ELSE
967 928 : lambda = rhs_vec%local_data(ir, ic)
968 : END IF
969 : END DO
970 : END DO
971 1856 : CALL para_env%sum(lambda)
972 19524 : CALL para_env%sum(charges)
973 : !
974 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
975 10690 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
976 :
977 1856 : CALL cp_fm_struct_release(vec_struct)
978 1856 : CALL cp_fm_release(rhs_vec)
979 :
980 1856 : te = m_walltime()
981 1856 : ftime = te - ti
982 1856 : CALL timestop(handle)
983 :
984 1856 : END SUBROUTINE mi_solver
985 :
986 : ! **************************************************************************************************
987 : !> \brief ...
988 : !> \param charges ...
989 : !> \param lambda ...
990 : !> \param eeq_energy ...
991 : !> \param eeq_mat ...
992 : !> \param particle_set ...
993 : !> \param kind_of ...
994 : !> \param cell ...
995 : !> \param chia ...
996 : !> \param gam ...
997 : !> \param gab ...
998 : !> \param qtot ...
999 : !> \param ewald_env ...
1000 : !> \param ewald_pw ...
1001 : !> \param eeq_sparam ...
1002 : !> \param ierror ...
1003 : !> \param iounit ...
1004 : ! **************************************************************************************************
1005 858 : SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1006 858 : kind_of, cell, chia, gam, gab, qtot, &
1007 : ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
1008 :
1009 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
1010 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
1011 : TYPE(cp_fm_type) :: eeq_mat
1012 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1013 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1014 : TYPE(cell_type), POINTER :: cell
1015 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1016 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
1017 : REAL(KIND=dp), INTENT(IN) :: qtot
1018 : TYPE(ewald_environment_type), POINTER :: ewald_env
1019 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1020 : TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
1021 : INTEGER, INTENT(OUT) :: ierror
1022 : INTEGER, OPTIONAL :: iounit
1023 :
1024 : CHARACTER(len=*), PARAMETER :: routineN = 'pbc_solver'
1025 :
1026 : INTEGER :: ewald_type, handle, i, iac, iar, ic, ikind, info, ir, iunit, iv, ix, iy, iz, &
1027 : jkind, max_diis, mdiis, natom, ncloc, ndiis, nkind, now, nrloc, ns, sdiis
1028 : INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1029 858 : INTEGER, DIMENSION(:), POINTER :: cind, rind
1030 : REAL(KIND=dp) :: ad, alpha, astep, deth, dr, eeqn, &
1031 : eps_diis, ftime, grc1, grc2, rcut, &
1032 : res, resin, rmax, te, ti
1033 858 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1034 858 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat, fvec, vmat, xvec
1035 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rijl, rj
1036 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1037 858 : REAL(KIND=dp), DIMENSION(:), POINTER :: rhs, rv0, xv0
1038 : TYPE(cp_fm_struct_type), POINTER :: mat_struct
1039 : TYPE(cp_fm_type) :: mmat, pmat
1040 : TYPE(mp_para_env_type), POINTER :: para_env
1041 :
1042 858 : CALL timeset(routineN, handle)
1043 858 : ti = m_walltime()
1044 :
1045 858 : ierror = 0
1046 :
1047 858 : iunit = -1
1048 858 : IF (PRESENT(iounit)) iunit = iounit
1049 :
1050 858 : natom = SIZE(particle_set)
1051 858 : nkind = SIZE(gam)
1052 : !
1053 858 : CALL get_cell(cell=cell, deth=deth)
1054 858 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1055 858 : ad = 2.0_dp*alpha*oorootpi
1056 858 : IF (ewald_type /= do_ewald_spme) THEN
1057 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1058 : END IF
1059 : !
1060 858 : rmax = 2.0_dp*rcut
1061 : ! max cells used
1062 858 : CALL get_cell(cell, h=hmat, periodic=periodic)
1063 858 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1064 858 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1065 858 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1066 858 : IF (periodic(1) == 0) ncell(1) = 0
1067 858 : IF (periodic(2) == 0) ncell(2) = 0
1068 858 : IF (periodic(3) == 0) ncell(3) = 0
1069 : !
1070 : CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1071 858 : chia, gam, gab, qtot, ftime)
1072 858 : IF (iunit > 0) THEN
1073 858 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1074 : END IF
1075 858 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1076 858 : CALL cp_fm_create(pmat, mat_struct)
1077 858 : CALL cp_fm_create(mmat, mat_struct)
1078 : !
1079 : ! response matrix
1080 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1081 858 : row_indices=rind, col_indices=cind)
1082 858 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1083 3668 : DO ir = 1, nrloc
1084 2810 : iar = rind(ir)
1085 2810 : ri = 0.0_dp
1086 2810 : IF (iar <= natom) THEN
1087 2381 : ikind = kind_of(iar)
1088 9524 : ri(1:3) = particle_set(iar)%r(1:3)
1089 : END IF
1090 27984 : DO ic = 1, ncloc
1091 24316 : iac = cind(ic)
1092 24316 : IF (iac > natom .AND. iar > natom) THEN
1093 429 : eeq_mat%local_data(ir, ic) = 0.0_dp
1094 429 : CYCLE
1095 23887 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1096 4762 : eeq_mat%local_data(ir, ic) = 1.0_dp
1097 4762 : CYCLE
1098 : END IF
1099 19125 : jkind = kind_of(iac)
1100 76500 : rj(1:3) = particle_set(iac)%r(1:3)
1101 76500 : rij(1:3) = ri(1:3) - rj(1:3)
1102 76500 : rij = pbc(rij, cell)
1103 155810 : DO ix = -ncell(1), ncell(1)
1104 1090125 : DO iy = -ncell(2), ncell(2)
1105 7630875 : DO iz = -ncell(3), ncell(3)
1106 26239500 : cvec = [ix, iy, iz]
1107 124637625 : rijl = rij + MATMUL(hmat, cvec)
1108 26239500 : dr = SQRT(SUM(rijl**2))
1109 6559875 : IF (dr > rmax) CYCLE
1110 1566705 : IF (iar == iac .AND. dr < 0.00001_dp) THEN
1111 2381 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1112 : ELSE
1113 1564324 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1114 : END IF
1115 2503830 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1116 : END DO
1117 : END DO
1118 : END DO
1119 : END DO
1120 : END DO
1121 : !
1122 : ! preconditioner
1123 : CALL cp_fm_get_info(pmat, nrow_local=nrloc, ncol_local=ncloc, &
1124 858 : row_indices=rind, col_indices=cind)
1125 858 : CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1126 3668 : DO ir = 1, nrloc
1127 2810 : iar = rind(ir)
1128 2810 : ri = 0.0_dp
1129 2810 : IF (iar <= natom) THEN
1130 2381 : ikind = kind_of(iar)
1131 9524 : ri(1:3) = particle_set(iar)%r(1:3)
1132 : END IF
1133 27984 : DO ic = 1, ncloc
1134 24316 : iac = cind(ic)
1135 24316 : IF (iac > natom .AND. iar > natom) THEN
1136 429 : pmat%local_data(ir, ic) = 0.0_dp
1137 429 : CYCLE
1138 23887 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1139 4762 : pmat%local_data(ir, ic) = 1.0_dp
1140 4762 : CYCLE
1141 : END IF
1142 19125 : jkind = kind_of(iac)
1143 76500 : rj(1:3) = particle_set(iac)%r(1:3)
1144 76500 : rij(1:3) = ri(1:3) - rj(1:3)
1145 76500 : rij = pbc(rij, cell)
1146 19125 : IF (iar == iac) THEN
1147 2381 : grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1148 : ELSE
1149 16744 : grc2 = erf(gab(ikind, jkind)*dr)/dr
1150 : END IF
1151 21935 : pmat%local_data(ir, ic) = grc2
1152 : END DO
1153 : END DO
1154 858 : CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1155 : ! preconditioner invers
1156 858 : CALL cp_fm_invert(pmat, mmat)
1157 : !
1158 : ! rhs
1159 858 : ns = natom + 1
1160 2574 : ALLOCATE (rhs(ns))
1161 5620 : rhs(1:natom) = chia(1:natom)
1162 858 : rhs(ns) = -qtot
1163 : !
1164 2574 : ALLOCATE (xv0(ns), rv0(ns))
1165 : ! initial guess
1166 5620 : xv0(1:natom) = charges(1:natom)
1167 858 : xv0(ns) = 0.0_dp
1168 : ! DIIS optimizer
1169 858 : max_diis = eeq_sparam%max_diis
1170 858 : mdiis = eeq_sparam%mdiis
1171 858 : sdiis = eeq_sparam%sdiis
1172 858 : eps_diis = eeq_sparam%eps_diis
1173 858 : astep = eeq_sparam%alpha
1174 6006 : ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1175 156330 : xvec = 0.0_dp; fvec = 0.0_dp
1176 7722 : ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1177 168168 : dmat = 0.0_dp; dvec = 0.0_dp
1178 858 : ndiis = 1
1179 858 : now = 1
1180 : CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1181 858 : cell, particle_set, xv0, rhs, rv0)
1182 6478 : resin = SQRT(SUM(rv0*rv0))
1183 : !
1184 8572 : DO iv = 1, max_diis
1185 75160 : res = SQRT(SUM(rv0*rv0))
1186 8572 : IF (res > 10._dp*resin) EXIT
1187 7842 : IF (res < eps_diis) EXIT
1188 : !
1189 7714 : now = MOD(iv - 1, mdiis) + 1
1190 7714 : ndiis = MIN(iv, mdiis)
1191 68682 : xvec(1:ns, now) = xv0(1:ns)
1192 68682 : fvec(1:ns, now) = rv0(1:ns)
1193 56506 : DO i = 1, ndiis
1194 477928 : vmat(now, i) = SUM(fvec(:, now)*fvec(:, i))
1195 56506 : vmat(i, now) = vmat(now, i)
1196 : END DO
1197 7714 : IF (ndiis < sdiis) THEN
1198 24156 : xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1199 : ELSE
1200 84084 : dvec = 0.0_dp
1201 6006 : dvec(ndiis + 1) = 1.0_dp
1202 485674 : dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1203 52236 : dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1204 52236 : dmat(1:ndiis, ndiis + 1) = 1.0_dp
1205 6006 : dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1206 6006 : CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1207 694618 : dvec(1:ndiis + 1) = MATMUL(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1208 514280 : xv0(1:ns) = MATMUL(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1209 576036 : xv0(1:ns) = xv0(1:ns) + MATMUL(fvec(1:ns, 1:ndiis), dvec(1:ndiis))
1210 : END IF
1211 : !
1212 : CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1213 8572 : cell, particle_set, xv0, rhs, rv0)
1214 : END DO
1215 5620 : charges(1:natom) = xv0(1:natom)
1216 858 : lambda = xv0(ns)
1217 858 : eeq_energy = eeqn
1218 858 : IF (res > eps_diis) ierror = 1
1219 : !
1220 858 : DEALLOCATE (xvec, fvec, bvec)
1221 858 : DEALLOCATE (vmat, dmat, dvec)
1222 858 : DEALLOCATE (xv0, rv0)
1223 858 : DEALLOCATE (rhs)
1224 858 : CALL cp_fm_release(pmat)
1225 858 : CALL cp_fm_release(mmat)
1226 :
1227 858 : te = m_walltime()
1228 858 : IF (iunit > 0) THEN
1229 858 : IF (ierror == 1) THEN
1230 730 : WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1231 : ELSE
1232 128 : WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1233 : END IF
1234 858 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1235 : END IF
1236 858 : CALL timestop(handle)
1237 :
1238 3432 : END SUBROUTINE pbc_solver
1239 :
1240 : ! **************************************************************************************************
1241 : !> \brief ...
1242 : !> \param charges ...
1243 : !> \param lambda ...
1244 : !> \param eeq_energy ...
1245 : !> \param eeq_mat ...
1246 : !> \param particle_set ...
1247 : !> \param kind_of ...
1248 : !> \param cell ...
1249 : !> \param chia ...
1250 : !> \param gam ...
1251 : !> \param gab ...
1252 : !> \param qtot ...
1253 : !> \param ewald_env ...
1254 : !> \param ewald_pw ...
1255 : !> \param iounit ...
1256 : ! **************************************************************************************************
1257 730 : SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1258 730 : kind_of, cell, chia, gam, gab, qtot, ewald_env, ewald_pw, iounit)
1259 :
1260 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
1261 : REAL(KIND=dp), INTENT(INOUT) :: lambda, eeq_energy
1262 : TYPE(cp_fm_type) :: eeq_mat
1263 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1264 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1265 : TYPE(cell_type), POINTER :: cell
1266 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1267 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gab
1268 : REAL(KIND=dp), INTENT(IN) :: qtot
1269 : TYPE(ewald_environment_type), POINTER :: ewald_env
1270 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1271 : INTEGER, INTENT(IN), OPTIONAL :: iounit
1272 :
1273 : CHARACTER(len=*), PARAMETER :: routineN = 'fpbc_solver'
1274 :
1275 : INTEGER :: ewald_type, handle, ia, iac, iar, ic, &
1276 : ikind, ir, iunit, ix, iy, iz, jkind, &
1277 : natom, ncloc, ncvloc, nkind, nrloc, &
1278 : nrvloc, ns
1279 : INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1280 730 : INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
1281 : REAL(KIND=dp) :: ad, alpha, deth, dr, grc1, rcut, rmax, &
1282 : te, ti, xr
1283 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rijl, rj
1284 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1285 730 : REAL(KIND=dp), DIMENSION(:), POINTER :: pval, xval
1286 : TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
1287 : TYPE(cp_fm_type) :: rhs_vec
1288 : TYPE(mp_para_env_type), POINTER :: para_env
1289 :
1290 730 : CALL timeset(routineN, handle)
1291 730 : ti = m_walltime()
1292 :
1293 730 : iunit = -1
1294 730 : IF (PRESENT(iounit)) iunit = iounit
1295 :
1296 730 : natom = SIZE(particle_set)
1297 730 : nkind = SIZE(gam)
1298 730 : ns = natom + 1
1299 : !
1300 730 : CALL get_cell(cell=cell, deth=deth)
1301 730 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1302 730 : ad = 2.0_dp*alpha*oorootpi
1303 730 : IF (ewald_type /= do_ewald_spme) THEN
1304 0 : CALL cp_abort(__LOCATION__, "Only SPME Ewald method available with EEQ.")
1305 : END IF
1306 : !
1307 730 : rmax = 2.0_dp*rcut
1308 : ! max cells used
1309 730 : CALL get_cell(cell, h=hmat, periodic=periodic)
1310 730 : ncell(1) = CEILING(rmax/plane_distance(1, 0, 0, cell))
1311 730 : ncell(2) = CEILING(rmax/plane_distance(0, 1, 0, cell))
1312 730 : ncell(3) = CEILING(rmax/plane_distance(0, 0, 1, cell))
1313 730 : IF (periodic(1) == 0) ncell(1) = 0
1314 730 : IF (periodic(2) == 0) ncell(2) = 0
1315 730 : IF (periodic(3) == 0) ncell(3) = 0
1316 : !
1317 730 : CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1318 730 : CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1319 : CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1320 730 : row_indices=rind, col_indices=cind)
1321 : CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1322 730 : nrow_global=ns, ncol_global=1)
1323 730 : CALL cp_fm_create(rhs_vec, vec_struct)
1324 : CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1325 730 : row_indices=rvind, col_indices=cvind)
1326 : ! response matrix
1327 3051 : DO ir = 1, nrloc
1328 2321 : iar = rind(ir)
1329 2321 : ri = 0.0_dp
1330 2321 : IF (iar <= natom) THEN
1331 1956 : ikind = kind_of(iar)
1332 7824 : ri(1:3) = particle_set(iar)%r(1:3)
1333 : END IF
1334 19210 : DO ic = 1, ncloc
1335 16159 : iac = cind(ic)
1336 16159 : IF (iac > natom .AND. iar > natom) THEN
1337 365 : eeq_mat%local_data(ir, ic) = 0.0_dp
1338 365 : CYCLE
1339 15794 : ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1340 3912 : eeq_mat%local_data(ir, ic) = 1.0_dp
1341 3912 : CYCLE
1342 : END IF
1343 11882 : jkind = kind_of(iac)
1344 47528 : rj(1:3) = particle_set(iac)%r(1:3)
1345 47528 : rij(1:3) = ri(1:3) - rj(1:3)
1346 47528 : rij = pbc(rij, cell)
1347 97377 : DO ix = -ncell(1), ncell(1)
1348 677274 : DO iy = -ncell(2), ncell(2)
1349 4740918 : DO iz = -ncell(3), ncell(3)
1350 16302104 : cvec = [ix, iy, iz]
1351 77434994 : rijl = rij + MATMUL(hmat, cvec)
1352 16302104 : dr = SQRT(SUM(rijl**2))
1353 4075526 : IF (dr > rmax) CYCLE
1354 972848 : IF (iar == iac .AND. dr < 0.0001_dp) THEN
1355 1956 : grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1356 : ELSE
1357 970892 : grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1358 : END IF
1359 1555066 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1360 : END DO
1361 : END DO
1362 : END DO
1363 : END DO
1364 : END DO
1365 : !
1366 2920 : ALLOCATE (xval(natom), pval(natom))
1367 4642 : DO ia = 1, natom
1368 27676 : xval = 0.0_dp
1369 3912 : xval(ia) = 1.0_dp
1370 3912 : CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1371 : !
1372 18480 : DO ir = 1, nrloc
1373 13838 : iar = rind(ir)
1374 13838 : IF (iar /= ia) CYCLE
1375 19706 : DO ic = 1, ncloc
1376 13838 : iac = cind(ic)
1377 13838 : IF (iac > natom) CYCLE
1378 27676 : eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1379 : END DO
1380 : END DO
1381 : END DO
1382 730 : DEALLOCATE (xval, pval)
1383 : !
1384 : ! set up rhs vector
1385 3051 : DO ir = 1, nrvloc
1386 2321 : iar = rvind(ir)
1387 5372 : DO ic = 1, ncvloc
1388 2321 : iac = cvind(ic)
1389 2321 : ia = MAX(iar, iac)
1390 2321 : IF (ia > natom) THEN
1391 365 : xr = qtot
1392 : ELSE
1393 1956 : xr = -chia(ia)
1394 : END IF
1395 4642 : rhs_vec%local_data(ir, ic) = xr
1396 : END DO
1397 : END DO
1398 : !
1399 730 : CALL cp_fm_solve(eeq_mat, rhs_vec)
1400 : !
1401 4642 : charges = 0.0_dp
1402 730 : lambda = 0.0_dp
1403 3051 : DO ir = 1, nrvloc
1404 2321 : iar = rvind(ir)
1405 5372 : DO ic = 1, ncvloc
1406 2321 : iac = cvind(ic)
1407 2321 : ia = MAX(iar, iac)
1408 4642 : IF (ia <= natom) THEN
1409 1956 : xr = rhs_vec%local_data(ir, ic)
1410 1956 : charges(ia) = xr
1411 : ELSE
1412 365 : lambda = rhs_vec%local_data(ir, ic)
1413 : END IF
1414 : END DO
1415 : END DO
1416 730 : CALL para_env%sum(lambda)
1417 8554 : CALL para_env%sum(charges)
1418 : !
1419 : ! energy: 0.5*(q^T.X - lambda*totalcharge)
1420 4642 : eeq_energy = 0.5*SUM(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1421 :
1422 730 : CALL cp_fm_struct_release(vec_struct)
1423 730 : CALL cp_fm_release(rhs_vec)
1424 :
1425 730 : te = m_walltime()
1426 730 : IF (iunit > 0) THEN
1427 730 : WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
1428 : END IF
1429 730 : CALL timestop(handle)
1430 :
1431 2920 : END SUBROUTINE fpbc_solver
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief ...
1435 : !> \param ewald_env ...
1436 : !> \param ewald_pw ...
1437 : !> \param cell ...
1438 : !> \param particle_set ...
1439 : !> \param charges ...
1440 : !> \param potential ...
1441 : ! **************************************************************************************************
1442 3912 : SUBROUTINE apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, potential)
1443 : TYPE(ewald_environment_type), POINTER :: ewald_env
1444 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1445 : TYPE(cell_type), POINTER :: cell
1446 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1447 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1448 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential
1449 :
1450 : TYPE(mp_para_env_type), POINTER :: para_env
1451 :
1452 3912 : CALL ewald_env_get(ewald_env, para_env=para_env)
1453 27676 : potential = 0.0_dp
1454 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
1455 3912 : particle_set, potential)
1456 51440 : CALL para_env%sum(potential)
1457 :
1458 3912 : END SUBROUTINE apply_potential
1459 :
1460 : ! **************************************************************************************************
1461 : !> \brief ...
1462 : !> \param eeqn ...
1463 : !> \param fm_mat ...
1464 : !> \param mmat ...
1465 : !> \param ewald_env ...
1466 : !> \param ewald_pw ...
1467 : !> \param cell ...
1468 : !> \param particle_set ...
1469 : !> \param charges ...
1470 : !> \param rhs ...
1471 : !> \param potential ...
1472 : ! **************************************************************************************************
1473 8572 : SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
1474 8572 : cell, particle_set, charges, rhs, potential)
1475 : REAL(KIND=dp), INTENT(INOUT) :: eeqn
1476 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat, mmat
1477 : TYPE(ewald_environment_type), POINTER :: ewald_env
1478 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1479 : TYPE(cell_type), POINTER :: cell
1480 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1481 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1482 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rhs
1483 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential
1484 :
1485 : INTEGER :: na, ns
1486 : REAL(KIND=dp) :: lambda
1487 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mvec
1488 : TYPE(mp_para_env_type), POINTER :: para_env
1489 :
1490 8572 : ns = SIZE(charges)
1491 8572 : na = ns - 1
1492 8572 : CALL ewald_env_get(ewald_env, para_env=para_env)
1493 75160 : potential = 0.0_dp
1494 : CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
1495 8572 : particle_set, potential(1:na))
1496 124604 : CALL para_env%sum(potential(1:na))
1497 8572 : CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
1498 133176 : eeqn = 0.5_dp*SUM(charges(1:na)*potential(1:na)) + SUM(charges(1:na)*rhs(1:na))
1499 75160 : potential(1:ns) = potential(1:ns) + rhs(1:ns)
1500 25716 : ALLOCATE (mvec(ns))
1501 8572 : CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
1502 66588 : lambda = -SUM(mvec(1:na))/REAL(na, KIND=dp)
1503 66588 : potential(1:na) = mvec(1:na) + lambda
1504 8572 : DEALLOCATE (mvec)
1505 :
1506 8572 : END SUBROUTINE get_energy_gradient
1507 :
1508 : ! **************************************************************************************************
1509 : !> \brief ...
1510 : !> \param qs_env ...
1511 : !> \param charges ...
1512 : !> \param ef_energy ...
1513 : ! **************************************************************************************************
1514 328 : SUBROUTINE eeq_efield_energy(qs_env, charges, ef_energy)
1515 : TYPE(qs_environment_type), POINTER :: qs_env
1516 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges
1517 : REAL(KIND=dp), INTENT(OUT) :: ef_energy
1518 :
1519 : COMPLEX(KIND=dp) :: zdeta
1520 : COMPLEX(KIND=dp), DIMENSION(3) :: zi
1521 : INTEGER :: ia, idir, natom
1522 : LOGICAL :: dfield
1523 : REAL(KIND=dp) :: kr, omega, q
1524 : REAL(KIND=dp), DIMENSION(3) :: ci, dfilter, fieldpol, fpolvec, kvec, &
1525 : qi, ria
1526 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1527 : TYPE(cell_type), POINTER :: cell
1528 : TYPE(dft_control_type), POINTER :: dft_control
1529 328 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1530 :
1531 328 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1532 328 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1533 :
1534 328 : IF (dft_control%apply_period_efield) THEN
1535 164 : dfield = dft_control%period_efield%displacement_field
1536 :
1537 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1538 0 : CPABORT("use of strength_list not implemented for eeq_efield_energy")
1539 : END IF
1540 :
1541 656 : fieldpol = dft_control%period_efield%polarisation
1542 1148 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1543 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
1544 2132 : hmat = cell%hmat(:, :)/twopi
1545 656 : DO idir = 1, 3
1546 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1547 656 : + fieldpol(3)*hmat(3, idir)
1548 : END DO
1549 :
1550 656 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1551 820 : DO ia = 1, natom
1552 656 : q = charges(ia)
1553 2624 : ria = particle_set(ia)%r
1554 2624 : ria = pbc(ria, cell)
1555 2788 : DO idir = 1, 3
1556 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
1557 7872 : kr = SUM(kvec(:)*ria(:))
1558 1968 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
1559 2624 : zi(idir) = zi(idir)*zdeta
1560 : END DO
1561 : END DO
1562 656 : qi = AIMAG(LOG(zi))
1563 164 : IF (dfield) THEN
1564 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1565 0 : omega = cell%deth
1566 0 : ci = MATMUL(hmat, qi)/omega
1567 0 : ef_energy = 0.0_dp
1568 0 : DO idir = 1, 3
1569 0 : ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
1570 : END DO
1571 0 : ef_energy = -0.25_dp*omega/twopi*ef_energy
1572 : ELSE
1573 656 : ef_energy = SUM(fpolvec(:)*qi(:))
1574 : END IF
1575 :
1576 164 : ELSE IF (dft_control%apply_efield) THEN
1577 :
1578 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1579 656 : dft_control%efield_fields(1)%efield%strength
1580 :
1581 164 : ef_energy = 0.0_dp
1582 820 : DO ia = 1, natom
1583 2624 : ria = particle_set(ia)%r
1584 2624 : ria = pbc(ria, cell)
1585 656 : q = charges(ia)
1586 2788 : ef_energy = ef_energy - q*SUM(fieldpol*ria)
1587 : END DO
1588 :
1589 : ELSE
1590 0 : CPABORT("apply field")
1591 : END IF
1592 :
1593 328 : END SUBROUTINE eeq_efield_energy
1594 :
1595 : ! **************************************************************************************************
1596 : !> \brief ...
1597 : !> \param qs_env ...
1598 : !> \param efr ...
1599 : ! **************************************************************************************************
1600 328 : SUBROUTINE eeq_efield_pot(qs_env, efr)
1601 : TYPE(qs_environment_type), POINTER :: qs_env
1602 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: efr
1603 :
1604 : INTEGER :: ia, idir, natom
1605 : LOGICAL :: dfield
1606 : REAL(KIND=dp) :: kr
1607 : REAL(KIND=dp), DIMENSION(3) :: fieldpol, fpolvec, kvec, ria
1608 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1609 : TYPE(cell_type), POINTER :: cell
1610 : TYPE(dft_control_type), POINTER :: dft_control
1611 328 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1612 :
1613 328 : CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1614 328 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1615 :
1616 328 : IF (dft_control%apply_period_efield) THEN
1617 164 : dfield = dft_control%period_efield%displacement_field
1618 :
1619 164 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1620 0 : CPABORT("use of strength_list not implemented for eeq_efield_pot")
1621 : END IF
1622 :
1623 656 : fieldpol = dft_control%period_efield%polarisation
1624 1148 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1625 656 : fieldpol = -fieldpol*dft_control%period_efield%strength
1626 2132 : hmat = cell%hmat(:, :)/twopi
1627 656 : DO idir = 1, 3
1628 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1629 656 : + fieldpol(3)*hmat(3, idir)
1630 : END DO
1631 :
1632 164 : IF (dfield) THEN
1633 : ! dE/dq depends on q, postpone calculation
1634 0 : efr = 0.0_dp
1635 : ELSE
1636 820 : efr = 0.0_dp
1637 820 : DO ia = 1, natom
1638 2624 : ria = particle_set(ia)%r
1639 2624 : ria = pbc(ria, cell)
1640 2788 : DO idir = 1, 3
1641 7872 : kvec(:) = twopi*cell%h_inv(idir, :)
1642 7872 : kr = SUM(kvec(:)*ria(:))
1643 2624 : efr(ia) = efr(ia) + kr*fpolvec(idir)
1644 : END DO
1645 : END DO
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 820 : DO ia = 1, natom
1654 2624 : ria = particle_set(ia)%r
1655 2624 : ria = pbc(ria, cell)
1656 2788 : efr(ia) = -SUM(fieldpol*ria)
1657 : END DO
1658 :
1659 : ELSE
1660 0 : CPABORT("apply field")
1661 : END IF
1662 :
1663 328 : END SUBROUTINE eeq_efield_pot
1664 :
1665 : ! **************************************************************************************************
1666 : !> \brief ...
1667 : !> \param charges ...
1668 : !> \param dft_control ...
1669 : !> \param particle_set ...
1670 : !> \param cell ...
1671 : !> \param efr ...
1672 : ! **************************************************************************************************
1673 0 : SUBROUTINE eeq_dfield_pot(charges, dft_control, particle_set, cell, efr)
1674 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges
1675 : TYPE(dft_control_type), POINTER :: dft_control
1676 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1677 : TYPE(cell_type), POINTER :: cell
1678 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: efr
1679 :
1680 : COMPLEX(KIND=dp) :: zdeta
1681 : COMPLEX(KIND=dp), DIMENSION(3) :: zi
1682 : INTEGER :: ia, idir, natom
1683 : REAL(KIND=dp) :: kr, omega, q
1684 : REAL(KIND=dp), DIMENSION(3) :: ci, dfilter, fieldpol, kvec, qi, ria
1685 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1686 :
1687 0 : natom = SIZE(particle_set)
1688 :
1689 0 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1690 0 : CPABORT("use of strength_list not implemented for eeq_dfield_pot")
1691 : END IF
1692 :
1693 0 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1694 0 : fieldpol = dft_control%period_efield%polarisation
1695 0 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1696 0 : fieldpol = -fieldpol*dft_control%period_efield%strength
1697 0 : hmat = cell%hmat(:, :)/twopi
1698 0 : omega = cell%deth
1699 : !
1700 0 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1701 0 : DO ia = 1, natom
1702 0 : q = charges(ia)
1703 0 : ria = particle_set(ia)%r
1704 0 : ria = pbc(ria, cell)
1705 0 : DO idir = 1, 3
1706 0 : kvec(:) = twopi*cell%h_inv(idir, :)
1707 0 : kr = SUM(kvec(:)*ria(:))
1708 0 : zdeta = CMPLX(COS(kr), SIN(kr), KIND=dp)**q
1709 0 : zi(idir) = zi(idir)*zdeta
1710 : END DO
1711 : END DO
1712 0 : qi = AIMAG(LOG(zi))
1713 0 : ci = MATMUL(hmat, qi)/omega
1714 0 : ci = dfilter*(fieldpol - 2._dp*twopi*ci)
1715 0 : DO ia = 1, natom
1716 0 : ria = particle_set(ia)%r
1717 0 : ria = pbc(ria, cell)
1718 0 : efr(ia) = efr(ia) - SUM(ci*ria)
1719 : END DO
1720 :
1721 0 : END SUBROUTINE eeq_dfield_pot
1722 :
1723 : ! **************************************************************************************************
1724 : !> \brief ...
1725 : !> \param qs_env ...
1726 : !> \param charges ...
1727 : !> \param qlag ...
1728 : ! **************************************************************************************************
1729 8 : SUBROUTINE eeq_efield_force_loc(qs_env, charges, qlag)
1730 : TYPE(qs_environment_type), POINTER :: qs_env
1731 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1732 :
1733 : INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1734 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1735 : REAL(KIND=dp) :: q
1736 : REAL(KIND=dp), DIMENSION(3) :: fieldpol
1737 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1738 : TYPE(cell_type), POINTER :: cell
1739 : TYPE(dft_control_type), POINTER :: dft_control
1740 : TYPE(distribution_1d_type), POINTER :: local_particles
1741 : TYPE(mp_para_env_type), POINTER :: para_env
1742 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1743 8 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1744 :
1745 : CALL get_qs_env(qs_env=qs_env, &
1746 : dft_control=dft_control, &
1747 : cell=cell, particle_set=particle_set, &
1748 : nkind=nkind, natom=natom, &
1749 : para_env=para_env, &
1750 8 : local_particles=local_particles)
1751 :
1752 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1753 32 : dft_control%efield_fields(1)%efield%strength
1754 :
1755 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1756 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1757 8 : CALL get_qs_env(qs_env=qs_env, force=force)
1758 :
1759 32 : DO ikind = 1, nkind
1760 152 : force(ikind)%efield = 0.0_dp
1761 40 : DO ia = 1, local_particles%n_el(ikind)
1762 16 : iatom = local_particles%list(ikind)%array(ia)
1763 16 : q = charges(iatom) - qlag(iatom)
1764 16 : atom_a = atom_of_kind(iatom)
1765 88 : force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
1766 : END DO
1767 288 : CALL para_env%sum(force(ikind)%efield)
1768 : END DO
1769 :
1770 16 : END SUBROUTINE eeq_efield_force_loc
1771 :
1772 : ! **************************************************************************************************
1773 : !> \brief ...
1774 : !> \param qs_env ...
1775 : !> \param charges ...
1776 : !> \param qlag ...
1777 : ! **************************************************************************************************
1778 8 : SUBROUTINE eeq_efield_force_periodic(qs_env, charges, qlag)
1779 : TYPE(qs_environment_type), POINTER :: qs_env
1780 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1781 :
1782 : INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1783 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1784 : LOGICAL :: dfield, use_virial
1785 : REAL(KIND=dp) :: q
1786 : REAL(KIND=dp), DIMENSION(3) :: fa, fieldpol, ria
1787 : REAL(KIND=dp), DIMENSION(3, 3) :: pve
1788 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1789 : TYPE(cell_type), POINTER :: cell
1790 : TYPE(dft_control_type), POINTER :: dft_control
1791 : TYPE(distribution_1d_type), POINTER :: local_particles
1792 : TYPE(mp_para_env_type), POINTER :: para_env
1793 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1794 8 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1795 : TYPE(virial_type), POINTER :: virial
1796 :
1797 : CALL get_qs_env(qs_env=qs_env, &
1798 : dft_control=dft_control, &
1799 : cell=cell, particle_set=particle_set, &
1800 : virial=virial, &
1801 : nkind=nkind, natom=natom, &
1802 : para_env=para_env, &
1803 8 : local_particles=local_particles)
1804 :
1805 8 : dfield = dft_control%period_efield%displacement_field
1806 8 : CPASSERT(.NOT. dfield)
1807 :
1808 8 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1809 0 : CPABORT("use of strength_list not implemented for eeq_efield_force_periodic")
1810 : END IF
1811 :
1812 32 : fieldpol = dft_control%period_efield%polarisation
1813 56 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
1814 32 : fieldpol = -fieldpol*dft_control%period_efield%strength
1815 :
1816 8 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1817 :
1818 8 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1819 8 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1820 8 : CALL get_qs_env(qs_env=qs_env, force=force)
1821 :
1822 8 : pve = 0.0_dp
1823 32 : DO ikind = 1, nkind
1824 152 : force(ikind)%efield = 0.0_dp
1825 40 : DO ia = 1, local_particles%n_el(ikind)
1826 16 : iatom = local_particles%list(ikind)%array(ia)
1827 16 : q = charges(iatom) - qlag(iatom)
1828 64 : fa(1:3) = q*fieldpol(1:3)
1829 16 : atom_a = atom_of_kind(iatom)
1830 64 : force(ikind)%efield(1:3, atom_a) = fa
1831 40 : IF (use_virial) THEN
1832 0 : ria = particle_set(ia)%r
1833 0 : ria = pbc(ria, cell)
1834 0 : CALL virial_pair_force(pve, -0.5_dp, fa, ria)
1835 0 : CALL virial_pair_force(pve, -0.5_dp, ria, fa)
1836 : END IF
1837 : END DO
1838 288 : CALL para_env%sum(force(ikind)%efield)
1839 : END DO
1840 104 : virial%pv_virial = virial%pv_virial + pve
1841 :
1842 16 : END SUBROUTINE eeq_efield_force_periodic
1843 :
1844 12012 : END MODULE eeq_method
|
