(git:01d6ff2)
Loading...
Searching...
No Matches
eeq_method.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
117CONTAINS
118
119! **************************************************************************************************
120!> \brief ...
121!> \param qs_env ...
122!> \param iounit ...
123!> \param print_level ...
124!> \param ext ...
125! **************************************************************************************************
126 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 REAL(kind=dp), DIMENSION(:), POINTER :: charges
135 TYPE(cell_type), POINTER :: cell
136 TYPE(eeq_solver_type) :: eeq_sparam
137 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
138
139 mark_used(print_level)
140
141 CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, cell=cell)
142 IF (ext) THEN
143 NULLIFY (charges)
144 CALL get_qs_env(qs_env, eeq=charges)
145 cpassert(ASSOCIATED(charges))
146 enshift_type = 0
147 ELSE
148 ALLOCATE (charges(natom))
149 ! enforce en shift method 1 (original/molecular)
150 ! method 2 from paper on PBC seems not to work
151 enshift_type = 1
152 !IF (ALL(cell%perd == 0)) enshift_type = 1
153 CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
154 END IF
155
156 IF (iounit > 0) THEN
157
158 IF (enshift_type == 0) THEN
159 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (External)"
160 ELSE IF (enshift_type == 1) THEN
161 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Molecules))"
162 ELSE IF (enshift_type == 2) THEN
163 WRITE (unit=iounit, fmt="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Crystals))"
164 ELSE
165 cpabort("Unknown enshift_type")
166 END IF
167 WRITE (unit=iounit, fmt="(/,T2,A)") &
168 "# Atom Element Kind Atomic Charge"
169
170 DO iatom = 1, natom
171 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
172 element_symbol=element_symbol, &
173 kind_number=ikind)
174 WRITE (unit=iounit, fmt="(T4,I8,T18,A2,I10,T43,F12.4)") &
175 iatom, element_symbol, ikind, charges(iatom)
176 END DO
177
178 END IF
179
180 IF (.NOT. ext) DEALLOCATE (charges)
181
182 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chia, cnumbers, efr, gam
209 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
210 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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
220 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
221 TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
222 print_section
223
224 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 dft_control=dft_control)
234 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
235
236 logger => cp_get_default_logger()
237 IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
238 iunit = cp_logger_get_default_unit_nr()
239 ELSE
240 iunit = -1
241 END IF
242
243 totalcharge = dft_control%charge
244
245 CALL get_cnumbers(qs_env, cnumbers, dcnum, .false.)
246
247 ! gamma[a,b]
248 ALLOCATE (gab(nkind, nkind), gam(nkind))
249 gab = 0.0_dp
250 gam = 0.0_dp
251 DO ikind = 1, nkind
252 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
253 CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
254 DO jkind = 1, nkind
255 CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
256 CALL get_eeq_data(zb, eeq_model, rad=alb)
257 !
258 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
259 !
260 END DO
261 END DO
262
263 ! Chi[a,a]
264 sgamma = 8.0_dp ! see D4 for periodic systems paper
265 esg = 1.0_dp + exp(sgamma)
266 ALLOCATE (chia(natom))
267 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
268 DO iatom = 1, natom
269 ikind = kind_of(iatom)
270 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
271 CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
272 !
273 IF (enshift_type == 1) THEN
274 scn = cnumbers(iatom)/sqrt(cnumbers(iatom) + 1.0e-14_dp)
275 ELSE IF (enshift_type == 2) THEN
276 scn = log(esg/(esg - cnumbers(iatom)))
277 ELSE
278 cpabort("Unknown enshift_type")
279 END IF
280 chia(iatom) = xi - kappa*scn
281 !
282 END DO
283
284 ! efield
285 IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
286 dft_control%apply_efield_field) THEN
287 ALLOCATE (efr(natom))
288 efr(1:natom) = 0.0_dp
289 CALL eeq_efield_pot(qs_env, efr)
290 chia(1:natom) = chia(1:natom) + efr(1:natom)
291 DEALLOCATE (efr)
292 END IF
293
294 CALL cnumber_release(cnumbers, dcnum, .false.)
295
296 CALL get_cell(cell, periodic=periodic)
297 do_ewald = .NOT. all(periodic == 0)
298 IF (do_ewald) THEN
299 ALLOCATE (ewald_env)
300 CALL ewald_env_create(ewald_env, para_env)
301 poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
302 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
303 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
304 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
305 CALL get_qs_env(qs_env, cell_ref=cell_ref)
306 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
307 silent=.true., pset="EEQ")
308 ALLOCATE (ewald_pw)
309 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 ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
316 !
317 CALL ewald_env_release(ewald_env)
318 CALL ewald_pw_release(ewald_pw)
319 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 totalcharge=totalcharge, iounit=iunit)
325 END IF
326
327 DEALLOCATE (gab, gam, chia)
328
329 CALL timestop(handle)
330
331 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
362 INTEGER, DIMENSION(3) :: periodic
363 LOGICAL :: do_ewald, use_virial
364 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: c_radius, cnumbers, gam, qlag
370 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 REAL(kind=dp), DIMENSION(:), POINTER :: chrgx, dchia
374 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 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 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
386 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
387 TYPE(mp_para_env_type), POINTER :: para_env
388 TYPE(neighbor_list_iterator_p_type), &
389 DIMENSION(:), POINTER :: nl_iterator
390 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
391 POINTER :: sab_ew
392 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
393 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
394 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 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 dft_control=dft_control)
412
413 logger => cp_get_default_logger()
414 IF (para_env%is_source() .AND. logger%iter_info%print_level >= medium_print_level) THEN
415 iunit = cp_logger_get_default_unit_nr()
416 ELSE
417 iunit = -1
418 END IF
419
420 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
421 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
422
423 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
424
425 totalcharge = dft_control%charge
426
427 CALL get_cnumbers(qs_env, cnumbers, dcnum, .true.)
428
429 ! gamma[a,b]
430 ALLOCATE (gab(nkind, nkind), gam(nkind))
431 gab = 0.0_dp
432 gam = 0.0_dp
433 DO ikind = 1, nkind
434 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
435 CALL get_eeq_data(za, eeq_model, eta=gam(ikind), rad=ala)
436 DO jkind = 1, nkind
437 CALL get_qs_kind(qs_kind_set(jkind), zatom=zb)
438 CALL get_eeq_data(zb, eeq_model, rad=alb)
439 !
440 gab(ikind, jkind) = sqrt(1._dp/(ala*ala + alb*alb))
441 !
442 END DO
443 END DO
444
445 ALLOCATE (qlag(natom))
446
447 CALL get_cell(cell, periodic=periodic)
448 do_ewald = .NOT. all(periodic == 0)
449 IF (do_ewald) THEN
450 ALLOCATE (ewald_env)
451 CALL ewald_env_create(ewald_env, para_env)
452 poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
453 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
454 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
455 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
456 CALL get_qs_env(qs_env, cell_ref=cell_ref)
457 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
458 silent=.true., pset="EEQ")
459 ALLOCATE (ewald_pw)
460 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 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 para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
470 END IF
471
472 sgamma = 8.0_dp ! see D4 for periodic systems paper
473 esg = 1.0_dp + exp(sgamma)
474 ALLOCATE (chrgx(natom), dchia(natom))
475 DO iatom = 1, natom
476 ikind = kind_of(iatom)
477 CALL get_qs_kind(qs_kind_set(ikind), zatom=za)
478 CALL get_eeq_data(za, eeq_model, chi=xi, kcn=kappa)
479 !
480 IF (response_only) THEN
481 ctot = -0.5_dp*qlag(iatom)
482 ELSE
483 ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
484 END IF
485 IF (enshift_type == 1) THEN
486 scn = sqrt(cnumbers(iatom)) + 1.0e-14_dp
487 dchia(iatom) = -ctot*kappa/scn
488 ELSE IF (enshift_type == 2) THEN
489 cn = cnumbers(iatom)
490 scn = 1.0_dp/(esg - cn)
491 dchia(iatom) = -ctot*kappa*scn
492 ELSE
493 cpabort("Unknown enshift_type")
494 END IF
495 END DO
496
497 ! Efield
498 IF (dft_control%apply_period_efield) THEN
499 CALL eeq_efield_force_periodic(qs_env, charges, qlag)
500 ELSE IF (dft_control%apply_efield) THEN
501 CALL eeq_efield_force_loc(qs_env, charges, qlag)
502 ELSE IF (dft_control%apply_efield_field) THEN
503 cpabort("apply field")
504 END IF
505
506 ! Forces from q*X
507 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
508 DO ikind = 1, nkind
509 DO ia = 1, local_particles%n_el(ikind)
510 iatom = local_particles%list(ikind)%array(ia)
511 DO i = 1, dcnum(iatom)%neighbors
512 katom = dcnum(iatom)%nlist(i)
513 rik = dcnum(iatom)%rik(:, i)
514 drk = sqrt(sum(rik(:)**2))
515 IF (drk > 1.e-3_dp) THEN
516 fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
517 gradient(:, iatom) = gradient(:, iatom) - fdik(:)
518 gradient(:, katom) = gradient(:, katom) + fdik(:)
519 IF (use_virial) THEN
520 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 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 molecule_set=molecule_set)
535 subcells = 2.0_dp
536 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
537 rcut = 2.0_dp*rcut
538 NULLIFY (sab_ew)
539 ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
540 c_radius(:) = rcut
541 default_present = .true.
542 ALLOCATE (atom2d(nkind))
543 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
544 molecule_set, .false., particle_set=particle_set)
545 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 subcells=subcells, operator_type="PP", nlname="sab_ew")
548 DEALLOCATE (c_radius, pair_radius, default_present)
549 CALL atom2d_cleanup(atom2d)
550 !
551 CALL neighbor_list_iterator_create(nl_iterator, sab_ew)
552 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
553 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
554 iatom=iatom, jatom=jatom, r=rij)
555 !
556 dr2 = sum(rij**2)
557 dr = sqrt(dr2)
558 IF (dr > rcut .OR. dr < 1.e-6_dp) cycle
559 fe = 1.0_dp
560 IF (iatom == jatom) fe = 0.5_dp
561 IF (response_only) THEN
562 qq = -qlag(iatom)*charges(jatom)
563 ELSE
564 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
565 END IF
566 gama = gab(ikind, jkind)
567 gam2 = gama*gama
568 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
569 - 2._dp*alpha*exp(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
570 IF (response_only) THEN
571 qq1 = -qlag(iatom)*charges(jatom)
572 qq2 = -qlag(jatom)*charges(iatom)
573 ELSE
574 qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
575 qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
576 END IF
577 fdik(:) = -qq1*grc*rij(:)/dr
578 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
579 gradient(:, jatom) = gradient(:, jatom) - fdik(:)
580 IF (use_virial) THEN
581 CALL virial_pair_force(stress, -fe, fdik, rij)
582 END IF
583 fdik(:) = qq2*grc*rij(:)/dr
584 gradient(:, iatom) = gradient(:, iatom) - fdik(:)
585 gradient(:, jatom) = gradient(:, jatom) + fdik(:)
586 IF (use_virial) THEN
587 CALL virial_pair_force(stress, fe, fdik, rij)
588 END IF
589 END DO
590 CALL neighbor_list_iterator_release(nl_iterator)
591 !
592 CALL release_neighbor_list_sets(sab_ew)
593 ELSE
594 DO ikind = 1, nkind
595 DO ia = 1, local_particles%n_el(ikind)
596 iatom = local_particles%list(ikind)%array(ia)
597 ri(1:3) = particle_set(iatom)%r(1:3)
598 DO jatom = 1, natom
599 IF (iatom == jatom) cycle
600 jkind = kind_of(jatom)
601 IF (response_only) THEN
602 qq = -qlag(iatom)*charges(jatom)
603 ELSE
604 qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
605 END IF
606 rj(1:3) = particle_set(jatom)%r(1:3)
607 rij(1:3) = ri(1:3) - rj(1:3)
608 rij = pbc(rij, cell)
609 dr2 = sum(rij**2)
610 dr = sqrt(dr2)
611 gama = gab(ikind, jkind)
612 gam2 = gama*gama
613 grc = 2._dp*gama*exp(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
614 fdik(:) = qq*grc*rij(:)/dr
615 gradient(:, iatom) = gradient(:, iatom) + fdik(:)
616 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 IF (do_ewald) THEN
624 ALLOCATE (epforce(3, natom))
625 epforce = 0.0_dp
626 IF (response_only) THEN
627 dchia(1:natom) = qlag(1:natom)
628 ELSE
629 dchia(1:natom) = -charges(1:natom) + qlag(1:natom)
630 END IF
631 chrgx(1:natom) = charges(1:natom)
632 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
633 particle_set, dchia, epforce)
634 dchia(1:natom) = charges(1:natom)
635 chrgx(1:natom) = qlag(1:natom)
636 CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
637 particle_set, dchia, epforce)
638 gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + epforce(1:3, 1:natom)
639 DEALLOCATE (epforce)
640
641 ! virial
642 IF (use_virial) THEN
643 chrgx(1:natom) = charges(1:natom) - qlag(1:natom)
644 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
645 stress = stress - pvir
646 chrgx(1:natom) = qlag(1:natom)
647 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
648 stress = stress + pvir
649 IF (response_only) THEN
650 chrgx(1:natom) = charges(1:natom)
651 CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
652 stress = stress + pvir
653 END IF
654 END IF
655 !
656 CALL ewald_env_release(ewald_env)
657 CALL ewald_pw_release(ewald_pw)
658 DEALLOCATE (ewald_env, ewald_pw)
659 END IF
660
661 CALL cnumber_release(cnumbers, dcnum, .true.)
662
663 DEALLOCATE (gab, gam, qlag, chrgx, dchia)
664
665 CALL timestop(handle)
666
667 END SUBROUTINE eeq_forces
668
669! **************************************************************************************************
670!> \brief ...
671!> \param qs_env ...
672!> \param cnumbers ...
673!> \param dcnum ...
674!> \param calculate_forces ...
675! **************************************************************************************************
676 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 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 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
693 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
694 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
695 POINTER :: sab_cn
696 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
697 TYPE(qs_dispersion_type), POINTER :: disp
698 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 molecule_set=molecule_set)
708 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
709
710 ! Check for dispersion_env and sab_cn needed for cnumbers
711 ALLOCATE (disp)
712 disp%k1 = 16.0_dp
713 disp%k2 = 4._dp/3._dp
714 disp%eps_cn = 1.e-6_dp
715 disp%max_elem = maxelem
716 ALLOCATE (disp%rcov(maxelem))
717 disp%rcov(1:maxelem) = bohr*disp%k2*rcov(1:maxelem)
718 subcells = 2.0_dp
719 ! Build the neighbor lists for the CN
720 NULLIFY (sab_cn)
721 ALLOCATE (c_radius(nkind), default_present(nkind), pair_radius(nkind, nkind))
722 c_radius(:) = 0.0_dp
723 default_present = .true.
724 DO ikind = 1, nkind
725 CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
726 c_radius(ikind) = 4._dp*rcov(za)*bohr
727 END DO
728 ALLOCATE (atom2d(nkind))
729 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
730 molecule_set, .false., particle_set=particle_set)
731 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 subcells=subcells, operator_type="PP", nlname="sab_cn")
734 disp%sab_cn => sab_cn
735 DEALLOCATE (c_radius, pair_radius, default_present)
736 CALL atom2d_cleanup(atom2d)
737
738 ! Calculate coordination numbers
739 CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces, disp_env=disp)
740
741 CALL qs_dispersion_release(disp)
742
743 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 SUBROUTINE eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, &
767 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 CALL timeset(routinen, handle)
796
797 do_ewald = .false.
798 IF (PRESENT(ewald)) do_ewald = ewald
799 !
800 qtot = 0.0_dp
801 IF (PRESENT(totalcharge)) qtot = totalcharge
802 !
803 iunit = -1
804 IF (PRESENT(iounit)) iunit = iounit
805
806 ! EEQ solver parameters
807 do_direct = eeq_sparam%direct
808 do_sparse = eeq_sparam%sparse
809
810 do_displ = .false.
811 IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
812 do_displ = dft_control%period_efield%displacement_field
813 END IF
814
815 ! sparse option NYA
816 IF (do_sparse) THEN
817 CALL cp_abort(__location__, "EEQ: Sparse option not yet available")
818 END IF
819 !
820 natom = SIZE(particle_set)
821 nkind = SIZE(gam)
822 ns = natom + 1
823 CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
824 nrow_global=ns, ncol_global=ns)
825 CALL cp_fm_create(eeq_mat, mat_struct)
826 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
827 !
828 IF (do_ewald) THEN
829 cpassert(PRESENT(ewald_env))
830 cpassert(PRESENT(ewald_pw))
831 IF (do_direct) THEN
832 IF (do_displ) THEN
833 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 ewald_env, ewald_pw, iounit)
838 END IF
839 ELSE
840 IF (do_displ) THEN
841 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 ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
847 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 ewald_env, ewald_pw, iounit)
852 END IF
853 END IF
854 END IF
855 IF (qtot /= 0._dp) THEN
856 CALL get_cell(cell=cell, deth=deth)
857 CALL ewald_env_get(ewald_env, alpha=alpha)
858 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 cell, chia, gam, gab, qtot, ftime)
863 IF (iounit > 0) THEN
864 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
865 END IF
866 END IF
867 CALL cp_fm_struct_release(mat_struct)
868 CALL cp_fm_release(eeq_mat)
869
870 CALL timestop(handle)
871
872 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 SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
890 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 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 CALL timeset(routinen, handle)
916 ti = m_walltime()
917
918 natom = SIZE(particle_set)
919 nkind = SIZE(gam)
920 !
921 ns = natom + 1
922 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 row_indices=rind, col_indices=cind)
925 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
926 nrow_global=ns, ncol_global=1)
927 CALL cp_fm_create(rhs_vec, vec_struct)
928 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
929 row_indices=rvind, col_indices=cvind)
930 !
931 ! set up matrix
932 CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
933 CALL cp_fm_set_all(rhs_vec, 0.0_dp)
934 DO ir = 1, nrloc
935 iar = rind(ir)
936 IF (iar > natom) cycle
937 ikind = kind_of(iar)
938 ri(1:3) = particle_set(iar)%r(1:3)
939 DO ic = 1, ncloc
940 iac = cind(ic)
941 IF (iac > natom) cycle
942 jkind = kind_of(iac)
943 rj(1:3) = particle_set(iac)%r(1:3)
944 IF (iar == iac) THEN
945 grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
946 ELSE
947 rij(1:3) = ri(1:3) - rj(1:3)
948 rij = pbc(rij, cell)
949 dr = sqrt(sum(rij**2))
950 grc = erf(gab(ikind, jkind)*dr)/dr
951 END IF
952 eeq_mat%local_data(ir, ic) = grc
953 END DO
954 END DO
955 ! set up rhs vector
956 DO ir = 1, nrvloc
957 iar = rvind(ir)
958 DO ic = 1, ncvloc
959 iac = cvind(ic)
960 ia = max(iar, iac)
961 IF (ia > natom) THEN
962 xr = qtot
963 ELSE
964 xr = -chia(ia)
965 END IF
966 rhs_vec%local_data(ir, ic) = xr
967 END DO
968 END DO
969 !
970 CALL cp_fm_solve(eeq_mat, rhs_vec)
971 !
972 charges = 0.0_dp
973 lambda = 0.0_dp
974 DO ir = 1, nrvloc
975 iar = rvind(ir)
976 DO ic = 1, ncvloc
977 iac = cvind(ic)
978 ia = max(iar, iac)
979 IF (ia <= natom) THEN
980 xr = rhs_vec%local_data(ir, ic)
981 charges(ia) = xr
982 ELSE
983 lambda = rhs_vec%local_data(ir, ic)
984 END IF
985 END DO
986 END DO
987 CALL para_env%sum(lambda)
988 CALL para_env%sum(charges)
989 !
990 ! energy: 0.5*(q^T.X - lambda*totalcharge)
991 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
992
993 CALL cp_fm_struct_release(vec_struct)
994 CALL cp_fm_release(rhs_vec)
995
996 te = m_walltime()
997 ftime = te - ti
998 CALL timestop(handle)
999
1000 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 SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1022 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1050 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 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 CALL timeset(routinen, handle)
1059 ti = m_walltime()
1060
1061 ierror = 0
1062
1063 iunit = -1
1064 IF (PRESENT(iounit)) iunit = iounit
1065
1066 natom = SIZE(particle_set)
1067 nkind = SIZE(gam)
1068 !
1069 CALL get_cell(cell=cell, deth=deth)
1070 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1071 ad = 2.0_dp*alpha*oorootpi
1072 IF (ewald_type /= do_ewald_spme) THEN
1073 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1074 END IF
1075 !
1076 rmax = 2.0_dp*rcut
1077 ! max cells used
1078 CALL get_cell(cell, h=hmat, periodic=periodic)
1079 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1080 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1081 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1082 IF (periodic(1) == 0) ncell(1) = 0
1083 IF (periodic(2) == 0) ncell(2) = 0
1084 IF (periodic(3) == 0) ncell(3) = 0
1085 !
1086 CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1087 chia, gam, gab, qtot, ftime)
1088 IF (iunit > 0) THEN
1089 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1090 END IF
1091 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1092 CALL cp_fm_create(pmat, mat_struct)
1093 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 row_indices=rind, col_indices=cind)
1098 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1099 DO ir = 1, nrloc
1100 iar = rind(ir)
1101 ri = 0.0_dp
1102 IF (iar <= natom) THEN
1103 ikind = kind_of(iar)
1104 ri(1:3) = particle_set(iar)%r(1:3)
1105 END IF
1106 DO ic = 1, ncloc
1107 iac = cind(ic)
1108 IF (iac > natom .AND. iar > natom) THEN
1109 eeq_mat%local_data(ir, ic) = 0.0_dp
1110 cycle
1111 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1112 eeq_mat%local_data(ir, ic) = 1.0_dp
1113 cycle
1114 END IF
1115 jkind = kind_of(iac)
1116 rj(1:3) = particle_set(iac)%r(1:3)
1117 rij(1:3) = ri(1:3) - rj(1:3)
1118 rij = pbc(rij, cell)
1119 DO ix = -ncell(1), ncell(1)
1120 DO iy = -ncell(2), ncell(2)
1121 DO iz = -ncell(3), ncell(3)
1122 cvec = [ix, iy, iz]
1123 rijl = rij + matmul(hmat, cvec)
1124 dr = sqrt(sum(rijl**2))
1125 IF (dr > rmax) cycle
1126 IF (iar == iac .AND. dr < 0.00001_dp) THEN
1127 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1128 ELSE
1129 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1130 END IF
1131 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 row_indices=rind, col_indices=cind)
1141 CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1142 DO ir = 1, nrloc
1143 iar = rind(ir)
1144 ri = 0.0_dp
1145 IF (iar <= natom) THEN
1146 ikind = kind_of(iar)
1147 ri(1:3) = particle_set(iar)%r(1:3)
1148 END IF
1149 DO ic = 1, ncloc
1150 iac = cind(ic)
1151 IF (iac > natom .AND. iar > natom) THEN
1152 pmat%local_data(ir, ic) = 0.0_dp
1153 cycle
1154 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1155 pmat%local_data(ir, ic) = 1.0_dp
1156 cycle
1157 END IF
1158 jkind = kind_of(iac)
1159 rj(1:3) = particle_set(iac)%r(1:3)
1160 rij(1:3) = ri(1:3) - rj(1:3)
1161 rij = pbc(rij, cell)
1162 IF (iar == iac) THEN
1163 grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1164 ELSE
1165 grc2 = erf(gab(ikind, jkind)*dr)/dr
1166 END IF
1167 pmat%local_data(ir, ic) = grc2
1168 END DO
1169 END DO
1170 CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1171 ! preconditioner invers
1172 CALL cp_fm_invert(pmat, mmat)
1173 !
1174 ! rhs
1175 ns = natom + 1
1176 ALLOCATE (rhs(ns))
1177 rhs(1:natom) = chia(1:natom)
1178 rhs(ns) = -qtot
1179 !
1180 ALLOCATE (xv0(ns), rv0(ns))
1181 ! initial guess
1182 xv0(1:natom) = charges(1:natom)
1183 xv0(ns) = 0.0_dp
1184 ! DIIS optimizer
1185 max_diis = eeq_sparam%max_diis
1186 mdiis = eeq_sparam%mdiis
1187 sdiis = eeq_sparam%sdiis
1188 eps_diis = eeq_sparam%eps_diis
1189 astep = eeq_sparam%alpha
1190 ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1191 xvec = 0.0_dp; fvec = 0.0_dp
1192 ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1193 dmat = 0.0_dp; dvec = 0.0_dp
1194 ndiis = 1
1195 now = 1
1196 CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1197 cell, particle_set, xv0, rhs, rv0)
1198 resin = sqrt(sum(rv0*rv0))
1199 !
1200 DO iv = 1, max_diis
1201 res = sqrt(sum(rv0*rv0))
1202 IF (res > 10._dp*resin) EXIT
1203 IF (res < eps_diis) EXIT
1204 !
1205 now = mod(iv - 1, mdiis) + 1
1206 ndiis = min(iv, mdiis)
1207 xvec(1:ns, now) = xv0(1:ns)
1208 fvec(1:ns, now) = rv0(1:ns)
1209 DO i = 1, ndiis
1210 vmat(now, i) = sum(fvec(:, now)*fvec(:, i))
1211 vmat(i, now) = vmat(now, i)
1212 END DO
1213 IF (ndiis < sdiis) THEN
1214 xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1215 ELSE
1216 dvec = 0.0_dp
1217 dvec(ndiis + 1) = 1.0_dp
1218 dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1219 dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1220 dmat(1:ndiis, ndiis + 1) = 1.0_dp
1221 dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1222 CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1223 dvec(1:ndiis + 1) = matmul(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1224 xv0(1:ns) = matmul(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1225 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 cell, particle_set, xv0, rhs, rv0)
1230 END DO
1231 charges(1:natom) = xv0(1:natom)
1232 lambda = xv0(ns)
1233 eeq_energy = eeqn
1234 IF (res > eps_diis) ierror = 1
1235 !
1236 DEALLOCATE (xvec, fvec, bvec)
1237 DEALLOCATE (vmat, dmat, dvec)
1238 DEALLOCATE (xv0, rv0)
1239 DEALLOCATE (rhs)
1240 CALL cp_fm_release(pmat)
1241 CALL cp_fm_release(mmat)
1242
1243 te = m_walltime()
1244 IF (iunit > 0) THEN
1245 IF (ierror == 1) THEN
1246 WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1247 ELSE
1248 WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1249 END IF
1250 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1251 END IF
1252 CALL timestop(handle)
1253
1254 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 SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1274 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 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 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 CALL timeset(routinen, handle)
1307 ti = m_walltime()
1308
1309 iunit = -1
1310 IF (PRESENT(iounit)) iunit = iounit
1311
1312 natom = SIZE(particle_set)
1313 nkind = SIZE(gam)
1314 ns = natom + 1
1315 !
1316 CALL get_cell(cell=cell, deth=deth)
1317 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1318 ad = 2.0_dp*alpha*oorootpi
1319 IF (ewald_type /= do_ewald_spme) THEN
1320 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1321 END IF
1322 !
1323 rmax = 2.0_dp*rcut
1324 ! max cells used
1325 CALL get_cell(cell, h=hmat, periodic=periodic)
1326 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1327 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1328 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1329 IF (periodic(1) == 0) ncell(1) = 0
1330 IF (periodic(2) == 0) ncell(2) = 0
1331 IF (periodic(3) == 0) ncell(3) = 0
1332 !
1333 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1334 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 row_indices=rind, col_indices=cind)
1337 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1338 nrow_global=ns, ncol_global=1)
1339 CALL cp_fm_create(rhs_vec, vec_struct)
1340 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1341 row_indices=rvind, col_indices=cvind)
1342 ! response matrix
1343 DO ir = 1, nrloc
1344 iar = rind(ir)
1345 ri = 0.0_dp
1346 IF (iar <= natom) THEN
1347 ikind = kind_of(iar)
1348 ri(1:3) = particle_set(iar)%r(1:3)
1349 END IF
1350 DO ic = 1, ncloc
1351 iac = cind(ic)
1352 IF (iac > natom .AND. iar > natom) THEN
1353 eeq_mat%local_data(ir, ic) = 0.0_dp
1354 cycle
1355 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1356 eeq_mat%local_data(ir, ic) = 1.0_dp
1357 cycle
1358 END IF
1359 jkind = kind_of(iac)
1360 rj(1:3) = particle_set(iac)%r(1:3)
1361 rij(1:3) = ri(1:3) - rj(1:3)
1362 rij = pbc(rij, cell)
1363 DO ix = -ncell(1), ncell(1)
1364 DO iy = -ncell(2), ncell(2)
1365 DO iz = -ncell(3), ncell(3)
1366 cvec = [ix, iy, iz]
1367 rijl = rij + matmul(hmat, cvec)
1368 dr = sqrt(sum(rijl**2))
1369 IF (dr > rmax) cycle
1370 IF (iar == iac .AND. dr < 0.0001_dp) THEN
1371 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1372 ELSE
1373 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1374 END IF
1375 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 ALLOCATE (xval(natom), pval(natom))
1383 DO ia = 1, natom
1384 xval = 0.0_dp
1385 xval(ia) = 1.0_dp
1386 CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1387 !
1388 DO ir = 1, nrloc
1389 iar = rind(ir)
1390 IF (iar /= ia) cycle
1391 DO ic = 1, ncloc
1392 iac = cind(ic)
1393 IF (iac > natom) cycle
1394 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1395 END DO
1396 END DO
1397 END DO
1398 DEALLOCATE (xval, pval)
1399 !
1400 ! set up rhs vector
1401 DO ir = 1, nrvloc
1402 iar = rvind(ir)
1403 DO ic = 1, ncvloc
1404 iac = cvind(ic)
1405 ia = max(iar, iac)
1406 IF (ia > natom) THEN
1407 xr = qtot
1408 ELSE
1409 xr = -chia(ia)
1410 END IF
1411 rhs_vec%local_data(ir, ic) = xr
1412 END DO
1413 END DO
1414 !
1415 CALL cp_fm_solve(eeq_mat, rhs_vec)
1416 !
1417 charges = 0.0_dp
1418 lambda = 0.0_dp
1419 DO ir = 1, nrvloc
1420 iar = rvind(ir)
1421 DO ic = 1, ncvloc
1422 iac = cvind(ic)
1423 ia = max(iar, iac)
1424 IF (ia <= natom) THEN
1425 xr = rhs_vec%local_data(ir, ic)
1426 charges(ia) = xr
1427 ELSE
1428 lambda = rhs_vec%local_data(ir, ic)
1429 END IF
1430 END DO
1431 END DO
1432 CALL para_env%sum(lambda)
1433 CALL para_env%sum(charges)
1434 !
1435 ! energy: 0.5*(q^T.X - lambda*totalcharge)
1436 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1437
1438 CALL cp_fm_struct_release(vec_struct)
1439 CALL cp_fm_release(rhs_vec)
1440
1441 te = m_walltime()
1442 IF (iunit > 0) THEN
1443 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
1444 END IF
1445 CALL timestop(handle)
1446
1447 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 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 CALL ewald_env_get(ewald_env, para_env=para_env)
1469 potential = 0.0_dp
1470 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
1471 particle_set, potential)
1472 CALL para_env%sum(potential)
1473
1474 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 SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
1490 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 ns = SIZE(charges)
1507 na = ns - 1
1508 CALL ewald_env_get(ewald_env, para_env=para_env)
1509 potential = 0.0_dp
1510 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
1511 particle_set, potential(1:na))
1512 CALL para_env%sum(potential(1:na))
1513 CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
1514 eeqn = 0.5_dp*sum(charges(1:na)*potential(1:na)) + sum(charges(1:na)*rhs(1:na))
1515 potential(1:ns) = potential(1:ns) + rhs(1:ns)
1516 ALLOCATE (mvec(ns))
1517 CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
1518 lambda = -sum(mvec(1:na))/real(na, kind=dp)
1519 potential(1:na) = mvec(1:na) + lambda
1520 DEALLOCATE (mvec)
1521
1522 END SUBROUTINE get_energy_gradient
1523
1524! **************************************************************************************************
1525!> \brief ...
1526!> \param qs_env ...
1527!> \param charges ...
1528!> \param ef_energy ...
1529! **************************************************************************************************
1530 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1546
1547 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1548 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1549
1550 IF (dft_control%apply_period_efield) THEN
1551 dfield = dft_control%period_efield%displacement_field
1552
1553 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1554 cpabort("use of strength_list not implemented for eeq_efield_energy")
1555 END IF
1556
1557 fieldpol = dft_control%period_efield%polarisation
1558 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1559 fieldpol = -fieldpol*dft_control%period_efield%strength
1560 hmat = cell%hmat(:, :)/twopi
1561 DO idir = 1, 3
1562 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1563 + fieldpol(3)*hmat(3, idir)
1564 END DO
1565
1566 zi(:) = cmplx(1._dp, 0._dp, dp)
1567 DO ia = 1, natom
1568 q = charges(ia)
1569 ria = particle_set(ia)%r
1570 ria = pbc(ria, cell)
1571 DO idir = 1, 3
1572 kvec(:) = twopi*cell%h_inv(idir, :)
1573 kr = sum(kvec(:)*ria(:))
1574 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
1575 zi(idir) = zi(idir)*zdeta
1576 END DO
1577 END DO
1578 qi = aimag(log(zi))
1579 IF (dfield) THEN
1580 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1581 omega = cell%deth
1582 ci = matmul(hmat, qi)/omega
1583 ef_energy = 0.0_dp
1584 DO idir = 1, 3
1585 ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
1586 END DO
1587 ef_energy = -0.25_dp*omega/twopi*ef_energy
1588 ELSE
1589 ef_energy = sum(fpolvec(:)*qi(:))
1590 END IF
1591
1592 ELSE IF (dft_control%apply_efield) THEN
1593
1594 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1595 dft_control%efield_fields(1)%efield%strength
1596
1597 ef_energy = 0.0_dp
1598 DO ia = 1, natom
1599 ria = particle_set(ia)%r
1600 ria = pbc(ria, cell)
1601 q = charges(ia)
1602 ef_energy = ef_energy - q*sum(fieldpol*ria)
1603 END DO
1604
1605 ELSE
1606 cpabort("apply field")
1607 END IF
1608
1609 END SUBROUTINE eeq_efield_energy
1610
1611! **************************************************************************************************
1612!> \brief ...
1613!> \param qs_env ...
1614!> \param efr ...
1615! **************************************************************************************************
1616 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1628
1629 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1630 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1631
1632 IF (dft_control%apply_period_efield) THEN
1633 dfield = dft_control%period_efield%displacement_field
1634
1635 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1636 cpabort("use of strength_list not implemented for eeq_efield_pot")
1637 END IF
1638
1639 fieldpol = dft_control%period_efield%polarisation
1640 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1641 fieldpol = -fieldpol*dft_control%period_efield%strength
1642 hmat = cell%hmat(:, :)/twopi
1643 DO idir = 1, 3
1644 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1645 + fieldpol(3)*hmat(3, idir)
1646 END DO
1647
1648 IF (dfield) THEN
1649 ! dE/dq depends on q, postpone calculation
1650 efr = 0.0_dp
1651 ELSE
1652 efr = 0.0_dp
1653 DO ia = 1, natom
1654 ria = particle_set(ia)%r
1655 ria = pbc(ria, cell)
1656 DO idir = 1, 3
1657 kvec(:) = twopi*cell%h_inv(idir, :)
1658 kr = sum(kvec(:)*ria(:))
1659 efr(ia) = efr(ia) + kr*fpolvec(idir)
1660 END DO
1661 END DO
1662 END IF
1663
1664 ELSE IF (dft_control%apply_efield) THEN
1665
1666 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1667 dft_control%efield_fields(1)%efield%strength
1668
1669 DO ia = 1, natom
1670 ria = particle_set(ia)%r
1671 ria = pbc(ria, cell)
1672 efr(ia) = -sum(fieldpol*ria)
1673 END DO
1674
1675 ELSE
1676 cpabort("apply field")
1677 END IF
1678
1679 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 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 natom = SIZE(particle_set)
1704
1705 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1706 cpabort("use of strength_list not implemented for eeq_dfield_pot")
1707 END IF
1708
1709 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1710 fieldpol = dft_control%period_efield%polarisation
1711 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1712 fieldpol = -fieldpol*dft_control%period_efield%strength
1713 hmat = cell%hmat(:, :)/twopi
1714 omega = cell%deth
1715 !
1716 zi(:) = cmplx(1._dp, 0._dp, dp)
1717 DO ia = 1, natom
1718 q = charges(ia)
1719 ria = particle_set(ia)%r
1720 ria = pbc(ria, cell)
1721 DO idir = 1, 3
1722 kvec(:) = twopi*cell%h_inv(idir, :)
1723 kr = sum(kvec(:)*ria(:))
1724 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
1725 zi(idir) = zi(idir)*zdeta
1726 END DO
1727 END DO
1728 qi = aimag(log(zi))
1729 ci = matmul(hmat, qi)/omega
1730 ci = dfilter*(fieldpol - 2._dp*twopi*ci)
1731 DO ia = 1, natom
1732 ria = particle_set(ia)%r
1733 ria = pbc(ria, cell)
1734 efr(ia) = efr(ia) - sum(ci*ria)
1735 END DO
1736
1737 END SUBROUTINE eeq_dfield_pot
1738
1739! **************************************************************************************************
1740!> \brief ...
1741!> \param qs_env ...
1742!> \param charges ...
1743!> \param qlag ...
1744! **************************************************************************************************
1745 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1751 REAL(kind=dp) :: q
1752 REAL(kind=dp), DIMENSION(3) :: fieldpol
1753 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1759 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 local_particles=local_particles)
1767
1768 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1769 dft_control%efield_fields(1)%efield%strength
1770
1771 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1772 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1773 CALL get_qs_env(qs_env=qs_env, force=force)
1774
1775 DO ikind = 1, nkind
1776 force(ikind)%efield = 0.0_dp
1777 DO ia = 1, local_particles%n_el(ikind)
1778 iatom = local_particles%list(ikind)%array(ia)
1779 q = charges(iatom) - qlag(iatom)
1780 atom_a = atom_of_kind(iatom)
1781 force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
1782 END DO
1783 CALL para_env%sum(force(ikind)%efield)
1784 END DO
1785
1786 END SUBROUTINE eeq_efield_force_loc
1787
1788! **************************************************************************************************
1789!> \brief ...
1790!> \param qs_env ...
1791!> \param charges ...
1792!> \param qlag ...
1793! **************************************************************************************************
1794 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 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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1810 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 local_particles=local_particles)
1820
1821 dfield = dft_control%period_efield%displacement_field
1822 cpassert(.NOT. dfield)
1823
1824 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1825 cpabort("use of strength_list not implemented for eeq_efield_force_periodic")
1826 END IF
1827
1828 fieldpol = dft_control%period_efield%polarisation
1829 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1830 fieldpol = -fieldpol*dft_control%period_efield%strength
1831
1832 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1833
1834 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1835 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1836 CALL get_qs_env(qs_env=qs_env, force=force)
1837
1838 pve = 0.0_dp
1839 DO ikind = 1, nkind
1840 force(ikind)%efield = 0.0_dp
1841 DO ia = 1, local_particles%n_el(ikind)
1842 iatom = local_particles%list(ikind)%array(ia)
1843 q = charges(iatom) - qlag(iatom)
1844 fa(1:3) = q*fieldpol(1:3)
1845 atom_a = atom_of_kind(iatom)
1846 force(ikind)%efield(1:3, atom_a) = fa
1847 IF (use_virial) THEN
1848 ria = particle_set(ia)%r
1849 ria = pbc(ria, cell)
1850 CALL virial_pair_force(pve, -0.5_dp, fa, ria)
1851 CALL virial_pair_force(pve, -0.5_dp, ria, fa)
1852 END IF
1853 END DO
1854 CALL para_env%sum(force(ikind)%efield)
1855 END DO
1856 virial%pv_virial = virial%pv_virial + pve
1857
1858 END SUBROUTINE eeq_efield_force_periodic
1859
1860END MODULE eeq_method
cell_types::pbc
Definition cell_types.F:97
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:200
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:257
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_solve
subroutine, public cp_fm_solve(matrix_a, general_a)
computes the the solution to A*b=A_general using lu decomposition
Definition cp_fm_basic_linalg.F:2401
cp_fm_basic_linalg::cp_fm_invert
subroutine, public cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix.
Definition cp_fm_basic_linalg.F:2082
cp_fm_basic_linalg::cp_fm_matvec
subroutine, public cp_fm_matvec(amat, xv, yv, alpha, beta)
Calculates yv = alpha*amat*xv + beta*yv where amat: fm matrix xv : vector replicated yv : vector repl...
Definition cp_fm_basic_linalg.F:3254
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:556
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:73
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
eeq_data
EEQ data from different sources.
Definition eeq_data.F:12
eeq_data::get_eeq_data
subroutine, public get_eeq_data(za, model, chi, eta, kcn, rad)
...
Definition eeq_data.F:240
eeq_input
Input definition and setup for EEQ model.
Definition eeq_input.F:12
eeq_method
Calculation of charge equilibration method.
Definition eeq_method.F:12
eeq_method::eeq_forces
subroutine, public eeq_forces(qs_env, charges, dcharges, gradient, stress, eeq_sparam, eeq_model, enshift_type, response_only)
...
Definition eeq_method.F:347
eeq_method::eeq_efield_energy
subroutine, public eeq_efield_energy(qs_env, charges, ef_energy)
...
Definition eeq_method.F:1531
eeq_method::eeq_efield_force_periodic
subroutine, public eeq_efield_force_periodic(qs_env, charges, qlag)
...
Definition eeq_method.F:1795
eeq_method::eeq_efield_pot
subroutine, public eeq_efield_pot(qs_env, efr)
...
Definition eeq_method.F:1617
eeq_method::eeq_charges
subroutine, public eeq_charges(qs_env, charges, eeq_sparam, eeq_model, enshift_type)
...
Definition eeq_method.F:193
eeq_method::eeq_print
subroutine, public eeq_print(qs_env, iounit, print_level, ext)
...
Definition eeq_method.F:127
eeq_method::eeq_efield_force_loc
subroutine, public eeq_efield_force_loc(qs_env, charges, qlag)
...
Definition eeq_method.F:1746
eeq_method::eeq_solver
subroutine, public eeq_solver(charges, lambda, eeq_energy, particle_set, kind_of, cell, chia, gam, gab, para_env, blacs_env, dft_control, eeq_sparam, totalcharge, ewald, ewald_env, ewald_pw, iounit)
...
Definition eeq_method.F:769
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_set
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
Definition ewald_environment_types.F:181
ewald_environment_types::ewald_env_create
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
Definition ewald_environment_types.F:222
ewald_environment_types::ewald_env_release
subroutine, public ewald_env_release(ewald_env)
releases the given ewald_env (see doc/ReferenceCounting.html)
Definition ewald_environment_types.F:239
ewald_environment_types::read_ewald_section_tb
subroutine, public read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset)
Purpose: read the EWALD section for TB methods.
Definition ewald_environment_types.F:389
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:124
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
ewald_pw_types::ewald_pw_release
subroutine, public ewald_pw_release(ewald_pw)
releases the memory used by the ewald_pw
Definition ewald_pw_types.F:102
ewald_pw_types::ewald_pw_create
subroutine, public ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
creates the structure ewald_pw_type
Definition ewald_pw_types.F:83
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:153
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::invmat
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition mathlib.F:551
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:76
qs_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
qs_dispersion_cnum::cnumber_release
subroutine, public cnumber_release(cnumbers, dcnum, derivatives)
...
Definition qs_dispersion_cnum.F:1426
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1364
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_dispersion_types::qs_dispersion_release
subroutine, public qs_dispersion_release(dispersion_env)
...
Definition qs_dispersion_types.F:143
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1095
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition qs_neighbor_lists.F:19
qs_neighbor_lists::atom2d_cleanup
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
Definition qs_neighbor_lists.F:132
qs_neighbor_lists::pair_radius_setup
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
Definition qs_neighbor_lists.F:1665
qs_neighbor_lists::build_neighbor_lists
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
Definition qs_neighbor_lists.F:1067
qs_neighbor_lists::atom2d_build
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Definition qs_neighbor_lists.F:174
spme
Calculate the electrostatic energy by the Smooth Particle Ewald method.
Definition spme.F:14
spme::spme_forces
subroutine, public spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, charges_b, forces_b)
Calculate the forces on particles B for the electrostatic interaction betrween particles A and B.
Definition spme.F:568
spme::spme_potential
subroutine, public spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, potential)
Calculate the electrostatic potential from particles A (charge A) at positions of particles B.
Definition spme.F:427
spme::spme_virial
subroutine, public spme_virial(ewald_env, ewald_pw, particle_set, box, mcharge, virial)
Internal Virial for 1/2 [rho||rho] (rho=mcharge)
Definition spme.F:719
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:631
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
eeq_input::eeq_solver_type
Definition eeq_input.F:27
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
molecule_types::molecule_type
Definition molecule_types.F:159
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_neighbor_lists::local_atoms_type
Definition qs_neighbor_lists.F:105
virial_types::virial_type
Definition virial_types.F:26