Octopus
ion_interaction.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2021 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
23 use atom_oct_m
24 use comm_oct_m
25 use debug_oct_m
26 use global_oct_m
33 use mpi_oct_m
36 use parser_oct_m
38 use ps_oct_m
40 use space_oct_m
43
44 implicit none
45
46 private
47 public :: &
55
57 real(real64) :: alpha
58 type(distributed_t) :: dist
59 end type ion_interaction_t
60
61 integer, parameter :: &
62 ION_COMPONENT_REAL = 1, &
66
67contains
68
69 subroutine ion_interaction_init(this, namespace, space, natoms)
70 type(ion_interaction_t), intent(out) :: this
71 type(namespace_t), intent(in) :: namespace
72 class(space_t), intent(in) :: space
73 integer, intent(in) :: natoms
74
75 push_sub(ion_interaction_init)
76
77 !%Variable EwaldAlpha
78 !%Type float
79 !%Default 0.21
80 !%Section Hamiltonian
81 !%Description
82 !% The value 'Alpha' that controls the splitting of the Coulomb
83 !% interaction in the Ewald sum used to calculation the ion-ion
84 !% interaction for periodic systems. This value affects the speed
85 !% of the calculation, normally users do not need to modify it.
86 !%End
87 call parse_variable(namespace, 'EwaldAlpha', 0.21_real64, this%alpha)
88
89 call distributed_nullify(this%dist, natoms)
90
91 if (space%periodic_dim == 1) then
92 call messages_write('For systems that are periodic in 1D, the interaction between', new_line = .true.)
93 call messages_write('ions is not implemented. This affects the calculation', new_line = .true.)
94 call messages_write('of total energy and forces, so both are zeroed.')
95 call messages_warning(namespace=namespace)
96 end if
97
99 end subroutine ion_interaction_init
100
101 subroutine ion_interaction_init_parallelization(this, natoms, mc)
102 type(ion_interaction_t), intent(inout) :: this
103 integer, intent(in) :: natoms
104 type(multicomm_t), intent(in) :: mc
105
107
108 !As the code below is not parallelized with any of k-point, states nor domain
109 !we can safely parallelize it over atoms
110 if (debug%info) then
111 call distributed_init(this%dist, natoms, mc%master_comm, "Ions")
112 else
113 call distributed_init(this%dist, natoms, mc%master_comm)
114 end if
118
119 subroutine ion_interaction_end(this)
120 type(ion_interaction_t), intent(inout) :: this
121
122 push_sub(ion_interaction_end)
123
124 this%alpha = -m_one
125
126 call distributed_end(this%dist)
127
128 pop_sub(ion_interaction_end)
129 end subroutine ion_interaction_end
130
135 subroutine ion_interaction_calculate(this, space, latt, atom, natoms, pos, lsize, energy, force, &
136 energy_components, force_components)
137 type(ion_interaction_t), intent(inout) :: this
138 class(space_t), intent(in) :: space
139 type(lattice_vectors_t), intent(in) :: latt
140 type(atom_t), intent(in) :: atom(:)
141 integer, intent(in) :: natoms
142 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
143 real(real64), intent(in) :: lsize(:)
144 real(real64), intent(out) :: energy
145 real(real64), intent(out) :: force(:, :)
146 real(real64), optional, intent(out) :: energy_components(:)
147 real(real64), optional, intent(out) :: force_components(:, :, :)
148
151 call profiling_in("ION_ION_INTERACTION")
152
153 if (present(energy_components)) then
154 assert(ubound(energy_components, dim = 1) == ion_num_components)
155 energy_components = m_zero
156 end if
157
158 if (present(force_components)) then
159 assert(all(ubound(force_components) == (/space%dim, natoms, ion_num_components/)))
160 force_components = m_zero
161 end if
163 if (space%is_periodic() .and. any_species_is_jellium_sphere(atom)) then
164 call messages_not_implemented('No periodic implementation of ion-ion energy for the jellium sphere')
165 end if
166
167 if (space%is_periodic()) then
168 if (all_species_are_jellium_slab(atom)) then
169 energy = jellium_slab_energy_periodic(space, atom, lsize)
170 force = 0._real64
171 else
172 call ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, energy_components, force_components)
173 end if
174 else
175 call ion_interaction_finite(this%dist, space, atom, pos, lsize, energy, force)
176 energy = energy + jellium_self_energy_finite(this%dist, latt, atom, lsize)
177 end if
178
179 call profiling_out("ION_ION_INTERACTION")
181
182 end subroutine ion_interaction_calculate
183
189 function jellium_slab_energy_periodic(space, atom, lsize) result(energy)
190 class(space_t), intent(in) :: space
191 type(atom_t), intent(in) :: atom(:)
192 real(real64), intent(in) :: lsize(:)
193 real(real64) :: energy
195 real(real64) :: area
196
197 ! Implementation assumes a single atom
198 assert(size(atom) == 1)
199 ! This is only allowed if periodic dim = 2. In that case the lattice volume is in fact an area.
200 assert(space%periodic_dim == 2)
201
202 select type(spec => atom(1)%species)
203 type is (jellium_slab_t)
204 area = lsize(1) * lsize(2) * m_four
205 energy = m_pi * spec%get_density(lsize) **2 * area * spec%thickness()**3 / m_three
206 class default
207 assert(.false.)
208 end select
209
211
225 function jellium_self_energy_finite(dist, latt, atom, lsize) result(energy)
226 type(distributed_t), intent(in) :: dist
227 type(lattice_vectors_t), intent(in) :: latt
228 type(atom_t), intent(in) :: atom(:)
229 real(real64), intent(in) :: lsize(:)
230 real(real64) :: energy
231
232 real(real64) :: zi
233 integer :: iatom
234 logical :: lattice_is_orthogonal
235 class(species_t), pointer :: spec
236
238
239 energy = 0._real64
240 lattice_is_orthogonal = .not. latt%nonorthogonal
241
242 do iatom = dist%start, dist%end
243 spec => atom(iatom)%species
244 zi = spec%get_zval()
245
246 select type(spec)
247 type is (jellium_sphere_t)
248 energy = energy + (m_three / m_five) * zi**2 / spec%radius()
249 ! The part depending on the simulation sphere is neglected
250
251 type is (jellium_slab_t)
252 ! Jellium slab energy only implemented for orthogonal cells.
253 ! One would need to replace (lsize(1) * lsize(2)) * spec%thickness()) with the triple product
254 assert(lattice_is_orthogonal)
255 energy = energy + m_pi * zi**2 / (m_four * lsize(1)*lsize(2)) * spec%thickness() / m_three
256 ! The part depending on the simulation box transverse dimension is neglected
257 end select
258 nullify(spec)
259 enddo
260
261 call comm_allreduce(dist%mpi_grp, energy)
262
264
265 end function jellium_self_energy_finite
266
268 subroutine ion_interaction_finite(dist, space, atom, pos, lsize, energy, force)
269 type(distributed_t), intent(in) :: dist
270 class(space_t), intent(in) :: space
271 type(atom_t), intent(in) :: atom(:)
272 real(real64), intent(in) :: pos(:,:)
273 real(real64), intent(in) :: lsize(:)
274 real(real64), intent(out) :: energy
275 real(real64), intent(out) :: force(:, :)
276
277 class(species_t), pointer :: species_i, species_j
278 real(real64) :: r(space%dim), f(space%dim)
279 real(real64) :: r_mag
280 real(real64) :: u_e
281 real(real64) :: zi, zj
282 integer :: iatom, jatom, natoms
283
284 push_sub(ion_interaction_finite)
285
286 natoms = size(atom)
287 energy = m_zero
288 force(1:space%dim, 1:natoms) = m_zero
289
290 do iatom = dist%start, dist%end
291 species_i => atom(iatom)%species
292 zi = species_i%get_zval()
293
294 do jatom = iatom + 1, natoms
295 species_j => atom(jatom)%species
296 zj = species_j%get_zval()
297
298 r = pos(:, iatom) - pos(:, jatom)
299 r_mag = norm2(r)
300 u_e = zi * zj / r_mag
301
302 energy = energy + u_e
303 f(1:space%dim) = (u_e / r_mag**2) * r(1:space%dim)
304 force(1:space%dim, iatom) = force(1:space%dim, iatom) + f(1:space%dim)
305 force(1:space%dim, jatom) = force(1:space%dim, jatom) - f(1:space%dim)
306 end do
307 end do
308
309 call comm_allreduce(dist%mpi_grp, energy)
310 call comm_allreduce(dist%mpi_grp, force)
311
312 nullify(species_i, species_j)
313
315
316 end subroutine ion_interaction_finite
317
319 subroutine ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, &
320 energy_components, force_components)
321 type(ion_interaction_t), intent(in) :: this
322 class(space_t), intent(in) :: space
323 type(lattice_vectors_t), intent(in) :: latt
324 type(atom_t), intent(in) :: atom(:)
325 integer, intent(in) :: natoms
326 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
327 real(real64), intent(out) :: energy
328 real(real64), intent(out) :: force(:, :)
329 real(real64), optional, intent(out) :: energy_components(:)
330 real(real64), optional, intent(out) :: force_components(:, :, :)
331
332 real(real64) :: ereal, efourier, epseudo, eself
333 real(real64) :: charge
334
336
337 energy = m_zero
338 force(1:space%dim, 1:natoms) = m_zero
339
340 call ewald_short(this%dist, space, latt, atom, pos, this%alpha, ereal, force)
341 if (present(force_components)) then
342 force_components(1:space%dim, 1:natoms, ion_component_real) = force(1:space%dim, 1:natoms)
343 end if
344
345 call ewald_self_interaction(this%dist, atom, this%alpha, eself, charge)
346
347 call profiling_in("EWALD_LONG")
348 select case (space%periodic_dim)
349 case (1)
350 ! Warning added in init routine, such that it is not displayed per SCF step
351 efourier = m_zero
352 ! Do not confuse the user and set to zero all the other components
353 ereal = m_zero
354 eself = m_zero
355 force = m_zero
356 epseudo = m_zero
357 case (2)
358 ! The energy contribution of the long range part of the pseudo is
359 ! not correctly accounted for in systems periodic in 1D or 2D, however
360 ! this term should not appear here anyway. See Issue #950.
361 epseudo = m_zero
362 call ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
363 case (3)
364 call ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
365 !TODO(Alex/Nicolas) Issue #950. Refactor: Move G=0 correction from ion-ion energy to pseudopotential energy
366 call pseudopotential_correction_3d(this%dist, latt, atom, charge, epseudo)
367 end select
368 call profiling_out("EWALD_LONG")
369
370 if (present(energy_components)) then
371 energy_components(ion_component_real) = ereal
372 energy_components(ion_component_self) = eself
373 energy_components(ion_component_fourier) = efourier
374 end if
375
376 if (present(force_components)) then
377 ! This is dependent on the order in which the force terms are computed
378 force_components(1:space%dim, 1:natoms, ion_component_fourier) = &
379 force(1:space%dim, 1:natoms) - force_components(1:space%dim, 1:natoms, ion_component_real)
380 end if
381
382 energy = ereal + efourier + eself + epseudo
383
385 end subroutine ion_interaction_periodic
386
407 subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
408 type(distributed_t), intent(in) :: dist
409 class(space_t), intent(in) :: space
410 type(lattice_vectors_t), intent(in) :: latt
411 type(atom_t), intent(in) :: atom(:)
412 real(real64), intent(in) :: pos(:, :)
413
414 real(real64), intent(in) :: alpha
415 real(real64), intent(out) :: ereal
416 real(real64), intent(inout) :: force(:, :)
417 ! Intent(inout) allows force contributions to be summed
418 integer :: iatom, jatom, icopy, natoms
419 real(real64) :: rnorm, xi(space%dim)
420 real(real64) :: force_real(space%dim)
421 real(real64) :: zi, zj
422 real(real64) :: erfc
423 real(real64) :: rcut
424 type(lattice_iterator_t) :: latt_iter
425
426 push_sub_with_profile(ewald_short)
427
428 ereal = m_zero
429 ! Origin of this value is not documented
430 rcut = 6.0_real64 / alpha
431 latt_iter = lattice_iterator_t(latt, rcut)
432 natoms = size(atom)
433
434 !$omp parallel default(shared) private(iatom, jatom, zi, zj, icopy, xi, rnorm, erfc, force_real) reduction(+:ereal, force)
435 do iatom = dist%start, dist%end
436 if (.not. atom(iatom)%species%represents_real_atom()) cycle
437 zi = atom(iatom)%species%get_zval()
438
439 ! Diagonal terms iatom == jatom for all cells, except T=(0,0,0)
440 !$omp do
441 do icopy = 1, latt_iter%n_cells
442 rnorm = norm2(latt_iter%get(icopy))
443 if (rnorm < r_min_atom_dist) cycle
444 if (rnorm > rcut) cycle
445 erfc = m_one - loct_erf(alpha * rnorm)
446 ereal = ereal + m_half * zi * zi * erfc /rnorm
447 enddo
448 !$omp end do
449
450 ! Upper triangle, for all replica cells
451 do jatom = iatom + 1, natoms
452 zj = atom(jatom)%species%get_zval()
453
454 !$omp do
455 do icopy = 1, latt_iter%n_cells
456 xi = pos(:, iatom) + latt_iter%get(icopy)
457 rnorm = norm2(xi - pos(:, jatom))
458 if (rnorm > rcut) cycle
459
460 erfc = m_one - loct_erf(alpha * rnorm)
461 force_real(:) = zj * zi * (xi - pos(:, jatom)) * &
462 (erfc / rnorm + m_two * alpha / sqrt(m_pi) *exp(-(alpha*rnorm)**2)) / rnorm**2
463
464 ! Factor 1/2 omitted as one is only summing over upper triangle
465 ereal = ereal + zj*zi*erfc/rnorm
466
467 ! Upper trianglar contribution
468 force(1:space%dim, jatom) = force(1:space%dim, jatom) - force_real
469
470 ! Lower triangular contribution
471 force(1:space%dim, iatom) = force(1:space%dim, iatom) + force_real
472 end do
473 !$omp end do
474
475 end do
476 end do
477 !$omp end parallel
478
479 call comm_allreduce(dist%mpi_grp, ereal)
480 call comm_allreduce(dist%mpi_grp, force)
481
482 pop_sub_with_profile(ewald_short)
483 end subroutine ewald_short
484
489 subroutine ewald_self_interaction(dist, atom, alpha, eself, charge)
490 type(distributed_t), intent(in) :: dist
491 type(atom_t), intent(in) :: atom(:)
492 real(real64), intent(in) :: alpha
493 real(real64), intent(out) :: eself
494 real(real64), intent(out) :: charge
495
496 integer :: iatom
497 real(real64) :: zi
498
499 push_sub(ewald_self_interaction)
501 eself = m_zero
502 charge = m_zero
503
504 do iatom = dist%start, dist%end
505 zi = atom(iatom)%species%get_zval()
506 charge = charge + zi
507 eself = eself - alpha / sqrt(m_pi) * zi**2
508 end do
509
510 call comm_allreduce(dist%mpi_grp, eself)
511 call comm_allreduce(dist%mpi_grp, charge)
512
514 end subroutine ewald_self_interaction
515
517 subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
518 type(ion_interaction_t), intent(in) :: this
519 class(space_t), intent(in) :: space
520 type(lattice_vectors_t), intent(in) :: latt
521 type(atom_t), intent(in) :: atom(:)
522 integer, intent(in) :: natoms
523 real(real64), intent(in) :: pos(:,:)
524 real(real64), intent(inout) :: efourier
525 real(real64), intent(inout) :: force(:, :)
526 real(real64), intent(in) :: charge
527
528 real(real64) :: rcut, gmax_squared
529 integer :: iatom
530 integer :: ix, iy, iz, isph
531 real(real64) :: gvec(3), gred(3), gg2, gx
532 real(real64) :: factor
533 complex(real64) :: sumatoms, tmp(3), aa
534
535 complex(real64), allocatable :: phase(:)
536
537 push_sub(ewald_long_3d)
538
539 assert(space%dim == 3)
540 assert(space%periodic_dim == 3)
541
542 ! And the long-range part, using an Ewald sum
543 safe_allocate(phase(1:natoms))
544
545 ! get a converged value for the cutoff in g
546 rcut = sqrt(minval(sum(latt%klattice**2, dim=1)))
547
548 ! 9.5 is a constant that controls the range separation
549 isph = ceiling(9.5_real64*this%alpha/rcut)
550
551 ! First the G = 0 term (charge was calculated previously)
552 efourier = -m_pi*charge**2/(m_two*this%alpha**2*latt%rcell_volume)
553
554 ! Getting the G-shell cutoff
555 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
556
557 do ix = -isph, isph
558 do iy = -isph, isph
559 do iz = -isph, isph
560
561 gred = [ix, iy, iz]
562 call kpoints_to_absolute(latt, gred, gvec)
563 gg2 = dot_product(gvec, gvec)
564
565 ! g=0 must be removed from the sum
566 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
567
568 gx = -0.25_real64*gg2/this%alpha**2
569
570 if (gx < -36.0_real64) cycle
571
572 factor = m_two*m_pi/latt%rcell_volume*exp(gx)/gg2
573
574 if (factor < epsilon(factor)) cycle
575
576 sumatoms = m_z0
577 !$omp parallel do private(iatom, gx, aa) reduction(+:sumatoms)
578 do iatom = 1, natoms
579 gx = sum(gvec*pos(:,iatom))
580 aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
581 phase(iatom) = aa
582 sumatoms = sumatoms + aa
583 end do
584
585 efourier = efourier + real(factor*sumatoms*conjg(sumatoms), real64)
586
587 do iatom = 1, natoms
588 tmp = m_zi*gvec*phase(iatom)
589 force(1:space%dim, iatom) = force(1:space%dim, iatom) - factor*real(conjg(tmp)*sumatoms + tmp*conjg(sumatoms), real64)
590
591 end do
592
593 end do
594 end do
595 end do
596
597 safe_deallocate_a(phase)
598
599 pop_sub(ewald_long_3d)
600
601 end subroutine ewald_long_3d
602
604 subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
605 type(ion_interaction_t), intent(in) :: this
606 class(space_t), intent(in) :: space
607 type(lattice_vectors_t), intent(in) :: latt
608 type(atom_t), intent(in) :: atom(:)
609 integer, intent(in) :: natoms
610 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
611 real(real64), intent(inout) :: efourier
612 real(real64), intent(inout) :: force(:, :)
613
614 real(real64) :: rcut, gmax_squared
615 integer :: iatom, jatom
616 integer :: ix, iy, ix_max, iy_max
617 real(real64) :: gvec(space%dim), gg2, gx, gg_abs
618 real(real64) :: factor,factor1,factor2, coeff
619 real(real64) :: dz_max, dz_ij, erfc1, erfc2, tmp_erf
620 real(real64), allocatable :: force_tmp(:,:)
621 real(real64), parameter :: tol = 1e-10_real64
622
623 push_sub(ewald_long_2d)
624
625 assert(space%periodic_dim == 2)
626 assert(space%dim == 2 .or. space%dim == 3)
627
628 ! And the long-range part, using an Ewald sum
629
630 ! Searching maximum distance
631 if (space%dim == 3) then
632 dz_max = m_zero
633 do iatom = 1, natoms
634 do jatom = iatom + 1, natoms
635 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
636 end do
637 end do
638 else
639 ! For a 2D system, all atoms are on the plane, so the distance is zero
640 dz_max = m_zero
641 end if
642
643 !get a converged value for the cutoff in g
644 rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
645 if (dz_max > tol) then
646 do
647 if (rcut * dz_max >= m_max_exp_arg) exit !Maximum double precision number
648 erfc1 = m_one - loct_erf(this%alpha*dz_max + m_half*rcut/this%alpha)
649 if (erfc1 * exp(rcut*dz_max) < 1.e-10_real64) exit
650 rcut = rcut * 1.414_real64
651 end do
652 end if
653
654 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
655 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
656
657 safe_allocate(force_tmp(1:space%dim, 1:natoms))
658 force_tmp = m_zero
659
660 ! First the G = 0 term
661 efourier = m_zero
662 factor = m_pi/latt%rcell_volume
663 !$omp parallel do private(jatom, dz_ij, tmp_erf, factor1, factor2) reduction(+:efourier,force_tmp) &
664 !$omp& collapse(2)
665 do iatom = this%dist%start, this%dist%end
666 do jatom = 1, natoms
667 ! efourier
668 if (space%dim == 3) then
669 dz_ij = pos(3, iatom) - pos(3, jatom)
670 else
671 dz_ij = m_zero
672 end if
673
674 tmp_erf = loct_erf(this%alpha*dz_ij)
675 factor1 = dz_ij*tmp_erf
676 factor2 = exp(-(this%alpha*dz_ij)**2)/(this%alpha*sqrt(m_pi))
677
678 efourier = efourier - factor &
679 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
680
681 ! force
682 if (iatom == jatom)cycle
683 if (abs(tmp_erf) < m_epsilon) cycle
684
685 if (space%dim == 3) then
686 force_tmp(3, iatom) = force_tmp(3, iatom) - (- m_two*factor) &
687 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * tmp_erf
688 end if
689
690 end do
691 end do
692
693 ! Getting the G-shell cutoff
694 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
695 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
696
697 !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, gx, dz_ij, erfc1, factor1, erfc2, factor2, coeff) &
698 !$omp& collapse(2) reduction(+:efourier, force_tmp)
699 do ix = -ix_max, ix_max
700 do iy = -iy_max, iy_max
701
702 gvec = ix*latt%klattice(:, 1) + iy*latt%klattice(:, 2)
703 gg2 = sum(gvec**2)
704
705 ! g=0 must be removed from the sum
706 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
707 gg_abs = sqrt(gg2)
708 factor = m_half*m_pi/(latt%rcell_volume*gg_abs)
709
710 do iatom = this%dist%start, this%dist%end
711 do jatom = iatom, natoms
712 ! efourier
713 gx = sum(gvec(1:2) * (pos(1:2, iatom) - pos(1:2, jatom)))
714 gx = gvec(1)*(pos(1, iatom) - pos(1, jatom)) + gvec(2)*(pos(2, iatom) - pos(2, jatom))
715 if (space%dim == 3) then
716 dz_ij = pos(3, iatom) - pos(3, jatom)
717 else
718 dz_ij = m_zero
719 end if
720
721 erfc1 = m_one - loct_erf(this%alpha*dz_ij + m_half*gg_abs/this%alpha)
722 if (abs(erfc1) > m_epsilon) then
723 factor1 = exp(gg_abs*dz_ij)*erfc1
724 else
725 factor1 = m_zero
726 end if
727 erfc2 = m_one - loct_erf(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)
728 if (abs(erfc2) > m_epsilon) then
729 factor2 = exp(-gg_abs*dz_ij)*erfc2
730 else
731 factor2 = m_zero
732 end if
733
734 if (iatom == jatom) then
735 coeff = m_one
736 else
737 coeff = m_two
738 end if
739
740 efourier = efourier &
741 + factor * coeff &
742 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
743 * cos(gx)* ( factor1 + factor2)
744
745 ! force
746 if (iatom == jatom) cycle
747
748 force_tmp(1:2, iatom) = force_tmp(1:2, iatom) &
749 + m_two * factor * gvec(1:2) &
750 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
751 *sin(gx)*(factor1 + factor2)
752
753 force_tmp(1:2, jatom) = force_tmp(1:2, jatom) &
754 - m_two * factor * gvec(1:2) &
755 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
756 *sin(gx)*(factor1 + factor2)
757
758 factor1 = gg_abs*erfc1 &
759 - m_two*this%alpha/sqrt(m_pi)*exp(-(this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
760 if (abs(factor1) > m_epsilon) then
761 factor1 = factor1*exp(gg_abs*dz_ij)
762 else
763 factor1 = m_zero
764 end if
765
766 factor2 = gg_abs*erfc2 &
767 - m_two*this%alpha/sqrt(m_pi)*exp(-(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
768 if (abs(factor2) > m_epsilon) then
769 factor2 = factor2*exp(-gg_abs*dz_ij)
770 else
771 factor2 = m_zero
772 end if
773
774 if (space%dim == 3) then
775 force_tmp(3, iatom) = force_tmp(3, iatom) &
776 - m_two*factor &
777 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
778 * cos(gx)* ( factor1 - factor2)
779 force_tmp(3, jatom) = force_tmp(3, jatom) &
780 + m_two*factor &
781 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
782 * cos(gx)* ( factor1 - factor2)
783 end if
784
785 end do
786 end do
787
788
789 end do
790 end do
791
792 call comm_allreduce(this%dist%mpi_grp, efourier)
793 call comm_allreduce(this%dist%mpi_grp, force_tmp)
794
795 force = force + force_tmp
796
797 safe_deallocate_a(force_tmp)
798
799 pop_sub(ewald_long_2d)
800 end subroutine ewald_long_2d
801
802 !TODO(Alex/Nicolas) Issue #950. Refactor: Move G=0 correction from ion-ion energy to pseudopotential energy
809 subroutine pseudopotential_correction_3d(dist, latt, atom, charge, epseudo)
810 type(distributed_t), intent(in) :: dist
811 type(lattice_vectors_t), intent(in) :: latt
812 type(atom_t), intent(in) :: atom(:)
813 real(real64), intent(out) :: epseudo
814
815 real(real64) :: zi
816 real(real64) :: charge
817 integer :: iatom
818
820
821 epseudo = m_zero
822 do iatom = dist%start, dist%end
823 select type(spec => atom(iatom)%species)
824 class is(pseudopotential_t)
825 zi = spec%get_zval()
826 epseudo = epseudo + m_pi *zi * &
827 (spec%ps%sigma_erf * sqrt(m_two))**2 / latt%rcell_volume * charge
828 end select
829 end do
830 call comm_allreduce(dist%mpi_grp, epseudo)
831
833
834 end subroutine pseudopotential_correction_3d
835
837 subroutine ion_interaction_stress(this, space, latt, atom, natoms, pos, stress_ii)
838 type(ion_interaction_t), intent(inout) :: this
839 class(space_t), intent(in) :: space
840 type(lattice_vectors_t), intent(in) :: latt
841 type(atom_t), intent(in) :: atom(:)
842 integer, intent(in) :: natoms
843 real(real64), intent(in) :: pos(:,:)
844 real(real64), intent(out) :: stress_ii(space%dim, space%dim)
845
846 real(real64) :: stress_short(1:space%dim, 1:space%dim), stress_Ewald(1:space%dim, 1:space%dim)
847
848 push_sub(ion_interaction_stress)
849
850 stress_ii = m_zero
851
852 ! Only implemented in the periodic case
853 assert(space%is_periodic())
854
855 ! Short range part in real space
856 call ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
857
858 ! Long range part in Fourier space
859 select case(space%periodic_dim)
860 case(3)
861 call ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
862 case(2)
863 call ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
864 case default
865 assert(.false.)
866 end select
867
868 stress_ii = stress_short + stress_ewald
869
871 end subroutine ion_interaction_stress
872
873 ! ---------------------------------------------------------
889
890 subroutine ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
891 type(ion_interaction_t), intent(inout) :: this
892 class(space_t), intent(in) :: space
893 type(lattice_vectors_t), intent(in) :: latt
894 type(atom_t), intent(in) :: atom(:)
895 integer, intent(in) :: natoms
896 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
897 real(real64), intent(out) :: stress_short(1:space%dim, 1:space%dim)
898
899 real(real64) :: xi(space%dim)
900 real(real64) :: r_ij, zi, zj, erfc, Hp, factor
901 integer :: iatom, jatom, icopy, idir, jdir
902 real(real64) :: alpha, rcut
903 type(lattice_iterator_t) :: latt_iter
904
906 call profiling_in("ION_ION_STRESS_SHORT")
907
908 ! Only implemented in the periodic case
909 assert(space%is_periodic())
910
911 alpha = this%alpha
912
913 ! See the code for the energy above to understand this parameter
914 rcut = 6.0_real64/alpha
915
916 ! the short-range part is calculated directly
917 stress_short = m_zero
918 latt_iter = lattice_iterator_t(latt, rcut)
919
920 do iatom = this%dist%start, this%dist%end
921 select type(spec => atom(iatom)%species)
922 class is(jellium_t)
923 cycle
924 end select
925 zi = atom(iatom)%species%get_zval()
926
927 do icopy = 1, latt_iter%n_cells
928 xi = pos(:, iatom) + latt_iter%get(icopy)
929
930 do jatom = 1, natoms
931 zj = atom(jatom)%species%get_zval()
932 r_ij = norm2(xi - pos(:, jatom))
933
934 if (r_ij < r_min_atom_dist) cycle
935
936 erfc = m_one - loct_erf(alpha*r_ij)
937 hp = -m_two/sqrt(m_pi)*exp(-(alpha*r_ij)**2) - erfc/(alpha*r_ij)
938 factor = m_half*zj*zi*alpha*hp
939 do idir = 1, space%periodic_dim
940 do jdir = 1, space%periodic_dim
941 stress_short(idir, jdir) = stress_short(idir, jdir) &
942 - factor*(xi(idir) - pos(idir, jatom))*(xi(jdir) - pos(jdir, jatom))/(r_ij**2)
943 end do
944 end do
945
946 end do
947 end do
948 end do
949
950 if (this%dist%parallel) then
951 call comm_allreduce(this%dist%mpi_grp, stress_short)
952 end if
953
954 stress_short = stress_short/latt%rcell_volume
955
956 call profiling_out("ION_ION_STRESS_SHORT")
957
959 end subroutine ion_interaction_stress_short
960
961
962
963 ! ---------------------------------------------------------
975 subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
976 type(ion_interaction_t), intent(inout) :: this
977 class(space_t), intent(in) :: space
978 type(lattice_vectors_t), intent(in) :: latt
979 type(atom_t), intent(in) :: atom(:)
980 integer, intent(in) :: natoms
981 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
982 real(real64), intent(out) :: stress_ewald(3, 3)
984 real(real64) :: zi, rcut, gmax_squared
985 integer :: iatom
986 integer :: ix, iy, iz, isph, idim, idir, jdir
987 real(real64) :: gred(3), gvec(3), gg2, gx
988 real(real64) :: factor, charge, charge_sq
989 complex(real64) :: sumatoms, aa
990
991 call profiling_in("STRESS_3D_EWALD")
992 push_sub(ewald_3d_stress)
993
994 ! Currently this is only implemented for 3D
995 assert(space%dim == 3)
996 assert(space%periodic_dim == 3) ! Not working for mixed periodicity
997 ! (klattice along the non-periodic directions is wrong)
998 ! Anyway gg/gg2 is not correct for mixed periodicity
999
1000 stress_ewald = m_zero
1001
1002 ! And the long-range part, using an Ewald sum
1003 charge = m_zero
1004 charge_sq = m_zero
1005 do iatom = 1, natoms
1006 zi = atom(iatom)%species%get_zval()
1007 charge = charge + zi
1008 charge_sq = charge_sq + zi**2
1009 end do
1010
1011 ! get a converged value for the cutoff in g
1012 rcut = huge(rcut)
1013 do idim = 1, space%periodic_dim
1014 rcut = min(rcut, sum(latt%klattice(1:space%periodic_dim, idim)**2))
1015 end do
1016
1017 rcut = sqrt(rcut)
1018
1019 isph = ceiling(9.5_real64*this%alpha/rcut)
1020
1021 ! Getting the G-shell cutoff
1022 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
1023
1024 do ix = -isph, isph
1025 do iy = -isph, isph
1026 do iz = -isph, isph
1027
1028 gred = [ix, iy, iz]
1029 call kpoints_to_absolute(latt, gred, gvec)
1030 gg2 = sum(gvec**2)
1031
1032 ! g=0 must be removed from the sum
1033 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1034
1035 gx = -0.25_real64*gg2/this%alpha**2
1036
1037 if (gx < -36.0_real64) cycle
1038
1039 factor = m_two*m_pi*exp(gx)/(latt%rcell_volume*gg2)
1040
1041 if (factor < epsilon(factor)) cycle
1042
1043 sumatoms = m_z0
1044
1045 do iatom = 1, natoms
1046 gx = sum(gvec*pos(:, iatom))
1047 aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
1048 sumatoms = sumatoms + aa
1049 end do
1050
1051 factor = factor*abs(sumatoms)**2
1052
1053 do idir = 1, 3
1054 do jdir = 1, 3
1055 stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
1056 - m_two*factor*gvec(idir)*gvec(jdir)/gg2*(0.25_real64*gg2/this%alpha**2+m_one)
1057 end do
1058 stress_ewald(idir, idir) = stress_ewald(idir, idir) + factor
1059 end do
1060
1061 end do
1062 end do
1063 end do
1064
1065
1066 ! The G = 0 term of the Ewald summation
1067 factor = m_half*m_pi*charge**2/(latt%rcell_volume*this%alpha**2)
1068 do idir = 1,3
1069 stress_ewald(idir,idir) = stress_ewald(idir,idir) - factor
1070 end do
1071
1072 stress_ewald = stress_ewald / latt%rcell_volume
1073
1074
1075 call profiling_out("STRESS_3D_EWALD")
1076 pop_sub(ewald_3d_stress)
1077
1078 end subroutine ewald_3d_stress
1079
1080 ! ---------------------------------------------------------
1093 subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1094 type(ion_interaction_t), intent(inout) :: this
1095 type(space_t), intent(in) :: space
1096 type(lattice_vectors_t), intent(in) :: latt
1097 type(atom_t), intent(in) :: atom(:)
1098 integer, intent(in) :: natoms
1099 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
1100 real(real64), intent(out) :: stress_Ewald(3, 3)
1101
1102 real(real64) :: rcut, efourier
1103 integer :: iatom, jatom, idir, jdir
1104 integer :: ix, iy, ix_max, iy_max
1105 real(real64) :: gvec(3), gred(3), gg2, cos_gx, gg_abs, gmax_squared
1106 real(real64) :: factor,factor1,factor2, coeff, e_ewald
1107 real(real64) :: dz_max, z_ij, erfc1, erfc2, diff(3)
1108 real(real64), parameter :: tol = 1e-10_real64
1109
1110 push_sub(ewald_2d_stress)
1111
1112 assert(space%periodic_dim == 2)
1113 assert(space%dim == 3)
1114
1115 stress_ewald = m_zero
1116
1117 ! Searching maximum distance
1118 dz_max = m_zero
1119 do iatom = 1, natoms
1120 do jatom = iatom + 1, natoms
1121 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
1122 end do
1123 end do
1124
1125 !get a converged value for the cutoff in g
1126 ! Note: to understand these numbers, one needs to look into the energy routine for Ewald 2D
1127 rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
1128 if (dz_max > tol) then ! Else the code here does not work properly
1129 do
1130 if (rcut * dz_max >= m_max_exp_arg) exit !Maximum double precision number
1131 erfc1 = m_one - loct_erf(this%alpha*dz_max + m_half*rcut/this%alpha)
1132 if (erfc1 * exp(rcut*dz_max) < tol) exit
1133 rcut = rcut * 1.414_real64
1134 end do
1135 end if
1136
1137 ! First the G = 0 term
1138 efourier = m_zero
1139 factor = m_pi/latt%rcell_volume
1140 !$omp parallel do private(jatom, z_ij, factor1, factor2) reduction(+:efourier) collapse(2)
1141 do iatom = 1, natoms
1142 do jatom = 1, natoms
1143 z_ij = pos(3, iatom) - pos(3, jatom)
1144
1145 factor1 = z_ij * loct_erf(this%alpha*z_ij)
1146 factor2 = exp(-(this%alpha*z_ij)**2)/(this%alpha*sqrt(m_pi))
1147
1148 efourier = efourier - factor &
1149 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
1150 end do
1151 end do
1152
1153 ! Adding the G=0 term
1154 do idir = 1, 2
1155 stress_ewald(idir, idir) = efourier
1156 end do
1157
1158 ! Getting the G-shell cutoff
1159 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
1160 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
1161 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
1162 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
1163
1164 !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, diff, cos_gx, z_ij, idir, jdir, erfc1, factor1) &
1165 !$omp& private(erfc2, factor2, coeff, e_ewald) &
1166 !$omp& collapse(2) reduction(+:stress_Ewald)
1167 do ix = -ix_max, ix_max
1168 do iy = -iy_max, iy_max
1169
1170 gred = [ix, iy, 0]
1171 call kpoints_to_absolute(latt, gred, gvec)
1172 gg2 = dot_product(gvec,gvec)
1173
1174 ! g=0 must be removed from the sum
1175 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1176
1177 gg_abs = sqrt(gg2)
1178 factor = m_fourth*m_pi/(latt%rcell_volume*this%alpha*gg2)
1179
1180 do iatom = 1, natoms
1181 do jatom = iatom, natoms
1182 diff = pos(:, iatom) - pos(:, jatom)
1183 cos_gx = cos(sum(gvec(1:2) * diff(1:2)))
1184 z_ij = diff(3)
1185
1186 factor1 = screening_function_2d(this%alpha, z_ij, gg_abs, erfc1)
1187 factor2 = screening_function_2d(this%alpha,-z_ij, gg_abs, erfc2)
1188
1189 if (iatom == jatom) then
1190 coeff = m_one
1191 else
1192 coeff = m_two
1193 end if
1194
1195 do idir = 1, 2
1196 do jdir = 1, 2
1197 stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
1198 - factor*gvec(idir)*gvec(jdir) * cos_gx * (factor1 + factor2) * coeff&
1199 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval()
1200 end do
1201 end do
1202
1203 if (abs(erfc1) > m_epsilon) then
1204 factor1 = exp(-gg_abs*z_ij)*erfc1
1205 else
1206 factor1 = m_zero
1207 end if
1208 if (abs(erfc2) > m_epsilon) then
1209 factor2 = exp(gg_abs*z_ij)*erfc2
1210 else
1211 factor2 = m_zero
1212 end if
1213
1214 e_ewald = m_half * m_pi/latt%rcell_volume * coeff &
1215 * atom(iatom)%species%get_zval() * atom(jatom)%species%get_zval() &
1216 * cos_gx / gg_abs * (factor1 + factor2)
1217
1218 do idir = 1, 2
1219 stress_ewald(idir, idir) = stress_ewald(idir, idir) + e_ewald
1220 end do
1221
1222 end do !jatom
1223 end do !iatom
1224 end do !iy
1225 end do !ix
1226
1227 !call comm_allreduce(this%dist%mpi_grp, stress_Ewald)
1228
1229 stress_ewald = stress_ewald / latt%rcell_volume
1230
1231 pop_sub(ewald_2d_stress)
1232 end subroutine ewald_2d_stress
1233
1234 ! ---------------------------------------------------------
1236 real(real64) function screening_function_2d(alpha, z_ij, gg_abs, erfc) result(factor)
1237 real(real64), intent(in) :: alpha
1238 real(real64), intent(in) :: z_ij
1239 real(real64), intent(in) :: gg_abs
1240 real(real64), intent(out) :: erfc
1241
1242 real(real64) :: arg
1243
1244 arg = -alpha*z_ij + m_half*gg_abs/alpha
1245 erfc = m_one - loct_erf(arg)
1246 factor = m_two*alpha*(m_one/gg_abs + z_ij)*erfc - m_two/sqrt(m_pi)*exp(-arg**2)
1247 factor = factor*exp(-gg_abs*z_ij)
1248
1249 end function screening_function_2d
1250
1251 ! ---------------------------------------------------------
1252
1253 subroutine ion_interaction_test(space, latt, atom, natoms, pos, lsize, &
1254 namespace, mc)
1255 class(space_t), intent(in) :: space
1256 type(lattice_vectors_t), intent(in) :: latt
1257 type(atom_t), intent(in) :: atom(:)
1258 integer, intent(in) :: natoms
1259 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
1260 real(real64), intent(in) :: lsize(:)
1261 type(namespace_t), intent(in) :: namespace
1262 type(multicomm_t), intent(in) :: mc
1263
1264 type(ion_interaction_t) :: ion_interaction
1265 real(real64) :: energy
1266 real(real64), allocatable :: force(:, :), force_components(:, :, :)
1267 real(real64) :: energy_components(1:ION_NUM_COMPONENTS)
1268 integer :: iatom, idir
1269
1270 push_sub(ion_interaction_test)
1271
1272 call ion_interaction_init(ion_interaction, namespace, space, natoms)
1273 call ion_interaction_init_parallelization(ion_interaction, natoms, mc)
1274
1275 safe_allocate(force(1:space%dim, 1:natoms))
1276 safe_allocate(force_components(1:space%dim, 1:natoms, 1:ion_num_components))
1277
1278 call ion_interaction_calculate(ion_interaction, space, latt, atom, natoms, pos, lsize, energy, force, &
1279 energy_components = energy_components, force_components = force_components)
1280
1281 call messages_write('Ionic energy =')
1282 call messages_write(energy, fmt = '(f20.10)')
1283 call messages_info(namespace=namespace)
1284
1285 call messages_write('Real space energy =')
1286 call messages_write(energy_components(ion_component_real), fmt = '(f20.10)')
1287 call messages_info(namespace=namespace)
1288
1289 call messages_write('Self energy =')
1290 call messages_write(energy_components(ion_component_self), fmt = '(f20.10)')
1291 call messages_info(namespace=namespace)
1292
1293 call messages_write('Fourier energy =')
1294 call messages_write(energy_components(ion_component_fourier), fmt = '(f20.10)')
1295 call messages_info(namespace=namespace)
1296
1297 call messages_info(namespace=namespace)
1298
1299 do iatom = 1, natoms
1300 call messages_write('Ionic force atom')
1301 call messages_write(iatom)
1302 call messages_write(' =')
1303 do idir = 1, space%dim
1304 call messages_write(force(idir, iatom), fmt = '(f20.10)')
1305 end do
1306 call messages_info(namespace=namespace)
1307
1308 call messages_write('Real space force atom')
1309 call messages_write(iatom)
1310 call messages_write(' =')
1311 do idir = 1, space%dim
1312 call messages_write(force_components(idir, iatom, ion_component_real), fmt = '(f20.10)')
1313 end do
1314 call messages_info(namespace=namespace)
1315
1316 call messages_write('Fourier space force atom')
1317 call messages_write(iatom)
1318 call messages_write(' =')
1319 do idir = 1, space%dim
1320 call messages_write(force_components(idir, iatom, ion_component_fourier), fmt = '(f20.10)')
1321 end do
1322 call messages_info(namespace=namespace)
1323
1324 call messages_info(namespace=namespace)
1325 end do
1326
1327 safe_deallocate_a(force)
1328 safe_deallocate_a(force_components)
1330 call ion_interaction_end(ion_interaction)
1331
1332 pop_sub(ion_interaction_test)
1333 end subroutine ion_interaction_test
1334
1335end module ion_interaction_oct_m
1336
1337!! Local Variables:
1338!! mode: f90
1339!! coding: utf-8
1340!! End:
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
pure logical function, public all_species_are_jellium_slab(atom)
Check if all species are jellium slab.
Definition: atom.F90:292
pure logical function, public any_species_is_jellium_sphere(atom)
Check if any species is a jellium sphere.
Definition: atom.F90:309
type(debug_t), save, public debug
Definition: debug.F90:156
subroutine, public distributed_end(this)
subroutine, public distributed_nullify(this, total)
subroutine, public distributed_init(this, total, comm, tag, scalapack_compat)
Distribute N instances across M processes of communicator comm
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_max_exp_arg
Definition: global.F90:208
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_fourth
Definition: global.F90:197
complex(real64), parameter, public m_z0
Definition: global.F90:198
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public r_min_atom_dist
Minimal distance between two distinguishable atoms.
Definition: global.F90:183
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
real(real64), parameter, public m_five
Definition: global.F90:193
real(real64) function screening_function_2d(alpha, z_ij, gg_abs, erfc)
Auxiliary function for the Ewald 2D stress.
subroutine, public ion_interaction_stress(this, space, latt, atom, natoms, pos, stress_ii)
Computes the contribution to the stress tensor the ion-ion energy.
subroutine, public ion_interaction_init_parallelization(this, natoms, mc)
integer, parameter ion_component_self
real(real64) function jellium_slab_energy_periodic(space, atom, lsize)
Electrostatic energy of a periodic jellium slab.
subroutine, public ion_interaction_test(space, latt, atom, natoms, pos, lsize, namespace, mc)
subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
In-Chul Yeh and Max L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
subroutine ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
Computes the short-range contribution to the stress tensor the ion-ion energy.
subroutine ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, energy_components, force_components)
Total Ewald electrostatic energy and forces, for 1D, 2D and 3D systems.
real(real64) function jellium_self_energy_finite(dist, latt, atom, lsize)
Electrostatic self-interaction for jellium instances, with orthogonal cells.
subroutine, public ion_interaction_init(this, namespace, space, natoms)
subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
Short range component of the Ewald electrostatic energy and force.
subroutine pseudopotential_correction_3d(dist, latt, atom, charge, epseudo)
G=0 component of Ewald energy arising from the pseudopotentials, for 3D systems.
subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
Computes the long-range part of the 3D Ewald summation.
integer, parameter ion_component_real
integer, parameter ion_num_components
subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
Computes the contribution to the stress tensor from the 3D Ewald sum.
integer, parameter ion_component_fourier
subroutine ion_interaction_finite(dist, space, atom, pos, lsize, energy, force)
Electrostatic Ewald energy and forces for finite systems.
subroutine, public ion_interaction_end(this)
subroutine, public ion_interaction_calculate(this, space, latt, atom, natoms, pos, lsize, energy, force, energy_components, force_components)
Top level routine for computing electrostatic energies and forces between ions.
subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
Computes the contribution to the stress tensor from the 2D Ewald sum.
subroutine ewald_self_interaction(dist, atom, alpha, eself, charge)
@ brief Ewald self-interaction energy
subroutine, public kpoints_to_absolute(latt, kin, kout)
Definition: kpoints.F90:1031
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
Definition: ps.F90:114
Distribution of N instances over mpi_grpsize processes, for the local rank mpi_grprank....
The following class implements a lattice iterator. It allows one to loop over all cells that are with...
An abstract class for species. Derived classes include jellium, all electron, and pseudopotential spe...
Definition: species.F90:143
int true(void)