I am simulating a system where a polymer is translocating from within a spherical cavity to outside the cavity. I am interested in calculating the vz of the atom as it is just about to exit the hole.
The atoms of the polymer leave in a linear manner one after the another, with atom index 1 leaving first and so on.
In order to calculate the vz of the atom just before it exits, I invoke the following commands:
(I have defined the spherical cavity previously as region spherical)
group inside_sphere dynamic polymer region spherical
compute v1 inside_sphere property/atom vz
variable vel_sphere equal c_v1[1]
fix v2 all print 1000 â${vel_sphere}â append sphere_velocity.txt screen no
I am assuming as the beads move, they are allocated into the dynamic region in ascending order(as atom 1 leaves first followed by the other atoms in a linear fashion). Therefore in every timestep the first index of the compute performed on the dynamic group should belong to the atom of lowest index that still hasnât left the cavity.
The input script does seem to work but the problem is that the velocity is always outputted as zero. Can this be rectified with a different fix or compute?
variable vel_sphere equal c_v1[1] does not access the âfirst element in the groupâ but the atom with atom ID 1. So it will always access the same atom and values of per-atom computes that are not members of the compute group are set to 0.
This is documented behavior.
The major problem of what you are doing is that you are âthinking globalâ while LAMMPS âoperates localâ, i.e. LAMMPS has distributed data and all operations that âglobalizeâ data are complex, require additional communication (and thus hinder effective parallelization). Thus any operation of the âfirst of a group kindâ are always difficult to get right, never the default operation, and rarely efficient.
The best way to output per-atom data that needs to conform to some condition is with a dump file and using dump_modify thresh to select which atoms to output. Please note that per-atom variables can do rather complex selections (but not really âfirst of a groupâ) through multiplication or addition of logical operation.
In general, I would output selected output into a dump and the post process that data. With the recently added âyamlâ dump style, loading dump file data into python and then processing it is easier than ever.
I was really hoping on being able to pull off this analysis in the input script itself and save the additional hassle of post processing. The globalized and localising mode of operating which you have explained was not immediately apparent to me from the documentation. Thank you for your suggestions!
If you consider post-processing a hassle, you are in the wrong business. More importantly, you are discarding the fact, that the effort to do the entire processing in a single input may be (much) more than just make a smart and simple selection of the data for post-processing and then write a simple tool for further analysis. Technically, with the forward and backward python interface in LAMMPS there are very few limits to what you can do. It is mostly a matter of your programming skills and determination.
There may be a rather convoluted way to approach what you are looking for,
you create an atom style variable that multiplies a rmask() or grmask() with the suitable arguments with the atom ID
you then use this function to generate chunks and use the compress yes option. Chunk ID 1 should then only contain the atom with the smallest atom ID that conforms to your region (and group) requirements
what options you have from here on to access the data you need from that info, you will have to dig out of the docs yourself.
Instead of manipulating the per atom vector, I have tried a simpler solution which utilizes the count function. In a 100 bead polymer if there were 38 beads found in the sphere, this would impy (100-38)=62 beads have exited. The index of the bead that just exited would be 62 and the index of the closest bead inside would be 63.
Hence the required index would be simply 101-count(inside_sphere)
I have tried to incorporate this logic into lammps using the following commands
:
group inside_sphere dynamic polymer region spherical
compute v1 inside_sphere property/atom vz
variable sphere_loc equal (101-count(inside_sphere))
variable vel_sphere equal c_v1[$${sphere_loc}] (I intended on using the single dollar sign instead of double, as the web formatting removes it)
fix v2 all print 1000 â${vel_sphere}â append sphere_velocity.txt screen no
However I am still not able to get the desired results as the outputted value is still zero.
the chunk/atom with compress would ideally work as it eliminates atom ids that do not belong to the chunk directly. However I am not able to find if there is a further fix/compute that operates on chunks and gives per atom data similar to the property/atom compute
This cannot work, since the â$â variable expansion happens immediately during preprocessing of the input file. Thus you get the current value at the time when that input line is processed.
Weâre back at the point that you are trying to use the LAMMPS script language as a full-featured script language, which it is not. I have now multiple times pointed out that this is best done with calling out to python in combination with using the LAMMPS python module to retrieve and process the data (taking into account parallel processing, too). And that leads us back to the point that you have to go through much more hoops and effort compared to post-processing.
You are missing the point of this approach. Due to the way the chunks are defined, there is only one atom per chunk. Furthermore, the chunk ID is a per-atom property, thus now you could, for example, use that in combination with dump_modify thresh to dump only the atoms with chunk ID 1, which should be the desired atom. You can also define atom style variables that have your desired per-atom property multiplied by a logical that compares the chunk ID to the value 1. Then the resulting variable will be zero for all atoms but the desired one and the âsumâ reduction would be the desired value of the selected atom. But be aware that reduction operations are very inefficient and interfere significantly with parallel efficiency.
In short, I think you are on the warpath to make your life more miserable than necessary because of a misconception of what is âeasyâ and what is âdifficultâ. And now you keep investing effort into something without success where you already could have the results in your hands for some time.
I am very thankful to you for your inputs again, could you please clarify if there is any way in using the instantaneous value of the variable in the array index c_v1[]. Is there an alternate syntax that can be used to effectively obtain c_v1[sphere_loc] instead of using the dollar sign.
With regards to the chunk ID, I am unable to see a way which allows you to label the atoms individually by their id or index, if i set it to molecule or type, all atoms belonging to the region are labelled as 1 and the rest as zero.
Suppose I did have a command that could do it, then I believe I may be able to access it through the chunk/property command which outputs the chunk ID prior to compression
I believe in the case of my system which is relatively simple linear translocation obtaining the instantaneous value of the starting index of the atom in the region(which can be obtained) and putting it the array index of the compute vector would solve my problem
No. I have explained multiple time why LAMMPS chooses not to provide such features as they are inherently very inefficient given the data model.
If you cannot follow what I am suggesting, then I cannot help you. I think this is straightforward with just a little insight into how âLAMMPS worksâ. Overall, my patience with people who prefer to ignore my suggestions and instead insist on doing things differently and then keep running into problems at every turn is limited and in your case now exhausted. Of course you are entirely at liberty to keep following your path. Just donât expect me to make any suggestions on this proverbial âdead horseâ.
I am very thankful to you for your suggestions and I do intend on following them. I only intended to obtain clarification from a pedagogical perspective on this procedure and to see if it is doable on lammps, which would aid me in my understanding lammps object manipulation.
I would like to express my sincere thanks for your patience and taking your time off towards resolving our queries and offering suggestions to help in our efforts.
Please note that this is not how it comes across. Rather than as holding tightly on to a preconceived notion and expectations of how things would have to work. If you want to learn how a software works that is not an effective strategy.
On the contrary, in my personal experience it always helped me to first get a sense of the âmindsetâ of a software (if you can describe it that way) and review documentation and examples and try to understand what they do how. That gives you a better sense of how things work than sticking to figuring out how to implement a specific workflow. Once you are familiar with the available workflows and data models, then it is much easier to construct inputs and design strategies to obtain a viable solution. In that process there are two major pitfalls: 1) premature optimization, i.e. going through a lot of hoops in trying to avoid (perceived) performance bottlenecks without confirming by some careful benchmarks whether those are relevant, 2) trying to solve all problems with the same tool (particularly when you are not yet very familiar with that tool). It is often more effective to use a combination of tools, each more optimized for a specific task (this is a core part of the so-called âunix philosophyâ where its power comes from the combination of many small tools that can just one task very well).
I sincerely appreciate your suggestions. I would like just comment on one inconsistency:
group inside_sphere dynamic polymer region spherical
compute v1 inside_sphere property/atom vz
variable sphere_loc equal (101-count(inside_sphere))
variable vel_sphere equal c_v1[$${sphere_loc}] (I intended on using the single dollar sign instead of double, as the web formatting removes it)
fix p1 all print 1000 â$${sphere_loc}â append sphere_pos.txt screen no (I intended on using the single dollar sign instead of double, as the web formatting removes it)
fix v2 all print 1000 â${vel_sphere}â append sphere_velocity.txt screen no
Executing this I notice that the updating of sphere_loc when called with the dollar sign is correctly updated when in the fix print command but the same variable fails to update when used as an index to a vector obtained from a compute and only sticks to the value it held during preprocessing almost like an internal style variable.
It would be enormously useful(amidst parallelization concerns) if you could allow for the updating of variables to also hold when called in the index of a vector as it may be very useful in other tasks.
This is documented behavior.
Please note that the dollar sign in fix print is inside quotation marks and thus not subject to expansion during pre-processing. So when you do c_v1[${sphere_loc}] the variable command never gets to see ${sphere_loc} while fix print does.
There are more cases where variables are accessed differently. Sometimes with a v_ prefix, sometimes not (when a variable is required and thus the v_ prefix to tell it apart from a reference to a compute or fix is not needed). Similar things happens to fix or compute references.
As I have mentioned multiple times, LAMMPSâ script language was never intended to be a fully formed script languages. and because most of these behaviors are not part of a single script interpreter, but rather implemented as local argument parsers in each individual command, some inconsistencies are to be expected, especially considering that LAMMPS is not authored by a single person and strictly reviewed by a committee or similar (we are far too few people with too little time for that). This is made worse by the fact that some parts of LAMMPS are very old and it is thus very difficult to make changes for consistency without breaking a lot of peopleâs inputs.
In recent times, and specifically thanks to additional available time due to the pandemic, we have started internal refactoring to increase code sharing and make behavior more consistent, but with a code base this big (> 1 million lines), it is a massive undertaking. Considering the rate at which new code is added, doesnât make this easier, even though we are now much more aware of potential inconsistencies and try to avoid issues ahead of time based on previous experiences.
I have replaced ${} with the v_ syntax for variables and I got it working now. I am once again very thankful to you for your time and patience on the matter.
as it may be very useful in other tasks.