So i am caculating bending force from mean curvature but after some time the forces are too large so they beak bonds. I have checked the force calculation and it seems correct. Is there something wrong the way I add these forces to lammps?

LAMMPS LAMMPS Development
1

This is my code to add forces to lammps

LAMMPS_NS::LAMMPS *lmp;
// custom argument vector for LAMMPS library
const char *lmpargv {“liblammps”, “-log”, “none”};
int lmpargc = sizeof(lmpargv)/sizeof(const char *);
// explicitly initialize MPI
MPI_Init(&argc, &argv);

// create LAMMPS instance
lmp = new LAMMPS_NS::LAMMPS(lmpargc, (char **)lmpargv, MPI_COMM_WORLD);
std::cout << "LAMMPS version ID: " << lmp->num_ver << std::endl;

writeInitialConfig(“init_config.txt”);
lmp->input->file(“infile.lmp”);

int natoms = static_cast<int>(lmp->atom->natoms);

int nlocal = lmp->atom->nlocal;
double **f = lmp->atom->f;
double **x = lmp->atom->x;
int *id = lmp->atom->tag;
for (int step = 0; step <= FTIME; ++step) 
{
    std::unordered_map<int, Vert*> id_to_vert;
    for (int i = 0; i < total_verts; ++i)
        id_to_vert[verts[i]->id] = verts[i];
    
    for (int i = 0; i < nlocal; ++i) 
    {
        int atom_id = id[i];
        Vert* v = id_to_vert[atom_id - 1];
        v->setCoord(x[i][0], x[i][1], x[i][2]);
        v->force[0] = v->force[1] = v->force[2] = 0.0;
    }


    if (step % OUTPUT_FREQ == 0) 
    {
        trials = 0;
        do
        {
            prob = udist();
            ran1 = (int) (udist() * total_edges) + 1;
            cout << ran1 << endl;
            if(prob<0.5)
            {
                flipSuccess = edgeFlip(edgef[ran1]);
                Nflip += flipSuccess;
            }
            else
            {
                flipSuccess = edgeFlip(edgeb[ran1]);
                Nflip += flipSuccess;
            }
            trials++;
        } while (!flipSuccess && trials <=MAX_TRIALS);
        
        if(flipSuccess)
        {
            cout <<"Vertices fliped : "<<vertices[0] <<" "<<vertices[1] <<" "<<vertices[2] <<" "<<vertices[3] <<endl;
            stringstream group_cmd;
            group_cmd << "group my_group id " << vertices[2] << " " << vertices[3];
            lammps_command(lmp, group_cmd.str().c_str());
            lammps_command(lmp, "delete_bonds my_group multi remove");
            lammps_command(lmp, "group my_group delete");

            stringstream bond_cmd;
            bond_cmd << "create_bonds single/bond " << 1 << " " << vertices[0] << " " << vertices[1];
            lammps_command(lmp, bond_cmd.str().c_str());
            lammps_command(lmp, "write_data restart.txt");
        }
    }

    lmp->input->one("run 0 pre yes post no");

    // Caculate bending Forces on each point

    for (int l = 0; l < total_verts; l++) 
    {
        if(verts[l]->boundary)continue;
        for (int j = 0; j < verts[l]->num_neighbors; j++) 
        {
            int n_id = verts[l]->vertex[j]->id;
            if(verts[n_id]->boundary)continue;
            bendingForce(n_id, l); 
        }
    }

    for (int i = 0; i < nlocal; ++i) 
    {
        int atom_id = id[i];
        Vert* v = id_to_vert[atom_id - 1];
        f[i][0] += v->force[0];
        f[i][1] += v->force[1];
        f[i][2] += v->force[2];
    }

    lmp->input->one("run 1 pre no post no"); 
    if (step % OUTPUT_FREQ == 0) outputVTKFile(step);
}

}
Also this is my input script
dimension 3
units lj
atom_style bond
boundary f f f

read_data init_config.txt extra/bond/per/atom 10 extra/special/per/atom 25
reset_timestep 0

bond_style harmonic
bond_coeff 1 500 1.2

pair_style lj/cut 1.12246
pair_coeff 1 1 1.0 1.0 1.12246
pair_modify shift yes
special_bonds lj 0.0 1.0 1.0

timestep 0.0005

neighbor 1 bin
neigh_modify every 1 delay 0 check yes
comm_modify cutoff 3.0

min_style cg
minimize 1e-6 1e-8 10000 100000

fix 1 all nvt temp 1.0 1.0 1.0

2

This cannot work. If you look at the flow of control of a step in LAMMPS, you’ll see that the forces are cleared at the beginning of the force computation, so your addition will be ignored.

The proper way to add forces during a run from an external code would be to use fix external.

Also, but I cannot be for certain without knowing more details and having a closer look (for which I don’t have the time right now), this doesn’t look like it would be properly compatible with running in parallel or even periodic boundaries.

3

What, exactly, are you trying to do here? (As a genuine question. Since your bendingForce method is opaque to us we cannot have any hope of making suggestions.)

And are you sure there is no way of turning your desired added force into either a bond_style, an angle_style, or even a fix addforce?