diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 622d9211a22..bf4add25617 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -1,6 +1,6 @@ file(GLOB EspressoCore_SRC - "*.cpp" - ) + "*.cpp" "virtual_sites/VirtualSitesRelative.cpp") + if( WITH_COVERAGE ) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -Og --coverage -fprofile-arcs -ftest-coverage") diff --git a/src/core/collision.cpp b/src/core/collision.cpp index c0e509f79c5..96cad48bb83 100755 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -26,7 +26,7 @@ #include "grid.hpp" #include "domain_decomposition.hpp" #include "rotation.hpp" -#include "virtual_sites_relative.hpp" +#include "virtual_sites/VirtualSitesRelative.hpp" using namespace std; diff --git a/src/core/collision.hpp b/src/core/collision.hpp index 1a415670a18..ce8e3cbad62 100755 --- a/src/core/collision.hpp +++ b/src/core/collision.hpp @@ -21,7 +21,6 @@ #include "particle_data.hpp" #include "interaction_data.hpp" -#include "virtual_sites_relative.hpp" #include "virtual_sites.hpp" #include "integrate.hpp" diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 0feda04fc9c..8bc9f6e89d7 100755 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -103,17 +103,7 @@ void force_calc() { } #endif -#if defined(VIRTUAL_SITES_RELATIVE) && defined(LB) - // This is on a workaround stage: - // When using virtual sites relative and LB at the same time, it is necessary - // to reassemble the cell lists after all position updates, also of virtual - // particles. - if ((lattice_switch & LATTICE_LB) && - cell_structure.type == CELL_STRUCTURE_DOMDEC && (!cell_structure.use_verlet_list) && - (std::dynamic_pointer_cast(virtual_sites()) != nullptr)) - cells_update_ghosts(); -#endif - + espressoSystemInterface.update(); #ifdef COLLISION_DETECTION diff --git a/src/core/initialize.cpp b/src/core/initialize.cpp index 9400a252917..58dbe5cd5f9 100755 --- a/src/core/initialize.cpp +++ b/src/core/initialize.cpp @@ -741,6 +741,4 @@ void on_ghost_flags_change() { }; #endif -// if (old_have_v != ghosts_have_v) -// cells_re_init(CELL_STRUCTURE_CURRENT); } diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index ce8a8afe4e4..159be8e7d4e 100755 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -26,7 +26,7 @@ #include "cells.hpp" #include "integrate.hpp" #include "initialize.hpp" -#include "virtual_sites_relative.hpp" +#include "virtual_sites.hpp" #include "npt.hpp" #include "p3m.hpp" #include "p3m-dipolar.hpp" @@ -122,11 +122,8 @@ void pressure_calc(double *result, double *result_t, double *result_nb, double * calc_long_range_virials(); -#ifdef VIRTUAL_SITES_RELATIVE - std::shared_ptr vsr =std::dynamic_pointer_cast(virtual_sites()); - if (vsr) { - vsr->pressure_and_stress_tensor_contribution(virials.vs_relative,p_tensor.vs_relative); - } +#ifdef VIRTUAL_SITES + virtual_sites()->pressure_and_stress_tensor_contribution(virials.virtual_sites,p_tensor.virtual_sites); #endif @@ -235,14 +232,14 @@ void init_virials(Observable_stat *stat) { // Determine number of contribution for different interaction types // bonded, nonbonded, coulomb, dipolar, rigid bodies - int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vsr; + int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vs; n_pre = 1; n_non_bonded = (n_particle_types*(n_particle_types+1))/2; n_coulomb = 0; n_dipolar = 0; - n_vsr=0; + n_vs=0; #ifdef ELECTROSTATICS switch (coulomb.method) { @@ -263,14 +260,13 @@ void init_virials(Observable_stat *stat) break; } #endif -#ifdef VIRTUAL_SITES_RELATIVE - // rigid bodies - n_vsr=1; +#ifdef VIRTUAL_SITES + n_vs=virtual_sites()->n_pressure_contribs(); #endif // Allocate memory for the data - obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vsr, 1); + obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vs, 1); stat->init_status = 0; } @@ -291,7 +287,7 @@ void init_p_tensor(Observable_stat *stat) { // Determine number of contribution for different interaction types // bonded, nonbonded, coulomb, dipolar, rigid bodies - int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vsr; + int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vs; n_pre = 1; @@ -299,7 +295,7 @@ void init_p_tensor(Observable_stat *stat) n_coulomb = 0; n_dipolar = 0; - n_vsr=0; + n_vs=0; #ifdef ELECTROSTATICS switch (coulomb.method) { @@ -319,12 +315,11 @@ void init_p_tensor(Observable_stat *stat) default: n_dipolar = 0; } #endif -#ifdef VIRTUAL_SITES_RELATIVE - // rigid bodies - n_vsr=1; +#ifdef VIRTUAL_SITES + n_vs=virtual_sites()->n_pressure_contribs(); #endif - obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vsr, 9); + obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vs, 9); stat->init_status = 0; } diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index b8e3dff03c3..a722925a018 100755 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -1184,12 +1184,12 @@ void analyze_activate(PartCfg &partCfg, int ind) { void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded, int n_non_bonded, int n_coulomb, int n_dipolar, - int n_vsr, int c_size) { + int n_vs, int c_size) { int i; // Number of doubles to store pressure in int total = c_size * (n_pre + n_bonded_ia + n_non_bonded + n_coulomb + - n_dipolar + n_vsr); + n_dipolar + n_vs); // Allocate mem for the double list stat->data.resize(total); @@ -1201,13 +1201,13 @@ void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded, stat->n_coulomb = n_coulomb; stat->n_dipolar = n_dipolar; stat->n_non_bonded = n_non_bonded; - stat->n_vs_relative = n_vsr; // virtual sites relative (rigid bodies) + stat->n_virtual_sites = n_vs; // Pointers to the start of different contributions stat->bonded = stat->data.e + c_size * n_pre; stat->non_bonded = stat->bonded + c_size * n_bonded_ia; stat->coulomb = stat->non_bonded + c_size * n_non_bonded; stat->dipolar = stat->coulomb + c_size * n_coulomb; - stat->vs_relative = stat->dipolar + c_size * n_dipolar; + stat->virtual_sites = stat->dipolar + c_size * n_dipolar; // Set all obseravables to zero for (i = 0; i < total; i++) diff --git a/src/core/statistics.hpp b/src/core/statistics.hpp index 0318235e432..0e682725d50 100755 --- a/src/core/statistics.hpp +++ b/src/core/statistics.hpp @@ -57,7 +57,7 @@ typedef struct { /** number of non bonded interactions */ int n_non_bonded; /** Number of virtual sites relative (rigid body) conributions */ - int n_vs_relative; + int n_virtual_sites; /** start of bonded interactions. Right after the special ones */ double *bonded; @@ -68,7 +68,7 @@ typedef struct { /** start of observables for coulomb interaction. */ double *dipolar; /** Start of observables for virtual sites relative (rigid bodies) */ - double *vs_relative; + double *virtual_sites; /** number of doubles per data item */ int chunk_size; @@ -437,7 +437,7 @@ void invalidate_obs(); void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded, int n_non_bonded, int n_coulomb, int n_dipolar, - int n_vsr, int chunk_size); + int n_vs, int chunk_size); void obsstat_realloc_and_clear_non_bonded(Observable_stat_non_bonded *stat_nb, int n_nonbonded, int chunk_size_nb); diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index f98b6004e22..4ae9f23e06e 100755 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -24,6 +24,7 @@ #include "initialize.hpp" #include "statistics.hpp" #include "integrate.hpp" +#include "rotation.hpp" #ifdef VIRTUAL_SITES namespace { @@ -41,4 +42,111 @@ void set_virtual_sites(std::shared_ptr v) { on_ghost_flags_change(); } +#ifdef VIRTUAL_SITES_RELATIVE +int vs_relate_to(int part_num, int relate_to) +{ + // Get the data for the particle we act on and the one we wnat to relate + // it to. + auto p_current = get_particle_data(part_num); + auto p_relate_to = get_particle_data(relate_to); + if (!p_current || !p_relate_to) { + runtimeErrorMsg() <<"Could not retrieve particle data for the given id"; + return ES_ERROR; + } + + // get the distance between the particles + double d[3]; + get_mi_vector(d, p_current->r.p,p_relate_to->r.p); + + + + // Check, if the distance between virtual and non-virtual particles is larger htan minimum global cutoff + // If so, warn user + double l=sqrt(sqrlen(d)); + if (l>min_global_cut) { + runtimeErrorMsg() << "Warning: The distance between virtual and non-virtual particle (" << l << ") is\nlarger than the minimum global cutoff (" << min_global_cut << "). This may lead to incorrect simulations\nunder certain conditions. Set the \"System()\" class property \"min_global_cut\" to\nincrease the minimum cutoff.\n"; + return ES_ERROR; + } + + // Now, calculate the quaternions which specify the angle between + // the director of the particel we relate to and the vector + // (paritlce_we_relate_to - this_particle) + double quat[4]; + // The vs_relative implemnation later obtains the direcotr by multiplying + // the quaternions representing the orientation of the real particle + // with those in the virtual particle. The re quulting quaternion is then + // converted to a director. + // Whe have quat_(real particle) *quat_(virtual particle) + // = quat_(obtained from desired director) + // Resolving this for the quat_(virtaul particle) + + //Normalize desired director + int i; + + // If the distance between real & virtual particle is 0 + // we just set the relative orientation to 1 0 0 0, as it is irrelevant but + // needs to be a valid quaternion + if (l!=0) + { + for (i=0;i<3;i++) + d[i]/=l; + + // Obtain quaternions from desired director + double quat_director[4]; + convert_quatu_to_quat(d, quat_director); + + // Define quat as described above: + double x=0; + for (i=0;i<4;i++) + x+=p_relate_to->r.quat[i]*p_relate_to->r.quat[i]; + + quat[0]=0; + for (i=0;i<4;i++) + quat[0] +=p_relate_to->r.quat[i]*quat_director[i]; + + quat[1] =-quat_director[0] *p_relate_to->r.quat[1] + +quat_director[1] *p_relate_to->r.quat[0] + +quat_director[2] *p_relate_to->r.quat[3] + -quat_director[3] *p_relate_to->r.quat[2]; + quat[2] =p_relate_to->r.quat[1] *quat_director[3] + + p_relate_to->r.quat[0] *quat_director[2] + - p_relate_to->r.quat[3] *quat_director[1] + - p_relate_to->r.quat[2] * quat_director[0]; + quat[3] =quat_director[3] *p_relate_to->r.quat[0] + - p_relate_to->r.quat[3] *quat_director[0] + + p_relate_to->r.quat[2] * quat_director[1] + - p_relate_to->r.quat[1] *quat_director[2]; + for (i=0;i<4;i++) + quat[i]/=x; + + + // Verify result + double qtemp[4]; + multiply_quaternions(p_relate_to->r.quat,quat,qtemp); + for (i=0;i<4;i++) + if (fabs(qtemp[i]-quat_director[i])>1E-9) + fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n", + i, qtemp[i], quat_director[i]); + } + else + { + quat[0]=1; + quat[1]=quat[2]=quat[3]=0; + } + + // Set the particle id of the particle we want to relate to, the distnace + // and the relative orientation + if (set_particle_vs_relative(part_num, relate_to, l, quat) == ES_ERROR) { + runtimeErrorMsg() << "setting the vs_relative attributes failed"; + return ES_ERROR; + } + set_particle_virtual(part_num,1); + + return ES_OK; +} + + + + +#endif #endif diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites/VirtualSites.hpp similarity index 68% rename from src/core/virtual_sites.hpp rename to src/core/virtual_sites/VirtualSites.hpp index 6a17a0ae162..4193b78c910 100755 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -18,8 +18,8 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . */ -#ifndef _VIRTUAL_SITES_HPP -#define _VIRTUAL_SITES_HPP +#ifndef VIRTUAL_SITES_VIRTUAL_SITES_HPP +#define VIRTUAL_SITES_VIRTUAL_SITES_HPP /** \file virtual_sites.hpp @@ -48,6 +48,10 @@ class VirtualSites { virtual void update(bool recalc_positions=true) const =0; /** Back-transfer forces (and torques) to non-virtual particles */ virtual void back_transfer_forces_and_torques() const =0; + /** @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& have_velocity() const { return m_have_velocity; }; @@ -62,26 +66,5 @@ class VirtualSites { }; - /** @brief Do-nothing virtual-sites implementation */ - class VirtualSitesOff : public VirtualSites { - /** @brief Update positions and/or velocities of virtual sites - - * Velocities are only updated update_velocities() return true - * @param recalc_positions can be used to skip the reculation of positions - */ - void update(bool recalc_positions=true) const override {}; - /** Back-transfer forces (and torques) to non-virtual particles */ - void back_transfer_forces_and_torques() const override {}; - /** @brief Is a ghost communication needed after position updates */ - bool need_ghost_comm_after_pos_update() const override { return false;} - /** Is a ghost comm needed before a velocity update */ - bool need_ghost_comm_before_vel_update() const override {return false;}; - /** Is a ghost comm needed before back_transfer */ - bool need_ghost_comm_before_back_transfer() const override {return false;}; - }; - -const std::shared_ptr& virtual_sites(); - -void set_virtual_sites(std::shared_ptr v); #endif #endif diff --git a/src/core/virtual_sites_relative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp similarity index 64% rename from src/core/virtual_sites_relative.cpp rename to src/core/virtual_sites/VirtualSitesRelative.cpp index 7c1beba26aa..10e16498747 100755 --- a/src/core/virtual_sites_relative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -18,9 +18,8 @@ */ #include "config.hpp" -#include "virtual_sites_relative.hpp" +#include "VirtualSitesRelative.hpp" #include "rotation.hpp" -#include "virtual_sites_relative.hpp" #include "errorhandling.hpp" #ifdef VIRTUAL_SITES_RELATIVE @@ -174,107 +173,6 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { } // Setup the virtual_sites_relative properties of a particle so that the given virtaul particle will follow the given real particle -int vs_relate_to(int part_num, int relate_to) -{ - // Get the data for the particle we act on and the one we wnat to relate - // it to. - auto p_current = get_particle_data(part_num); - auto p_relate_to = get_particle_data(relate_to); - if (!p_current || !p_relate_to) { - runtimeErrorMsg() <<"Could not retrieve particle data for the given id"; - return ES_ERROR; - } - - // get the distance between the particles - double d[3]; - get_mi_vector(d, p_current->r.p,p_relate_to->r.p); - - - - // Check, if the distance between virtual and non-virtual particles is larger htan minimum global cutoff - // If so, warn user - double l=sqrt(sqrlen(d)); - if (l>min_global_cut) { - runtimeErrorMsg() << "Warning: The distance between virtual and non-virtual particle (" << l << ") is\nlarger than the minimum global cutoff (" << min_global_cut << "). This may lead to incorrect simulations\nunder certain conditions. Set the \"System()\" class property \"min_global_cut\" to\nincrease the minimum cutoff.\n"; - return ES_ERROR; - } - - // Now, calculate the quaternions which specify the angle between - // the director of the particel we relate to and the vector - // (paritlce_we_relate_to - this_particle) - double quat[4]; - // The vs_relative implemnation later obtains the direcotr by multiplying - // the quaternions representing the orientation of the real particle - // with those in the virtual particle. The re quulting quaternion is then - // converted to a director. - // Whe have quat_(real particle) *quat_(virtual particle) - // = quat_(obtained from desired director) - // Resolving this for the quat_(virtaul particle) - - //Normalize desired director - int i; - - // If the distance between real & virtual particle is 0 - // we just set the relative orientation to 1 0 0 0, as it is irrelevant but - // needs to be a valid quaternion - if (l!=0) - { - for (i=0;i<3;i++) - d[i]/=l; - - // Obtain quaternions from desired director - double quat_director[4]; - convert_quatu_to_quat(d, quat_director); - - // Define quat as described above: - double x=0; - for (i=0;i<4;i++) - x+=p_relate_to->r.quat[i]*p_relate_to->r.quat[i]; - - quat[0]=0; - for (i=0;i<4;i++) - quat[0] +=p_relate_to->r.quat[i]*quat_director[i]; - - quat[1] =-quat_director[0] *p_relate_to->r.quat[1] - +quat_director[1] *p_relate_to->r.quat[0] - +quat_director[2] *p_relate_to->r.quat[3] - -quat_director[3] *p_relate_to->r.quat[2]; - quat[2] =p_relate_to->r.quat[1] *quat_director[3] - + p_relate_to->r.quat[0] *quat_director[2] - - p_relate_to->r.quat[3] *quat_director[1] - - p_relate_to->r.quat[2] * quat_director[0]; - quat[3] =quat_director[3] *p_relate_to->r.quat[0] - - p_relate_to->r.quat[3] *quat_director[0] - + p_relate_to->r.quat[2] * quat_director[1] - - p_relate_to->r.quat[1] *quat_director[2]; - for (i=0;i<4;i++) - quat[i]/=x; - - - // Verify result - double qtemp[4]; - multiply_quaternions(p_relate_to->r.quat,quat,qtemp); - for (i=0;i<4;i++) - if (fabs(qtemp[i]-quat_director[i])>1E-9) - fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n", - i, qtemp[i], quat_director[i]); - } - else - { - quat[0]=1; - quat[1]=quat[2]=quat[3]=0; - } - - // Set the particle id of the particle we want to relate to, the distnace - // and the relative orientation - if (set_particle_vs_relative(part_num, relate_to, l, quat) == ES_ERROR) { - runtimeErrorMsg() << "setting the vs_relative attributes failed"; - return ES_ERROR; - } - set_particle_virtual(part_num,1); - - return ES_OK; -} diff --git a/src/core/virtual_sites_relative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp similarity index 90% rename from src/core/virtual_sites_relative.hpp rename to src/core/virtual_sites/VirtualSitesRelative.hpp index 929b425bf3b..bd56c4803ce 100755 --- a/src/core/virtual_sites_relative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -17,9 +17,10 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . */ -#ifndef _VIRTUAL_SITES_RELATIVE_HPP -#define _VIRTUAL_SITES_RELATIVE_HPP +#ifndef VIRTUAL_SITES_VIRTUAL_SITES_RELATIVE_HPP +#define VIRTUAL_SITES_VIRTUAL_SITES_RELATIVE_HPP +#include "config.hpp" #ifdef VIRTUAL_SITES_RELATIVE #include "virtual_sites.hpp" @@ -44,14 +45,14 @@ /** Is a ghost comm needed before a velocity update */ bool need_ghost_comm_before_vel_update() const override {return (n_nodes>1) && have_velocity();}; bool need_ghost_comm_before_back_transfer() const override {return true;}; - void pressure_and_stress_tensor_contribution(double* pressure, double* stress_tensor) const; + int n_pressure_contribs() const override {return 1;}; + void pressure_and_stress_tensor_contribution(double* pressure, double* stress_tensor) const override; private: void update_pos(Particle& p) const; void update_vel(Particle& p) const; }; -int vs_relate_to(int part_num, int relate_to); #endif #endif diff --git a/src/core/virtual_sites_com.cpp b/src/core/virtual_sites_com.cpp deleted file mode 100755 index e1fdb48907a..00000000000 --- a/src/core/virtual_sites_com.cpp +++ /dev/null @@ -1,238 +0,0 @@ -/* - Copyright (C) 2010,2011,2012,2013,2014,2015,2016 The ESPResSo project - Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - Max-Planck-Institute for Polymer Research, Theory Group - - This file is part of ESPResSo. - - ESPResSo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ESPResSo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#include "virtual_sites_com.hpp" - -#ifdef VIRTUAL_SITES_COM - -#include "virtual_sites.hpp" -#include "cells.hpp" -#include "topology.hpp" -#include "forces.hpp" -#include "partCfg_global.hpp" -#include "integrate.hpp" - -// forward declarations -void calc_mol_vel(Particle *p_com,double v_com[3]); -void calc_mol_pos(Particle *p_com,double r_com[3]); -void put_mol_force_on_parts(Particle *p_com); - -void update_mol_vel_particle(Particle *p){ - int j; - double v_com[3]; - if (p->p.isVirtual) { - calc_mol_vel(p,v_com); - for (j=0;j<3;j++){ - p->m.v[j] = v_com[j]; - } - } -} - - -void update_mol_pos_particle(Particle *p){ - int j; - double r_com[3]; - if (p->p.isVirtual) { - calc_mol_pos(p,r_com); - for (j=0;j<3;j++){ - p->r.p[j] = r_com[j]; - } - } -} - -void distribute_mol_force() { - for (auto &p : local_cells.particles()) { - if (p.p.isVirtual) { - if (p.p.isVirtual) { - if (sqrlen(p.f.f) != 0) { - put_mol_force_on_parts(&p); - } - } - } -} - -void calc_mol_vel(Particle *p_com,double v_com[3]){ - int i,j,mol_id; - double M=0; - Particle *p; -#ifdef VIRTUAL_SITES_DEBUG - int count=0; -#endif - for (i=0;i<3;i++){ - v_com[i]=0.0; - } - mol_id=p_com->p.mol_id; - for (i=0;ip.isVirtual) continue; - for (j=0;j<3;j++){ - v_com[j] += (*p).p.mass*p->m.v[j]; - } - M+=(*p).p.mass; -#ifdef VIRTUAL_SITES_DEBUG - count++; -#endif - } - for (j=0;j<3;j++){ - v_com[j] /= M; - } -#ifdef VIRTUAL_SITES_DEBUG - if (count!=topology[mol_id].part.n-1){ - runtimeErrorMsg() <<"There is more than one COM in calc_mol_vel! mol_id= " << mol_id << "\n"; - return; - } -#endif -} - -/* this is a local version of center of mass, because ghosts don't have image boxes*/ -/* but p_com is a real particle */ -void calc_mol_pos(Particle *p_com,double r_com[3]){ - int i,j,mol_id; - double M=0; - double vec12[3]; - Particle *p; -#ifdef VIRTUAL_SITES_DEBUG - int count=0; -#endif - for (i=0;i<3;i++){ - r_com[i]=0.0; - } - mol_id=p_com->p.mol_id; - for (i=0;ip.isVirtual) continue; - get_mi_vector(vec12,p->r.p, p_com->r.p); - for (j=0;j<3;j++){ - r_com[j] += (*p).p.mass*vec12[j]; - } - M+=(*p).p.mass; -#ifdef VIRTUAL_SITES_DEBUG - count++; -#endif - } - for (j=0;j<3;j++){ - r_com[j] /= M; - r_com[j] += p_com->r.p[j]; - } -#ifdef VIRTUAL_SITES_DEBUG - if (count!=topology[mol_id].part.n-1){ - runtimeErrorMsg() <<"There is more than one COM in calc_mol_pos! mol_id= " << mol_id << "\n"; - return; - } -#endif -} - -void put_mol_force_on_parts(Particle *p_com){ - int i,j,mol_id; - Particle *p; - double force[3],M; -#ifdef VIRTUAL_SITES_DEBUG - int count=0; -#endif - mol_id=p_com->p.mol_id; - for (i=0;i<3;i++){ - force[i]=p_com->f.f[i]; - p_com->f.f[i]=0.0; - } -#ifdef MASS - M=0; - for (i=0;ip.isVirtual) continue; - M+=(*p).p.mass; - } -#else - M=topology[mol_id].part.n-1; -#endif - for (i=0;ip.isVirtual) { - for (j=0;j<3;j++){ - p->f.f[j]+=(*p).p.mass*force[j]/M; - } -#ifdef VIRTUAL_SITES_DEBUG - count++; -#endif - } - } -#ifdef VIRTUAL_SITES_DEBUG - if (count!=topology[mol_id].part.n-1){ - runtimeErrorMsg() <<"There is more than one COM input_mol_force_on_parts! mol_id= " << mol_id << "\n"; - return; - } -#endif -} - -Particle *get_mol_com_particle(Particle *calling_p){ - int mol_id; - int i; - Particle *p; - - mol_id=calling_p->p.mol_id; - - if (mol_id < 0) { - runtimeErrorMsg() <<"Particle does not have a mol id! pnr= " << calling_p->p.identity << "\n"; - return nullptr; - } - for (i=0;ip.isVirtual) { - return p; - } - } - - runtimeErrorMsg() <<"No com found in get_mol_com_particleParticle does not exist in put_mol_force_on_parts! pnr= " << calling_p->p.identity << "\n"; - return nullptr; - - return calling_p; -} - - -#endif diff --git a/src/core/virtual_sites_com.hpp b/src/core/virtual_sites_com.hpp deleted file mode 100755 index a022658ad45..00000000000 --- a/src/core/virtual_sites_com.hpp +++ /dev/null @@ -1,46 +0,0 @@ -/* - Copyright (C) 2010,2011,2012,2013,2014,2015,2016 The ESPResSo project - Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - Max-Planck-Institute for Polymer Research, Theory Group - - This file is part of ESPResSo. - - ESPResSo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ESPResSo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#ifndef _VIRTUAL_SITES_COM_H -#define _VIRTUAL_SITES_COM_H - -#include "config.hpp" - -#ifdef VIRTUAL_SITES_COM -#include "particle_data.hpp" - -// The following three functions have to be provided by all implementations -// of virtual sites -// Update the vel/pos of the given virtual particle as defined by the real -// particles in the same molecule -void update_mol_pos_particle(Particle *); -void update_mol_vel_particle(Particle *); - -// Distribute forces that have accumulated on virtual particles to the -// associated real particles -void distribute_mol_force(); - -// Gets the (first) virtual particle of the same molecule as the given (real) -// particle -Particle *get_mol_com_particle(Particle *calling_p); - -#endif - -#endif diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index c5a2db304f4..6e3fb47ceeb 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -367,10 +367,7 @@ class Analysis(object): * "nonbonded_inter" type_i, type_j", nonboned pressure between short ranged forces between type i and j and different mol_ids * "coulomb", Maxwell stress, how it is calculated depends on the method * "dipolar", TODO - * "vs_relative", In case of rigid body rotation, virial contribution from torques is not included. - The pressure contribution for rigid bodies constructed by means of the - VIRTUAL\_SITES\_RELATIVE mechanism is included. On the other hand, the - pressure contribution for rigid bonds is not included. + * "virtual_sites", Stress contribution due to virtual sites """ v_comp=int(v_comp) @@ -448,8 +445,12 @@ class Analysis(object): p["dipolar"] = total_dipolar # virtual sites - IF VIRTUAL_SITES_RELATIVE == 1: - p["vs_relative"] = c_analyze.total_pressure.vs_relative[0] + IF VIRTUAL_SITES == 1: + p_vs=0. + for i in range(c_analyze.total_pressure.n_virtual_sites): + p_vs += c_analyze.total_pressure.virtual_sites[i] + p["virtual_sites",i] = c_analyze.total_pressure.virtual_sites[0] + p["virtual_sites"] = p_vs return p @@ -471,10 +472,7 @@ class Analysis(object): * "nonbonded_inter type_i", type_j, nonboned stress tensor between short ranged forces between type i and j and different mol_ids * "coulomb", Maxwell stress tensor, how it is calculated depends on the method * "dipolar", TODO - * "vs_relative", In case of rigid body rotation, virial contribution from torques is not included. - The stress_tensor contribution for rigid bodies constructed by means of the - VIRTUAL\_SITES\_RELATIVE mechanism is included. On the other hand, the - pressure contribution for rigid bonds is not included. + * "virtual_sites", Stress tensor contribution for virtual sites """ v_comp=int(v_comp) @@ -566,8 +564,13 @@ class Analysis(object): # virtual sites IF VIRTUAL_SITES_RELATIVE == 1: - p["vs_relative"] = np.reshape(create_nparray_from_double_array( - c_analyze.total_p_tensor.vs_relative, 9), (3, 3)) + total_vs = np.zeros(9) + for i in range(c_analyze.total_p_tensor.n_virtual_sites): + p["virtual_sites", i] = np.reshape( + create_nparray_from_double_array( + c_analyze.total_p_tensor.virtual_sites +9*i, 9), (3, 3) ) + total_virtual_sites += p["virtual_sites", i] + p["virtual_sites"] = total_virtual_sites return p diff --git a/src/python/espressomd/c_analyze.pxd b/src/python/espressomd/c_analyze.pxd index 332ce8d3c42..1887947b6af 100755 --- a/src/python/espressomd/c_analyze.pxd +++ b/src/python/espressomd/c_analyze.pxd @@ -50,11 +50,12 @@ cdef extern from "statistics.hpp": int n_coulomb int n_dipolar int n_non_bonded + int n_virtual_sites double * bonded double * non_bonded double * coulomb double * dipolar - double * vs_relative + double * virtual_sites cdef extern from "statistics.hpp": ctypedef struct Observable_stat_non_bonded: diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 24eb1c39b8b..9e599a0a036 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -206,7 +206,7 @@ cdef extern from "particle_data.hpp": bool particle_exists(int part) -cdef extern from "virtual_sites_relative.hpp": +cdef extern from "virtual_sites.hpp": IF VIRTUAL_SITES_RELATIVE == 1: int vs_relate_to(int part_num, int relate_to) int set_particle_vs_relative(int part, int vs_relative_to, double vs_distance, double * vs_quat) diff --git a/src/script_interface/virtual_sites/ActiveVirtualSitesHandle.hpp b/src/script_interface/virtual_sites/ActiveVirtualSitesHandle.hpp index 950c971acc8..9d6378038b7 100644 --- a/src/script_interface/virtual_sites/ActiveVirtualSitesHandle.hpp +++ b/src/script_interface/virtual_sites/ActiveVirtualSitesHandle.hpp @@ -41,12 +41,7 @@ class ActiveVirtualSitesHandle : public AutoParameters { [this](Variant const &value) { m_active_implementation = get_value>(value); - if (m_active_implementation) { - ::set_virtual_sites(m_active_implementation->virtual_sites()); - } - else { - throw std::runtime_error("the implementation has to be an instance of a sub-class of VirtualSites"); - }; + ::set_virtual_sites(m_active_implementation->virtual_sites()); }, [this]() { return (this->m_active_implementation != nullptr) ? this->m_active_implementation->id() : ObjectId(); diff --git a/src/script_interface/virtual_sites/VirtualSitesOff.hpp b/src/script_interface/virtual_sites/VirtualSitesOff.hpp index da8854f0b4c..489456bb85e 100644 --- a/src/script_interface/virtual_sites/VirtualSitesOff.hpp +++ b/src/script_interface/virtual_sites/VirtualSitesOff.hpp @@ -24,7 +24,7 @@ #include "config.hpp" #include "VirtualSites.hpp" -#include "core/virtual_sites.hpp" +#include "core/virtual_sites/VirtualSitesOff.hpp" #ifdef VIRTUAL_SITES namespace ScriptInterface { diff --git a/src/script_interface/virtual_sites/VirtualSitesRelative.hpp b/src/script_interface/virtual_sites/VirtualSitesRelative.hpp index b6acdaffa13..fc0b8fe2678 100644 --- a/src/script_interface/virtual_sites/VirtualSitesRelative.hpp +++ b/src/script_interface/virtual_sites/VirtualSitesRelative.hpp @@ -24,7 +24,7 @@ #include "config.hpp" #include "VirtualSites.hpp" -#include "core/virtual_sites_relative.hpp" +#include "core/virtual_sites/VirtualSitesRelative.hpp" #ifdef VIRTUAL_SITES_RELATIVE namespace ScriptInterface {