marrink-lab / vermouth-martinize

Describe and apply transformation on molecular structures and topologies
Apache License 2.0
84 stars 37 forks source link

interaction removal #599

Closed csbrasnett closed 1 month ago

csbrasnett commented 1 month ago

This was partially addressed by #483 but I don't think completely. I've written a small script to rewrite topologies to aid visualisation.

One issue we've had is here:

https://github.com/csbrasnett/martini_vis/blob/86d7036635c65f6bf7cec2a54a31d6f868838cb9/vis_top_writer.py#L78-L83

ie. for the removal of constraints. This should looks through constraints, and if there's meta eg. ifdef, then the interaction gets removed. However, this doesn't seem to happen at all consistently. Comparing the constraints before and after for an example molecule, some seem to have been removed, some are still there, and some have had the associated atoms changed.

csbrasnett commented 1 month ago

To give an example, if I run:

a = copy.copy(block.interactions['constraints'])
for i in a:
    print(i)
for bond in block.interactions['constraints']:
    if bond.meta.get('ifndef'):
        print(bond)
        block.remove_interaction('constraints', bond.atoms)
print('\n')
for bond in block.interactions['constraints']:
    if bond.meta.get('ifndef'):
        print(bond)

I get the following output:

Interaction(atoms=[15, 16], parameters=['1', '0.335'], meta={'ifndef': 'FLEXIBLE'})                                                                                                                                
Interaction(atoms=[16, 19], parameters=['1', '0.412'], meta={'ifndef': 'FLEXIBLE'})                                                                                                                                
Interaction(atoms=[18, 19], parameters=['1', '0.293'], meta={'ifndef': 'FLEXIBLE'})                                                                                                                                
Interaction(atoms=[15, 18], parameters=['1', '0.404'], meta={'ifndef': 'FLEXIBLE'})                                                                                                                                
Interaction(atoms=[16, 18], parameters=['1', '0.470'], meta={'ifndef': 'FLEXIBLE'})                                                                                                                                
...

Interaction(atoms=[15, 16], parameters=['1', '0.335'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[18, 19], parameters=['1', '0.293'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[16, 18], parameters=['1', '0.470'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[30, 31], parameters=['1', '0.341'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[40, 42], parameters=['1', '0.300'], meta={'ifndef': 'FLEXIBLE'})
...

Interaction(atoms=[16, 19], parameters=['1', '0.412'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[15, 18], parameters=['1', '0.404'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[24, 25], parameters=['1', '0.270'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[40, 41], parameters=['1', '0.320'], meta={'ifndef': 'FLEXIBLE'})
Interaction(atoms=[41, 42], parameters=['1', '0.270'], meta={'ifndef': 'FLEXIBLE'})
...

so, the first block is what should be there to begin, and I've checked it against what's present in the input itp. The second block should be the same as the first, because it just goes through the interaction and finds anything with the meta 'ifndef', which is all of them in this case. But! It's only finding every other one for some reason. Then the third block is the leftovers from the ones that weren't found first time around, which is strange because the loop to find them is no different.

I've worked out a workaround for now, but it'd be nice to work out why this isn't working as expected.

pckroon commented 1 month ago

I think this might be a funny side effect of iterating over something that you're modifying. Try the following?

a = copy.copy(block.interactions['constraints'])
for i in a:
    print(i)
for bond in list(block.interactions['constraints']):
    if bond.meta.get('ifndef'):
        print(bond)
        block.remove_interaction('constraints', bond.atoms)
print('\n')
for bond in block.interactions['constraints']:
    if bond.meta.get('ifndef'):
        print(bond)
csbrasnett commented 1 month ago

ahhh thank you v much that's sorted it!