Skip to content

Commit

Permalink
refactor creation of bonded particles (2)
Browse files Browse the repository at this point in the history 
reduce code duplication
  • Loading branch information
@danielque
danielque committed Aug 9, 2018
1 parent cd88c96 commit be85da6
Show file tree
Hide file tree
Showing 4 changed files with 65 additions and 116 deletions.
86 changes: 37 additions & 49 deletions src/fix_template_multiplespheres.cpp
View file
Expand Up @@ -76,8 +76,7 @@ FixTemplateMultiplespheres::FixTemplateMultiplespheres(LAMMPS *lmp, int narg, ch
r_sphere = new double[nspheres];
atom_type_sphere = 0;

bonded_implicit = false;
bonded_explicit = false;
bonded = false;
// number of partners and partner array
np = NULL;
p = NULL;
Expand Down Expand Up @@ -174,24 +173,34 @@ FixTemplateMultiplespheres::FixTemplateMultiplespheres(LAMMPS *lmp, int narg, ch
if (narg < iarg+2)
error->fix_error(FLERR,this,"not enough arguments for 'bonded'");

++iarg;
// Note: Implicitly creates the bonds between the atoms based on the neighbour list and a bonding distance
if(strcmp(arg[iarg+1],"yes/implicit") == 0)
bonded_implicit = true;

// Note: Exlicitly creates bonds between the atoms without using any bonding distnace equation and neighbor lists.
if(strcmp(arg[iarg],"yes") == 0 || strcmp(arg[iarg],"yes/implicit") == 0)
{
++iarg;
bonded = true;
}
// Note: Explicitly creates bonds between the atoms without using any bonding distance equation and neighbor lists.
// This bonding scheme has been implemented for the specific application. User-defined partner list is required
// for creating the bonds.
else if(strcmp(arg[iarg+1],"yes/explicit") == 0)
else if(strcmp(arg[iarg],"yes/explicit") == 0)
{
if(narg<iarg+3) {error->fix_error(FLERR,this,"Partner list is required for explicit bonds!");}
if(narg<iarg+2)
error->fix_error(FLERR,this,"Partner list is required for explicit bonds!");

if(strcmp(arg[iarg+2],"nbond_pairs") == 0){
int nbond_pairs = atoi(arg[iarg+3]);
++iarg;
if(strcmp(arg[iarg],"nbond_pairs") == 0)
{
++iarg;
int nbond_pairs = atoi(arg[iarg]);

if (narg < iarg + 3 + 2*nbond_pairs) error->fix_error(FLERR,this,"not enough arguments: define the partner list here!");
if (narg < iarg + 2*nbond_pairs + 1)
{
error->fix_error(FLERR,this,"not enough arguments: define the partner list here!");
}
else
{
int counter = 0;
++iarg;
// number of partners and partner array
np = new int[nspheres]();
p = new int*[nspheres];
Expand All @@ -200,9 +209,8 @@ FixTemplateMultiplespheres::FixTemplateMultiplespheres(LAMMPS *lmp, int narg, ch

for(int i = 0; i < nbond_pairs; ++i)
{
int ipartner = atoi(arg[iarg+4+counter++])-1;
int jpartner = atoi(arg[iarg+4+counter++])-1;

int ipartner = atoi(arg[iarg++])-1;
int jpartner = atoi(arg[iarg++])-1;
{
p[ipartner][np[ipartner]] = jpartner;
p[jpartner][np[jpartner]] = ipartner;
Expand All @@ -212,18 +220,20 @@ FixTemplateMultiplespheres::FixTemplateMultiplespheres(LAMMPS *lmp, int narg, ch
}
}
}
bonded_explicit = true;
}
bonded = true;
}


}
//Note: This switch can be used to create a multiplesphere particle without any bonds.
else if(strcmp(arg[iarg+1],"no") == 0){
bonded_implicit = false;
bonded_explicit = false;}
else if(strcmp(arg[iarg],"no") == 0)
{
++iarg;
bonded = false;
}
else
{
error->fix_error(FLERR,this,"expecting 'yes/implicit' or 'yes/explicit' or 'no' after 'bonded'");
iarg+=2;
}
hasargs = true;
}
else if(strcmp(style,"particletemplate/multiplespheres") == 0)
error->fix_error(FLERR,this,"unknown keyword");
Expand Down Expand Up @@ -270,7 +280,7 @@ void FixTemplateMultiplespheres::post_create()
calc_bounding_sphere();
calc_center_of_mass();

if((bonded_implicit || bonded_explicit) && !fix_bond_random_id)
if(bonded && !fix_bond_random_id)
{
fix_bond_random_id = static_cast<FixPropertyAtom*>(modify->find_fix_property("bond_random_id","property/atom","scalar",0,0,this->style,false));

Expand Down Expand Up @@ -552,7 +562,7 @@ void FixTemplateMultiplespheres::randomize_ptilist(int n_random,int distribution

pti->groupbit = groupbit | distribution_groupbit; //NP also contains insert_groupbit

if(bonded_implicit || bonded_explicit)
if(bonded)
{
if (!pti->fix_property)
{
Expand All @@ -579,29 +589,7 @@ void FixTemplateMultiplespheres::randomize_ptilist(int n_random,int distribution

void FixTemplateMultiplespheres::finalize_insertion()
{
if(bonded_implicit)
{
// check each pti since we don't know which of them have actually been inserted
for(int i = 0; i < n_pti_max; ++i)
{
ParticleToInsert *pti = pti_list[i];

// only need to create bonds if fix_property is fix_bond_random_id
if(pti->fix_property && pti->fix_property[0] == fix_bond_random_id)
{
// only need to create bonds for pti created this time step
// timestep is integer part of fix_property_value
// Note: it's possble that not each pti gets inserted, the pti itself
// does the final check if bond creation is necessary
bigint insertionstep = static_cast<bigint>(pti->fix_property_value[0][0]);
if(insertionstep == update->ntimestep)
pti->create_bonds_implicit();
}
}
}


if(bonded_explicit)
if(bonded)
{
// check each pti since we don't know which of them have actually been inserted
for(int i = 0; i < n_pti_max; ++i)
Expand All @@ -616,7 +604,7 @@ void FixTemplateMultiplespheres::finalize_insertion()
// does the final check if bond creation is necessary
bigint insertionstep = static_cast<bigint>(pti->fix_property_value[0][0]);
if(insertionstep == update->ntimestep)
pti->create_bonds_explicit(np, p);
pti->create_bonds(np, p);
}
}
}
Expand Down
5 changes: 2 additions & 3 deletions src/fix_template_multiplespheres.h
View file
Expand Up @@ -48,7 +48,7 @@ class FixTemplateMultiplespheres : public FixTemplateSphere {
int maxtype();
int mintype();
int number_spheres();
bool is_bonded() const { return bonded_implicit || bonded_explicit; }
bool is_bonded() const { return bonded; }

// multi insertion
virtual void init_ptilist(int);
Expand Down Expand Up @@ -103,8 +103,7 @@ class FixTemplateMultiplespheres : public FixTemplateSphere {
// number of tries for mc
int ntry;

bool bonded_implicit;
bool bonded_explicit;
bool bonded;

class FixPropertyAtom *fix_bond_random_id;
};
Expand Down
85 changes: 23 additions & 62 deletions src/particleToInsert.cpp
View file
Expand Up @@ -160,16 +160,10 @@ int ParticleToInsert::insert()

/* ---------------------------------------------------------------------- */

int ParticleToInsert::create_bonds_implicit()
int ParticleToInsert::create_bond_partners(int*& npartner, int**&partner)
{
if(nspheres == 1 || !needs_bonding || local_start < 0)
return 0;

needs_bonding = false; // reset in case pti gets reused

// find bond partners
int *npartner = new int[nspheres](); // convert to member variable to be set from fix template
int **partner = new int*[nspheres]; // convert to member variable to be set from fix template
npartner = new int[nspheres]();
partner = new int*[nspheres];

for(int i = 0; i < nspheres; ++i)
partner[i] = new int[nspheres-1]();
Expand Down Expand Up @@ -201,70 +195,27 @@ int ParticleToInsert::create_bonds_implicit()
partner[j][npartner[j]] = i;
npartner[i]++;
npartner[j]++;
create_bonds = 1;
create_bonds = 2;
}
}
}

int ncreate = 0;

if(create_bonds)
{
// create bonds
int **_bond_type = atom->bond_type;
int **_bond_atom = atom->bond_atom;
int *num_bond = atom->num_bond;
int newton_bond = force->newton_bond;
int n_bondhist = atom->n_bondhist;
double ***bond_hist = atom->bond_hist;

// Note: atoms are created with dummy tag = 0, but
// actual tags must be available at this point,
// i.e. atom->tag_extend() must have been called
for(int i = 0; i < nspheres; ++i)
{
if (npartner[i] == 0) continue;

for(int k = 0; k < npartner[i]; ++k)
{
const int j = partner[i][k];
if (!newton_bond || i < j)
{
const int ilocal = local_start + i;

if (num_bond[ilocal] == atom->bond_per_atom)
{
error->one(FLERR,"New bond exceeded bonds per atom in fix bond/create");
}

_bond_type[ilocal][num_bond[ilocal]] = bond_type;
_bond_atom[ilocal][num_bond[ilocal]] = atom->tag[local_start+j];

// reset history
for (int ih = 0; ih < n_bondhist; ++ih)
{
bond_hist[ilocal][num_bond[ilocal]][ih] = 0.;
}
num_bond[ilocal]++;
}
return create_bonds;
}

if(i < j)
++ncreate;
}
}
}
/* ---------------------------------------------------------------------- */

void ParticleToInsert::destroy_bond_partners(int* npartner, int** partner)
{
for(int i = 0; i < nspheres; ++i)
delete [] partner[i];
delete [] partner;
delete [] npartner;

return ncreate;
}

/* ---------------------------------------------------------------------- */

int ParticleToInsert::create_bonds_explicit(int *npartner, int **partner)
int ParticleToInsert::create_bonds(int *npartner, int **partner)
{
if(nspheres == 1 || !needs_bonding || local_start < 0)
return 0;
Expand All @@ -273,7 +224,12 @@ int ParticleToInsert::create_bonds_explicit(int *npartner, int **partner)

int create_bonds = 1;

int ncreate = 0;
if(!npartner && !partner)
{
create_bonds = create_bond_partners(npartner, partner);
}

int ncreated = 0;

if(create_bonds)
{
Expand Down Expand Up @@ -316,12 +272,17 @@ int ParticleToInsert::create_bonds_explicit(int *npartner, int **partner)
}

if(i < j)
++ncreate;
++ncreated;
}
}
}

return ncreate;
if(create_bonds != 1) // create_bond_partners allocated memory
{
destroy_bond_partners(npartner, partner);
}

return ncreated;
}

/* ---------------------------------------------------------------------- */
Expand Down
5 changes: 3 additions & 2 deletions src/particleToInsert.h
View file
Expand Up @@ -89,11 +89,12 @@ namespace LAMMPS_NS {
virtual int set_x_v_omega(double *,double *,double *,double *);

virtual void scale_pti(double r_scale);
int create_bonds_implicit();
int create_bonds_explicit(int *npartner, int **partner);
int create_bonds(int *npartner=NULL, int **partner=NULL);
private:
int local_start;
bool needs_bonding;
int create_bond_partners(int *&npartner, int **&partner);
void destroy_bond_partners(int *npartner, int **partner);
};

}
Expand Down

0 comments on commit be85da6

Please sign in to comment.