Issue with Superquadric Curvature

Submitted by Narayani on Tue, 01/26/2021 - 15:16

Hello!
I have been trying to use the gaussian curvature for the force calculations between two superquadric particles. To do so , I have used the command :

pair_style gran model hertz/stiffness  surface superquadric gaussianCurvature yes tangential_damping off

However , the radius that it is using for the forces calculations is the same as volume equivalent radius. Similarly I tried with meanCurvature yes and with CurvatureLimitFactor as more too , but even then the volume equivalent radius is being used...

It says in the description that the radius of curvature is limited by the volume equivalent radius, but the dimensions of the body in question are such that the gaussian curvature is lesser than the volume equivalent curvature.

How do I make use of the gaussian curvature ?

I am new to the world of LIGGGHTS , please excuse me if this is silly mistake.

Thanks in advanced
Narayani

mschramm | Tue, 01/26/2021 - 19:25

Hello,
how are you currently tracking to see which radius is used?
Looking at surface_model_superquadric.h, ri and rj are calculated by taking the cbrt of the product of the shape values. However, ri and rj are only used to calculate the contact point.

The two methods are set to trigger which method "calc_curvature_coefficient" is used to obtain coefficients. These coefficients are then used to calculate the effective radius.
You could add some debug print statements to check that things are behaving as expected.

Narayani | Wed, 01/27/2021 - 10:43

Thank you Matt for your quick reply.
I am trying to find the contact forces between two ellipsoids . The forces between the particles are modeled as hertzian forces ... so i am using the hertz formula to calculate the forces and the stress energy by hand and i am comparing with the result from LIGGGHTS. They match only if i consider the volume equivalent radius.
Like I mentioned before I am using following command for specifying the nature of the forces as:
pair_style gran model hertz/stiffness surface superquadric curvatureLimitFactor more gaussianCurvature yes tangential_damping off

and I am calculating the stress stored with compute st all stress/atom virial .
You mentioned about using some print statements.... which files would i need to look at?
The calc_curvature_coefficient which i found in surface_model_superquadric.h is in the following lines
if(curvatureLimitFactor > 0.0) {
int curvature_radius_method = meanCurvature ? 0 : 1;
double koefi = particle_i.calc_curvature_coefficient(curvature_radius_method, contact_point_i); //mean curvature coefficient
double koefj = particle_j.calc_curvature_coefficient(curvature_radius_method, contact_point_j); //mean curvature coefficient
Like the check here for mean curvature coefficient , there is no such check for gaussian curvature in the file surface_model_superquadric.h
I tried putting the print statements to check if the mean and gaussian curvatures are being printed, but haven't got any print outputs.

mschramm | Wed, 01/27/2021 - 16:32

The reason that there is no check for the gaussian method is because it can only be one or the other. When the function "calc_curvature_coefficient" is called, a curvature_radius_method value of 0 will use the meanCurvature method, whereas a curvature_radius_method value of 1 will use the ggaussianCurvature method.

Where did you place the print commands?
Looking at your pair_style command and comparing it to the LIGGGHTS' example, try rearranging it to.
pair_style gran model hertz/stiffness tangential_damping off surface superquadric curvatureLimitFactor more gaussianCurvature yes

Narayani | Wed, 01/27/2021 - 16:46

Hello Matt!
Thanks a ton for the suggestion to have a print debug...
My supervisor suggested that i use a non zero value instead of more and it worked! So now the command looks like
pair_style gran model hertz/stiffness surface superquadric curvatureLimitFactor 10 gaussianCurvature yes tangential_damping off

perhaps this was what was intended in the documentation.
Narayani