[lammps-users] Harmonic potential as a pair-wise potential

Hi LAMMPS Users

Does anybody know how to implement harmonic bonded potential as a pairwise potential with a cut-off. Do i have to modify the code for this??

Thanks

Best Regards
Sreekanth

If you want a pair style that will implement a harmonic spring which turns
off at some cutoff distance, you'll need to implement it yourself as
a new pair style, or tabulate it and use pair_style table.

Steve

Thanks Steve. I am a beginner with LAMMPS. Can you please suggest me any useful reference to implement new pair style and any helpful suggestions.

Thanks

Regards
Sreekanth

doc/Section_modify.html has instructions on this

Steve

Hi Steve

I have changed the following lines in Void PairLJCut::compute(int eflag, int vflag) to modify pair_style lj/cut for harmonic non-bonded interactions. Now, epsilon would correspond to K and sigma would correspond to r0 (equilibrium length). Please let me know if this is sufficient. In my simulation, i also want to use lj/cut along with harmonic non-bonded interactions. How can i give this modification a new name to use it as command in my input script.

rsq = delxdelx + delydely + delz*delz;
// Harmonic_non-bonded 03/10/10
r = sqrt(rsq);
// Harmonic_non-bonded 03/10/10
jtype = type[j];

if (rsq < cutsq[itype][jtype]) {
// r2inv = 1.0/rsq;
// r6inv = r2invr2invr2inv;
// forcelj = r6inv * (lj1[itype][jtype]r6inv - lj2[itype][jtype]);
// fpair = factor_lj
forcelj*r2inv;

dr = r-sigma[itype][jtype]; // Harmonic_non-bonded
rk = epsilon[itype][jtype]dr; // Harmoin_non-bonded
fpair = -2.0
rk/r; // Harmonic_non-bonded

Please comment.

Thanks

Regards
Sreekanth

For the latter part of my question regarding the change of name, I have added to files named pair_harmonic_cut.cpp and pair_harmonic_cut.h. I modified existing the pair_lj_cut.cpp and pair_lj_cut.h files accordingly. So now, i can use it as

pair_style harmonic/cut 1.5
pair_coeff * * k r0

please let me know

Thanks

Regards
Sreekanth

Looks fine to me. I presume you are debugging it.

Steve