Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Separate short-range forces by signature #4768

Merged
merged 1 commit into from
Aug 8, 2023
Merged

Separate short-range forces by signature #4768

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 11 additions & 4 deletions src/core/constraints/ShapeBasedConstraint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,11 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,

if (dist > 0) {
outer_normal_vec = -dist_vec / dist;
pf = calc_non_bonded_pair_force(p, part_rep, ia_params, dist_vec, dist,
get_ptr(coulomb_kernel));
pf = calc_central_radial_force(p, part_rep, ia_params, dist_vec, dist) +
calc_central_radial_charge_force(p, part_rep, ia_params, dist_vec,
dist, get_ptr(coulomb_kernel)) +
calc_non_central_force(p, part_rep, ia_params, dist_vec, dist);

#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_force =
Expand All @@ -103,8 +106,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
#endif
} else if (m_penetrable && (dist <= 0)) {
if ((!m_only_positive) && (dist < 0)) {
pf = calc_non_bonded_pair_force(p, part_rep, ia_params, dist_vec, -dist,
get_ptr(coulomb_kernel));
pf =
calc_central_radial_force(p, part_rep, ia_params, dist_vec, -dist) +
calc_central_radial_charge_force(p, part_rep, ia_params, dist_vec,
-dist, get_ptr(coulomb_kernel)) +
calc_non_central_force(p, part_rep, ia_params, dist_vec, -dist);

#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_force = dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,
Expand Down
50 changes: 38 additions & 12 deletions src/core/forces_inline.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,10 +74,11 @@

#include <tuple>

inline ParticleForce calc_non_bonded_pair_force(
Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
Utils::Vector3d const &d, double const dist,
Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel) {
inline ParticleForce calc_central_radial_force(Particle const &p1,
Particle const &p2,
IA_parameters const &ia_params,
Utils::Vector3d const &d,
double const dist) {

ParticleForce pf{};
double force_factor = 0;
Expand Down Expand Up @@ -133,19 +134,39 @@ inline ParticleForce calc_non_bonded_pair_force(
#ifdef LJCOS2
force_factor += ljcos2_pair_force_factor(ia_params, dist);
#endif
/* Thole damping */
#ifdef THOLE
pf.f += thole_pair_force(p1, p2, ia_params, d, dist, coulomb_kernel);
#endif
/* tabulated */
#ifdef TABULATED
force_factor += tabulated_pair_force_factor(ia_params, dist);
#endif
pf.f += force_factor * d;
return pf;
}

inline ParticleForce calc_central_radial_charge_force(
Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
Utils::Vector3d const &d, double const dist,
Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel) {

ParticleForce pf{};
/* Thole damping */
#ifdef THOLE
pf.f += thole_pair_force(p1, p2, ia_params, d, dist, coulomb_kernel);
#endif

return pf;
}

inline ParticleForce calc_non_central_force(Particle const &p1,
Particle const &p2,
IA_parameters const &ia_params,
Utils::Vector3d const &d,
double const dist) {

ParticleForce pf{};
/* Gay-Berne */
#ifdef GAY_BERNE
pf += gb_pair_force(p1.quat(), p2.quat(), ia_params, d, dist);
#endif
pf.f += force_factor * d;
return pf;
}

Expand Down Expand Up @@ -189,10 +210,15 @@ inline void add_non_bonded_pair_force(

if (dist < ia_params.max_cut) {
#ifdef EXCLUSIONS
if (do_nonbonded(p1, p2))
if (do_nonbonded(p1, p2)) {
#endif
pf += calc_central_radial_force(p1, p2, ia_params, d, dist);
pf += calc_central_radial_charge_force(p1, p2, ia_params, d, dist,
coulomb_kernel);
pf += calc_non_central_force(p1, p2, ia_params, d, dist);
#ifdef EXCLUSIONS
}
#endif
pf += calc_non_bonded_pair_force(p1, p2, ia_params, d, dist,
coulomb_kernel);
}

/***********************************************/
Expand Down
7 changes: 5 additions & 2 deletions src/core/pressure_inline.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,11 @@ inline void add_non_bonded_pair_virials(
#endif
{
auto const &ia_params = get_ia_param(p1.type(), p2.type());
auto const force =
calc_non_bonded_pair_force(p1, p2, ia_params, d, dist, kernel_forces).f;
auto const force = calc_central_radial_force(p1, p2, ia_params, d, dist).f +
calc_central_radial_charge_force(p1, p2, ia_params, d,
dist, kernel_forces)
.f +
calc_non_central_force(p1, p2, ia_params, d, dist).f;
auto const stress = Utils::tensor_product(d, force);

auto const type1 = p1.mol_id();
Expand Down