Dear, Dr Axel,
I’ve been trying to modify the bond_harmonic.cpp to change r0 in the original formula to the distance of the initial state of the atom.
I download the LAMMPS source with git, release branch, on ubuntu 20.04, compile lammps by make.
In the example I ran,the initial configuration of the atom was generated by the matlab program,read into the input file with read_data.
I add the following lines of code to compute() in bond_harmonic.cpp:
And then it comes up ‘segmentation fault’. There is no error if you comment out this code.
I used GDB and valgrind to check for bugs, which showed that reading x0 was invaild.
I don’t know how to fix this problem.
And if I want to get the initial state of the atom how do I get it?
How can I achieve my goal?
If you check atom.h you will see that the x0 array is used by the PERI package – I don’t know what it stores, but LAMMPS does not store each atom’s initial position by default. (It would take more memory, for a piece of data rarely used in most MD calculations.)
In your position I would look at fix spring/self which also stores initial coordinates. Once you understand what it is doing you should be able to implement something similar for your fix.
This is a very bad idea for too many reasons to list them all here.
Please just use fix spring/self, which was written to do exactly what you want to achieve.
Otherwise, you would have to add immobile dummy particles that have no pairwise interactions with each other and other particles (you can exclude them from the neighbor list to achieve that) at the initial location of the original atoms and then define bonds between those dummy atoms and their originals. Then you can use bond_style harmonic for those bonds. However, that is just implementing with more complexity and overhead what fix spring/self already does.
@akohlmey@srtee
I think @ZhaoyanHuang is looking to change r0 from a per-bond-type quantity to a per-bond quantity, defined as the initial bone distance, based on his code.
@ZhaoyanHuang, if that’s what you want, i think it will be nearly impossible to implement. How would you keep it consistent between runs? How would you store it in a data file or a restart file? And more importantly, why do you want to do this?
I think we all have misunderstood what the shown code example does. If the forces between bonded atoms are computed from their initial positions, then we actually have a constant per-atom force. As I said before, modifying the bond potential is not a good idea. The main problem for that is for the given application that atoms may migrate between subdomains, but the original positions remain in place. But regardless of how this is set up or what the intentions are, I strongly believe that there is no need to modify the source code, but instead make creative use of the existing features.
So if the forces of the bonds should be the forces between the atoms at the initial positions, that this means, that those forces will not change over the course of the trajectory and it is no longer required to walk the bond list, but just use the total forces for only those bonded contribution. That means, they could be stored and then applied with fix addforce.
The simplest way to do this would be to set up the initial geometry, define only the bonded force styles (use the corresponding “zero” styles for everything else), set up three custom dump files that only output the atom-ID and the fx, or fy, or fz force component and do a run 0 command. Those files need to be stripped off their header parts and can then be used to define atomfile style variables that will then contain the total added force in x-, y-, or z-direction to each atom based on the initial positions. Then those variables can be applied to the system with fix addforce on second input that sets up the same system but this time with the intended interactions.
A more sophisticated version would do the same thing within the same input, where the bond force is stored with fix store/state without writing it out and then the pair/bond/angle/dihedral or whatever interactions are changed to their proper settings instead of zero force dummies and then, atom style variables are defined for the three force components that read the force from fix store/state and apply it with fix addforce.
My goal is to implement the bond potential in one of my advisor’s papers, which is to change r0 to the initial distance of the atom:
E = K(r_ij - L0(ri,rj))^2, where L0 = |x0,ij| is the initial length of the bond.
He wrote the code bond_harmonic1.cpp, but it got wrong, it’s in his code to compute L0 the way I wrote it up here. So I went with his code and changed bond_harmonic.cpp.
I’m a programmer, and I don’t know much about molecular dynamics, so I might have made a stupid mistake.
Thank you for your suggestions, I will continue to revise
That is not the correct description, you want to set r0 to the initial length of the bond. That is also what it says in the text that you quoted. Those differences matter, when you are asking people.
you don’t know enough about molecular dynamics force fields to know what kind of data you need to access and how and particularly how to test it
you don’t know enough about LAMMPS, its data model and parallelization scheme to do what you want to do correctly
So to resolve these problems, you will need to find a suitable collaborator (or get the corresponding help from your advisor) that can either train you or will do the implementation for you. It is going to be next to impossible to resolve this through an online forum, unless somebody here feels this has enough hack value and does the implementation for you.
Luckily for you (perhaps), I do feel that this has “hack value” and also could be useful to other LAMMPS users, if done correctly. Please find attached below two files that implement a new bond style bond/restrain, that should implement what you described. It is similar to bond style harmonic, but the bond_coeff command only takes one argument, the force constant. This is lightly tested and should be compatible with parallel execution and restarts (it worked for a very simple test case).
If you compare the code to bond style harmonic, you should see that the amount of changes required are fairly minimal, since the storage of the initial state has been “outsourced” to an internal instance of fix property/atom, so all the nasty work moving the data with atoms as they migrate between subdomains and do the forward communications to update its data to be available with ghost atoms is done by that fix code and does not have to be repeated. Yet, there are a few subtleties that are important so that restarting and using multiple run commands works correctly.
A more efficient implementation would store the initial length of the bonds along with the bond atoms instead of recomputing it over and over again, but that would require implementing a new atom style and that would be too much hack for too little value.