How to determine the maximum contact force

richti83's picture
Submitted by richti83 on Tue, 08/20/2013 - 08:59

I've a question about computing the maximum contact force
I know I can get the pair-force with
compute 2 all pair/gran/local force
But How to calculate the magnitude ?
(I tried variable fmag atom sqrt(c_2[1]^2+c_2[2]^2+c_2[3]^2) but got ERROR: Mismatched compute in variable formula)

Next problem is how to find the maximum of fmag, I tried
variable fmax equal max(fmag)
but got "ERROR: Invalid special function in variable formula"

One solution seems to be compute reduce

compute maxX all reduce max c_2[1]
compute maxY all reduce max c_2[2]
compute maxZ all reduce max c_2[3]
variable fmag equal sqrt(c_maxX^2+c_maxY^2+c_maxZ^2)

BUT this would sum the pair-forces of different atoms I think

finaly I added a 4th column to compute_pair_gran_local

if(fflag)
{
array[ipair][n++] = fx;
array[ipair][n++] = fy;
array[ipair][n++] = fz;
array[ipair][n++] = sqrt(fx*fx+fy*fy+fz*fz);
}

and reduced this on scriptlevel with

compute 2 all pair/gran/local force
compute 3 all wall/gran/local force
compute fmagP all reduce max c_2[4]
compute fmagW all reduce max c_3[4]

I'm wondering if there is any simpler solution ..
Thanks,
Christian.

PaulWinkler's picture

PaulWinkler | Sun, 08/25/2013 - 19:48

Hi Christian,

yes, if you add the maximal force components you will get an unrealistic fmag, since the compute reduce special function calculate the results independed from other per atom vectors.

Dealing with variables in liggghts can be very hilarious. There are many reasons why your

variable fmag atom sqrt(c_2[1]^2+c_2[2]^2+c_2[3]^2)

can fail. If you just need the variable at just one timestep, you may test a work around with fix store/state, that will eliminate lots of errors in variable handling. First you have to store your values in a fix at the moment the fix is executed:

fix c2 all store/state 0 c_2[1] c_2[2] c_2[3]

Then you should be able to run:

variable fmag atom sqrt(f_c2[1]^2+f_c2[2]^2+f_c2[3]^2)

without problems. The values will not hold the current value if the force vectors change! You should use a value > 0 if you need to update the result frequently (e.g. dump steps).
Since there are a lot of reasons, associated to the delayed evaluation of equal-style variables in combination with the 'variable used between current runs' – issue, you can try evaluating the variable frequently with the thermo command. But this depends on the structure of the script and won't work without further testing.

Best wishes,

Paul

anandmds's picture

anandmds | Sun, 10/02/2016 - 13:52

Hey, I want to use the data computed by " compute pair/gran/local & wall/gran/local" in my output. When I used it directly in atom style variable, I was getting the error " Mismatched compute in variable formula". So I used " fix store/state" as suggested in the previous post to store the values in the current timestep and then use in the atom style variable to hold it and then finally dump it using the dump command. Now the error that has propped is "ERROR: Fix store/state compute does not calculate per-atom values (../fix_store_state.cpp:285)". Please help me to solve this issue. I need to use the data from pair/gran/local and wall/gran/local. How do i go about doing this ?

Thanks & Regards

Anand M.

Anand.M
M.S. Mechanical Engineering,
I.I.T. Madras