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
409 subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
410 type(distributed_t), intent(in) :: dist
411 class(space_t), intent(in) :: space
412 type(lattice_vectors_t), intent(in) :: latt
413 type(atom_t), intent(in) :: atom(:)
414 real(real64), intent(in) :: pos(:, :)
415
416 real(real64), intent(in) :: alpha
417 real(real64), intent(out) :: ereal
418 real(real64), intent(inout) :: force(:, :)
419 ! Intent(inout) allows force contributions to be summed
420 integer :: iatom, jatom, icopy, natoms
421 real(real64) :: rnorm, xi(space%dim)
422 real(real64) :: force_real(space%dim)
423 real(real64) :: zi, zj
424 real(real64) :: erfc
425 real(real64) :: rcut
426 type(lattice_iterator_t) :: latt_iter
427 real(real64) :: charge, coeff
428
429 push_sub_with_profile(ewald_short)
430
431 ereal = m_zero
432 ! Origin of this value is not documented
433 rcut = 6.0_real64 / alpha
434 latt_iter = lattice_iterator_t(latt, rcut)
435 natoms = size(atom)
436
437 charge = m_zero
438 do iatom = dist%start, dist%end
439 if (.not. atom(iatom)%species%represents_real_atom()) cycle
440 zi = atom(iatom)%species%get_zval()
441 charge = charge + zi**2
442 end do
443
444 ! Diagonal terms iatom == jatom for all cells, except T=(0,0,0)
445 ! Note: Only half of the copies are needed, by symmetries
446 do icopy = 1, latt_iter%n_cells
447 rnorm = norm2(latt_iter%get(icopy))
448 if (rnorm < r_min_atom_dist) cycle
449 if (rnorm > rcut) cycle
450 erfc = loct_erfc(alpha * rnorm)
451 ereal = ereal + m_half * charge * erfc /rnorm
452 end do
453
454 coeff = m_two * alpha / sqrt(m_pi)
455
456 !$omp parallel default(shared) &
457 !$omp& private(iatom, jatom, zi, zj, icopy, xi, rnorm, erfc, force_real, charge) &
458 !$omp& reduction(+: ereal, force)
459 do iatom = dist%start, dist%end
460 if (.not. atom(iatom)%species%represents_real_atom()) cycle
461 zi = atom(iatom)%species%get_zval()
462
463 ! Upper triangle, for all replica cells
464 do jatom = iatom + 1, natoms
465 zj = atom(jatom)%species%get_zval()
466
467 charge = zi*zj
468
469 !$omp do
470 do icopy = 1, latt_iter%n_cells
471 xi = pos(:, iatom) + latt_iter%get(icopy)
472 rnorm = norm2(xi - pos(:, jatom))
473 if (rnorm > rcut) cycle
474
475 erfc = loct_erfc(alpha * rnorm) / rnorm
476
477 ! Factor 1/2 omitted as one is only summing over upper triangle
478 ereal = ereal + charge * erfc
479
480 force_real(:) = charge * (xi - pos(:, jatom)) * &
481 (erfc + coeff *exp(-(alpha*rnorm)**2)) / rnorm**2
482
483 ! Upper trianglar contribution
484 force(1:space%dim, jatom) = force(1:space%dim, jatom) - force_real
485
486 ! Lower triangular contribution
487 force(1:space%dim, iatom) = force(1:space%dim, iatom) + force_real
488 end do
489 !$omp end do
490
491 end do
492 end do
493 !$omp end parallel
494
495 call comm_allreduce(dist%mpi_grp, ereal)
496 call comm_allreduce(dist%mpi_grp, force)
497
498 pop_sub_with_profile(ewald_short)
499 end subroutine ewald_short
500
505 subroutine ewald_self_interaction(dist, atom, alpha, eself, charge)
506 type(distributed_t), intent(in) :: dist
507 type(atom_t), intent(in) :: atom(:)
508 real(real64), intent(in) :: alpha
509 real(real64), intent(out) :: eself
510 real(real64), intent(out) :: charge
511
512 integer :: iatom
513 real(real64) :: zi
514
515 push_sub(ewald_self_interaction)
516
517 eself = m_zero
518 charge = m_zero
519
520 do iatom = dist%start, dist%end
521 zi = atom(iatom)%species%get_zval()
522 charge = charge + zi
523 eself = eself - alpha / sqrt(m_pi) * zi**2
524 end do
525
526 call comm_allreduce(dist%mpi_grp, eself)
527 call comm_allreduce(dist%mpi_grp, charge)
528
530 end subroutine ewald_self_interaction
531
533 subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
534 type(ion_interaction_t), intent(in) :: this
535 class(space_t), intent(in) :: space
536 type(lattice_vectors_t), intent(in) :: latt
537 type(atom_t), intent(in) :: atom(:)
538 integer, intent(in) :: natoms
539 real(real64), intent(in) :: pos(:,:)
540 real(real64), intent(inout) :: efourier
541 real(real64), intent(inout) :: force(:, :)
542 real(real64), intent(in) :: charge
543
544 real(real64) :: rcut, gmax_squared
545 integer :: iatom
546 integer :: ix, iy, iz, isph
547 real(real64) :: gvec(3), gred(3), gg2, gx
548 real(real64) :: factor
549 complex(real64) :: sumatoms, tmp(3), aa
550
551 complex(real64), allocatable :: phase(:)
552
553 push_sub(ewald_long_3d)
554
555 assert(space%dim == 3)
556 assert(space%periodic_dim == 3)
557
558 ! And the long-range part, using an Ewald sum
559 safe_allocate(phase(1:natoms))
560
561 ! get a converged value for the cutoff in g
562 rcut = sqrt(minval(sum(latt%klattice**2, dim=1)))
563
564 ! 9.5 is a constant that controls the range separation
565 isph = ceiling(9.5_real64*this%alpha/rcut)
566
567 ! First the G = 0 term (charge was calculated previously)
568 efourier = -m_pi*charge**2/(m_two*this%alpha**2*latt%rcell_volume)
569
570 ! Getting the G-shell cutoff
571 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
572
573 do ix = -isph, isph
574 do iy = -isph, isph
575 do iz = -isph, isph
576
577 ! Exploit k <-> -k symmetry
578 ! Only process one half of reciprocal space.
579 ! g=0 must also be removed from the sum
580 if (ix < 0) cycle
581 if (ix == 0 .and. iy < 0) cycle
582 if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
583
584 gred = [ix, iy, iz]
585 call kpoints_to_absolute(latt, gred, gvec)
586 gg2 = dot_product(gvec, gvec)
587
588 if (gg2 > gmax_squared*1.001_real64) cycle
589
590 gx = -0.25_real64*gg2/this%alpha**2
591
592 if (gx < -36.0_real64) cycle
593
594 ! We have used the k-> -k symmetry, hence the factor 4
595 factor = m_four*m_pi/latt%rcell_volume*exp(gx)/gg2
596
597 if (factor < epsilon(factor)) cycle
599 sumatoms = m_z0
600 !$omp parallel do private(iatom, gx, aa) reduction(+:sumatoms)
601 do iatom = 1, natoms
602 gx = sum(gvec*pos(:,iatom))
603 aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
604 phase(iatom) = aa
605 sumatoms = sumatoms + aa
606 end do
607
608 efourier = efourier + factor * real(sumatoms*conjg(sumatoms), real64)
609
610 do iatom = 1, natoms
611 tmp = m_zi*gvec*phase(iatom)
612 force(1:space%dim, iatom) = force(1:space%dim, iatom) - factor*real(conjg(tmp)*sumatoms + tmp*conjg(sumatoms), real64)
613
614 end do
615
616 end do
617 end do
618 end do
619
620 safe_deallocate_a(phase)
621
622 pop_sub(ewald_long_3d)
623
624 end subroutine ewald_long_3d
625
629 subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
630 type(ion_interaction_t), intent(in) :: this
631 class(space_t), intent(in) :: space
632 type(lattice_vectors_t), intent(in) :: latt
633 type(atom_t), intent(in) :: atom(:)
634 integer, intent(in) :: natoms
635 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
636 real(real64), intent(inout) :: efourier
637 real(real64), intent(inout) :: force(:, :)
638
639 real(real64) :: rcut, gmax_squared
640 integer :: iatom, jatom
641 integer :: ix, iy, ix_max, iy_max
642 real(real64) :: gvec(space%dim), gg2, gx, gg_abs
643 real(real64) :: factor,factor1,factor2, coeff
644 real(real64) :: dz_max, dz_ij, erfc1, erfc2, tmp_erf
645 real(real64), allocatable :: force_tmp(:,:)
646 real(real64), parameter :: tol = 1e-10_real64
647
648 push_sub(ewald_long_2d)
649
650 assert(space%periodic_dim == 2)
651 assert(space%dim == 2 .or. space%dim == 3)
652
653 ! And the long-range part, using an Ewald sum
654
655 ! Searching maximum distance
656 if (space%dim == 3) then
657 dz_max = m_zero
658 do iatom = 1, natoms
659 do jatom = iatom + 1, natoms
660 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
661 end do
662 end do
663 else
664 ! For a 2D system, all atoms are on the plane, so the distance is zero
665 dz_max = m_zero
666 end if
667
668 !get a converged value for the cutoff in g
669 rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
670 if (dz_max > tol) then
671 do
672 if (rcut * dz_max >= m_max_exp_arg) exit !Maximum double precision number
673 erfc1 = m_one - loct_erf(this%alpha*dz_max + m_half*rcut/this%alpha)
674 if (erfc1 * exp(rcut*dz_max) < 1.e-10_real64) exit
675 rcut = rcut * 1.414_real64
676 end do
677 end if
678
679 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
680 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
681
682 safe_allocate(force_tmp(1:space%dim, 1:natoms))
683 force_tmp = m_zero
684
685 ! First the G = 0 term
686 efourier = m_zero
687 factor = m_pi/latt%rcell_volume
688 !$omp parallel do private(jatom, dz_ij, tmp_erf, factor1, factor2) reduction(+:efourier,force_tmp) &
689 !$omp& collapse(2)
690 do iatom = this%dist%start, this%dist%end
691 do jatom = 1, natoms
692 ! efourier
693 if (space%dim == 3) then
694 dz_ij = pos(3, iatom) - pos(3, jatom)
695 else
696 dz_ij = m_zero
697 end if
698
699 tmp_erf = loct_erf(this%alpha*dz_ij)
700 factor1 = dz_ij*tmp_erf
701 factor2 = exp(-(this%alpha*dz_ij)**2)/(this%alpha*sqrt(m_pi))
702
703 efourier = efourier - factor &
704 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
705
706 ! force
707 if (iatom == jatom)cycle
708 if (abs(tmp_erf) < m_epsilon) cycle
709
710 if (space%dim == 3) then
711 force_tmp(3, iatom) = force_tmp(3, iatom) - (- m_two*factor) &
712 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * tmp_erf
713 end if
714
715 end do
716 end do
717
718 ! Getting the G-shell cutoff
719 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
720 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
721
722 !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, gx, dz_ij, erfc1, factor1, erfc2, factor2, coeff) &
723 !$omp& collapse(2) reduction(+:efourier, force_tmp)
724 do ix = -ix_max, ix_max
725 do iy = -iy_max, iy_max
726
727 gvec = ix*latt%klattice(:, 1) + iy*latt%klattice(:, 2)
728 gg2 = sum(gvec**2)
729
730 ! g=0 must be removed from the sum
731 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
732 gg_abs = sqrt(gg2)
733 factor = m_half*m_pi/(latt%rcell_volume*gg_abs)
734
735 do iatom = this%dist%start, this%dist%end
736 do jatom = iatom, natoms
737 ! efourier
738 gx = sum(gvec(1:2) * (pos(1:2, iatom) - pos(1:2, jatom)))
739 gx = gvec(1)*(pos(1, iatom) - pos(1, jatom)) + gvec(2)*(pos(2, iatom) - pos(2, jatom))
740 if (space%dim == 3) then
741 dz_ij = pos(3, iatom) - pos(3, jatom)
742 else
743 dz_ij = m_zero
744 end if
745
746 erfc1 = m_one - loct_erf(this%alpha*dz_ij + m_half*gg_abs/this%alpha)
747 if (abs(erfc1) > m_epsilon) then
748 factor1 = exp(gg_abs*dz_ij)*erfc1
749 else
750 factor1 = m_zero
751 end if
752 erfc2 = m_one - loct_erf(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)
753 if (abs(erfc2) > m_epsilon) then
754 factor2 = exp(-gg_abs*dz_ij)*erfc2
755 else
756 factor2 = m_zero
757 end if
758
759 if (iatom == jatom) then
760 coeff = m_one
761 else
762 coeff = m_two
763 end if
764
765 efourier = efourier &
766 + factor * coeff &
767 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
768 * cos(gx)* ( factor1 + factor2)
769
770 ! force
771 if (iatom == jatom) cycle
772
773 force_tmp(1:2, iatom) = force_tmp(1:2, iatom) &
774 + m_two * factor * gvec(1:2) &
775 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
776 *sin(gx)*(factor1 + factor2)
777
778 force_tmp(1:2, jatom) = force_tmp(1:2, jatom) &
779 - m_two * factor * gvec(1:2) &
780 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
781 *sin(gx)*(factor1 + factor2)
782
783 factor1 = gg_abs*erfc1 &
784 - m_two*this%alpha/sqrt(m_pi)*exp(-(this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
785 if (abs(factor1) > m_epsilon) then
786 factor1 = factor1*exp(gg_abs*dz_ij)
787 else
788 factor1 = m_zero
789 end if
790
791 factor2 = gg_abs*erfc2 &
792 - m_two*this%alpha/sqrt(m_pi)*exp(-(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
793 if (abs(factor2) > m_epsilon) then
794 factor2 = factor2*exp(-gg_abs*dz_ij)
795 else
796 factor2 = m_zero
797 end if
798
799 if (space%dim == 3) then
800 force_tmp(3, iatom) = force_tmp(3, iatom) &
801 - m_two*factor &
802 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
803 * cos(gx)* ( factor1 - factor2)
804 force_tmp(3, jatom) = force_tmp(3, jatom) &
805 + m_two*factor &
806 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
807 * cos(gx)* ( factor1 - factor2)
808 end if
809
810 end do
811 end do
812
813
814 end do
815 end do
816
817 call comm_allreduce(this%dist%mpi_grp, efourier)
818 call comm_allreduce(this%dist%mpi_grp, force_tmp)
819
820 force = force + force_tmp
821
822 safe_deallocate_a(force_tmp)
823
824 pop_sub(ewald_long_2d)
825 end subroutine ewald_long_2d
826
827 !TODO(Alex/Nicolas) Issue #950. Refactor: Move G=0 correction from ion-ion energy to pseudopotential energy
834 subroutine pseudopotential_correction_3d(dist, latt, atom, charge, epseudo)
835 type(distributed_t), intent(in) :: dist
836 type(lattice_vectors_t), intent(in) :: latt
837 type(atom_t), intent(in) :: atom(:)
838 real(real64), intent(out) :: epseudo
839
840 real(real64) :: zi
841 real(real64) :: charge
842 integer :: iatom
843
845
846 epseudo = m_zero
847 do iatom = dist%start, dist%end
848 select type(spec => atom(iatom)%species)
849 class is(pseudopotential_t)
850 zi = spec%get_zval()
851 epseudo = epseudo + m_pi *zi * &
852 (spec%ps%sigma_erf * sqrt(m_two))**2 / latt%rcell_volume * charge
853 end select
854 end do
855 call comm_allreduce(dist%mpi_grp, epseudo)
856
858
859 end subroutine pseudopotential_correction_3d
860
862 subroutine ion_interaction_stress(this, space, latt, atom, natoms, pos, stress_ii)
863 type(ion_interaction_t), intent(inout) :: this
864 class(space_t), intent(in) :: space
865 type(lattice_vectors_t), intent(in) :: latt
866 type(atom_t), intent(in) :: atom(:)
867 integer, intent(in) :: natoms
868 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
869 real(real64), intent(out) :: stress_ii(space%dim, space%dim)
870
871 real(real64) :: stress_short(1:space%dim, 1:space%dim), stress_Ewald(1:space%dim, 1:space%dim)
872
873 push_sub(ion_interaction_stress)
874
875 stress_ii = m_zero
876
877 ! Only implemented in the periodic case
878 assert(space%is_periodic())
879
880 ! Short range part in real space
881 call ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
882
883 ! Long range part in Fourier space
884 select case(space%periodic_dim)
885 case(3)
886 call ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
887 case(2)
888 call ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
889 case default
890 assert(.false.)
891 end select
892
893 stress_ii = stress_short + stress_ewald
894
896 end subroutine ion_interaction_stress
897
898 ! ---------------------------------------------------------
916
917 subroutine ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
918 type(ion_interaction_t), intent(inout) :: this
919 class(space_t), intent(in) :: space
920 type(lattice_vectors_t), intent(in) :: latt
921 type(atom_t), intent(in) :: atom(:)
922 integer, intent(in) :: natoms
923 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
924 real(real64), intent(out) :: stress_short(1:space%dim, 1:space%dim)
925
926 real(real64) :: xi(space%dim)
927 real(real64) :: r_ij, zi, zj, erfc, Hp, factor
928 integer :: iatom, jatom, icopy, idir, jdir
929 real(real64) :: alpha, rcut
930 type(lattice_iterator_t) :: latt_iter
931
933 call profiling_in("ION_ION_STRESS_SHORT")
934
935 ! Only implemented in the periodic case
936 assert(space%is_periodic())
937
938 alpha = this%alpha
939
940 ! See the code for the energy above to understand this parameter
941 rcut = 6.0_real64/alpha
942
943 ! the short-range part is calculated directly
944 stress_short = m_zero
945 latt_iter = lattice_iterator_t(latt, rcut)
946
947 do iatom = this%dist%start, this%dist%end
948 select type(spec => atom(iatom)%species)
949 class is(jellium_t)
950 cycle
951 end select
952 zi = atom(iatom)%species%get_zval()
953
954 do icopy = 1, latt_iter%n_cells
955 xi = pos(:, iatom) + latt_iter%get(icopy)
956
957 do jatom = 1, natoms
958 zj = atom(jatom)%species%get_zval()
959 r_ij = norm2(xi - pos(:, jatom))
960
961 if (r_ij < r_min_atom_dist) cycle
962
963 erfc = loct_erfc(alpha*r_ij)
964 hp = -m_two/sqrt(m_pi)*exp(-(alpha*r_ij)**2) - erfc/(alpha*r_ij)
965 factor = m_half*zj*zi*alpha*hp
966 do idir = 1, space%periodic_dim
967 do jdir = 1, space%periodic_dim
968 stress_short(idir, jdir) = stress_short(idir, jdir) &
969 - factor*(xi(idir) - pos(idir, jatom))*(xi(jdir) - pos(jdir, jatom))/(r_ij**2)
970 end do
971 end do
972
973 end do
974 end do
975 end do
976
977 if (this%dist%parallel) then
978 call comm_allreduce(this%dist%mpi_grp, stress_short)
979 end if
980
981 stress_short = stress_short/latt%rcell_volume
982
983 call profiling_out("ION_ION_STRESS_SHORT")
984
986 end subroutine ion_interaction_stress_short
987
988
989
990 ! ---------------------------------------------------------
1005 subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1006 type(ion_interaction_t), intent(inout) :: this
1007 class(space_t), intent(in) :: space
1008 type(lattice_vectors_t), intent(in) :: latt
1009 type(atom_t), intent(in) :: atom(:)
1010 integer, intent(in) :: natoms
1011 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
1012 real(real64), intent(out) :: stress_Ewald(3, 3)
1013
1014 real(real64) :: zi, rcut, gmax_squared
1015 integer :: iatom
1016 integer :: ix, iy, iz, isph, idim, idir, jdir
1017 real(real64) :: gred(3), gvec(3), gg2, gx
1018 real(real64) :: factor, charge, charge_sq, off_diagonal_weight
1019 complex(real64) :: sumatoms, aa
1020
1021 call profiling_in("STRESS_3D_EWALD")
1022 push_sub(ewald_3d_stress)
1023
1024 ! Currently this is only implemented for 3D
1025 assert(space%dim == 3)
1026 assert(space%periodic_dim == 3) ! Not working for mixed periodicity
1027 ! (klattice along the non-periodic directions is wrong)
1028 ! Anyway gg/gg2 is not correct for mixed periodicity
1029
1030 stress_ewald = m_zero
1031
1032 ! And the long-range part, using an Ewald sum
1033 charge = m_zero
1034 charge_sq = m_zero
1035 do iatom = 1, natoms
1036 zi = atom(iatom)%species%get_zval()
1037 charge = charge + zi
1038 charge_sq = charge_sq + zi**2
1039 end do
1040
1041 ! get a converged value for the cutoff in g
1042 rcut = huge(rcut)
1043 do idim = 1, space%periodic_dim
1044 rcut = min(rcut, sum(latt%klattice(1:space%periodic_dim, idim)**2))
1045 end do
1046
1047 rcut = sqrt(rcut)
1048
1049 isph = ceiling(9.5_real64*this%alpha/rcut)
1050
1051 ! Getting the G-shell cutoff
1052 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
1053
1054 do ix = -isph, isph
1055 do iy = -isph, isph
1056 do iz = -isph, isph
1057
1058 ! Exploit k <-> -k symmetry
1059 ! Only process one half of reciprocal space.
1060 ! g=0 must also be removed from the sum
1061 if (ix < 0) cycle
1062 if (ix == 0 .and. iy < 0) cycle
1063 if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
1064
1065 gred = [ix, iy, iz]
1066 call kpoints_to_absolute(latt, gred, gvec)
1067 gg2 = sum(gvec**2)
1068
1069 ! g=0 must be removed from the sum
1070 if (gg2 > gmax_squared*1.001_real64) cycle
1071
1072 gx = -0.25_real64*gg2/this%alpha**2
1073
1074 if (gx < -36.0_real64) cycle
1075
1076 ! We have used the k-> -k symmetry, hence the factor 4
1077 factor = m_four*m_pi*exp(gx)/(latt%rcell_volume*gg2)
1078
1079 if (factor < epsilon(factor)) cycle
1080
1081 sumatoms = m_z0
1082
1083 do iatom = 1, natoms
1084 gx = sum(gvec*pos(:, iatom))
1085 aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
1086 sumatoms = sumatoms + aa
1087 end do
1088
1089 factor = factor*abs(sumatoms)**2
1090 off_diagonal_weight = - m_two*factor/gg2*(0.25_real64*gg2/this%alpha**2+m_one)
1091
1092 do idir = 1, 3
1093 do jdir = 1, 3
1094 stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
1095 + gvec(idir) * gvec(jdir) * off_diagonal_weight
1096 end do
1097 stress_ewald(idir, idir) = stress_ewald(idir, idir) + factor
1098 end do
1099
1100 end do
1101 end do
1102 end do
1103
1104
1105 ! The G = 0 term of the Ewald summation
1106 factor = m_half*m_pi*charge**2/(latt%rcell_volume*this%alpha**2)
1107 do idir = 1,3
1108 stress_ewald(idir,idir) = stress_ewald(idir,idir) - factor
1109 end do
1110
1111 stress_ewald = stress_ewald / latt%rcell_volume
1112
1113
1114 call profiling_out("STRESS_3D_EWALD")
1115 pop_sub(ewald_3d_stress)
1116
1117 end subroutine ewald_3d_stress
1118
1119 ! ---------------------------------------------------------
1135 subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1136 type(ion_interaction_t), intent(inout) :: this
1137 type(space_t), intent(in) :: space
1138 type(lattice_vectors_t), intent(in) :: latt
1139 type(atom_t), intent(in) :: atom(:)
1140 integer, intent(in) :: natoms
1141 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
1142 real(real64), intent(out) :: stress_Ewald(3, 3)
1143
1144 real(real64) :: rcut, efourier
1145 integer :: iatom, jatom, idir, jdir
1146 integer :: ix, iy, ix_max, iy_max
1147 real(real64) :: gvec(3), gred(3), gg2, cos_gx, gg_abs, gmax_squared
1148 real(real64) :: factor,factor1,factor2, coeff, e_ewald
1149 real(real64) :: dz_max, z_ij, erfc1, erfc2, diff(3)
1150 real(real64), parameter :: tol = 1e-10_real64
1151
1152 push_sub(ewald_2d_stress)
1153
1154 assert(space%periodic_dim == 2)
1155 assert(space%dim == 3)
1156
1157 stress_ewald = m_zero
1158
1159 ! Searching maximum distance
1160 dz_max = m_zero
1161 do iatom = 1, natoms
1162 do jatom = iatom + 1, natoms
1163 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
1164 end do
1165 end do
1166
1167 !get a converged value for the cutoff in g
1168 ! Note: to understand these numbers, one needs to look into the energy routine for Ewald 2D
1169 rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
1170 if (dz_max > tol) then ! Else the code here does not work properly
1171 do
1172 if (rcut * dz_max >= m_max_exp_arg) exit !Maximum double precision number
1173 erfc1 = m_one - loct_erf(this%alpha*dz_max + m_half*rcut/this%alpha)
1174 if (erfc1 * exp(rcut*dz_max) < tol) exit
1175 rcut = rcut * 1.414_real64
1176 end do
1177 end if
1178
1179 ! First the G = 0 term
1180 efourier = m_zero
1181 factor = m_pi/latt%rcell_volume
1182 !$omp parallel do private(jatom, z_ij, factor1, factor2) reduction(+:efourier) collapse(2)
1183 do iatom = 1, natoms
1184 do jatom = 1, natoms
1185 z_ij = pos(3, iatom) - pos(3, jatom)
1186
1187 factor1 = z_ij * loct_erf(this%alpha*z_ij)
1188 factor2 = exp(-(this%alpha*z_ij)**2)/(this%alpha*sqrt(m_pi))
1189
1190 efourier = efourier - factor &
1191 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
1192 end do
1193 end do
1194
1195 ! Adding the G=0 term
1196 do idir = 1, 2
1197 stress_ewald(idir, idir) = efourier
1198 end do
1199
1200 ! Getting the G-shell cutoff
1201 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
1202 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
1203 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
1204 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
1205
1206 !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, diff, cos_gx, z_ij, idir, jdir, erfc1, factor1) &
1207 !$omp& private(erfc2, factor2, coeff, e_ewald) &
1208 !$omp& collapse(2) reduction(+:stress_Ewald)
1209 do ix = -ix_max, ix_max
1210 do iy = -iy_max, iy_max
1211
1212 gred = [ix, iy, 0]
1213 call kpoints_to_absolute(latt, gred, gvec)
1214 gg2 = dot_product(gvec,gvec)
1215
1216 ! g=0 must be removed from the sum
1217 if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1218
1219 gg_abs = sqrt(gg2)
1220 factor = m_fourth*m_pi/(latt%rcell_volume*this%alpha*gg2)
1221
1222 do iatom = 1, natoms
1223 do jatom = iatom, natoms
1224 diff = pos(:, iatom) - pos(:, jatom)
1225 cos_gx = cos(sum(gvec(1:2) * diff(1:2)))
1226 z_ij = diff(3)
1227
1228 factor1 = screening_function_2d(this%alpha, z_ij, gg_abs, erfc1)
1229 factor2 = screening_function_2d(this%alpha,-z_ij, gg_abs, erfc2)
1230
1231 if (iatom == jatom) then
1232 coeff = m_one
1233 else
1234 coeff = m_two
1235 end if
1236
1237 do idir = 1, 2
1238 do jdir = 1, 2
1239 stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
1240 - factor*gvec(idir)*gvec(jdir) * cos_gx * (factor1 + factor2) * coeff&
1241 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval()
1242 end do
1243 end do
1244
1245 if (abs(erfc1) > m_epsilon) then
1246 factor1 = exp(-gg_abs*z_ij)*erfc1
1247 else
1248 factor1 = m_zero
1249 end if
1250 if (abs(erfc2) > m_epsilon) then
1251 factor2 = exp(gg_abs*z_ij)*erfc2
1252 else
1253 factor2 = m_zero
1254 end if
1255
1256 e_ewald = m_half * m_pi/latt%rcell_volume * coeff &
1257 * atom(iatom)%species%get_zval() * atom(jatom)%species%get_zval() &
1258 * cos_gx / gg_abs * (factor1 + factor2)
1259
1260 do idir = 1, 2
1261 stress_ewald(idir, idir) = stress_ewald(idir, idir) + e_ewald
1262 end do
1263
1264 end do !jatom
1265 end do !iatom
1266 end do !iy
1267 end do !ix
1268
1269 !call comm_allreduce(this%dist%mpi_grp, stress_Ewald)
1270
1271 stress_ewald = stress_ewald / latt%rcell_volume
1272
1273 pop_sub(ewald_2d_stress)
1274 end subroutine ewald_2d_stress
1275
1276 ! ---------------------------------------------------------
1278 real(real64) function screening_function_2d(alpha, z_ij, gg_abs, erfc) result(factor)
1279 real(real64), intent(in) :: alpha
1280 real(real64), intent(in) :: z_ij
1281 real(real64), intent(in) :: gg_abs
1282 real(real64), intent(out) :: erfc
1283
1284 real(real64) :: arg
1285
1286 arg = -alpha*z_ij + m_half*gg_abs/alpha
1287 erfc = m_one - loct_erf(arg)
1288 factor = m_two*alpha*(m_one/gg_abs + z_ij)*erfc - m_two/sqrt(m_pi)*exp(-arg**2)
1289 factor = factor*exp(-gg_abs*z_ij)
1290
1291 end function screening_function_2d
1292
1293 ! ---------------------------------------------------------
1294
1295 subroutine ion_interaction_test(space, latt, atom, natoms, pos, lsize, &
1296 namespace, mc)
1297 class(space_t), intent(in) :: space
1298 type(lattice_vectors_t), intent(in) :: latt
1299 type(atom_t), intent(in) :: atom(:)
1300 integer, intent(in) :: natoms
1301 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
1302 real(real64), intent(in) :: lsize(:)
1303 type(namespace_t), intent(in) :: namespace
1304 type(multicomm_t), intent(in) :: mc
1305
1306 type(ion_interaction_t) :: ion_interaction
1307 real(real64) :: energy
1308 real(real64), allocatable :: force(:, :), force_components(:, :, :)
1309 real(real64) :: energy_components(1:ION_NUM_COMPONENTS)
1310 integer :: iatom, idir
1311
1312 push_sub(ion_interaction_test)
1313
1314 call ion_interaction_init(ion_interaction, namespace, space, natoms)
1315 call ion_interaction_init_parallelization(ion_interaction, natoms, mc)
1316
1317 safe_allocate(force(1:space%dim, 1:natoms))
1318 safe_allocate(force_components(1:space%dim, 1:natoms, 1:ion_num_components))
1319
1320 call ion_interaction_calculate(ion_interaction, space, latt, atom, natoms, pos, lsize, energy, force, &
1321 energy_components = energy_components, force_components = force_components)
1322
1323 call messages_write('Ionic energy =')
1324 call messages_write(energy, fmt = '(f20.10)')
1325 call messages_info(namespace=namespace)
1326
1327 call messages_write('Real space energy =')
1328 call messages_write(energy_components(ion_component_real), fmt = '(f20.10)')
1329 call messages_info(namespace=namespace)
1330
1331 call messages_write('Self energy =')
1332 call messages_write(energy_components(ion_component_self), fmt = '(f20.10)')
1333 call messages_info(namespace=namespace)
1334
1335 call messages_write('Fourier energy =')
1336 call messages_write(energy_components(ion_component_fourier), fmt = '(f20.10)')
1337 call messages_info(namespace=namespace)
1338
1339 call messages_info(namespace=namespace)
1340
1341 do iatom = 1, natoms
1342 call messages_write('Ionic force atom')
1343 call messages_write(iatom)
1344 call messages_write(' =')
1345 do idir = 1, space%dim
1346 call messages_write(force(idir, iatom), fmt = '(f20.10)')
1347 end do
1348 call messages_info(namespace=namespace)
1349
1350 call messages_write('Real space force atom')
1351 call messages_write(iatom)
1352 call messages_write(' =')
1353 do idir = 1, space%dim
1354 call messages_write(force_components(idir, iatom, ion_component_real), fmt = '(f20.10)')
1355 end do
1356 call messages_info(namespace=namespace)
1357
1358 call messages_write('Fourier space force atom')
1359 call messages_write(iatom)
1360 call messages_write(' =')
1361 do idir = 1, space%dim
1362 call messages_write(force_components(idir, iatom, ion_component_fourier), fmt = '(f20.10)')
1363 end do
1364 call messages_info(namespace=namespace)
1365
1366 call messages_info(namespace=namespace)
1367 end do
1368
1369 safe_deallocate_a(force)
1370 safe_deallocate_a(force_components)
1372 call ion_interaction_end(ion_interaction)
1373
1374 pop_sub(ion_interaction_test)
1375 end subroutine ion_interaction_test
1376
1377end module ion_interaction_oct_m
1378
1379!! Local Variables:
1380!! mode: f90
1381!! coding: utf-8
1382!! 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:293
pure logical function, public any_species_is_jellium_sphere(atom)
Check if any species is a jellium sphere.
Definition: atom.F90:310
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)
Computes the long-range part of the 2D Ewald summation.
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
static double f(double w, void *p)
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)