diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index 7f542dee264..7dbc6cfd7b5 100644 --- a/doc/sphinx/particles.rst +++ b/doc/sphinx/particles.rst @@ -309,12 +309,11 @@ To switch the active scheme, the attribute :attr:`espressomd.system.System.virtu from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative system = espressomd.System() - system.virtual_sites = VirtualSitesRelative(have_velocity=True, have_quaternion=False) + system.virtual_sites = VirtualSitesRelative(have_quaternion=False) # or system.virtual_sites = VirtualSitesOff() By default, :class:`espressomd.virtual_sites.VirtualSitesOff` is selected. This means that virtual particles are not touched during integration. -The ``have_velocity`` parameter determines whether or not the velocity of virtual sites is calculated, which carries a performance cost. The ``have_quaternion`` parameter determines whether the quaternion of the virtual particle is updated (useful in combination with the :attr:`espressomd.particle_data.ParticleHandle.vs_quat` property of the virtual particle which defines the orientation of the virtual particle in the body fixed frame of the related real particle. diff --git a/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb b/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb index d10ff376b34..a4035b00950 100644 --- a/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb +++ b/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb @@ -294,7 +294,7 @@ "outputs": [], "source": [ "# Select the desired implementation for virtual sites\n", - "system.virtual_sites = VirtualSitesRelative(have_velocity=True)\n", + "system.virtual_sites = VirtualSitesRelative()\n", "# Setting min_global_cut is necessary when there is no interaction defined with a range larger than\n", "# the colloid such that the virtual particles are able to communicate their forces to the real particle\n", "# at the center of the colloid\n", diff --git a/samples/drude_bmimpf6.py b/samples/drude_bmimpf6.py index b0714d7c752..86cc66ddc9a 100644 --- a/samples/drude_bmimpf6.py +++ b/samples/drude_bmimpf6.py @@ -75,7 +75,7 @@ print("\n-->Ion pairs:", n_ionpairs, "Box size:", box_l) system = espressomd.System(box_l=[box_l, box_l, box_l]) -system.virtual_sites = VirtualSitesRelative(have_velocity=True) +system.virtual_sites = VirtualSitesRelative() if args.visu: d_scale = 0.988 * 0.5 diff --git a/samples/rigid_body.py b/samples/rigid_body.py index 841e0f771db..b3b70cd9ac2 100644 --- a/samples/rigid_body.py +++ b/samples/rigid_body.py @@ -32,8 +32,7 @@ system = espressomd.System(box_l=[10.0] * 3) -system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_velocity=True) +system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() system.time_step = 0.01 system.thermostat.set_langevin(kT=1.0, gamma=20.0, seed=42) diff --git a/src/core/event.cpp b/src/core/event.cpp index 1c202e3374a..5df85f482e2 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -423,7 +423,7 @@ unsigned global_ghost_flags() { void update_dependent_particles() { #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); cells_update_ghosts(global_ghost_flags()); #endif diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index e01fe2f0563..74f3fe72a4c 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -320,6 +320,9 @@ static bool is_poststorable(GhostCommunication const &ghost_comm) { } void ghost_communicator(GhostCommunicator *gcr, unsigned int data_parts) { + if (GHOSTTRANS_NONE == data_parts) + return; + static CommBuf send_buffer, recv_buffer; for (auto it = gcr->comm.begin(); it != gcr->comm.end(); ++it) { diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 9b39347faed..7989c8312c2 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -182,7 +182,7 @@ int integrate(int n_steps, int reuse_forces) { lb_lbcoupling_deactivate(); #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); #endif // Communication step: distribute ghost positions @@ -239,7 +239,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); #endif // Communication step: distribute ghost positions @@ -292,8 +292,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - // VIRTUAL_SITES update vel - virtual_sites()->update(false); // Recalc positions = false + virtual_sites()->update(); #endif /* verlet list statistics */ diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp index 902fe3c2505..dcf246398b4 100644 --- a/src/core/integrators/brownian_inline.hpp +++ b/src/core/integrators/brownian_inline.hpp @@ -507,11 +507,8 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, const ParticleRange &particles) { for (auto &p : particles) { // Don't propagate translational degrees of freedom of vs -#ifdef VIRTUAL_SITES extern bool thermo_virtual; - if (!(p.p.is_virtual) or thermo_virtual) -#endif - { + if (!(p.p.is_virtual) or thermo_virtual) { p.r.p += bd_drag(brownian.gamma, p, time_step); p.m.v = bd_drag_vel(brownian.gamma, p); p.r.p += bd_random_walk(brownian, p, time_step); diff --git a/src/core/integrators/velocity_verlet_inline.hpp b/src/core/integrators/velocity_verlet_inline.hpp index 7d51915f621..8dafc9a1e50 100644 --- a/src/core/integrators/velocity_verlet_inline.hpp +++ b/src/core/integrators/velocity_verlet_inline.hpp @@ -67,11 +67,10 @@ inline void velocity_verlet_propagate_vel_final(const ParticleRange &particles) { for (auto &p : particles) { -#ifdef VIRTUAL_SITES // Virtual sites are not propagated during integration if (p.p.is_virtual) continue; -#endif + for (int j = 0; j < 3; j++) { if (!(p.p.ext_flag & COORD_FIXED(j))) { /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */ diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 70b4c40e3a1..1b63fc88c9b 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -29,6 +29,7 @@ #include "npt.hpp" #include "pressure_inline.hpp" #include "virtual_sites.hpp" +#include #include "short_range_loop.hpp" @@ -127,8 +128,12 @@ void pressure_calc(double *result, double *result_t, double *result_nb, calc_long_range_virials(cell_structure.local_cells().particles()); #ifdef VIRTUAL_SITES - virtual_sites()->pressure_and_stress_tensor_contribution( - virials.virtual_sites, p_tensor.virtual_sites); + { + auto const vs_stress = virtual_sites()->stress_tensor(); + + *virials.virtual_sites += trace(vs_stress); + boost::copy(flatten(vs_stress), p_tensor.virtual_sites); + } #endif for (int n = 1; n < virials.data.n; n++) @@ -190,7 +195,7 @@ void init_virials(Observable_stat *stat) { Dipole::pressure_n(); #endif #ifdef VIRTUAL_SITES - n_vs = virtual_sites()->n_pressure_contribs(); + n_vs = 1; #endif // Allocate memory for the data @@ -228,7 +233,7 @@ void init_p_tensor(Observable_stat *stat) { auto const n_dipolar = 0; #endif #ifdef VIRTUAL_SITES - n_vs = virtual_sites()->n_pressure_contribs(); + n_vs = 1; #endif obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded, diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 4d6a3d9c1ba..fb4a8ee0a58 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -73,17 +73,17 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, } #ifdef ELECTROSTATICS - /* real space Coulomb */ - auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist); + { + /* real space Coulomb */ + auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist); - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; + } } - } - for (int i = 0; i < 3; i++) { - virials.coulomb[0] += p_coulomb[i][i]; + virials.coulomb[0] += trace(p_coulomb); } #endif /*ifdef ELECTROSTATICS */ diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index 8aab54a63ad..3bb3dd719a7 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -26,6 +26,8 @@ #include "virtual_sites/VirtualSites.hpp" +#include + /** @brief get active virtual sites implementation */ const std::shared_ptr &virtual_sites(); diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index 454b416af09..9e63d826668 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -33,42 +33,33 @@ */ #ifdef VIRTUAL_SITES -#include +#include /** @brief Base class for virtual sites implementations */ class VirtualSites { public: - VirtualSites() : m_have_velocity(true), m_have_quaternion(false){}; - /** @brief Update positions and/or velocities of virtual sites. - * Velocities are only updated if get_have_velocity() returns true. - * @param recalc_positions Skip the recalculation of positions if false. + VirtualSites() = default; + virtual ~VirtualSites() = default; + + /** + * @brief Update positions and velocities of virtual sites. */ - virtual void update(bool recalc_positions) const = 0; + virtual void update() const {} /** Back-transfer forces (and torques) to non-virtual particles. */ - virtual void back_transfer_forces_and_torques() const = 0; + virtual void back_transfer_forces_and_torques() const {} /** @brief Called after force calculation (and before rattle/shake) */ virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; - /** @brief Number of pressure contributions */ - virtual int n_pressure_contribs() const { return 0; }; /** @brief Pressure contribution. */ - virtual void - pressure_and_stress_tensor_contribution(double *pressure, - double *stress_tensor) const {}; - /** @brief Enable/disable velocity calculations for vs. */ - void set_have_velocity(const bool &v) { m_have_velocity = v; }; - const bool &get_have_velocity() const { return m_have_velocity; }; + virtual Utils::Matrix stress_tensor() const { return {}; }; /** @brief Enable/disable quaternion calculations for vs.*/ void set_have_quaternion(const bool &have_quaternion) { m_have_quaternion = have_quaternion; }; bool get_have_quaternion() const { return m_have_quaternion; }; - virtual ~VirtualSites() = default; - private: - bool m_have_velocity; - bool m_have_quaternion; + bool m_have_quaternion = false; }; #endif diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp index 3eb8973a5a4..c22a66687f1 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp @@ -28,8 +28,6 @@ * instantaneously transferred to the fluid */ class VirtualSitesInertialessTracers : public VirtualSites { - void update(bool recalc_positions) const override{}; - void back_transfer_forces_and_torques() const override{}; void after_force_calc() override; void after_lb_propagation() override; }; diff --git a/src/core/virtual_sites/VirtualSitesOff.hpp b/src/core/virtual_sites/VirtualSitesOff.hpp index ba4ddc8959b..75918ba616f 100644 --- a/src/core/virtual_sites/VirtualSitesOff.hpp +++ b/src/core/virtual_sites/VirtualSitesOff.hpp @@ -24,10 +24,7 @@ #include "VirtualSites.hpp" /** @brief Do-nothing virtual-sites implementation */ -class VirtualSitesOff : public VirtualSites { - void update(bool) const override{}; - void back_transfer_forces_and_torques() const override{}; -}; +class VirtualSitesOff : public VirtualSites {}; #endif #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 5dc59824aa8..660e404e718 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -18,72 +18,41 @@ */ #include "VirtualSitesRelative.hpp" -#include "forces_inline.hpp" -#include -#include #ifdef VIRTUAL_SITES_RELATIVE #include "cells.hpp" -#include "config.hpp" -#include "errorhandling.hpp" +#include "forces_inline.hpp" #include "grid.hpp" #include "integrate.hpp" +#include "particle_data.hpp" #include "rotation.hpp" -void VirtualSitesRelative::update(bool recalc_positions) const { - // Ghost update logic - if (n_nodes > 1) { - auto const data_parts = - (recalc_positions ? GHOSTTRANS_POSITION : GHOSTTRANS_NONE) | - (get_have_velocity() ? (GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM) - : GHOSTTRANS_NONE); - - if (recalc_positions or get_have_velocity()) { - ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); - } - } - for (auto &p : cell_structure.local_cells().particles()) { - if (!p.p.is_virtual) - continue; - - if (recalc_positions) - update_pos(p); - - if (get_have_velocity()) - update_vel(p); - - if (get_have_quaternion()) - update_virtual_particle_quaternion(p); - } -} +#include +#include -void VirtualSitesRelative::update_virtual_particle_quaternion( - Particle &p) const { - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - if (!p_real) { - throw std::runtime_error( - "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"); - } - p.r.quat = multiply_quaternions(p_real->r.quat, p.p.vs_relative.quat); +namespace { +/** + * @brief Orientation of the virtual site. + * @param p_ref Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Orientation quaternion of the virtual site. + */ +Utils::Vector4d +orientation(Particle const *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return multiply_quaternions(p_ref->r.quat, vs_rel.quat); } -void VirtualSitesRelative::update_pos(Particle &p) const { - // First obtain the real particle responsible for this virtual particle: - // Find the 1st real particle in the topology for the virtual particle's - // mol_id - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - // Check, if a real particle was found - if (!p_real) { - runtimeErrorMsg() - << "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"; - return; - } - +/** + * @brief Vector pointing from the real particle + * to the virtual site. + * + * @return Relative distance. + */ +Utils::Vector3d connection_vector( + Particle const *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { // Calculate the quaternion defining the orientation of the vector connecting // the virtual site and the real particle // This is obtained, by multiplying the quaternion representing the director @@ -91,40 +60,116 @@ void VirtualSitesRelative::update_pos(Particle &p) const { // specifies the relative orientation. auto const director = Utils::convert_quaternion_to_director( - Utils::multiply_quaternions(p_real->r.quat, - p.p.vs_relative.rel_orientation)) + Utils::multiply_quaternions(p_ref->r.quat, vs_rel.rel_orientation)) .normalize(); - auto const new_pos = p_real->r.p + director * p.p.vs_relative.distance; - /* The shift has to respect periodic boundaries: if the reference particles - * is not in the same image box, we potentially avoid to shift to the other - * side of the box. */ - auto const shift = get_mi_vector(new_pos, p.r.p, box_geo); - p.r.p += shift; - - if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) - set_resort_particles(Cells::RESORT_LOCAL); + return vs_rel.distance * director; } -void VirtualSitesRelative::update_vel(Particle &p) const { - // First obtain the real particle responsible for this virtual particle: - Particle *p_real = get_local_particle_data(p.p.vs_relative.to_particle_id); - // Check, if a real particle was found - if (!p_real) { - runtimeErrorMsg() - << "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"; - return; - } +/** + * @brief Position of the virtual site + * @param p_ref Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Position of the virtual site. + */ +Utils::Vector3d +position(Particle const *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return p_ref->r.p + connection_vector(p_ref, vs_rel); +} - auto const d = get_mi_vector(p.r.p, p_real->r.p, box_geo); +/** + * @brief Velocity of the virtual site + * @param p_ref Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Velocity of the virtual site. + */ +Utils::Vector3d +velocity(const Particle *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + auto const d = connection_vector(p_ref, vs_rel); // Get omega of real particle in space-fixed frame - Utils::Vector3d omega_space_frame = - convert_vector_body_to_space(*p_real, p_real->m.omega); + auto const omega_space_frame = + convert_vector_body_to_space(*p_ref, p_ref->m.omega); // Obtain velocity from v=v_real particle + omega_real_particle \times // director - p.m.v = vector_product(omega_space_frame, d) + p_real->m.v; + return vector_product(omega_space_frame, d) + p_ref->m.v; +} + +/** + * @brief Get reference particle. + * + * @param vs_rel Parameters to get the reference particle for. + * @return Pointer to reference particle. + */ +Particle *get_reference_particle( + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + auto p_ref = get_local_particle_data(vs_rel.to_particle_id); + if (!p_ref) { + throw std::runtime_error("No real particle associated with virtual site."); + } + + return p_ref; +} + +/** + * @brief Contraint force on the real particle. + * + * Calculates the force exerted by the constraint on the + * reference particle. + */ +ParticleForce constraint_force( + const ParticleForce &f, const Particle *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return {f.f, + vector_product(connection_vector(p_ref, vs_rel), f.f) + f.torque}; +} + +/** + * @brief Constraint force to hold the particles at its prescribed position. + * + * @param f Force on the virtual site. + * @param p_ref Reference particle. + * @param vs_rel Parameters. + * @return Constraint force. + */ +auto constraint_stress( + const Utils::Vector3d &f, const Particle *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + /* The constraint force is minus the force on the particle, make it force + * free. The counter force is translated by the connection vector to the + * real particle, hence the virial stress is */ + return tensor_product(connection_vector(p_ref, vs_rel), -f); +} +} // namespace + +void VirtualSitesRelative::update() const { + // Ghost update logic + auto const data_parts = GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM; + + ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); + + for (auto &p : cell_structure.local_cells().particles()) { + if (!p.p.is_virtual) + continue; + + const Particle *p_ref = get_reference_particle(p.p.vs_relative); + + auto const new_pos = position(p_ref, p.p.vs_relative); + /* The shift has to respect periodic boundaries: if the reference + * particles is not in the same image box, we potentially avoid to shift + * to the other side of the box. */ + p.r.p += get_mi_vector(new_pos, p.r.p, box_geo); + + p.m.v = velocity(p_ref, p.p.vs_relative); + + if (get_have_quaternion()) + p.r.quat = orientation(p_ref, p.p.vs_relative); + + if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) + set_resort_particles(Cells::RESORT_LOCAL); + } // namespace } // Distribute forces that have accumulated on virtual particles to the @@ -139,53 +184,28 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // We only care about virtual particles if (p.p.is_virtual) { // First obtain the real particle responsible for this virtual particle: - Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - - // The rules for transferring forces are: - // F_realParticle +=F_virtualParticle - // T_realParticle +=f_realParticle \times - // (r_virtualParticle-r_realParticle) + Particle *p_ref = get_reference_particle(p.p.vs_relative); // Add forces and torques - p_real->f.torque += - vector_product(get_mi_vector(p.r.p, p_real->r.p, box_geo), p.f.f) + - p.f.torque; - p_real->f.f += p.f.f; + p_ref->f += constraint_force(p.f, p_ref, p.p.vs_relative); } } } // Rigid body contribution to scalar pressure and stress tensor -void VirtualSitesRelative::pressure_and_stress_tensor_contribution( - double *pressure, double *stress_tensor) const { - // Division by 3 volume is somewhere else. (pressure.cpp after all pressure - // calculations) Iterate over all the particles in the local cells +Utils::Matrix VirtualSitesRelative::stress_tensor() const { + Utils::Matrix stress_tensor = {}; for (auto &p : cell_structure.local_cells().particles()) { if (!p.p.is_virtual) continue; - update_pos(p); - // First obtain the real particle responsible for this virtual particle: - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - - // Get distance vector pointing from real to virtual particle, respecting - // periodic boundary i - // conditions - auto const d = get_mi_vector(p_real->r.p, p.r.p, box_geo); - - // Stress tensor contribution - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - stress_tensor[k * 3 + l] += p.f.f[k] * d[l]; - - // Pressure = 1/3 trace of stress tensor - // but the 1/3 is applied somewhere else. - *pressure += p.f.f * d; + const Particle *p_ref = get_reference_particle(p.p.vs_relative); + + stress_tensor += constraint_stress(p.f.f, p_ref, p.p.vs_relative); } -} + return stress_tensor; +} #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 975cecee5de..cfe4e512928 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -23,38 +23,20 @@ #include "config.hpp" #ifdef VIRTUAL_SITES_RELATIVE -#include "Particle.hpp" -#include "communication.hpp" -#include "virtual_sites.hpp" +#include "VirtualSites.hpp" + +#include /** @brief Virtual sites implementation for rigid bodies */ class VirtualSitesRelative : public VirtualSites { public: VirtualSitesRelative() = default; /** @copydoc VirtualSites::update */ - void update(bool recalc_positions) const override; + void update() const override; /** @copydoc VirtualSites::back_transfer_forces_and_torques */ void back_transfer_forces_and_torques() const override; - /** @copydoc VirtualSites::n_pressure_contribs */ - int n_pressure_contribs() const override { return 1; }; - /** @copydoc VirtualSites::pressure_and_stress_tensor_contribution */ - void - pressure_and_stress_tensor_contribution(double *pressure, - double *stress_tensor) const override; - -private: - /** Update the position of the given virtual particle as defined by the real - * particles in the same molecule - */ - void update_pos(Particle &p) const; - /** Update the velocity of the given virtual particle as defined by the real - * particles in the same molecule - */ - void update_vel(Particle &p) const; - /** @brief Update the orientation of the virtual particles with respect - * to the real particle. - */ - void update_virtual_particle_quaternion(Particle &p) const; + /** @copydoc VirtualSites::stress_tensor */ + Utils::Matrix stress_tensor() const override; }; #endif diff --git a/src/python/espressomd/virtual_sites.py b/src/python/espressomd/virtual_sites.py index ac2e721fe8f..2987a054d8a 100644 --- a/src/python/espressomd/virtual_sites.py +++ b/src/python/espressomd/virtual_sites.py @@ -58,14 +58,5 @@ class VirtualSitesRelative(ScriptInterfaceHelper): """Virtual sites implementation placing virtual sites relative to other particles. See :ref:`Rigid arrangements of particles` for details. - Attributes can be set on the instance or passed to the constructor as - keyword arguments. - - Attributes - ---------- - have_velocity : :obj:`bool` - Determines whether the velocity of the virtual sites is calculated. - This carries a performance cost. - """ _so_name = "VirtualSites::VirtualSitesRelative" diff --git a/src/script_interface/virtual_sites/VirtualSites.hpp b/src/script_interface/virtual_sites/VirtualSites.hpp index 5bccbc3bd80..228b4ffff82 100644 --- a/src/script_interface/virtual_sites/VirtualSites.hpp +++ b/src/script_interface/virtual_sites/VirtualSites.hpp @@ -34,12 +34,7 @@ class VirtualSites : public AutoParameters { public: VirtualSites() { add_parameters( - {{"have_velocity", - [this](const Variant &v) { - virtual_sites()->set_have_velocity(get_value(v)); - }, - [this]() { return virtual_sites()->get_have_velocity(); }}, - {"have_quaternion", + {{"have_quaternion", [this](const Variant &v) { virtual_sites()->set_have_quaternion(get_value(v)); }, diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index 8e0467b09a9..bdebb05c061 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -155,6 +155,42 @@ using Vector3i = VectorXi<3>; template using Matrix = Vector, N>; +/** + * @brief Trace of a matrix. + * + * Returns the sum of the diagonal elements + * of a square matrix. + * + * @tparam T Arithmetic type + * @tparam N Matrix dimension + * @param m Input matrix + * @return Trace of matrix. + */ +template T trace(Matrix const &m) { + auto tr = T{}; + for (size_t i = 0; i < N; i++) + tr += m[i][i]; + + return tr; +} + +/** + * @brief Flatten a matrix to a linear vector. + * + * @param m Input Matrix + * @return Flat vector with elements of the matrix. + */ +template +Vector flatten(Matrix const &m) { + Vector ret; + + for (size_t i = 0; i < N; i++) + for (size_t j = 0; j < M; j++) + ret[i * M + j] = m[j][i]; + + return ret; +} + namespace detail { template auto binary_op(Vector const &a, Vector const &b, Op op) { diff --git a/src/utils/tests/Vector_test.cpp b/src/utils/tests/Vector_test.cpp index 421bad9f825..689faaff662 100644 --- a/src/utils/tests/Vector_test.cpp +++ b/src/utils/tests/Vector_test.cpp @@ -381,3 +381,19 @@ BOOST_AUTO_TEST_CASE(diag_matrix) { for (int j = 0; j < 3; j++) BOOST_CHECK_EQUAL(result[i][j], (i == j) ? v[i] : 0); } + +BOOST_AUTO_TEST_CASE(trace_) { + auto const A = Utils::Matrix{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}; + auto const result = trace(A); + auto const expected = A[0][0] + A[1][1] + A[2][2]; + + BOOST_CHECK_EQUAL(expected, result); +} + +BOOST_AUTO_TEST_CASE(flatten_) { + auto const A = Utils::Matrix{{1, 2}, {3, 4}}; + auto const result = flatten(A); + auto const expected = Utils::Vector{1, 3, 2, 4}; + + BOOST_CHECK(result == expected); +} diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index c2c9b7253ab..7d5c524af99 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -161,7 +161,7 @@ if espressomd.has_features(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']): system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_velocity=True, have_quaternion=True) + have_quaternion=True) system.part[1].vs_auto_relate_to(0) if espressomd.has_features(['LENNARD_JONES']) and 'LJ' in modes: diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index b7a125c7ed7..35c72e6f038 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -77,9 +77,8 @@ def test_aa_method_switching(self): self.assertIsInstance(self.system.virtual_sites, VirtualSitesOff) # Switch implementation - self.system.virtual_sites = VirtualSitesRelative(have_velocity=False) + self.system.virtual_sites = VirtualSitesRelative() self.assertIsInstance(self.system.virtual_sites, VirtualSitesRelative) - self.assertFalse(self.system.virtual_sites.have_velocity) def test_vs_quat(self): self.system.part.clear() @@ -128,7 +127,7 @@ def test_vs_quat(self): def test_pos_vel_forces(self): system = self.system system.cell_system.skin = 0.3 - system.virtual_sites = VirtualSitesRelative(have_velocity=True) + system.virtual_sites = VirtualSitesRelative() system.box_l = [10, 10, 10] system.part.clear() system.time_step = 0.004 @@ -203,21 +202,12 @@ def test_pos_vel_forces(self): # Check self.assertLessEqual(np.linalg.norm(t_exp - t), 1E-6) - # Check virtual sites without velocity - system.virtual_sites.have_velocity = False - - # Velocity should not change - v2 = np.copy(system.part[2].v) - system.part[1].v = [17, -13.5, 2] - system.integrator.run(0, recalc_forces=True) - np.testing.assert_array_equal(v2, np.copy(system.part[2].v)) - def run_test_lj(self): """This fills the system with vs-based dumbells, adds a lj potential, integrates and verifies forces. This is to make sure that no pairs get lost or are outdated in the short range loop""" system = self.system - system.virtual_sites = VirtualSitesRelative(have_velocity=True) + system.virtual_sites = VirtualSitesRelative() # Parameters n = 90 phi = 0.6 @@ -320,10 +310,12 @@ def test_zz_stress_tensor(self): system.time_step = 0.01 system.cell_system.skin = 0.1 system.min_global_cut = 0.2 - # Should not have one if vs are turned off + # Should not have a pressure system.virtual_sites = VirtualSitesOff() - self.assertNotIn("virtual_sites", system.analysis.pressure()) - self.assertNotIn("virtual_sites", system.analysis.stress_tensor()) + stress_vs = system.analysis.stress_tensor()["virtual_sites", 0] + p_vs = system.analysis.pressure()["virtual_sites", 0] + np.testing.assert_allclose(stress_vs, 0., atol=1e-10) + np.testing.assert_allclose(p_vs, 0., atol=1e-10) # vs relative contrib system.virtual_sites = VirtualSitesRelative() diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index 4ebab12362a..9328e7cf179 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -64,7 +64,6 @@ def test_aa_method_switching(self): self.system.virtual_sites = VirtualSitesInertialessTracers() self.assertIsInstance( self.system.virtual_sites, VirtualSitesInertialessTracers) - self.assertTrue(self.system.virtual_sites.have_velocity) def test_advection(self): self.reset_lb(ext_force_density=[0.1, 0, 0])