Hello,

In the peridynamics PMB pair style source (pair_peri_pmb.cpp), when the scalar force state is calculated and added to the total force on the “node” it is first made negative:

pair_peri_pmb.cpp, line 242

“”"

stretch = dr / r0[i][jj];

rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;

if (r > 0.0) fbond = -(rk/r);

“”"

Here stretch is the strain (amount stretched from equilibrium), kspring is the spring constant, fbond is the bond force, and r is the magnitude of the total stretch

My question: why is fbond the negative of rk? Is this simply Hooke’s law: F=-kx?

Thanks

Hello,

In the peridynamics PMB pair style source (pair_peri_pmb.cpp), when the

scalar force state is calculated and added to the total force on the "node"

it is first made negative:

pair_peri_pmb.cpp, line 242

"""

stretch = dr / r0[i][jj];

rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;

if (r > 0.0) fbond = -(rk/r);

"""

Here stretch is the strain (amount stretched from equilibrium), kspring is

the spring constant, fbond is the bond force, and r is the magnitude of the

total stretch

My question: why is fbond the negative of rk? Is this simply Hooke's law:

F=-kx?

this question is rather pointless without also stating, whether this

is f_AB or f_BA. remember newton's third law?

axel.

Axel,

I understand you point. I should have been more specific but I’m trying to learn a pretty large code base so cut me a little slack

So, to be more specific, here’s the full code:

“”"

for (i = 0; i < nlocal; i++) {

xtmp = x[i][0];

ytmp = x[i][1];

ztmp = x[i][2];

itype = type[i];

jnum = npartner[i];

s0_new[i] = DBL_MAX;

first = true;

for (jj = 0; jj < jnum; jj++) {

if (partner[i][jj] == 0) continue;

j = atom->map(partner[i][jj]);

// check if lost a partner without first breaking bond

if (j < 0) {

partner[i][jj] = 0;

continue;

}

// compute force density, add to PD equation of motion

delx = xtmp - x[j][0];

dely = ytmp - x[j][1];

delz = ztmp - x[j][2];

if (periodic) domain->minimum_image(delx,dely,delz);

rsq = delx*delx + dely*dely + delz*delz;

jtype = type[j];

delta = cut[itype][jtype];

r = sqrt(rsq);

dr = r - r0[i][jj];

// avoid roundoff errors

if (fabs(dr) < 2.2204e-016) dr = 0.0;

// scale vfrac[j] if particle j near the horizon

if ((fabs(r0[i][jj] - delta)) <= half_lc)

vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +

(1.0 + ((delta - half_lc)/(2*half_lc) ) );

else vfrac_scale = 1.0;

stretch = dr / r0[i][jj];

rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;

if (r > 0.0) fbond = -(rk/r);

else fbond = 0.0;

f[i][0] += delx*fbond;*

f[i][1] += delyfbond;

f[i][2] += delz*fbond;

“”"

So, F=-k*displacement?

Axel,

I understand you point. I should have been more specific but I'm trying to

learn a pretty large code base so cut me a little slack

if you would be looking at a complex piece of code (and LAMMPS has

those, too), yes, i would be willing to cut you some more slack, but

this one is not a complex one.

So, to be more specific, here's the full code:

[...]

So, F=-k*displacement?

this isn't really more specific. which displacement of what are you

talking about here?

as my (final) advice, here is what i can read from the code, purely

based on looking at the code and some educated guesses. please note,

that my understanding of the peridynamics method itself is extremely

limited. however, this code seems rather straightforward:

for each atom i, you loop over its "partners", i.e. atoms that were

neighbors of i at the beginning of the calculation. if a partner

doesn't exist anymore, you continue with the next.

then you compute the (current) distance r between atoms i and j and

compute the difference dr between the current distance and the

distance at the beginning (r0). from this you compute the relative

extension or compression, stretch, and multiply it with the force

constant. finally you project on the x, y and z components via delx/r

dely/r and delz/r. the /r term is applied ahead of time for obvious

computational efficiency. as for the sign of fbond, you must keep in

mind in mind that dr can be both positive or negative depending on

whether the distance between atoms i and j has been reduced or grown.

from here on, i strongly suggest you take a piece of paper and make a

sketch of the situation and see what the sign of the (restoring) force

should be.

reading code correctly requires more than just matching individual

lines with a formula, that you recognize from somewhere.

axel.