I’m trying to compute a pairwise quantity that is dependent on it’s value in the previous timestep and current time step. Precisely, I would like to compute the work done by friction at a contact,
W=integral(F_{friction}.ds).
I know my ‘ds’ which is the transverse velocities multiplied by dt. I also know F_{friction} at the present timestep. I would like to store the W variable in an atom-type variable as this is how it was included by a colleague. The problem is right now I think he is computing the integral using:
W_new=W_old+F_{friction}.ds.
I don’t think this a good method of performing numerical integration. Therefore I want to try and use the trapezoidal rule to perform integration. My concern is if I were to do this, I would have to use F_{friction} from not only the current timestep but the previous timestep as well. I believe shear vector within compute is where I store my history dependent quantities.
So my question is, how can I access the frictional forces from the previous timestep by storing it in the shear vector. I tried to naively change
shear = &allshear[8*jj];
within ::compute by increasing 8 by the desired number of variables I want to store. I did manage to compile this code but I ended up getting a segmentation fault when running the simulation. I feel like I must edit a few other files but I am kind of lost. Not the best C++ programmer out there so it isn’t obvious to me what I must edit.
What you are trying to do is not as straightforward as you think and there is likely no way to do it unless you improve your skills in reading, understanding, and writing C++ code.
The main question for now is: is this F_{friction} term a per-pair or a per-atom property? That determines what you need to modify.
In the first case you could indeed expand the storage used for storing shearing data, but that is maintained by an internal fix implemented in src/fix_neigh_history.cpp/.h
In the second case, you can add an array to the pair style and need to grow it from within and also implement the necessary pack/unpack routines and adjust the settings so that the data is migrated with the atoms when they migrate between MPI ranks, i.e. subdomains. Then you may also need to implement a forward communication in order to communicate the properties to ghost atoms.
for all of this there is no simple and generic step-by-step instruction of a “do this, then that” kind, since details are different in different cases. we are trying to expand the programmer guide section of the manual to include some more high-level information, but in the end any kind of manipulations and modifications of this kind require a sufficient understanding of some of the lower level LAMMPS code and that in turn is best gained from reading existing code and making sense of it. We can answer specific questions on specific details, but to give you a more complete tutorial would significantly more work than just implementing what you are looking for.
Thanks as always! I am learning c++ through LAMMPS but it is a very steep learning curve having to modify it although i’ve done it in the past using your guidance
So yes, F_{friction} is a pairwise quantity. I’m using an old version of LAMMPS (Nov 2016) and looking at the src files I do not see fix_neigh_history.cpp but I do see a fix_shear_history.cpp. I think this is what I must modify as I do see arrays only allocating for 8 terms in shear. I bumped that upto 11 in pair_gran but not in shear_history so I guess I will look at modifying this code.
Do you think more modifications will be needed in other files to allocate more memory?
I have already expressed in a previous discussion that I consider your approach to modify code based on a very old version of LAMMPS and backporting features from more recent versions a very bad idea. We have moved on and improved LAMMPS quite a bit in recent years to make it easier to modify and maintain.
Thus I have no interest in providing assistance with versions of LAMMPS that are this old and having to deal with all of their problems and shortcomings we have eliminated since.
You also are inheriting all the bugs that we have fixed since (and there were quite a few that were discovered and fixed(!) thanks to more rigorous testing and static code analysis) and are missing on the improved APIs and documentation that only apply to the current code. That you do this with limited C++ experience is making your life even harder by doing something that is harder than what you are trying to avoid doing.
That said, open source is all about the liberty to modify it as you see fit and thus there is nothing stopping you to continue on your path. Just don’t expect input from me.
In addition to the issues Axel has mentioned, I am not sure if trapezoidal integration is in fact that much more accurate than a simple sum over time steps, especially as to be worth the effort put in.
Have you tried dumping out the needed raw quantities for, say, 100 time steps, and comparing your trapezoidal integration algorithm to the current results? If the accuracy is good enough, it’s good enough.
