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