How to reject a path if not in ensamble

@dwhswenson

HI @dwhswenson
I had this issue in TIS simulation (if you remember) where gas bubble was changing from sphere to cylinder in some paths. I have figured out a way to find if a frame is a sphere or cylinder. My goal now is to check at every frame if the gas is in sphere phase. If not, I want to reject that path (mc move) and sample again

My question is
1- Is it an issue with detail balance to reject such a move ?
2- What is the best way to do it. Currently I am doing is if it is Cylinder then I return a value of CV that is outside the sampling region (like CV = ∞ ) . This will reject the path (I think so, still checking ).
What would be a better way to do it in OPS ?

To rephrase the problem a little more generally, I believe the issue is that there are certain configurations (i.e., a volume in phase space) that you want to forbid your trajectories from visiting. You have the function to identify those configurations, so it's just a matter of telling OPS not to visit those areas.

Is it an issue with detail balance to reject such a move ?

I think this is fine; the main danger is whether the underlying physics changes if you're near the cylinder bubble state, not just in it. Otherwise, you're simply disallowing that region of phase space because it isn't physically valid. We do this when a move generates a NaN (discussion in #570, especially #570 (comment)). The cylinder-bubble is analogous to forbidden values of x in that example.

What is the best way to do it. Currently I am doing is if it is Cylinder then I return a value of CV that is outside the sampling region (like CV = ∞ ) . This will reject the path

I don't think this will work. If your state is defined as cv >= boundary, then cv=float('inf') will be in that state.

What would be a better way to do it in OPS ?

I drafted up a really long and convoluted suggestion here, then I came across something in the code that I implemented specifically to solve this kind of problem! 😆

The key is the intersect_with parameter for VolumeInterfaceSets. See this notebook from the additional examples for details. Basically, the idea is:

Create a CV for your cylinder bubble. Let's say it returns +1 if cylinder, -1 if not.
Create a CVDefinedVolume for the cylinder bubble, e.g., with lambda_min=0, lambda_max=float('inf').
Use the kwarg intersect_with=~cylinder_bubble_vol when creating your volume interface set (where the ~ says to only allow things outside this; you could always create allowed_bubble_vol; it shouldn't make any real difference in performance).
To avoid having your initial/final states having the cylinder bubble (not that it should cause a problem, but formally the correct thing to do), use the & operator to intersect with your existing initial/final state volumes: new_initial = initial & ~cylinder_bubble_vol.

Hey David,
Sounds good. although a bit stuck at creating network. I get the following error
AttributeError: 'function' object has no attribute 'name'

Will be glad if you can have a quick look at the script
github_issue_TIS.zip

Quick look: I do the following

cylinder_bubble_vol= paths.CVDefinedVolume(cylinder_check, lambda_min=0, lambda_max=float('inf')).named("cylinder")
interface_1 = paths.VolumeInterfaceSet(MCG, float("-inf"), [5,6,7, 8, 9 ,10,11], intersect_with=~cylinder_bubble_vol)

Looks (at a brief glance) like cyclinder_check is a function -- you didn't wrap it in an OPS collective variable (e.g., CoordinateFunctionCV or MDTrajFunctionCV). The input to CVDefinedVolume needs to be an OPS CV.

I'd recommend getting the CV working on a template snapshot before trying to run your whole script; making sure it has the right input/output types can take a little work (needs to input an OPS snapshot, and output a single float/int).

Oh, My bad. It works now. Will get back if there is any issue

@bolhuis

Hi David.
I used the above definition to sample paths and recently got some weird accepted paths.

I have a function cylinder_check that returns +1 if a frame is in cylinder and -1 if a frame is in sphere bubble shape. I define the states and interfaces as below.

cylinder = paths.MDTrajFunctionCV("cylinder", cylinder_check, template.topology)
cylinder_bubble_vol= paths.CVDefinedVolume(cylinder, lambda_min=0, lambda_max=float('inf')).named("cylinder")
interface_1 = paths.VolumeInterfaceSet(MCG, float("-inf"), [2,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100],intersect_with=~cylinder_bubble_vol)
crystal_new  =crystal & ~cylinder_bubble_vol
liquid_new = liquid & ~cylinder_bubble_vol

When I was sampling interface 60, I got paths like following.
1- Clearly the path did not cross interface 60 (ie lambda > 60)
2- On visual inspection, it converted to "Cylinder" when it reaches state A (marked by black arrow).

Clearly the logic for ensamble definition is not right. @bolhuis suggested me 2 ways define the ensemble :

(In (A & not cylinder)
Out (A || B || cylinder)
In ((A || B) & not cylinder ) ) and exists Out (60 || cylinder).

or

(In (A & not cylinder)
In ( not (A || B) & not cylinder)
In ((A || B) & not cylinder ) ) and exists In ( MCG> 60 & not cylinder)

My question is : What is the best way to include this definition? Will I need to use SequentialEnsemble as follows? I need to run RETIS simulation and using SequentialEnsemble to build RETIS ensamble will require thinking. If there is any reference notebook, please do share

vol = paths.CVDefinedVolume(MCG, float("-inf"),60)
paths.SequentialEnsemble([
    paths.LengthEnsemble(1) & paths.AllInXEnsemble(liquid_new) & paths.AllInXEnsemble(~cylinder_bubble_vol),
    paths.OptionalEnsemble(paths.AllOutXEnsemble(liquid_new | crystal_new | ~cylinder_bubble_vol)),
    paths.LengthEnsemble(1) & paths.AllInXEnsemble(liquid_new | crystal_new & ~cylinder_bubble_vol ),  
    paths.LengthEnsemble(1) & AllOutXEnsemble(vol)
])

I have attached my log file just incase for reference
cylinder_stopping.zip

I think there's a reasonably straightforward way to hack this, as long as you're using the MISTIS network. The trick is to use the strict_sampling flag in MISTISNetwork (see #509). Essentially, strict_sampling says that the A->B sampling transition should reject paths that visit any state in the network other than A and B.

More mathematically, strict_sampling can be used to give you something like normal_TIS_interface & AllOut(cylinder), which is, I think, what you really want.

So I think you can accomplish what you want with the following:

Add a fake transition involving the cylinder volume.
Use strict_sampling=True to forbid visiting that volume.
Be careful with the move scheme -- you'll need to forbid sampling of ensembles associated with the fake transition. This part requires careful attention, and I'm not going to go into the details here.

Here's a quick example of what this could look like:

import openpathsampling as paths
cv_x = paths.FunctionCV("cv_x", lambda snap: snap.xyz[0][0])
is_evil = paths.FunctionCV("is_evil", lambda snap: snap.xyz[0][1])
forbidden = paths.CVDefinedVolume(is_evil, lambda_min=0, lambda_max=float('inf'))
state_A = paths.CVDefinedVolume(cv_x, float('-inf'), 2)
state_B = paths.CVDefinedVolume(cv_x, 10, float('inf'))
interfaces = paths.VolumeInterfaceSet(cv_x, float('-inf'), [2, 4, 6])
fake_interfaces = paths.InterfaceSet(volumes=[paths.FullVolume()])
network = paths.MISTISNetwork([
    (state_A, interfaces, state_B),
    (paths.EmptyVolume(), fake_interfaces, forbidden)
], strict_sampling=True)

I use paths.EmptyVolume() and paths.FullVolume() here because (1) they ensure that no trajectory is valid (the transition must start somewhere that doesn't exist, and transition out of everywhere that does exist) and (2) they get automatically ignored by by volume combinations (vol_A | EmptyVolume() is vol_A -- it doesn't even create a new object).

Here's an example of the ensembles this outputs:

>>> print(network.sampling_ensembles[0])
(
  (
    [
    (
      x[t] in {x|cv_x(x) in [-inf, 2]} for all t
    )
    and
    (
      len(x) = 1
    ),
    {x[t] in (not ({x|cv_x(x) in [-inf, 2]} or {x|cv_x(x) in [10, inf]})) for all t} (OPTIONAL),
    (
      x[t] in ({x|cv_x(x) in [-inf, 2]} or {x|cv_x(x) in [10, inf]}) for all t
    )
    and
    (
      len(x) = 1
    )
    ]
  )
  and
  (
    exists t such that x[t] in (not {x|cv_x(x) in [-inf, 2]})
  )
)
and
(
  x[t] in (not {x|is_evil(x) in [0, inf]}) for all t
)

This is just the normal TIS ensemble, with the extra requirement that you never visit the forbidden state. I think that is what you want, right?

Hi David,

Yes this is what I want to achieve and your example makes sense however doesn't work right away for RETIS.
For RETIS, I have an initial trajectory that is all in "Sphere" phase (so no cylinder frame is present). It looks like this

When I incorporate fake interface as above, I get stuck in a loop that cant find a suitable initial trajectory. (This trajectory used to work fine before. I also tried strict_sampling=False and still get same issue)

My script (attached is full script)

crystal = (paths.CVDefinedVolume(MCG, lambda_min=25, lambda_max=350)
        )
liquid = (paths.CVDefinedVolume(MCG,lambda_min = 1, lambda_max = 2.1)
         )

cylinder = paths.MDTrajFunctionCV("cylinder", cylinder_check, template.topology)
cylinder_bubble_vol= paths.CVDefinedVolume(cylinder, lambda_min=0, lambda_max=float('inf')).named("cylinder")


interface_1 = paths.VolumeInterfaceSet(MCG, float("-inf"), [2.9,3.9,4.9,5.9,6.9,7.9])
fake_interfaces = paths.InterfaceSet(volumes=[paths.FullVolume()])


forbidden = cylinder_bubble_vol

#network = paths.MISTISNetwork([(liquid_new, interface_1 ,crystal_new)])
network = paths.MISTISNetwork([
    (liquid, interface_1, crystal),
    (paths.EmptyVolume(), fake_interfaces, forbidden)
], strict_sampling=True)

scheme = paths.DefaultScheme(network, engine = engine)

Attached is the detail log file .

finished generating: still missing 2 samples
Missing ensembles:
*  [B->cylinder 0 strict]
*  [B MIS minus]
No extra ensembles.

Not sure if my implementation is wrong or I am missing something
cylinder_stopping_2.zip

Don't use DefaultScheme. This is what I meant when I said you need to be careful with the move scheme. You need to specify the ensembles for each move type.

Finding the correct minus ensemble is a little tricky. The rest is easy:

allowed_ensembles = network.sampling_transitions[0].ensembles
minus_ens = [ens for ens in network.minus_ensembles if ens.state_vol == state_A][0]
scheme = paths.MoveScheme(network)
scheme.append([
    paths.strategies.OneWayShootingStrategy(ensembles=allowed_ensembles),
    paths.strategies.PathReversalStrategy(ensembles=allowed_ensembles),
    paths.strategies.NearestNeighborRepExStrategy(ensembles=allowed_ensembles),
    paths.strategies.MinusMoveStrategy(ensembles=minus_ens),
    paths.strategies.OrganizeByMoveGroupStrategy()
])

After that, you can see the list of ensembles that are expected to be used and the ones that are not expected to be used:

>>> scheme.find_used_ensembles()
{<openpathsampling.ensemble.IntersectionEnsemble at 0x115ad9630>,
 <openpathsampling.ensemble.IntersectionEnsemble at 0x115ad95f8>,
 <openpathsampling.ensemble.IntersectionEnsemble at 0x115ad95c0>,
 <openpathsampling.ensemble.MinusInterfaceEnsemble at 0x115ae9dd8>}
>>> scheme.find_unused_ensembles()
{<openpathsampling.ensemble.IntersectionEnsemble at 0x115ae9eb8>,
 <openpathsampling.ensemble.MinusInterfaceEnsemble at 0x115af88d0>}

As expected, we're not using the ensemble for the fake transition, nor the minus interface that is automatically associated with that transition (you can print the ensemble definitions to prove it.) The check_initial_conditions function only checks that there are trajectories for the ensembles that are used -- it exists to ensure that the move scheme won't crash on startup. By not making movers for those ensembles, you don't need trajectories for them either.

By the way, this is also how you really should do TIS when you do it with each interface as a separate run. Always use the same network, but use different move schemes. If you do it that way, all the default analysis works with very little effort.

HI David, Thanks for explaining the custom-made scheme. As of now it works and I will get back incase there is an issue

PS: For someone reading this in future, add engine = engine to minus and one way shooting strategy to make the above code block run MD simulation.

scheme.append([
    paths.strategies.OneWayShootingStrategy(ensembles=allowed_ensembles, selector=paths.UniformSelector(),engine = engine),
    paths.strategies.NearestNeighborRepExStrategy(ensembles=allowed_ensembles),
    paths.strategies.MinusMoveStrategy(ensembles=minus_ens, engine=engine),
    paths.strategies.OrganizeByMoveGroupStrategy()
])

Navigation Menu

Repository navigation

How to reject a path if not in ensamble #854

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

How to reject a path if not in ensamble #854

Description

Activity

dwhswenson commented on Sep 19, 2019

arjunwadhawan commented on Sep 20, 2019

dwhswenson commented on Sep 20, 2019

arjunwadhawan commented on Sep 21, 2019

arjunwadhawan commented on Nov 26, 2019

dwhswenson commented on Nov 26, 2019

arjunwadhawan commented on Nov 26, 2019

dwhswenson commented on Nov 26, 2019

arjunwadhawan commented on Nov 28, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions