simtk.org

I have a system with a COM force between two groups that is in a separate force group, 1. All other forces are in group 0. To report only forces from group 0, one would call `system.getState(..., groups={0})`. My question is if this can be done using a Reporter? From the manual example, one would naively expand the return tuple to `(steps, False, False, True, False, False, False, {0})`.

A related question to center-of-mass forces: The `CustomCentroidBondForce::usesPeriodicBoundaryConditions()` function returns false. Does this mean that the `distance()` function does not return the minimum image distance?

You just need to call simulation.context.getState() yourself. Normally the simulation gets a state for you and passes it to the reporter, but that's just an optimization. It's so that if multiple reporters are all getting called at the same time, only one state needs to be created. But instead, you can have describeNextReport() return false for all the arguments (since it won't use any information from the state passed to it), then retrieve a state yourself.

Does this mean that the `distance()` function does not return the minimum image distance?

Correct.

Peter

Peter, thanks for your quick response.

A few more questions:

• 1. Is there a way to get the MI distance between groups? I would call state.getPeriodicBoxVectors() and use the side lengths in a custom distance function and manually calculate COMs for g1, g2 etc. also from a custom function. A periodicdistance() function as in i.e. CustomExternalForce would seem useful.
  
  2. Does the mass center calculation respect periodic boundaries? I.e. if a collection of particles is wrapped around the boundaries, one would normally not want the COM to be in the middle of the box. Something like this:

  Code: Select all

  def centerOfMass(positions, box):
      """ Calculates the geometric center taking into account periodic boundaries
      
      More here: https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions
      """
      theta=np.divide(positions, box).astype(np.float) * 2*np.pi
      x1=np.array( [np.cos(theta[:,0]).mean(), np.cos(theta[:,1]).mean(), np.cos(theta[:,2]).mean()] )
      x2=np.array( [np.sin(theta[:,0]).mean(), np.sin(theta[:,1]).mean(), np.sin(theta[:,2]).mean()] )
      return box * (np.arctan2(-x1,-x2)+np.pi) / (2*np.pi)

(my questions arise from writing a force pulling scheme for PMF calculations)