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