-
Notifications
You must be signed in to change notification settings - Fork 34
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.
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
Added dipole-dipole interactions. #136
Conversation
if (type == "qpotential") | ||
selfenergy_prefactor = 0.5; | ||
|
||
selfenergy_ion_prefactor = 0.0; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
document in docs/energy.md
src/energy.cpp
Outdated
@@ -344,6 +397,7 @@ Hamiltonian::Hamiltonian(Space &spc, const json &j) { | |||
using namespace Potential; | |||
|
|||
typedef CombinedPairPotential<CoulombGalore, LennardJones> CoulombLJ; | |||
typedef CombinedPairPotential<DipoleDipoleGalore, LennardJones> Stockmayer; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is this needed? Just use nonbonded_exact
?
src/multipole.h
Outdated
/** | ||
* @brief Help-function for `sfQpotential` using order 3. | ||
*/ | ||
inline void _DipoleDipoleQ2Help_3(double &a3, double &b3, double q) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
function names should be "camelCase", i.e. dipoleDipoleQ2Help3
template <class Titer> double monopoleMoment(Titer begin, Titer end) { | ||
double z = 0; | ||
for (auto it = begin; it != end; ++it) | ||
z += it->charge; | ||
return z; | ||
} //!< Calculates dipole moment vector for a set of particles | ||
|
||
/** |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
nice!
@@ -82,21 +82,21 @@ void Faunus::Potential::CoulombGalore::sfQpotential(const Faunus::json &j) { | |||
order = j.value("order", 300); | |||
table = sf.generate([&](double q) { return qPochhammerSymbol(q, 1, order); }, 0, 1); | |||
calcDielectric = [&](double M2V) { return 1 + 3 * M2V; }; | |||
selfenergy_prefactor = 0.5; | |||
selfenergy_prefactor = -1.0; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
document self energies in md files?
This PR relates to issue #135. |
@@ -58,7 +58,7 @@ void EwaldData::update(const Point &box) { | |||
factor *= 2; | |||
double dkz2 = double(kz * kz); | |||
Point kv = 2 * pc::pi * Point(kx / L.x(), ky / L.y(), kz / L.z()); | |||
double k2 = kv.dot(kv); | |||
double k2 = kv.dot(kv) + kappa2; // last term is only for Yukawa-Ewald |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
where is Yuwawa ewald?
src/potentials.cpp
Outdated
} | ||
|
||
catch (std::exception &e) { | ||
std::cerr << "DipoleDipoleGalore error: " << e.what(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
use instead throw std::runtime_error(name + ": " + e.what());
src/potentials.cpp
Outdated
j["type"] = type; | ||
if (type == "wolf" || type == "fennell" || type == "ewald") | ||
j["alpha"] = alpha; | ||
if (type == "qpotential" || type == "q2potential") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
elseif
Particle a = atoms[0]; | ||
Particle b = atoms[1]; | ||
Point r = {2, 0, 0}; | ||
CHECK(u(a, a, r) == Approx(dipoledipole(a, a, r))); // interaction between two parallell dipoles, directed parallell to their seperation |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this checks if FunctorPotential gives same results as dipoledipole. Add test that compares with precalculated values? (probably for some less idealized dipole orientations)
src/potentials.h
Outdated
using doctest::Approx; | ||
|
||
json j = R"({ "atomlist" : [ | ||
{"A": { "mu":[3.0,0.0,0.0] }}, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
according to documentation, mu
should be a unit vector
Please check that:
make all pyfaunus usagetips unittests test
passes with no errorsthere are no compiler warnings
no temporary files (.o, .pqr, state etc) are included in the commit
naming policies for classes, functions, and variables adhere to:
source files use soft-indendation of four. To setup automatic
code formatting of proposed changes using
clang-format
, installthe provided git commit hook: