(git:e5b1968)
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_displ, 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 do_displ = .false.
795 IF (dft_control%apply_period_efield .AND. ASSOCIATED(dft_control%period_efield)) THEN
796 do_displ = dft_control%period_efield%displacement_field
797 END IF
798
799 ! sparse option NYA
800 IF (do_sparse) THEN
801 CALL cp_abort(__location__, "EEQ: Sparse option not yet available")
802 END IF
803 !
804 natom = SIZE(particle_set)
805 nkind = SIZE(gam)
806 ns = natom + 1
807 CALL cp_fm_struct_create(mat_struct, context=blacs_env, para_env=para_env, &
808 nrow_global=ns, ncol_global=ns)
809 CALL cp_fm_create(eeq_mat, mat_struct)
810 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
811 !
812 IF (do_ewald) THEN
813 cpassert(PRESENT(ewald_env))
814 cpassert(PRESENT(ewald_pw))
815 IF (do_direct) THEN
816 IF (do_displ) THEN
817 cpabort("NYA")
818 ELSE
819 CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
820 kind_of, cell, chia, gam, gab, qtot, &
821 ewald_env, ewald_pw, iounit)
822 END IF
823 ELSE
824 IF (do_displ) THEN
825 cpabort("NYA")
826 ELSE
827 ierror = 0
828 CALL pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
829 kind_of, cell, chia, gam, gab, qtot, &
830 ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
831 IF (ierror /= 0) THEN
832 ! backup to non-iterative method
833 CALL fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
834 kind_of, cell, chia, gam, gab, qtot, &
835 ewald_env, ewald_pw, iounit)
836 END IF
837 END IF
838 END IF
839 IF (qtot /= 0._dp) THEN
840 CALL get_cell(cell=cell, deth=deth)
841 CALL ewald_env_get(ewald_env, alpha=alpha)
842 eeq_energy = eeq_energy - 0.5_dp*qtot**2/alpha**2/deth
843 END IF
844 ELSE
845 CALL mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, &
846 cell, chia, gam, gab, qtot, ftime)
847 IF (iounit > 0) THEN
848 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Molecular solver time[s]", ftime
849 END IF
850 END IF
851 CALL cp_fm_struct_release(mat_struct)
852 CALL cp_fm_release(eeq_mat)
853
854 CALL timestop(handle)
855
856 END SUBROUTINE eeq_solver
857
858! **************************************************************************************************
859!> \brief ...
860!> \param charges ...
861!> \param lambda ...
862!> \param eeq_energy ...
863!> \param eeq_mat ...
864!> \param particle_set ...
865!> \param kind_of ...
866!> \param cell ...
867!> \param chia ...
868!> \param gam ...
869!> \param gab ...
870!> \param qtot ...
871!> \param ftime ...
872! **************************************************************************************************
873 SUBROUTINE mi_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, kind_of, cell, &
874 chia, gam, gab, qtot, ftime)
875
876 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
877 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
878 TYPE(cp_fm_type) :: eeq_mat
879 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
880 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
881 TYPE(cell_type), POINTER :: cell
882 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
883 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
884 REAL(kind=dp), INTENT(IN) :: qtot
885 REAL(kind=dp), INTENT(OUT) :: ftime
886
887 CHARACTER(len=*), PARAMETER :: routinen = 'mi_solver'
888
889 INTEGER :: handle, ia, iac, iar, ic, ikind, ir, &
890 jkind, natom, ncloc, ncvloc, nkind, &
891 nrloc, nrvloc, ns
892 INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
893 REAL(kind=dp) :: dr, grc, te, ti, xr
894 REAL(kind=dp), DIMENSION(3) :: ri, rij, rj
895 TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
896 TYPE(cp_fm_type) :: rhs_vec
897 TYPE(mp_para_env_type), POINTER :: para_env
898
899 CALL timeset(routinen, handle)
900 ti = m_walltime()
901
902 natom = SIZE(particle_set)
903 nkind = SIZE(gam)
904 !
905 ns = natom + 1
906 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
907 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
908 row_indices=rind, col_indices=cind)
909 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
910 nrow_global=ns, ncol_global=1)
911 CALL cp_fm_create(rhs_vec, vec_struct)
912 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
913 row_indices=rvind, col_indices=cvind)
914 !
915 ! set up matrix
916 CALL cp_fm_set_all(eeq_mat, 1.0_dp, 0.0_dp)
917 CALL cp_fm_set_all(rhs_vec, 0.0_dp)
918 DO ir = 1, nrloc
919 iar = rind(ir)
920 IF (iar > natom) cycle
921 ikind = kind_of(iar)
922 ri(1:3) = particle_set(iar)%r(1:3)
923 DO ic = 1, ncloc
924 iac = cind(ic)
925 IF (iac > natom) cycle
926 jkind = kind_of(iac)
927 rj(1:3) = particle_set(iac)%r(1:3)
928 IF (iar == iac) THEN
929 grc = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
930 ELSE
931 rij(1:3) = ri(1:3) - rj(1:3)
932 rij = pbc(rij, cell)
933 dr = sqrt(sum(rij**2))
934 grc = erf(gab(ikind, jkind)*dr)/dr
935 END IF
936 eeq_mat%local_data(ir, ic) = grc
937 END DO
938 END DO
939 ! set up rhs vector
940 DO ir = 1, nrvloc
941 iar = rvind(ir)
942 DO ic = 1, ncvloc
943 iac = cvind(ic)
944 ia = max(iar, iac)
945 IF (ia > natom) THEN
946 xr = qtot
947 ELSE
948 xr = -chia(ia)
949 END IF
950 rhs_vec%local_data(ir, ic) = xr
951 END DO
952 END DO
953 !
954 CALL cp_fm_solve(eeq_mat, rhs_vec)
955 !
956 charges = 0.0_dp
957 lambda = 0.0_dp
958 DO ir = 1, nrvloc
959 iar = rvind(ir)
960 DO ic = 1, ncvloc
961 iac = cvind(ic)
962 ia = max(iar, iac)
963 IF (ia <= natom) THEN
964 xr = rhs_vec%local_data(ir, ic)
965 charges(ia) = xr
966 ELSE
967 lambda = rhs_vec%local_data(ir, ic)
968 END IF
969 END DO
970 END DO
971 CALL para_env%sum(lambda)
972 CALL para_env%sum(charges)
973 !
974 ! energy: 0.5*(q^T.X - lambda*totalcharge)
975 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
976
977 CALL cp_fm_struct_release(vec_struct)
978 CALL cp_fm_release(rhs_vec)
979
980 te = m_walltime()
981 ftime = te - ti
982 CALL timestop(handle)
983
984 END SUBROUTINE mi_solver
985
986! **************************************************************************************************
987!> \brief ...
988!> \param charges ...
989!> \param lambda ...
990!> \param eeq_energy ...
991!> \param eeq_mat ...
992!> \param particle_set ...
993!> \param kind_of ...
994!> \param cell ...
995!> \param chia ...
996!> \param gam ...
997!> \param gab ...
998!> \param qtot ...
999!> \param ewald_env ...
1000!> \param ewald_pw ...
1001!> \param eeq_sparam ...
1002!> \param ierror ...
1003!> \param iounit ...
1004! **************************************************************************************************
1005 SUBROUTINE pbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1006 kind_of, cell, chia, gam, gab, qtot, &
1007 ewald_env, ewald_pw, eeq_sparam, ierror, iounit)
1008
1009 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
1010 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
1011 TYPE(cp_fm_type) :: eeq_mat
1012 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1013 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1014 TYPE(cell_type), POINTER :: cell
1015 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1016 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1017 REAL(kind=dp), INTENT(IN) :: qtot
1018 TYPE(ewald_environment_type), POINTER :: ewald_env
1019 TYPE(ewald_pw_type), POINTER :: ewald_pw
1020 TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
1021 INTEGER, INTENT(OUT) :: ierror
1022 INTEGER, OPTIONAL :: iounit
1023
1024 CHARACTER(len=*), PARAMETER :: routinen = 'pbc_solver'
1025
1026 INTEGER :: ewald_type, handle, i, iac, iar, ic, ikind, info, ir, iunit, iv, ix, iy, iz, &
1027 jkind, max_diis, mdiis, natom, ncloc, ndiis, nkind, now, nrloc, ns, sdiis
1028 INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1029 INTEGER, DIMENSION(:), POINTER :: cind, rind
1030 REAL(kind=dp) :: ad, alpha, astep, deth, dr, eeqn, &
1031 eps_diis, ftime, grc1, grc2, rcut, &
1032 res, resin, rmax, te, ti
1033 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bvec, dvec
1034 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat, fvec, vmat, xvec
1035 REAL(kind=dp), DIMENSION(3) :: ri, rij, rijl, rj
1036 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1037 REAL(kind=dp), DIMENSION(:), POINTER :: rhs, rv0, xv0
1038 TYPE(cp_fm_struct_type), POINTER :: mat_struct
1039 TYPE(cp_fm_type) :: mmat, pmat
1040 TYPE(mp_para_env_type), POINTER :: para_env
1041
1042 CALL timeset(routinen, handle)
1043 ti = m_walltime()
1044
1045 ierror = 0
1046
1047 iunit = -1
1048 IF (PRESENT(iounit)) iunit = iounit
1049
1050 natom = SIZE(particle_set)
1051 nkind = SIZE(gam)
1052 !
1053 CALL get_cell(cell=cell, deth=deth)
1054 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1055 ad = 2.0_dp*alpha*oorootpi
1056 IF (ewald_type /= do_ewald_spme) THEN
1057 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1058 END IF
1059 !
1060 rmax = 2.0_dp*rcut
1061 ! max cells used
1062 CALL get_cell(cell, h=hmat, periodic=periodic)
1063 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1064 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1065 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1066 IF (periodic(1) == 0) ncell(1) = 0
1067 IF (periodic(2) == 0) ncell(2) = 0
1068 IF (periodic(3) == 0) ncell(3) = 0
1069 !
1070 CALL mi_solver(charges, lambda, eeqn, eeq_mat, particle_set, kind_of, cell, &
1071 chia, gam, gab, qtot, ftime)
1072 IF (iunit > 0) THEN
1073 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC guess time[s]", ftime
1074 END IF
1075 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1076 CALL cp_fm_create(pmat, mat_struct)
1077 CALL cp_fm_create(mmat, mat_struct)
1078 !
1079 ! response matrix
1080 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1081 row_indices=rind, col_indices=cind)
1082 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1083 DO ir = 1, nrloc
1084 iar = rind(ir)
1085 ri = 0.0_dp
1086 IF (iar <= natom) THEN
1087 ikind = kind_of(iar)
1088 ri(1:3) = particle_set(iar)%r(1:3)
1089 END IF
1090 DO ic = 1, ncloc
1091 iac = cind(ic)
1092 IF (iac > natom .AND. iar > natom) THEN
1093 eeq_mat%local_data(ir, ic) = 0.0_dp
1094 cycle
1095 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1096 eeq_mat%local_data(ir, ic) = 1.0_dp
1097 cycle
1098 END IF
1099 jkind = kind_of(iac)
1100 rj(1:3) = particle_set(iac)%r(1:3)
1101 rij(1:3) = ri(1:3) - rj(1:3)
1102 rij = pbc(rij, cell)
1103 DO ix = -ncell(1), ncell(1)
1104 DO iy = -ncell(2), ncell(2)
1105 DO iz = -ncell(3), ncell(3)
1106 cvec = [ix, iy, iz]
1107 rijl = rij + matmul(hmat, cvec)
1108 dr = sqrt(sum(rijl**2))
1109 IF (dr > rmax) cycle
1110 IF (iar == iac .AND. dr < 0.00001_dp) THEN
1111 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1112 ELSE
1113 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1114 END IF
1115 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1116 END DO
1117 END DO
1118 END DO
1119 END DO
1120 END DO
1121 !
1122 ! preconditioner
1123 CALL cp_fm_get_info(pmat, nrow_local=nrloc, ncol_local=ncloc, &
1124 row_indices=rind, col_indices=cind)
1125 CALL cp_fm_set_all(pmat, 0.0_dp, 0.0_dp)
1126 DO ir = 1, nrloc
1127 iar = rind(ir)
1128 ri = 0.0_dp
1129 IF (iar <= natom) THEN
1130 ikind = kind_of(iar)
1131 ri(1:3) = particle_set(iar)%r(1:3)
1132 END IF
1133 DO ic = 1, ncloc
1134 iac = cind(ic)
1135 IF (iac > natom .AND. iar > natom) THEN
1136 pmat%local_data(ir, ic) = 0.0_dp
1137 cycle
1138 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1139 pmat%local_data(ir, ic) = 1.0_dp
1140 cycle
1141 END IF
1142 jkind = kind_of(iac)
1143 rj(1:3) = particle_set(iac)%r(1:3)
1144 rij(1:3) = ri(1:3) - rj(1:3)
1145 rij = pbc(rij, cell)
1146 IF (iar == iac) THEN
1147 grc2 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi
1148 ELSE
1149 grc2 = erf(gab(ikind, jkind)*dr)/dr
1150 END IF
1151 pmat%local_data(ir, ic) = grc2
1152 END DO
1153 END DO
1154 CALL cp_fm_set_all(mmat, 0.0_dp, 0.0_dp)
1155 ! preconditioner invers
1156 CALL cp_fm_invert(pmat, mmat)
1157 !
1158 ! rhs
1159 ns = natom + 1
1160 ALLOCATE (rhs(ns))
1161 rhs(1:natom) = chia(1:natom)
1162 rhs(ns) = -qtot
1163 !
1164 ALLOCATE (xv0(ns), rv0(ns))
1165 ! initial guess
1166 xv0(1:natom) = charges(1:natom)
1167 xv0(ns) = 0.0_dp
1168 ! DIIS optimizer
1169 max_diis = eeq_sparam%max_diis
1170 mdiis = eeq_sparam%mdiis
1171 sdiis = eeq_sparam%sdiis
1172 eps_diis = eeq_sparam%eps_diis
1173 astep = eeq_sparam%alpha
1174 ALLOCATE (xvec(ns, mdiis), fvec(ns, mdiis), bvec(ns))
1175 xvec = 0.0_dp; fvec = 0.0_dp
1176 ALLOCATE (vmat(mdiis, mdiis), dmat(mdiis + 1, mdiis + 1), dvec(mdiis + 1))
1177 dmat = 0.0_dp; dvec = 0.0_dp
1178 ndiis = 1
1179 now = 1
1180 CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1181 cell, particle_set, xv0, rhs, rv0)
1182 resin = sqrt(sum(rv0*rv0))
1183 !
1184 DO iv = 1, max_diis
1185 res = sqrt(sum(rv0*rv0))
1186 IF (res > 10._dp*resin) EXIT
1187 IF (res < eps_diis) EXIT
1188 !
1189 now = mod(iv - 1, mdiis) + 1
1190 ndiis = min(iv, mdiis)
1191 xvec(1:ns, now) = xv0(1:ns)
1192 fvec(1:ns, now) = rv0(1:ns)
1193 DO i = 1, ndiis
1194 vmat(now, i) = sum(fvec(:, now)*fvec(:, i))
1195 vmat(i, now) = vmat(now, i)
1196 END DO
1197 IF (ndiis < sdiis) THEN
1198 xv0(1:ns) = xv0(1:ns) + astep*rv0(1:ns)
1199 ELSE
1200 dvec = 0.0_dp
1201 dvec(ndiis + 1) = 1.0_dp
1202 dmat(1:ndiis, 1:ndiis) = vmat(1:ndiis, 1:ndiis)
1203 dmat(ndiis + 1, 1:ndiis) = 1.0_dp
1204 dmat(1:ndiis, ndiis + 1) = 1.0_dp
1205 dmat(ndiis + 1, ndiis + 1) = 0.0_dp
1206 CALL invmat(dmat(1:ndiis + 1, 1:ndiis + 1), info)
1207 dvec(1:ndiis + 1) = matmul(dmat(1:ndiis + 1, 1:ndiis + 1), dvec(1:ndiis + 1))
1208 xv0(1:ns) = matmul(xvec(1:ns, 1:ndiis), dvec(1:ndiis))
1209 xv0(1:ns) = xv0(1:ns) + matmul(fvec(1:ns, 1:ndiis), dvec(1:ndiis))
1210 END IF
1211 !
1212 CALL get_energy_gradient(eeqn, eeq_mat, mmat, ewald_env, ewald_pw, &
1213 cell, particle_set, xv0, rhs, rv0)
1214 END DO
1215 charges(1:natom) = xv0(1:natom)
1216 lambda = xv0(ns)
1217 eeq_energy = eeqn
1218 IF (res > eps_diis) ierror = 1
1219 !
1220 DEALLOCATE (xvec, fvec, bvec)
1221 DEALLOCATE (vmat, dmat, dvec)
1222 DEALLOCATE (xv0, rv0)
1223 DEALLOCATE (rhs)
1224 CALL cp_fm_release(pmat)
1225 CALL cp_fm_release(mmat)
1226
1227 te = m_walltime()
1228 IF (iunit > 0) THEN
1229 IF (ierror == 1) THEN
1230 WRITE (iunit, '(A)') " EEQ| PBC solver failed to converge "
1231 ELSE
1232 WRITE (iunit, '(A,T50,I4,T61,E20.5)') " EEQ| PBC solver: iterations/accuracy ", iv, res
1233 END IF
1234 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Iterative PBC solver: time[s]", te - ti
1235 END IF
1236 CALL timestop(handle)
1237
1238 END SUBROUTINE pbc_solver
1239
1240! **************************************************************************************************
1241!> \brief ...
1242!> \param charges ...
1243!> \param lambda ...
1244!> \param eeq_energy ...
1245!> \param eeq_mat ...
1246!> \param particle_set ...
1247!> \param kind_of ...
1248!> \param cell ...
1249!> \param chia ...
1250!> \param gam ...
1251!> \param gab ...
1252!> \param qtot ...
1253!> \param ewald_env ...
1254!> \param ewald_pw ...
1255!> \param iounit ...
1256! **************************************************************************************************
1257 SUBROUTINE fpbc_solver(charges, lambda, eeq_energy, eeq_mat, particle_set, &
1258 kind_of, cell, chia, gam, gab, qtot, ewald_env, ewald_pw, iounit)
1259
1260 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: charges
1261 REAL(kind=dp), INTENT(INOUT) :: lambda, eeq_energy
1262 TYPE(cp_fm_type) :: eeq_mat
1263 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1264 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1265 TYPE(cell_type), POINTER :: cell
1266 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: chia, gam
1267 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gab
1268 REAL(kind=dp), INTENT(IN) :: qtot
1269 TYPE(ewald_environment_type), POINTER :: ewald_env
1270 TYPE(ewald_pw_type), POINTER :: ewald_pw
1271 INTEGER, INTENT(IN), OPTIONAL :: iounit
1272
1273 CHARACTER(len=*), PARAMETER :: routinen = 'fpbc_solver'
1274
1275 INTEGER :: ewald_type, handle, ia, iac, iar, ic, &
1276 ikind, ir, iunit, ix, iy, iz, jkind, &
1277 natom, ncloc, ncvloc, nkind, nrloc, &
1278 nrvloc, ns
1279 INTEGER, DIMENSION(3) :: cvec, ncell, periodic
1280 INTEGER, DIMENSION(:), POINTER :: cind, cvind, rind, rvind
1281 REAL(kind=dp) :: ad, alpha, deth, dr, grc1, rcut, rmax, &
1282 te, ti, xr
1283 REAL(kind=dp), DIMENSION(3) :: ri, rij, rijl, rj
1284 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1285 REAL(kind=dp), DIMENSION(:), POINTER :: pval, xval
1286 TYPE(cp_fm_struct_type), POINTER :: mat_struct, vec_struct
1287 TYPE(cp_fm_type) :: rhs_vec
1288 TYPE(mp_para_env_type), POINTER :: para_env
1289
1290 CALL timeset(routinen, handle)
1291 ti = m_walltime()
1292
1293 iunit = -1
1294 IF (PRESENT(iounit)) iunit = iounit
1295
1296 natom = SIZE(particle_set)
1297 nkind = SIZE(gam)
1298 ns = natom + 1
1299 !
1300 CALL get_cell(cell=cell, deth=deth)
1301 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut, ewald_type=ewald_type)
1302 ad = 2.0_dp*alpha*oorootpi
1303 IF (ewald_type /= do_ewald_spme) THEN
1304 CALL cp_abort(__location__, "Only SPME Ewald method available with EEQ.")
1305 END IF
1306 !
1307 rmax = 2.0_dp*rcut
1308 ! max cells used
1309 CALL get_cell(cell, h=hmat, periodic=periodic)
1310 ncell(1) = ceiling(rmax/plane_distance(1, 0, 0, cell))
1311 ncell(2) = ceiling(rmax/plane_distance(0, 1, 0, cell))
1312 ncell(3) = ceiling(rmax/plane_distance(0, 0, 1, cell))
1313 IF (periodic(1) == 0) ncell(1) = 0
1314 IF (periodic(2) == 0) ncell(2) = 0
1315 IF (periodic(3) == 0) ncell(3) = 0
1316 !
1317 CALL cp_fm_get_info(eeq_mat, matrix_struct=mat_struct, para_env=para_env)
1318 CALL cp_fm_set_all(eeq_mat, 0.0_dp, 0.0_dp)
1319 CALL cp_fm_get_info(eeq_mat, nrow_local=nrloc, ncol_local=ncloc, &
1320 row_indices=rind, col_indices=cind)
1321 CALL cp_fm_struct_create(vec_struct, template_fmstruct=mat_struct, &
1322 nrow_global=ns, ncol_global=1)
1323 CALL cp_fm_create(rhs_vec, vec_struct)
1324 CALL cp_fm_get_info(rhs_vec, nrow_local=nrvloc, ncol_local=ncvloc, &
1325 row_indices=rvind, col_indices=cvind)
1326 ! response matrix
1327 DO ir = 1, nrloc
1328 iar = rind(ir)
1329 ri = 0.0_dp
1330 IF (iar <= natom) THEN
1331 ikind = kind_of(iar)
1332 ri(1:3) = particle_set(iar)%r(1:3)
1333 END IF
1334 DO ic = 1, ncloc
1335 iac = cind(ic)
1336 IF (iac > natom .AND. iar > natom) THEN
1337 eeq_mat%local_data(ir, ic) = 0.0_dp
1338 cycle
1339 ELSE IF ((iac > natom) .OR. (iar > natom)) THEN
1340 eeq_mat%local_data(ir, ic) = 1.0_dp
1341 cycle
1342 END IF
1343 jkind = kind_of(iac)
1344 rj(1:3) = particle_set(iac)%r(1:3)
1345 rij(1:3) = ri(1:3) - rj(1:3)
1346 rij = pbc(rij, cell)
1347 DO ix = -ncell(1), ncell(1)
1348 DO iy = -ncell(2), ncell(2)
1349 DO iz = -ncell(3), ncell(3)
1350 cvec = [ix, iy, iz]
1351 rijl = rij + matmul(hmat, cvec)
1352 dr = sqrt(sum(rijl**2))
1353 IF (dr > rmax) cycle
1354 IF (iar == iac .AND. dr < 0.0001_dp) THEN
1355 grc1 = gam(ikind) + 2.0_dp*gab(ikind, ikind)*oorootpi - ad
1356 ELSE
1357 grc1 = erf(gab(ikind, jkind)*dr)/dr - erf(alpha*dr)/dr
1358 END IF
1359 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + grc1
1360 END DO
1361 END DO
1362 END DO
1363 END DO
1364 END DO
1365 !
1366 ALLOCATE (xval(natom), pval(natom))
1367 DO ia = 1, natom
1368 xval = 0.0_dp
1369 xval(ia) = 1.0_dp
1370 CALL apply_potential(ewald_env, ewald_pw, cell, particle_set, xval, pval)
1371 !
1372 DO ir = 1, nrloc
1373 iar = rind(ir)
1374 IF (iar /= ia) cycle
1375 DO ic = 1, ncloc
1376 iac = cind(ic)
1377 IF (iac > natom) cycle
1378 eeq_mat%local_data(ir, ic) = eeq_mat%local_data(ir, ic) + pval(iac)
1379 END DO
1380 END DO
1381 END DO
1382 DEALLOCATE (xval, pval)
1383 !
1384 ! set up rhs vector
1385 DO ir = 1, nrvloc
1386 iar = rvind(ir)
1387 DO ic = 1, ncvloc
1388 iac = cvind(ic)
1389 ia = max(iar, iac)
1390 IF (ia > natom) THEN
1391 xr = qtot
1392 ELSE
1393 xr = -chia(ia)
1394 END IF
1395 rhs_vec%local_data(ir, ic) = xr
1396 END DO
1397 END DO
1398 !
1399 CALL cp_fm_solve(eeq_mat, rhs_vec)
1400 !
1401 charges = 0.0_dp
1402 lambda = 0.0_dp
1403 DO ir = 1, nrvloc
1404 iar = rvind(ir)
1405 DO ic = 1, ncvloc
1406 iac = cvind(ic)
1407 ia = max(iar, iac)
1408 IF (ia <= natom) THEN
1409 xr = rhs_vec%local_data(ir, ic)
1410 charges(ia) = xr
1411 ELSE
1412 lambda = rhs_vec%local_data(ir, ic)
1413 END IF
1414 END DO
1415 END DO
1416 CALL para_env%sum(lambda)
1417 CALL para_env%sum(charges)
1418 !
1419 ! energy: 0.5*(q^T.X - lambda*totalcharge)
1420 eeq_energy = 0.5*sum(charges(1:natom)*chia(1:natom)) - 0.5_dp*lambda*qtot
1421
1422 CALL cp_fm_struct_release(vec_struct)
1423 CALL cp_fm_release(rhs_vec)
1424
1425 te = m_walltime()
1426 IF (iunit > 0) THEN
1427 WRITE (iunit, '(A,T67,F14.3)') " EEQ| Direct PBC solver: time[s]", te - ti
1428 END IF
1429 CALL timestop(handle)
1430
1431 END SUBROUTINE fpbc_solver
1432
1433! **************************************************************************************************
1434!> \brief ...
1435!> \param ewald_env ...
1436!> \param ewald_pw ...
1437!> \param cell ...
1438!> \param particle_set ...
1439!> \param charges ...
1440!> \param potential ...
1441! **************************************************************************************************
1442 SUBROUTINE apply_potential(ewald_env, ewald_pw, cell, particle_set, charges, potential)
1443 TYPE(ewald_environment_type), POINTER :: ewald_env
1444 TYPE(ewald_pw_type), POINTER :: ewald_pw
1445 TYPE(cell_type), POINTER :: cell
1446 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1447 REAL(kind=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1448 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: potential
1449
1450 TYPE(mp_para_env_type), POINTER :: para_env
1451
1452 CALL ewald_env_get(ewald_env, para_env=para_env)
1453 potential = 0.0_dp
1454 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges, &
1455 particle_set, potential)
1456 CALL para_env%sum(potential)
1457
1458 END SUBROUTINE apply_potential
1459
1460! **************************************************************************************************
1461!> \brief ...
1462!> \param eeqn ...
1463!> \param fm_mat ...
1464!> \param mmat ...
1465!> \param ewald_env ...
1466!> \param ewald_pw ...
1467!> \param cell ...
1468!> \param particle_set ...
1469!> \param charges ...
1470!> \param rhs ...
1471!> \param potential ...
1472! **************************************************************************************************
1473 SUBROUTINE get_energy_gradient(eeqn, fm_mat, mmat, ewald_env, ewald_pw, &
1474 cell, particle_set, charges, rhs, potential)
1475 REAL(kind=dp), INTENT(INOUT) :: eeqn
1476 TYPE(cp_fm_type), INTENT(IN) :: fm_mat, mmat
1477 TYPE(ewald_environment_type), POINTER :: ewald_env
1478 TYPE(ewald_pw_type), POINTER :: ewald_pw
1479 TYPE(cell_type), POINTER :: cell
1480 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1481 REAL(kind=dp), DIMENSION(:), INTENT(IN), POINTER :: charges
1482 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rhs
1483 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: potential
1484
1485 INTEGER :: na, ns
1486 REAL(kind=dp) :: lambda
1487 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mvec
1488 TYPE(mp_para_env_type), POINTER :: para_env
1489
1490 ns = SIZE(charges)
1491 na = ns - 1
1492 CALL ewald_env_get(ewald_env, para_env=para_env)
1493 potential = 0.0_dp
1494 CALL spme_potential(ewald_env, ewald_pw, cell, particle_set, charges(1:na), &
1495 particle_set, potential(1:na))
1496 CALL para_env%sum(potential(1:na))
1497 CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
1498 eeqn = 0.5_dp*sum(charges(1:na)*potential(1:na)) + sum(charges(1:na)*rhs(1:na))
1499 potential(1:ns) = potential(1:ns) + rhs(1:ns)
1500 ALLOCATE (mvec(ns))
1501 CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
1502 lambda = -sum(mvec(1:na))/real(na, kind=dp)
1503 potential(1:na) = mvec(1:na) + lambda
1504 DEALLOCATE (mvec)
1505
1506 END SUBROUTINE get_energy_gradient
1507
1508! **************************************************************************************************
1509!> \brief ...
1510!> \param qs_env ...
1511!> \param charges ...
1512!> \param ef_energy ...
1513! **************************************************************************************************
1514 SUBROUTINE eeq_efield_energy(qs_env, charges, ef_energy)
1515 TYPE(qs_environment_type), POINTER :: qs_env
1516 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges
1517 REAL(kind=dp), INTENT(OUT) :: ef_energy
1518
1519 COMPLEX(KIND=dp) :: zdeta
1520 COMPLEX(KIND=dp), DIMENSION(3) :: zi
1521 INTEGER :: ia, idir, natom
1522 LOGICAL :: dfield
1523 REAL(kind=dp) :: kr, omega, q
1524 REAL(kind=dp), DIMENSION(3) :: ci, dfilter, fieldpol, fpolvec, kvec, &
1525 qi, ria
1526 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1527 TYPE(cell_type), POINTER :: cell
1528 TYPE(dft_control_type), POINTER :: dft_control
1529 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1530
1531 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1532 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1533
1534 IF (dft_control%apply_period_efield) THEN
1535 dfield = dft_control%period_efield%displacement_field
1536
1537 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1538 cpabort("use of strength_list not implemented for eeq_efield_energy")
1539 END IF
1540
1541 fieldpol = dft_control%period_efield%polarisation
1542 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1543 fieldpol = -fieldpol*dft_control%period_efield%strength
1544 hmat = cell%hmat(:, :)/twopi
1545 DO idir = 1, 3
1546 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1547 + fieldpol(3)*hmat(3, idir)
1548 END DO
1549
1550 zi(:) = cmplx(1._dp, 0._dp, dp)
1551 DO ia = 1, natom
1552 q = charges(ia)
1553 ria = particle_set(ia)%r
1554 ria = pbc(ria, cell)
1555 DO idir = 1, 3
1556 kvec(:) = twopi*cell%h_inv(idir, :)
1557 kr = sum(kvec(:)*ria(:))
1558 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
1559 zi(idir) = zi(idir)*zdeta
1560 END DO
1561 END DO
1562 qi = aimag(log(zi))
1563 IF (dfield) THEN
1564 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1565 omega = cell%deth
1566 ci = matmul(hmat, qi)/omega
1567 ef_energy = 0.0_dp
1568 DO idir = 1, 3
1569 ef_energy = ef_energy + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi*ci(idir))**2
1570 END DO
1571 ef_energy = -0.25_dp*omega/twopi*ef_energy
1572 ELSE
1573 ef_energy = sum(fpolvec(:)*qi(:))
1574 END IF
1575
1576 ELSE IF (dft_control%apply_efield) THEN
1577
1578 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1579 dft_control%efield_fields(1)%efield%strength
1580
1581 ef_energy = 0.0_dp
1582 DO ia = 1, natom
1583 ria = particle_set(ia)%r
1584 ria = pbc(ria, cell)
1585 q = charges(ia)
1586 ef_energy = ef_energy - q*sum(fieldpol*ria)
1587 END DO
1588
1589 ELSE
1590 cpabort("apply field")
1591 END IF
1592
1593 END SUBROUTINE eeq_efield_energy
1594
1595! **************************************************************************************************
1596!> \brief ...
1597!> \param qs_env ...
1598!> \param efr ...
1599! **************************************************************************************************
1600 SUBROUTINE eeq_efield_pot(qs_env, efr)
1601 TYPE(qs_environment_type), POINTER :: qs_env
1602 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: efr
1603
1604 INTEGER :: ia, idir, natom
1605 LOGICAL :: dfield
1606 REAL(kind=dp) :: kr
1607 REAL(kind=dp), DIMENSION(3) :: fieldpol, fpolvec, kvec, ria
1608 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1609 TYPE(cell_type), POINTER :: cell
1610 TYPE(dft_control_type), POINTER :: dft_control
1611 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1612
1613 CALL get_qs_env(qs_env, natom=natom, dft_control=dft_control)
1614 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1615
1616 IF (dft_control%apply_period_efield) THEN
1617 dfield = dft_control%period_efield%displacement_field
1618
1619 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1620 cpabort("use of strength_list not implemented for eeq_efield_pot")
1621 END IF
1622
1623 fieldpol = dft_control%period_efield%polarisation
1624 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1625 fieldpol = -fieldpol*dft_control%period_efield%strength
1626 hmat = cell%hmat(:, :)/twopi
1627 DO idir = 1, 3
1628 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) &
1629 + fieldpol(3)*hmat(3, idir)
1630 END DO
1631
1632 IF (dfield) THEN
1633 ! dE/dq depends on q, postpone calculation
1634 efr = 0.0_dp
1635 ELSE
1636 efr = 0.0_dp
1637 DO ia = 1, natom
1638 ria = particle_set(ia)%r
1639 ria = pbc(ria, cell)
1640 DO idir = 1, 3
1641 kvec(:) = twopi*cell%h_inv(idir, :)
1642 kr = sum(kvec(:)*ria(:))
1643 efr(ia) = efr(ia) + kr*fpolvec(idir)
1644 END DO
1645 END DO
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 DO ia = 1, natom
1654 ria = particle_set(ia)%r
1655 ria = pbc(ria, cell)
1656 efr(ia) = -sum(fieldpol*ria)
1657 END DO
1658
1659 ELSE
1660 cpabort("apply field")
1661 END IF
1662
1663 END SUBROUTINE eeq_efield_pot
1664
1665! **************************************************************************************************
1666!> \brief ...
1667!> \param charges ...
1668!> \param dft_control ...
1669!> \param particle_set ...
1670!> \param cell ...
1671!> \param efr ...
1672! **************************************************************************************************
1673 SUBROUTINE eeq_dfield_pot(charges, dft_control, particle_set, cell, efr)
1674 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges
1675 TYPE(dft_control_type), POINTER :: dft_control
1676 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
1677 TYPE(cell_type), POINTER :: cell
1678 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: efr
1679
1680 COMPLEX(KIND=dp) :: zdeta
1681 COMPLEX(KIND=dp), DIMENSION(3) :: zi
1682 INTEGER :: ia, idir, natom
1683 REAL(kind=dp) :: kr, omega, q
1684 REAL(kind=dp), DIMENSION(3) :: ci, dfilter, fieldpol, kvec, qi, ria
1685 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1686
1687 natom = SIZE(particle_set)
1688
1689 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1690 cpabort("use of strength_list not implemented for eeq_dfield_pot")
1691 END IF
1692
1693 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
1694 fieldpol = dft_control%period_efield%polarisation
1695 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1696 fieldpol = -fieldpol*dft_control%period_efield%strength
1697 hmat = cell%hmat(:, :)/twopi
1698 omega = cell%deth
1699 !
1700 zi(:) = cmplx(1._dp, 0._dp, dp)
1701 DO ia = 1, natom
1702 q = charges(ia)
1703 ria = particle_set(ia)%r
1704 ria = pbc(ria, cell)
1705 DO idir = 1, 3
1706 kvec(:) = twopi*cell%h_inv(idir, :)
1707 kr = sum(kvec(:)*ria(:))
1708 zdeta = cmplx(cos(kr), sin(kr), kind=dp)**q
1709 zi(idir) = zi(idir)*zdeta
1710 END DO
1711 END DO
1712 qi = aimag(log(zi))
1713 ci = matmul(hmat, qi)/omega
1714 ci = dfilter*(fieldpol - 2._dp*twopi*ci)
1715 DO ia = 1, natom
1716 ria = particle_set(ia)%r
1717 ria = pbc(ria, cell)
1718 efr(ia) = efr(ia) - sum(ci*ria)
1719 END DO
1720
1721 END SUBROUTINE eeq_dfield_pot
1722
1723! **************************************************************************************************
1724!> \brief ...
1725!> \param qs_env ...
1726!> \param charges ...
1727!> \param qlag ...
1728! **************************************************************************************************
1729 SUBROUTINE eeq_efield_force_loc(qs_env, charges, qlag)
1730 TYPE(qs_environment_type), POINTER :: qs_env
1731 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1732
1733 INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1734 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1735 REAL(kind=dp) :: q
1736 REAL(kind=dp), DIMENSION(3) :: fieldpol
1737 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1738 TYPE(cell_type), POINTER :: cell
1739 TYPE(dft_control_type), POINTER :: dft_control
1740 TYPE(distribution_1d_type), POINTER :: local_particles
1741 TYPE(mp_para_env_type), POINTER :: para_env
1742 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1743 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1744
1745 CALL get_qs_env(qs_env=qs_env, &
1746 dft_control=dft_control, &
1747 cell=cell, particle_set=particle_set, &
1748 nkind=nkind, natom=natom, &
1749 para_env=para_env, &
1750 local_particles=local_particles)
1751
1752 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
1753 dft_control%efield_fields(1)%efield%strength
1754
1755 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1756 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1757 CALL get_qs_env(qs_env=qs_env, force=force)
1758
1759 DO ikind = 1, nkind
1760 force(ikind)%efield = 0.0_dp
1761 DO ia = 1, local_particles%n_el(ikind)
1762 iatom = local_particles%list(ikind)%array(ia)
1763 q = charges(iatom) - qlag(iatom)
1764 atom_a = atom_of_kind(iatom)
1765 force(ikind)%efield(1:3, atom_a) = -q*fieldpol(1:3)
1766 END DO
1767 CALL para_env%sum(force(ikind)%efield)
1768 END DO
1769
1770 END SUBROUTINE eeq_efield_force_loc
1771
1772! **************************************************************************************************
1773!> \brief ...
1774!> \param qs_env ...
1775!> \param charges ...
1776!> \param qlag ...
1777! **************************************************************************************************
1778 SUBROUTINE eeq_efield_force_periodic(qs_env, charges, qlag)
1779 TYPE(qs_environment_type), POINTER :: qs_env
1780 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges, qlag
1781
1782 INTEGER :: atom_a, ia, iatom, ikind, natom, nkind
1783 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1784 LOGICAL :: dfield, use_virial
1785 REAL(kind=dp) :: q
1786 REAL(kind=dp), DIMENSION(3) :: fa, fieldpol, ria
1787 REAL(kind=dp), DIMENSION(3, 3) :: pve
1788 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1789 TYPE(cell_type), POINTER :: cell
1790 TYPE(dft_control_type), POINTER :: dft_control
1791 TYPE(distribution_1d_type), POINTER :: local_particles
1792 TYPE(mp_para_env_type), POINTER :: para_env
1793 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1794 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1795 TYPE(virial_type), POINTER :: virial
1796
1797 CALL get_qs_env(qs_env=qs_env, &
1798 dft_control=dft_control, &
1799 cell=cell, particle_set=particle_set, &
1800 virial=virial, &
1801 nkind=nkind, natom=natom, &
1802 para_env=para_env, &
1803 local_particles=local_particles)
1804
1805 dfield = dft_control%period_efield%displacement_field
1806 cpassert(.NOT. dfield)
1807
1808 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
1809 cpabort("use of strength_list not implemented for eeq_efield_force_periodic")
1810 END IF
1811
1812 fieldpol = dft_control%period_efield%polarisation
1813 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
1814 fieldpol = -fieldpol*dft_control%period_efield%strength
1815
1816 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1817
1818 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1819 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
1820 CALL get_qs_env(qs_env=qs_env, force=force)
1821
1822 pve = 0.0_dp
1823 DO ikind = 1, nkind
1824 force(ikind)%efield = 0.0_dp
1825 DO ia = 1, local_particles%n_el(ikind)
1826 iatom = local_particles%list(ikind)%array(ia)
1827 q = charges(iatom) - qlag(iatom)
1828 fa(1:3) = q*fieldpol(1:3)
1829 atom_a = atom_of_kind(iatom)
1830 force(ikind)%efield(1:3, atom_a) = fa
1831 IF (use_virial) THEN
1832 ria = particle_set(ia)%r
1833 ria = pbc(ria, cell)
1834 CALL virial_pair_force(pve, -0.5_dp, fa, ria)
1835 CALL virial_pair_force(pve, -0.5_dp, ria, fa)
1836 END IF
1837 END DO
1838 CALL para_env%sum(force(ikind)%efield)
1839 END DO
1840 virial%pv_virial = virial%pv_virial + pve
1841
1842 END SUBROUTINE eeq_efield_force_periodic
1843
1844END 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:3285
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:1515
eeq_method::eeq_efield_force_periodic
subroutine, public eeq_efield_force_periodic(qs_env, charges, qlag)
...
Definition eeq_method.F:1779
eeq_method::eeq_efield_pot
subroutine, public eeq_efield_pot(qs_env, efr)
...
Definition eeq_method.F:1601
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:1730
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, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:518
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:1629
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:1031
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:586
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:219
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