Dear all ,
I am trying to understand AQS(athermal quasistatic shear deformation) protocol by replicating the results of the following paper (10.1103/PhysRevLett.93.016001) to understand the process of athermal quasi-static shear deformation. In the following paper they are using a modified version of the formula used in pair_style harmonic/cut. I have made changes in pair_harmonic_cut.cpp to modify the formula. I am currently using lammps-patch_7Feb2024_update1 version of lammps. I have attached the edited pair_harmonic_cut.cpp file and pasted the snippet of the part of code where I have done the changes.
const double rsum = radius[itype] + radius[jtype];
const double r = sqrt(rsq)/rsum;
double PairHarmonicCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double &fforce)
{
double *radius = atom->radius;
if (rsq >= cutsq[itype][jtype]) {
fforce = 0.0;
return 0.0;
}
const double rsum = radius[itype] + radius[jtype];
const double r = sqrt(rsq)/rsum;
void PairHarmonicCut::born_matrix(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double &dupair,
double &du2pair)
{
double *radius = atom->radius;
double rsum = radius[itype] + radius[jtype];
double r = sqrt(rsq)/rsum;
After applying these changes, I am trying to minimise the initial configuration using conjugate gradient method and the simulation is stopping with a message lineseach alpha is zero, instead of stopping with a force tolerance criteria. I have attached the .cpp file and input file.
Can anyone help me understand the mistake I am doing to update the .cpp file which is causing the issue during energy minimisation.
in.data_file (1.2 KB)
pair_harmonic_cut.cpp (10.3 KB)