openpathsampling / openpathsampling

An open source Python framework for transition interface and path sampling calculations.
http://openpathsampling.org
MIT License
99 stars 47 forks source link

Allow for periodic collective variables #308

Open jhprinz opened 8 years ago

jhprinz commented 8 years ago

Maybe adding a flag, that a cv can be periodic might make it easier to analyze these later. it could somehow replace the PeriodicVolume or parts of it.

It would also be useful for plotting trajectories in phi,psi space. A flip in the periodic position could be easily detected and so could avoid having these flip lines in the phi-psi plots for trajectories.

dwhswenson commented 8 years ago

PeriodicVolume is designed so that the user doesn't need to know what the bounds of the periodic domain are (only the size of the period, and mostly, not even that). The whole point is to remove any concern about the details about the domain of which a periodic CV is implemented. Let's keep that as a target: maybe the user needs to know that an angle has a periodicity of 360 degrees, but the user should never need to worry about whether the internal implementation returns angles between 0 and 360 or between -180 and 180. In fact, if the user wants to display trajectories using bounds -90 to 270, that should also be easy. What I describe below should make that possible, with minimal effort.

Your primary concern is about the time series having large hops (this is the reason you see jump lines when a trajectory wraps), so the issue is not with the CV (which deals with only one snapshot at a single time). Rather, it is with the analysis (which takes a whole time series). When you only know one frame, you can never know whether the previous frame was in this periodic image or a neighboring one. It's when you draw a line between two periodic images that there are problems. So it won't be fixed by a new class of CVs.

The way to fix this is in the analysis, and is based on changes of greater than $L/2$ (where $L$ is the period). My suggestion: add an option to CVs to define a period (default to None), so that all analysis routines have knowledge of the period for the CV. Creating a separate PeriodicCV would require checking the class; adding an attribute to the CV means you just check that attribute.

For analysis purposes, I would suggest creating a function that unwraps according to a given CV. I'd call it CollectiveVariable.unwrap(self, time_series) and basically, it returns a list by adding/subtracting cv.period every time abs(time_series[i-1] - time_series[i]) > cv.period * 0.5. (Remove the abs to figure out add/subtract.) If a CV isn't periodic, or if the whole trajectory stays in one image, cv.unwrap(cv(traj)) == cv(traj). One important point: it alway trusts that the first number is in the desired periodic image, and just builds a continuous trajectory from there. (See below on drawing routines; this is also why I say it should take a time series of numbers, not a trajectory.)

This will work for normal analysis routines, because it gives you the continuous trajectory. However, it won't work for the drawing routines, where you actually do want to wrap the trajectory into a single periodic image. The cheap and easy way to do it is to modify to drawing routines to split the trajectory into segments, with different segments any time there's a jump of $L/2$ or greater. But in this case, you won't see the line indicating that you've gone out of the visualized periodic image (it will just stop at the last frame in that image; no line extending to the edge). To go out of that, the best thing to do is probably to use the unwrapping function above on versions of the trajectory shifted by a periodic image in each direction: so for a plot of 2 periodic CVs, you'd draw 5 total lines. Then you'll just have to set the viewing bounds to avoid showing much replication.

Trajectory plots of 1 periodic CV are best done with that unwrap function directly, although an approach similar to 2D could be used. Trajectory plots of 2 periodic CVs can be done as described in the previous paragraph. Trajectory plots of 3 periodic CVs gets harder to visualize: I think something analogous to the 2D case can be done, but I'm not 100% sure that it will look good. Plots of 4 or more periodic CVs probably require psychedelic drugs to visualize, so they are beyond the scope of this project.