# Function in LAMMPS C++ code to calculate/call the original positions of each atoms

I am trying to use simple harmonic springs to keep their geometry in crystalline material. It looks like an atom will interact with their neighbour atoms by spring within the cut-off distance r_cut-off via this equation: E = k(ij)(r(ij) -r_0(ij))^2 within cut-off distance r_cut-off, where k(ij) is stiffness constant , and r(ij), r(ij)_0 are the current distance and original distance, respectively, between atoms i and atom j. Base on the last discussion in LAMMPS USER forum http://lammps.sandia.gov/threads/msg11009.html in which advised that it should be started from modifying lj/cut force field and create a new one (harmonic.cpp & harmonic.h). The authors proposed that the k, r_0, r_cut-off input in this new harmonic/cut was respect with epsilon, sigma, and r_cut-off for lj/cut. However, r_0 need to be input each an i-j pairwise, how I can deal with many different values of r_0 for different pairs of i-j automatically instead of input this coefficient manually. For example, for a simple cubic with size a, an atom at a certain edge can interact with other atoms. The distance between an atom and its first neighbour is a, but with the second neighbour is sqrt(2)*a, …. Thus, r_0, thought, should be calculated from equilibrium state ( initial step) for each pair of i-j interaction. Can anyone let me know which function in LAMMPS C++ code can help me calculate or call the original positions of each atoms ?

I really appreciate your help !

Thi

Why don't you simply use fix spring/self?

Axel.

Because he needs a pair interaction.

If your system is as simple as a cubic one why not defining the bonded atoms manually via the bond command and then replicate to make it as large as desired. Shouldnt that work?
Carlos

If your system is as simple as a cubic one why not defining the bonded atoms
manually via the bond command and then replicate to make it as large as
desired. Shouldnt that work?

nope. replicate cannot handle bonds across periodic boundaries, as
would be needed in this case, cubic or otherwise.

axel.

That’s one for me to keep in mind.
Carlos

Because he needs a pair interaction.

who says so? the description said "to keep the geometry in a crystal".
there is a much simpler way to solve this. rather than explaining in a
long-winded tirade that LAMMPS doesn't work as expected and that it
would be much more helpful to first familiarize yourself with the
overall code design before asking such a questions, i thought it might
be more productive to suggest a simpler alternative. (...and i get
less nastygrams about being so nasty to beginners ;-)).

that being said, there is already a method to have arbitrary harmonic
pairwise interactions via the "list" pair style. all that would be
required is to write a little script/program that lists all pairs that
should be handled, computes the equilibrium distance and then prints
out the individual parameters. i still believe that this is likely not
a good solution for the indicated problem, and was hoping that my
question would result in prompting the necessary explanation why
either approach would be suitable or not.

axel.

Because he needs a pair interaction.

who says so? the description said "to keep the geometry in a crystal".

Yeah, maybe. He asked for it and he only really knows. To tether atoms
down to a crystal spring self works well.
A simple pair interaction like he asked for would be trivial to set up
with pair style table.
No programming necessary, just tabulate energy and force in a text file.
Daniel

Many thanks to you Axel !
Thi