Goal

• The goal is, to have a virtual particle at the center of mass of a bunch of other particles (say, a molecule).
• Forces acting on this center of mass are distributed back onto the constituents
• Note that the pattern is opposite to VS Relative: the position/velocity/mass of the virtaul site at center of mass is calculated from the positions/velocities/masses of many particles.

Design choices

• the implementation is optimized for a small number of objects for which a center of mass particle is placed, but theses objects consist of many particles

• the mol_id particle property is used to distinguish between the different objects/molecules

• A map of (mol_id, virtual_stie_particle_id) pairs is stored in the class implementing the algorithm. Obtaining this informaiton on the fly would require a second pass over the particles and a global communication.

• The center of mass is calculated using unfolded positions. I.e. if the user creates a molecule wrapping across the periodic boundaries, they need to make sure to set this correctly:

  p1.pos = syste.box_l-dx; p2.pos= system.box_l +dx rather than p2.pos=dx/2

Steps

Core class

• in src/core/virtual_sites/VirtualSitesCenterOfMass.cpp/hpp

• data member std::map<int,int> m_vs_pid_for_mol_mid

• implement void VirtualSitesCenterOfMass::update()

  • define a struc type ComINfo holidng members for total mass and weighted positoin sum., both default constructed to zero.

  • Initialize a `std::map<int, ComInfo> com_by_mol_id.

  • Add entries for all mol_ids contained in the m_vs_pid_by_mol_id

  • Loop over all particles cell_sturcuter().local_particles()

    • if the mol_id of a particle is in the com_by_mol_id map, add itts mass to te ComInfo::total_mass and the UNFOLDED position times mass to thei wieghted_sum

  • loop over all (mol_id, ComINfo) pairs

    • skip the particle if it is virtual.
    • calculate the center of mass from wieghed_sum/total_mass
    • if Es runs on more than one MPI rank, through an exception. Later, a reduction step can be added to make this parallel.
    • Look up the virtual site particle id for the mol id from the m_vs_pid_by_mol_id map
    • get the virtual site by its particle id cell_srucutre().get_local_particle()
    • if a nullptr is returend, through (will have to e be changed in a parlallel implementation)
    • set the particle position to the cetner of mass position, and the mass to the recorded total mass
    • Fold the particle positoin (and image box) box_geo.fold_position()

• Implement the VirtualStiesCenterOfMass::backtransfer_forces_and_torques() funciont
  to be contineud