From 7b78ce885c39713386ad24a0fddb8f599a8d8f9a Mon Sep 17 00:00:00 2001 From: bjornstenqvist Date: Mon, 22 Jun 2020 11:06:06 +0200 Subject: [PATCH 1/5] Add reciprocal ewald force --- src/energy.cpp | 57 ++++++++++++++++++++++++++++++++++++++------------ src/energy.h | 1 + 2 files changed, 45 insertions(+), 13 deletions(-) diff --git a/src/energy.cpp b/src/energy.cpp index a78a173f2..747f437f4 100644 --- a/src/energy.cpp +++ b/src/energy.cpp @@ -375,6 +375,37 @@ double Ewald::energy(Change &change) { return u; } +void Ewald::force(std::vector &force) { + assert(force.size() == spc.p.size()); + double volume = spc.geo.getVolume(); + auto force_iter = force.begin(); + *force_iter = {0.0, 0.0, 0.0}; + + // Surface force contribution + Point total_dipole_moment = {0.0, 0.0, 0.0}; + for (auto &p : spc.p) { + total_dipole_moment += p.pos * p.charge;// + dipoles[i]; + } + + for (auto &p : spc.p) { + for (size_t k = 0; k < data.k_vectors.cols(); k++) { + std::complex Q = data.Q_ion[k] + data.Q_dipole[k]; + double kDotR = data.k_vectors.col(k).dot(p.pos); + double coskDotR = std::cos(kDotR); + double sinkDotR = std::sin(kDotR); + std::complex expKri(coskDotR,sinkDotR); + std::complex qmu(0.0, p.charge); + std::complex repart = expKri * qmu * std::conj(Q); + (*force_iter) += std::real(repart) * data.k_vectors.col(k) * data.Aks[k]; + } + (*force_iter) += (total_dipole_moment * p.charge) / (2.0*data.surface_dielectric_constant + 1.0); + (*force_iter) *= -4.0 * pc::pi / volume; + + + force_iter++; + } +} + /** * @todo Implement a sync() function in EwaldData to selectively copy information */ @@ -382,9 +413,9 @@ void Ewald::sync(Energybase *energybase_pointer, Change &change) { auto other = dynamic_cast(energybase_pointer); assert(other); if (other->key == OLD) { - old_groups = - &(other->spc - .groups); // give NEW access to OLD space for optimized updates + old_groups = + &(other->spc + .groups); // give NEW access to OLD space for optimized updates } // hard-coded sync; should be expanded when dipolar ewald is supported @@ -579,7 +610,7 @@ double Bonded::energy(Change &change) { int offset = std::distance(spc.p.begin(), spc.groups[group.index].begin()); // add an offset to the group atom indices to get the absolute indices std::transform(group.atoms.begin(), group.atoms.end(), std::back_inserter(atoms_ndx), - [offset](int i) { return i + offset; }); + [offset](int i) { return i + offset; }); energy += sum_energy(intra_group, atoms_ndx); } } @@ -702,7 +733,7 @@ Hamiltonian::Hamiltonian(Space &spc, const json &j) { #ifdef ENABLE_MPI emplace_back(it.value(), spc); #else - emplace_back(it.value(), spc); + emplace_back(it.value(), spc); #endif #if defined ENABLE_FREESASA else if (it.key() == "sasa") @@ -771,19 +802,19 @@ SASAEnergy::SASAEnergy(Space &spc, double cosolute_concentration, double probe_r SASAEnergy::SASAEnergy(const json &j, Space &spc) : SASAEnergy(spc, j.value("molarity", 0.0) * 1.0_molar, j.value("radius", 1.4) * 1.0_angstrom) {} -void SASAEnergy::updatePositions([[gnu::unused]] const ParticleVector &p) { - assert(p.size() == spc.positions().size()); - positions.clear(); - for(auto pos: spc.positions()) { - auto xyz = pos.data(); - positions.insert(positions.end(), xyz, xyz+3); + void SASAEnergy::updatePositions([[gnu::unused]] const ParticleVector &p) { + assert(p.size() == spc.positions().size()); + positions.clear(); + for(auto pos: spc.positions()) { + auto xyz = pos.data(); + positions.insert(positions.end(), xyz, xyz+3); + } } -} void SASAEnergy::updateRadii(const ParticleVector &p) { radii.resize(p.size()); std::transform(p.begin(), p.end(), radii.begin(), - [](auto &a) { return atoms[a.id].sigma * 0.5; }); + [](auto &a) { return atoms[a.id].sigma * 0.5; }); } void SASAEnergy::updateSASA(const ParticleVector &p, const Change &) { diff --git a/src/energy.h b/src/energy.h index 8d9715e34..b6b79ccfc 100644 --- a/src/energy.h +++ b/src/energy.h @@ -205,6 +205,7 @@ class Ewald : public Energybase { Change &) override; //!< Called after a move is rejected/accepted //! as well as before simulation void to_json(json &) const override; + void force(std::vector &) override; // update forces on all particles }; class Isobaric : public Energybase { From 75da0c3b10d9bdff189c4ea76fe336a345322e34 Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Mon, 22 Jun 2020 15:02:25 +0200 Subject: [PATCH 2/5] Update --- src/energy.cpp | 52 +++++++++++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/src/energy.cpp b/src/energy.cpp index 747f437f4..e58850faf 100644 --- a/src/energy.cpp +++ b/src/energy.cpp @@ -375,34 +375,38 @@ double Ewald::energy(Change &change) { return u; } -void Ewald::force(std::vector &force) { - assert(force.size() == spc.p.size()); - double volume = spc.geo.getVolume(); - auto force_iter = force.begin(); - *force_iter = {0.0, 0.0, 0.0}; +/** + * @param forces Destination force vector + * + * Calculate forces from reciprocal space. Note that + * the destination force vector will *not* be zeroed + * before addition. + * + * @todo Dipole contribution incomplete + */ +void Ewald::force(std::vector &forces) { + assert(forces.size() == spc.p.size()); + const double volume = spc.geo.getVolume(); + auto force = forces.begin(); - // Surface force contribution + // Surface contribution Point total_dipole_moment = {0.0, 0.0, 0.0}; - for (auto &p : spc.p) { - total_dipole_moment += p.pos * p.charge;// + dipoles[i]; - } - - for (auto &p : spc.p) { - for (size_t k = 0; k < data.k_vectors.cols(); k++) { - std::complex Q = data.Q_ion[k] + data.Q_dipole[k]; - double kDotR = data.k_vectors.col(k).dot(p.pos); - double coskDotR = std::cos(kDotR); - double sinkDotR = std::sin(kDotR); - std::complex expKri(coskDotR,sinkDotR); - std::complex qmu(0.0, p.charge); + for (auto &particle : spc.p) { + total_dipole_moment += particle.pos * particle.charge; // + dipoles[i]; + } + + for (auto &particle : spc.p) { // loop over particles + (*force) = total_dipole_moment * particle.charge / (2.0 * data.surface_dielectric_constant + 1.0); + for (size_t i = 0; i < data.k_vectors.cols(); i++) { // loop over k vectors + std::complex Q = data.Q_ion[i] + data.Q_dipole[i]; + double k_dot_r = data.k_vectors.col(i).dot(particle.pos); + std::complex expKri(std::cos(k_dot_r), std::sin(k_dot_r)); + std::complex qmu(0.0, particle.charge); std::complex repart = expKri * qmu * std::conj(Q); - (*force_iter) += std::real(repart) * data.k_vectors.col(k) * data.Aks[k]; + (*force) += std::real(repart) * data.k_vectors.col(i) * data.Aks[i]; } - (*force_iter) += (total_dipole_moment * p.charge) / (2.0*data.surface_dielectric_constant + 1.0); - (*force_iter) *= -4.0 * pc::pi / volume; - - - force_iter++; + (*force) *= -4.0 * pc::pi / volume; + force++; } } From d5ff7002cada5dded659dfa2757eca2942173939 Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Mon, 22 Jun 2020 15:08:55 +0200 Subject: [PATCH 3/5] Revert indentation --- src/energy.cpp | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/energy.cpp b/src/energy.cpp index e58850faf..e6ea11531 100644 --- a/src/energy.cpp +++ b/src/energy.cpp @@ -417,9 +417,7 @@ void Ewald::sync(Energybase *energybase_pointer, Change &change) { auto other = dynamic_cast(energybase_pointer); assert(other); if (other->key == OLD) { - old_groups = - &(other->spc - .groups); // give NEW access to OLD space for optimized updates + old_groups = &(other->spc.groups); // give NEW access to OLD space for optimized updates } // hard-coded sync; should be expanded when dipolar ewald is supported @@ -614,7 +612,7 @@ double Bonded::energy(Change &change) { int offset = std::distance(spc.p.begin(), spc.groups[group.index].begin()); // add an offset to the group atom indices to get the absolute indices std::transform(group.atoms.begin(), group.atoms.end(), std::back_inserter(atoms_ndx), - [offset](int i) { return i + offset; }); + [offset](int i) { return i + offset; }); energy += sum_energy(intra_group, atoms_ndx); } } @@ -737,7 +735,7 @@ Hamiltonian::Hamiltonian(Space &spc, const json &j) { #ifdef ENABLE_MPI emplace_back(it.value(), spc); #else - emplace_back(it.value(), spc); + emplace_back(it.value(), spc); #endif #if defined ENABLE_FREESASA else if (it.key() == "sasa") @@ -806,19 +804,18 @@ SASAEnergy::SASAEnergy(Space &spc, double cosolute_concentration, double probe_r SASAEnergy::SASAEnergy(const json &j, Space &spc) : SASAEnergy(spc, j.value("molarity", 0.0) * 1.0_molar, j.value("radius", 1.4) * 1.0_angstrom) {} - void SASAEnergy::updatePositions([[gnu::unused]] const ParticleVector &p) { - assert(p.size() == spc.positions().size()); - positions.clear(); - for(auto pos: spc.positions()) { - auto xyz = pos.data(); - positions.insert(positions.end(), xyz, xyz+3); - } +void SASAEnergy::updatePositions([[gnu::unused]] const ParticleVector &p) { + assert(p.size() == spc.positions().size()); + positions.clear(); + for (auto pos : spc.positions()) { + auto xyz = pos.data(); + positions.insert(positions.end(), xyz, xyz + 3); } +} void SASAEnergy::updateRadii(const ParticleVector &p) { radii.resize(p.size()); - std::transform(p.begin(), p.end(), radii.begin(), - [](auto &a) { return atoms[a.id].sigma * 0.5; }); + std::transform(p.begin(), p.end(), radii.begin(), [](auto &a) { return atoms[a.id].sigma * 0.5; }); } void SASAEnergy::updateSASA(const ParticleVector &p, const Change &) { From fc132d10694a4ebdc9676bea9beef371c02d0ad7 Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Mon, 22 Jun 2020 15:14:47 +0200 Subject: [PATCH 4/5] Revert indentation --- src/energy.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/energy.cpp b/src/energy.cpp index e6ea11531..2de84805c 100644 --- a/src/energy.cpp +++ b/src/energy.cpp @@ -417,7 +417,9 @@ void Ewald::sync(Energybase *energybase_pointer, Change &change) { auto other = dynamic_cast(energybase_pointer); assert(other); if (other->key == OLD) { - old_groups = &(other->spc.groups); // give NEW access to OLD space for optimized updates + old_groups = + &(other->spc + .groups); // give NEW access to OLD space for optimized updates } // hard-coded sync; should be expanded when dipolar ewald is supported @@ -807,15 +809,16 @@ SASAEnergy::SASAEnergy(const json &j, Space &spc) void SASAEnergy::updatePositions([[gnu::unused]] const ParticleVector &p) { assert(p.size() == spc.positions().size()); positions.clear(); - for (auto pos : spc.positions()) { + for(auto pos: spc.positions()) { auto xyz = pos.data(); - positions.insert(positions.end(), xyz, xyz + 3); + positions.insert(positions.end(), xyz, xyz+3); } } void SASAEnergy::updateRadii(const ParticleVector &p) { radii.resize(p.size()); - std::transform(p.begin(), p.end(), radii.begin(), [](auto &a) { return atoms[a.id].sigma * 0.5; }); + std::transform(p.begin(), p.end(), radii.begin(), + [](auto &a) { return atoms[a.id].sigma * 0.5; }); } void SASAEnergy::updateSASA(const ParticleVector &p, const Change &) { From 7be8af067dedd404c3c590609751037ef8db3b01 Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Mon, 22 Jun 2020 15:41:20 +0200 Subject: [PATCH 5/5] Include dipoles --- src/energy.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/energy.cpp b/src/energy.cpp index 2de84805c..1bf591ce4 100644 --- a/src/energy.cpp +++ b/src/energy.cpp @@ -381,32 +381,33 @@ double Ewald::energy(Change &change) { * Calculate forces from reciprocal space. Note that * the destination force vector will *not* be zeroed * before addition. - * - * @todo Dipole contribution incomplete */ void Ewald::force(std::vector &forces) { assert(forces.size() == spc.p.size()); const double volume = spc.geo.getVolume(); - auto force = forces.begin(); // Surface contribution Point total_dipole_moment = {0.0, 0.0, 0.0}; for (auto &particle : spc.p) { - total_dipole_moment += particle.pos * particle.charge; // + dipoles[i]; + auto mu = particle.hasExtension() ? particle.getExt().mu * particle.getExt().mulen : Point(0, 0, 0); + total_dipole_moment += particle.pos * particle.charge + mu; } + auto force = forces.begin(); // iterator to force vector on first particle + for (auto &particle : spc.p) { // loop over particles (*force) = total_dipole_moment * particle.charge / (2.0 * data.surface_dielectric_constant + 1.0); + double mu_scalar = particle.hasExtension() ? particle.getExt().mulen : 0.0; + std::complex qmu(mu_scalar, particle.charge); for (size_t i = 0; i < data.k_vectors.cols(); i++) { // loop over k vectors std::complex Q = data.Q_ion[i] + data.Q_dipole[i]; double k_dot_r = data.k_vectors.col(i).dot(particle.pos); std::complex expKri(std::cos(k_dot_r), std::sin(k_dot_r)); - std::complex qmu(0.0, particle.charge); std::complex repart = expKri * qmu * std::conj(Q); (*force) += std::real(repart) * data.k_vectors.col(i) * data.Aks[i]; } - (*force) *= -4.0 * pc::pi / volume; - force++; + (*force) *= -4.0 * pc::pi / volume * data.bjerrum_length; // to units of kT/Angstrom^2 + force++; // advance to next force vector } }