58 type(distributed_t) :: dist
61 integer,
parameter :: &
62 ION_COMPONENT_REAL = 1, &
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
87 call parse_variable(namespace,
'EwaldAlpha', 0.21_real64, this%alpha)
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.')
102 type(ion_interaction_t),
intent(inout) :: this
103 integer,
intent(in) :: natoms
104 type(multicomm_t),
intent(in) :: mc
120 type(ion_interaction_t),
intent(inout) :: this
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(:, :, :)
153 if (
present(energy_components))
then
155 energy_components =
m_zero
158 if (
present(force_components))
then
167 if (space%is_periodic())
then
172 call ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, energy_components, force_components)
190 class(
space_t),
intent(in) :: space
191 type(
atom_t),
intent(in) :: atom(:)
192 real(real64),
intent(in) :: lsize(:)
193 real(real64) :: energy
198 assert(
size(atom) == 1)
200 assert(space%periodic_dim == 2)
202 select type(spec => atom(1)%species)
204 area = lsize(1) * lsize(2) *
m_four
205 energy =
m_pi * spec%get_density(lsize) **2 * area * spec%thickness()**3 /
m_three
229 real(real64),
intent(in) :: lsize(:)
230 real(real64) :: energy
234 logical :: lattice_is_orthogonal
240 lattice_is_orthogonal = .not. latt%nonorthogonal
242 do iatom = dist%start, dist%end
243 spec => atom(iatom)%species
254 assert(lattice_is_orthogonal)
255 energy = energy +
m_pi * zi**2 / (
m_four * lsize(1)*lsize(2)) * spec%thickness() /
m_three
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(:, :)
277 class(
species_t),
pointer :: species_i, species_j
278 real(real64) :: r(space%dim), f(space%dim)
279 real(real64) :: r_mag
281 real(real64) :: zi, zj
282 integer :: iatom, jatom, natoms
288 force(1:space%dim, 1:natoms) =
m_zero
290 do iatom = dist%start, dist%end
291 species_i => atom(iatom)%species
292 zi = species_i%get_zval()
294 do jatom = iatom + 1, natoms
295 species_j => atom(jatom)%species
296 zj = species_j%get_zval()
298 r = pos(:, iatom) - pos(:, jatom)
300 u_e = zi * zj / r_mag
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)
312 nullify(species_i, species_j)
320 energy_components, force_components)
322 class(
space_t),
intent(in) :: space
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(:, :, :)
332 real(real64) :: ereal, efourier, epseudo, eself
333 real(real64) :: charge
338 force(1:space%dim, 1:natoms) =
m_zero
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)
348 select case (space%periodic_dim)
362 call ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
364 call ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
370 if (
present(energy_components))
then
371 energy_components(ion_component_real) = ereal
376 if (
present(force_components))
then
379 force(1:space%dim, 1:natoms) - force_components(1:space%dim, 1:natoms, ion_component_real)
382 energy = ereal + efourier + eself + epseudo
407 subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
409 class(
space_t),
intent(in) :: space
411 type(
atom_t),
intent(in) :: atom(:)
412 real(real64),
intent(in) :: pos(:, :)
414 real(real64),
intent(in) :: alpha
415 real(real64),
intent(out) :: ereal
416 real(real64),
intent(inout) :: force(:, :)
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
430 rcut = 6.0_real64 / alpha
435 do iatom = dist%start, dist%end
436 if (.not. atom(iatom)%species%represents_real_atom()) cycle
437 zi = atom(iatom)%species%get_zval()
441 do icopy = 1, latt_iter%n_cells
442 rnorm = norm2(latt_iter%get(icopy))
444 if (rnorm > rcut) cycle
446 ereal = ereal +
m_half * zi * zi * erfc /rnorm
451 do jatom = iatom + 1, natoms
452 zj = atom(jatom)%species%get_zval()
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
461 force_real(:) = zj * zi * (xi - pos(:, jatom)) * &
462 (erfc / rnorm +
m_two * alpha /
sqrt(
m_pi) *
exp(-(alpha*rnorm)**2)) / rnorm**2
465 ereal = ereal + zj*zi*erfc/rnorm
468 force(1:space%dim, jatom) = force(1:space%dim, jatom) - force_real
471 force(1:space%dim, iatom) = force(1:space%dim, iatom) + force_real
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
504 do iatom = dist%start, dist%end
505 zi = atom(iatom)%species%get_zval()
507 eself = eself - alpha /
sqrt(
m_pi) * zi**2
517 subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
519 class(
space_t),
intent(in) :: space
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
528 real(real64) :: rcut, gmax_squared
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
535 complex(real64),
allocatable :: phase(:)
539 assert(space%dim == 3)
540 assert(space%periodic_dim == 3)
543 safe_allocate(phase(1:natoms))
546 rcut =
sqrt(minval(sum(latt%klattice**2, dim=1)))
549 isph = ceiling(9.5_real64*this%alpha/rcut)
552 efourier = -
m_pi*charge**2/(
m_two*this%alpha**2*latt%rcell_volume)
555 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
563 gg2 = dot_product(gvec, gvec)
566 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
568 gx = -0.25_real64*gg2/this%alpha**2
570 if (gx < -36.0_real64) cycle
574 if (factor < epsilon(factor)) cycle
579 gx = sum(gvec*pos(:,iatom))
580 aa = atom(iatom)%species%get_zval()*cmplx(
cos(gx),
sin(gx), real64)
582 sumatoms = sumatoms + aa
585 efourier = efourier + real(factor*sumatoms*conjg(sumatoms), real64)
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)
597 safe_deallocate_a(phase)
604 subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
606 class(
space_t),
intent(in) :: space
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(:, :)
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
625 assert(space%periodic_dim == 2)
626 assert(space%dim == 2 .or. space%dim == 3)
631 if (space%dim == 3)
then
634 do jatom = iatom + 1, natoms
635 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
644 rcut =
m_two*this%alpha*4.6_real64 +
m_two*this%alpha**2*dz_max
645 if (dz_max > tol)
then
649 if (erfc1 *
exp(rcut*dz_max) < 1.e-10_real64)
exit
650 rcut = rcut * 1.414_real64
654 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
655 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
657 safe_allocate(force_tmp(1:space%dim, 1:natoms))
662 factor =
m_pi/latt%rcell_volume
665 do iatom = this%dist%start, this%dist%end
668 if (space%dim == 3)
then
669 dz_ij = pos(3, iatom) - pos(3, jatom)
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))
678 efourier = efourier - factor &
679 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
682 if (iatom == jatom)cycle
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
694 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
695 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
699 do ix = -ix_max, ix_max
700 do iy = -iy_max, iy_max
702 gvec = ix*latt%klattice(:, 1) + iy*latt%klattice(:, 2)
706 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
708 factor =
m_half*
m_pi/(latt%rcell_volume*gg_abs)
710 do iatom = this%dist%start, this%dist%end
711 do jatom = iatom, natoms
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)
723 factor1 =
exp(gg_abs*dz_ij)*erfc1
729 factor2 =
exp(-gg_abs*dz_ij)*erfc2
734 if (iatom == jatom)
then
740 efourier = efourier &
742 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
743 *
cos(gx)* ( factor1 + factor2)
746 if (iatom == jatom) cycle
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)
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)
758 factor1 = gg_abs*erfc1 &
761 factor1 = factor1*
exp(gg_abs*dz_ij)
766 factor2 = gg_abs*erfc2 &
769 factor2 = factor2*
exp(-gg_abs*dz_ij)
774 if (space%dim == 3)
then
775 force_tmp(3, iatom) = force_tmp(3, iatom) &
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) &
781 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
782 *
cos(gx)* ( factor1 - factor2)
795 force = force + force_tmp
797 safe_deallocate_a(force_tmp)
812 type(
atom_t),
intent(in) :: atom(:)
813 real(real64),
intent(out) :: epseudo
816 real(real64) :: charge
822 do iatom = dist%start, dist%end
823 select type(spec => atom(iatom)%species)
826 epseudo = epseudo +
m_pi *zi * &
827 (spec%ps%sigma_erf *
sqrt(
m_two))**2 / latt%rcell_volume * charge
839 class(
space_t),
intent(in) :: space
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)
846 real(real64) :: stress_short(1:space%dim, 1:space%dim), stress_Ewald(1:space%dim, 1:space%dim)
853 assert(space%is_periodic())
859 select case(space%periodic_dim)
861 call ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
863 call ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
868 stress_ii = stress_short + stress_ewald
892 class(
space_t),
intent(in) :: space
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)
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
909 assert(space%is_periodic())
914 rcut = 6.0_real64/alpha
920 do iatom = this%dist%start, this%dist%end
921 select type(spec => atom(iatom)%species)
925 zi = atom(iatom)%species%get_zval()
927 do icopy = 1, latt_iter%n_cells
928 xi = pos(:, iatom) + latt_iter%get(icopy)
931 zj = atom(jatom)%species%get_zval()
932 r_ij = norm2(xi - pos(:, jatom))
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)
950 if (this%dist%parallel)
then
954 stress_short = stress_short/latt%rcell_volume
975 subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
977 class(
space_t),
intent(in) :: space
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
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
995 assert(space%dim == 3)
996 assert(space%periodic_dim == 3)
1005 do iatom = 1, natoms
1006 zi = atom(iatom)%species%get_zval()
1007 charge = charge + zi
1008 charge_sq = charge_sq + zi**2
1013 do idim = 1, space%periodic_dim
1014 rcut = min(rcut, sum(latt%klattice(1:space%periodic_dim, idim)**2))
1019 isph = ceiling(9.5_real64*this%alpha/rcut)
1022 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
1033 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1035 gx = -0.25_real64*gg2/this%alpha**2
1037 if (gx < -36.0_real64) cycle
1041 if (factor < epsilon(factor)) cycle
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
1051 factor = factor*abs(sumatoms)**2
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)
1058 stress_ewald(idir, idir) = stress_ewald(idir, idir) + factor
1067 factor =
m_half*
m_pi*charge**2/(latt%rcell_volume*this%alpha**2)
1069 stress_ewald(idir,idir) = stress_ewald(idir,idir) - factor
1072 stress_ewald = stress_ewald / latt%rcell_volume
1093 subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1095 type(
space_t),
intent(in) :: space
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)
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
1112 assert(space%periodic_dim == 2)
1113 assert(space%dim == 3)
1119 do iatom = 1, natoms
1120 do jatom = iatom + 1, natoms
1121 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
1127 rcut =
m_two*this%alpha*4.6_real64 +
m_two*this%alpha**2*dz_max
1128 if (dz_max > tol)
then
1132 if (erfc1 *
exp(rcut*dz_max) < tol)
exit
1133 rcut = rcut * 1.414_real64
1139 factor =
m_pi/latt%rcell_volume
1141 do iatom = 1, natoms
1142 do jatom = 1, natoms
1143 z_ij = pos(3, iatom) - pos(3, jatom)
1145 factor1 = z_ij *
loct_erf(this%alpha*z_ij)
1146 factor2 =
exp(-(this%alpha*z_ij)**2)/(this%alpha*
sqrt(
m_pi))
1148 efourier = efourier - factor &
1149 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
1155 stress_ewald(idir, idir) = efourier
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))
1167 do ix = -ix_max, ix_max
1168 do iy = -iy_max, iy_max
1172 gg2 = dot_product(gvec,gvec)
1175 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1178 factor =
m_fourth*
m_pi/(latt%rcell_volume*this%alpha*gg2)
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)))
1189 if (iatom == jatom)
then
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()
1204 factor1 =
exp(-gg_abs*z_ij)*erfc1
1209 factor2 =
exp(gg_abs*z_ij)*erfc2
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)
1219 stress_ewald(idir, idir) = stress_ewald(idir, idir) + e_ewald
1229 stress_ewald = stress_ewald / latt%rcell_volume
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
1244 arg = -alpha*z_ij +
m_half*gg_abs/alpha
1247 factor = factor*
exp(-gg_abs*z_ij)
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
1265 real(real64) :: energy
1266 real(real64),
allocatable :: force(:, :), force_components(:, :, :)
1267 real(real64) :: energy_components(1:ION_NUM_COMPONENTS)
1268 integer :: iatom, idir
1275 safe_allocate(force(1:space%dim, 1:natoms))
1276 safe_allocate(force_components(1:space%dim, 1:natoms, 1:ion_num_components))
1279 energy_components = energy_components, force_components = force_components)
1281 call messages_write(
'Ionic energy =')
1282 call messages_write(energy, fmt =
'(f20.10)')
1283 call messages_info(namespace=namespace)
1285 call messages_write(
'Real space energy =')
1287 call messages_info(namespace=namespace)
1289 call messages_write(
'Self energy =')
1291 call messages_info(namespace=namespace)
1293 call messages_write(
'Fourier energy =')
1295 call messages_info(namespace=namespace)
1297 call messages_info(namespace=namespace)
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)')
1306 call messages_info(namespace=namespace)
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)')
1314 call messages_info(namespace=namespace)
1316 call messages_write(
'Fourier space force atom')
1317 call messages_write(iatom)
1318 call messages_write(
' =')
1319 do idir = 1, space%dim
1322 call messages_info(namespace=namespace)
1324 call messages_info(namespace=namespace)
1327 safe_deallocate_a(force)
1328 safe_deallocate_a(force_components)
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.
pure logical function, public any_species_is_jellium_sphere(atom)
Check if any species is a jellium sphere.
type(debug_t), save, public debug
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
real(real64), parameter, public m_max_exp_arg
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_fourth
complex(real64), parameter, public m_z0
complex(real64), parameter, public m_zi
real(real64), parameter, public r_min_atom_dist
Minimal distance between two distinguishable atoms.
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_three
real(real64), parameter, public m_five
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)
subroutine, public messages_not_implemented(feature, namespace)
subroutine, public messages_warning(no_lines, all_nodes, namespace)
This module handles the communicators for the various parallelization strategies.
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
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...