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
409 subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
411 class(
space_t),
intent(in) :: space
413 type(
atom_t),
intent(in) :: atom(:)
414 real(real64),
intent(in) :: pos(:, :)
416 real(real64),
intent(in) :: alpha
417 real(real64),
intent(out) :: ereal
418 real(real64),
intent(inout) :: force(:, :)
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
427 real(real64) :: charge, coeff
433 rcut = 6.0_real64 / alpha
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
446 do icopy = 1, latt_iter%n_cells
447 rnorm = norm2(latt_iter%get(icopy))
449 if (rnorm > rcut) cycle
451 ereal = ereal +
m_half * charge * erfc /rnorm
459 do iatom = dist%start, dist%end
460 if (.not. atom(iatom)%species%represents_real_atom()) cycle
461 zi = atom(iatom)%species%get_zval()
464 do jatom = iatom + 1, natoms
465 zj = atom(jatom)%species%get_zval()
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
478 ereal = ereal + charge * erfc
480 force_real(:) = charge * (xi - pos(:, jatom)) * &
481 (erfc + coeff *
exp(-(alpha*rnorm)**2)) / rnorm**2
484 force(1:space%dim, jatom) = force(1:space%dim, jatom) - force_real
487 force(1:space%dim, iatom) = force(1:space%dim, iatom) + force_real
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
520 do iatom = dist%start, dist%end
521 zi = atom(iatom)%species%get_zval()
523 eself = eself - alpha /
sqrt(
m_pi) * zi**2
533 subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
535 class(
space_t),
intent(in) :: space
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
544 real(real64) :: rcut, gmax_squared
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
551 complex(real64),
allocatable :: phase(:)
555 assert(space%dim == 3)
556 assert(space%periodic_dim == 3)
559 safe_allocate(phase(1:natoms))
562 rcut =
sqrt(minval(sum(latt%klattice**2, dim=1)))
565 isph = ceiling(9.5_real64*this%alpha/rcut)
568 efourier = -
m_pi*charge**2/(
m_two*this%alpha**2*latt%rcell_volume)
571 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
581 if (ix == 0 .and. iy < 0) cycle
582 if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
586 gg2 = dot_product(gvec, gvec)
588 if (gg2 > gmax_squared*1.001_real64) cycle
590 gx = -0.25_real64*gg2/this%alpha**2
592 if (gx < -36.0_real64) cycle
597 if (factor < epsilon(factor)) cycle
602 gx = sum(gvec*pos(:,iatom))
603 aa = atom(iatom)%species%get_zval()*cmplx(
cos(gx),
sin(gx), real64)
605 sumatoms = sumatoms + aa
608 efourier = efourier + factor * real(sumatoms*conjg(sumatoms), real64)
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)
620 safe_deallocate_a(phase)
629 subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
631 class(
space_t),
intent(in) :: space
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(:, :)
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
650 assert(space%periodic_dim == 2)
651 assert(space%dim == 2 .or. space%dim == 3)
656 if (space%dim == 3)
then
659 do jatom = iatom + 1, natoms
660 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
669 rcut =
m_two*this%alpha*4.6_real64 +
m_two*this%alpha**2*dz_max
670 if (dz_max > tol)
then
674 if (erfc1 *
exp(rcut*dz_max) < 1.e-10_real64)
exit
675 rcut = rcut * 1.414_real64
679 ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
680 iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
682 safe_allocate(force_tmp(1:space%dim, 1:natoms))
687 factor =
m_pi/latt%rcell_volume
690 do iatom = this%dist%start, this%dist%end
693 if (space%dim == 3)
then
694 dz_ij = pos(3, iatom) - pos(3, jatom)
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))
703 efourier = efourier - factor &
704 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
707 if (iatom == jatom)cycle
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
719 gmax_squared = sum(ix_max*latt%klattice(:, 1)**2)
720 gmax_squared = min(gmax_squared, sum(iy_max*latt%klattice(:, 2)**2))
724 do ix = -ix_max, ix_max
725 do iy = -iy_max, iy_max
727 gvec = ix*latt%klattice(:, 1) + iy*latt%klattice(:, 2)
731 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
733 factor =
m_half*
m_pi/(latt%rcell_volume*gg_abs)
735 do iatom = this%dist%start, this%dist%end
736 do jatom = iatom, natoms
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)
748 factor1 =
exp(gg_abs*dz_ij)*erfc1
754 factor2 =
exp(-gg_abs*dz_ij)*erfc2
759 if (iatom == jatom)
then
765 efourier = efourier &
767 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
768 *
cos(gx)* ( factor1 + factor2)
771 if (iatom == jatom) cycle
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)
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)
783 factor1 = gg_abs*erfc1 &
786 factor1 = factor1*
exp(gg_abs*dz_ij)
791 factor2 = gg_abs*erfc2 &
794 factor2 = factor2*
exp(-gg_abs*dz_ij)
799 if (space%dim == 3)
then
800 force_tmp(3, iatom) = force_tmp(3, iatom) &
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) &
806 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
807 *
cos(gx)* ( factor1 - factor2)
820 force = force + force_tmp
822 safe_deallocate_a(force_tmp)
837 type(
atom_t),
intent(in) :: atom(:)
838 real(real64),
intent(out) :: epseudo
841 real(real64) :: charge
847 do iatom = dist%start, dist%end
848 select type(spec => atom(iatom)%species)
851 epseudo = epseudo +
m_pi *zi * &
852 (spec%ps%sigma_erf *
sqrt(
m_two))**2 / latt%rcell_volume * charge
864 class(
space_t),
intent(in) :: space
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)
871 real(real64) :: stress_short(1:space%dim, 1:space%dim), stress_Ewald(1:space%dim, 1:space%dim)
878 assert(space%is_periodic())
884 select case(space%periodic_dim)
886 call ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
888 call ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
893 stress_ii = stress_short + stress_ewald
919 class(
space_t),
intent(in) :: space
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)
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
936 assert(space%is_periodic())
941 rcut = 6.0_real64/alpha
947 do iatom = this%dist%start, this%dist%end
948 select type(spec => atom(iatom)%species)
952 zi = atom(iatom)%species%get_zval()
954 do icopy = 1, latt_iter%n_cells
955 xi = pos(:, iatom) + latt_iter%get(icopy)
958 zj = atom(jatom)%species%get_zval()
959 r_ij = norm2(xi - pos(:, jatom))
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)
977 if (this%dist%parallel)
then
981 stress_short = stress_short/latt%rcell_volume
1005 subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1007 class(
space_t),
intent(in) :: space
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)
1014 real(real64) :: zi, rcut, gmax_squared
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
1025 assert(space%dim == 3)
1026 assert(space%periodic_dim == 3)
1035 do iatom = 1, natoms
1036 zi = atom(iatom)%species%get_zval()
1037 charge = charge + zi
1038 charge_sq = charge_sq + zi**2
1043 do idim = 1, space%periodic_dim
1044 rcut = min(rcut, sum(latt%klattice(1:space%periodic_dim, idim)**2))
1049 isph = ceiling(9.5_real64*this%alpha/rcut)
1052 gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
1062 if (ix == 0 .and. iy < 0) cycle
1063 if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
1070 if (gg2 > gmax_squared*1.001_real64) cycle
1072 gx = -0.25_real64*gg2/this%alpha**2
1074 if (gx < -36.0_real64) cycle
1079 if (factor < epsilon(factor)) cycle
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
1089 factor = factor*abs(sumatoms)**2
1090 off_diagonal_weight = -
m_two*factor/gg2*(0.25_real64*gg2/this%alpha**2+
m_one)
1094 stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
1095 + gvec(idir) * gvec(jdir) * off_diagonal_weight
1097 stress_ewald(idir, idir) = stress_ewald(idir, idir) + factor
1106 factor =
m_half*
m_pi*charge**2/(latt%rcell_volume*this%alpha**2)
1108 stress_ewald(idir,idir) = stress_ewald(idir,idir) - factor
1111 stress_ewald = stress_ewald / latt%rcell_volume
1135 subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
1137 type(
space_t),
intent(in) :: space
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)
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
1154 assert(space%periodic_dim == 2)
1155 assert(space%dim == 3)
1161 do iatom = 1, natoms
1162 do jatom = iatom + 1, natoms
1163 dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
1169 rcut =
m_two*this%alpha*4.6_real64 +
m_two*this%alpha**2*dz_max
1170 if (dz_max > tol)
then
1174 if (erfc1 *
exp(rcut*dz_max) < tol)
exit
1175 rcut = rcut * 1.414_real64
1181 factor =
m_pi/latt%rcell_volume
1183 do iatom = 1, natoms
1184 do jatom = 1, natoms
1185 z_ij = pos(3, iatom) - pos(3, jatom)
1187 factor1 = z_ij *
loct_erf(this%alpha*z_ij)
1188 factor2 =
exp(-(this%alpha*z_ij)**2)/(this%alpha*
sqrt(
m_pi))
1190 efourier = efourier - factor &
1191 * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
1197 stress_ewald(idir, idir) = efourier
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))
1209 do ix = -ix_max, ix_max
1210 do iy = -iy_max, iy_max
1214 gg2 = dot_product(gvec,gvec)
1217 if (gg2 <
m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
1220 factor =
m_fourth*
m_pi/(latt%rcell_volume*this%alpha*gg2)
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)))
1231 if (iatom == jatom)
then
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()
1246 factor1 =
exp(-gg_abs*z_ij)*erfc1
1251 factor2 =
exp(gg_abs*z_ij)*erfc2
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)
1261 stress_ewald(idir, idir) = stress_ewald(idir, idir) + e_ewald
1271 stress_ewald = stress_ewald / latt%rcell_volume
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
1286 arg = -alpha*z_ij +
m_half*gg_abs/alpha
1289 factor = factor*
exp(-gg_abs*z_ij)
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
1307 real(real64) :: energy
1308 real(real64),
allocatable :: force(:, :), force_components(:, :, :)
1309 real(real64) :: energy_components(1:ION_NUM_COMPONENTS)
1310 integer :: iatom, idir
1317 safe_allocate(force(1:space%dim, 1:natoms))
1318 safe_allocate(force_components(1:space%dim, 1:natoms, 1:ion_num_components))
1321 energy_components = energy_components, force_components = force_components)
1323 call messages_write(
'Ionic energy =')
1324 call messages_write(energy, fmt =
'(f20.10)')
1325 call messages_info(namespace=namespace)
1327 call messages_write(
'Real space energy =')
1329 call messages_info(namespace=namespace)
1331 call messages_write(
'Self energy =')
1333 call messages_info(namespace=namespace)
1335 call messages_write(
'Fourier energy =')
1337 call messages_info(namespace=namespace)
1339 call messages_info(namespace=namespace)
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)')
1348 call messages_info(namespace=namespace)
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)')
1356 call messages_info(namespace=namespace)
1358 call messages_write(
'Fourier space force atom')
1359 call messages_write(iatom)
1360 call messages_write(
' =')
1361 do idir = 1, space%dim
1364 call messages_info(namespace=namespace)
1366 call messages_info(namespace=namespace)
1369 safe_deallocate_a(force)
1370 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)
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)
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.
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...