Re-introduce virtual sites at center of mass

[RudolfWeeber]

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

[RudolfWeeber]

Fixup

Please assign not only the center of mass position/velocity to the virtual site, but also the total mass of the non-virtual particles making up the object

Backtransfer

the backtransfer funciton should do the follow

• Go over all paritcles in the system

• If the partice's molid is in m_vs_pid_for_mol_id map:

  • Retreive the virtual site for that mol_id (i.e., the pid from the map)
  • Get the force from the virtual site
  • Assign the virutal site's force to the non-virtual particle with the prefactor m_{nonvirtual}/m{Total}, where m_total is the mass of the entire object (sotred on the virtual site, see fixup)

Script interface

Add a new virtual site class in src/script_interface/virtual_sites. Look for VirtualStiesRelative in the script interface and use as template

@jngrad Where can an exammple for the script interface stuff for the std::map<int,int> be stolen?