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