How to get forces (and other quantities) for a set of tags

[trvsst]

Hi all,

I am trying to get forces and energies for several subsets of particles in my simulation. The concrete code is somewhat long, so I will delineate how is done in a general way, which hopefully will be useful for others as well. I am using hoomd 4.1.0 (and somewhere in the code gsd 3.1.1, but that should not matter)

I create a list of filters with hoomd.filter.Tags

list_filters=[hoomd.filter.Tags([0,1,2]), hoomd.filter.Tags([3,4,5])]

then I create a logger

logger_local = hoomd.logging.Logger()
logger_local.add(lj, quantities=['energies', 'forces'])

and finally, I create a new gsd file for each filter

list_writer=2*[0]
for ind in range(2):
fl = 'my_file' + str(ind)
list_writer[ind] = hoomd.write.GSD(filename=fl, trigger=trig, mode='ab', filter=list_filter[ind])
list_writer[ind].logger = logger_local

This is it.

I would expect that running the simulation this way, my_file0 and my_file1 would contain the energies and forces for particles [0,1,2] and [3,4,5] respectively. What I find instead is that both files contain the energies and forces for all the particles in the system, as if I had set filter=hoomd.filter.All() in hoomd.write.GSD.

Is there something incorrect in the logic of this setup?

[joaander]

The filter attribute of GSD applies to the particles section of the GSD file. GSD writes all logger provided quantities exactly as they are provided. To log a subset of particle energies, provide the expected array to the logger:

logger[('energies')] = (lambda: lj.energies[0,1,2], 'sequence')

[trvsst]

Thanks a lot Josh,

But then, I do not need a filter to get the energies and forces for a given set of particles?

My logic was that logger.add would compute the local energy and forces

If I create the filter
list_filters=[hoomd.filter.Tags([0,1,2]), hoomd.filter.Tags([3,4,5])]

then filter=list_filer[ind] would return the corresponding quantities for the given particles on each filter
list_writer[ind] = hoomd.write.GSD(filename=fl, trigger=trig, mode='ab', filter=list_filter[ind])

I see that my logic is incorrect because from what you say, if gsd receives the entire set of energies and forces it will output as it received it, but then if filter=list_filter[ind] gives you the same as filter=hoomd.filter.All(),
what is the purpose of filter in hoomd.write.GSD, it it cannot modify what receives from logger?
isn't list_filer[ind] is acting on particles?

[trvsst]

Also, given a filter object, is it possible to transform it back to a set of tags, such as, for example

np_array = filter_object.tags

and np.array is an array with all the tags that belong to that filter object?

[joaander]

list_writer[ind] = hoomd.write.GSD(filename=fl, trigger=trig, mode='ab', filter=list_filter[ind]) should not give you exactly the same output as All(). If your filtered GSD file contains particles/position, particles/orientation, particles/... for all particles (ignoring the filter), the please submit a bug report with a complete script (and minimal) that demonstrates the problem.

If your intention is to store only logged quantities, then use the Null filter with write.GSD or use write.HDF5Log.

The __call__ method of a filter returns the tags that match the filter. However, this method is only required to return the matches on the local MPI rank. Feel free to submit a pull request that adds a global_tags property which evaluates the set union across all ranks.

[trvsst]

Josh,

I include two files: It is mostly boilerplate code except for a few lines where I provided a lot of comments inside the file. Still, it is significant amount of code to place inside this message, so I provide a link instead. (It is not allowed to attach *.py files).

The first file runs a simple lj simulation: if at the top of the file you set
method_write=True implements filters as discussed, and it does not work: gives same as filter=hoomd.filter.All()
method_write=False implements the sequence as you suggested, and it does what I intended with method_write=True
https://www.dropbox.com/scl/fi/rf8v55u72ickorw2s1qhp/analyze_filter_lj.py?rlkey=pjhf5lf9itftms4hq2b0r2i1x&dl=0

This is the file I used to analyze the files, does not need to be modified whether you use method_write=True or False
https://www.dropbox.com/scl/fi/j71g2zzzpe3d3not1grx8/test_filter_lj.py?rlkey=eqmdbe37kiiuqc1ftquu674c3&dl=0

Does the code with method_write=True works as you expected or it is a bug? If it is a bug, what would you expect that the code described in method_write=True should give?

[joaander]

No, this is not a bug. Your test checks the shape of fl_data[ind_fl][ind].log[val], not the shape of fl_data[ind_fl][ind].particles.position.

[trvsst]

Thanks a lot Josh,

This was really very helpful. Thanks a lot. Now I think I understand the logic.

I do have a couple of questions, where I have spent a lot of time and cannot figure it out. I have prepared these examples. This is much closer to what I am trying to implement. The code does the following:

With a simple lj creates a number of filters (the example is set to 10) each of them contains energies and forces for 3 particles in sequence (see line 70). Then I store the energies and forces for each filter in a separate file, lines 70-117 with lots of comments.

Here are the problems:

question 1: If I replace the commented lines 92:93 with 94:95 in test_filter_lj (basically, using getattr instead of writing lj.(energies, forces) explicitly) the attributes are written incorrectly into the gsd file, energies are replaced with the forces, I get two force attributes instead of energies, forces.

question 2: Even when I run the code with explicit attributes, it does not work as I expected. I generate the 10 files, and when I analyze the 10 files I get the exact same result in all of them. If there is a trivial error, I apologize, but I have spend a ton of time analyzing it and I honestly think there is no trivial error. You can run the code and use analyze_filter_lj to see what there is inside those 10 gsd files, the code outputs only the last configuration. You can see that the result in all files are identical.

https://www.dropbox.com/scl/fi/j71g2zzzpe3d3not1grx8/test_filter_lj.py?rlkey=eqmdbe37kiiuqc1ftquu674c3&dl=0

https://www.dropbox.com/scl/fi/rf8v55u72ickorw2s1qhp/analyze_filter_lj.py?rlkey=pjhf5lf9itftms4hq2b0r2i1x&dl=0

[joaander]

You change the value of tags in the loop. The last value will be kept and the previous lambdas will reference it. This is discussed at length on stack overflow. Personally, I would use only 1 log file or create a lightweight class to store a copy of tags and lj to create the filtered logged quantities.

[trvsst]

Thanks a lot! I am sorry this last question was not a strict hoomd question, this is the first time in a very long time that I used lambda functions... It would have taken me a lot of hours to find this one, so thanks for saving the day.
Everything works as expected and now that I see the solution is obvious!

[joaander]

I posted a new answer that shows how to do this without lambda for any future users that read this.

[joaander]

The filter attribute of GSD applies to the particles section of the GSD file. GSD writes all logger provided quantities exactly as they are provided. To log a subset of particle energies, provide the expected array to the logger:

import copy

class LogSubset:
    def __init__(self, lj, tags):
        self.lj = lj
        self.tags = copy.copy(tags)
    
     @property
     def energies(self):
        return self.lj.energies[self.tags]


log_subset = LogSubset(lj, [0, 1, 2])
logger[('energies')] = (log_subset, 'energies', 'sequence')

[trvsst]

Thanks a lot Josh. This is very nice!

[trvsst]

others may find this useful,

this snippet implements the discussion above using lamda functions in a way that overcome the problem discussed above

# number of filters
n_filters = 10

# you can write here other quantities
log_quantities = ['energies','forces']

# create n filters, each of them containing three particles only
list_fl=n_filters*[0]
for ind in range(n_filters):
    list_fl[ind]=hoomd.filter.Tags([3*ind,3*ind+1,3*ind+2])

# define a logger to store the local quantities  one for each filter
l_logs = n_filters*[0]
for ind in range(n_filters):
        tags=list_fl[ind].tags
        l_logs[ind]= hoomd.logging.Logger()
        for qtes in log_quantities:
            l_logs[ind][qtes] = (lambda x=qtes,y=tags: getattr(lj,x)[y[:]], 'sequence')

# save the results for  each filter in a gsd file
for ind in range(n_filters):
    fl = 'my_file' + str(ind)+'.gsd'
    list_writer[ind] = hoomd.write.GSD(filename=fl, trigger=trig, mode='wb', filter=hoomd.filter.All())
    list_writer[ind].logger = l_logs[ind]