removing interactions (again)

[csbrasnett]

I think this is a slight rehash of #482 but it really doesn't seem to be working.

In short, to get cg_bonds working, I'm trying to write a script that takes an input .top file, reads the itps contained, removes all the interactions apart from the bonds, and write some 'visualisation' itps, and a corresponding .top with the full paths. (cg_bonds doesn't like virtual sites, position restraints, etc., so they need to be taken out)

script goes like:

from vermouth.gmx import read_itp, write_molecule_itp
from vermouth.forcefield import ForceField
import os

if __name__ == '__main__':
    #get files in the present directory
    files = os.listdir(os.getcwd())
    
    #get the topology file
    top = os.getcwd() + '/topol.top'
    
    #read the topology file, find the header lines which are not core martini itps
    with open(top) as f:
        topol_lines = f.readlines()
    
    topol_itps = [a.split('"')[1] for a in topol_lines if "#include" in a and "martini" not in a]
    
    #for each molecule in the system, read in the itp
    d = {}
    for i,j in enumerate(topol_itps):
        with open(j) as f:
            d[i] = f.readlines()
    
    #read the molecules into the forcefield
    ff = ForceField('martini3001')
    
    for i in d.keys():
        read_itp(d[i], ff)
    
    #for each of the molecules
    for i in ff.blocks:

        #get all interactions
        all_interactions = list(ff.blocks[i].interactions.keys())
        print(all_interactions)
        #for any interaction that isn't a bond
        for j in all_interactions:
            if j != 'bonds':
                # get the interaction
                ints = ff.blocks[i].interactions[j]
                # and remove it
                for k in ints:
                    ff.blocks[i].remove_interaction(j, k.atoms)
    
        #write out the molecule    
        mol_out = ff.blocks[i].to_molecule()
        mol_out.meta['moltype'] = i+'_vis'

        with open(i+'_vis.itp', 'w') as outfile:
            write_molecule_itp(mol_out, outfile=outfile)

However, it doesn't actually seem to modify the topology in any way. The mysterious thing is if I interactively use remove_interaction(etc), then it does remove the interaction specified. But it can't do it in a loop for every one in one go. So I'm not sure what's going on. I can't tell as well if removing all the interactions within a directive causes issues as well - would there be a way of removing a whole directive in one go, rather than iterating through the whole thing?

python 3.8.7
vermouth==0.9.6

files in the directory:

topol.top
martini_v3.0.0.itp
martini_v3.0.0_solvents_v1.itp
martini_v3.0.0_ions_v1.itp
molecule_0.itp

[pckroon]

...
    for molname, block in ff.blocks.items():

        for interaction_type in block.interactions:
            if interaction_type != 'bonds':
                del block.interactions[interaction_type]
    
        #write out the molecule    
        mol_out = block.to_molecule()
        mol_out.meta['moltype'] = molname+'_vis'

        with open(molname+'_vis.itp', 'w') as outfile:
            write_molecule_itp(mol_out, outfile=outfile)

should also work (answering your second question). I'm not sure why your solution wouldn't work though.

[csbrasnett]

Thanks! That seems to work with

...
for interaction_type in list(block.interactions):
...

otherwise I get a RuntimeError: dictionary changed size during iteration. But yes, no idea why the first solution doesn't work.

[pckroon]

Oh right, good catch.