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