Commit e5b86910 authored by Stan Moore's avatar Stan Moore
Browse files

Merge branch 'master' of https://github.com/lammps/lammps into kk_snap

parents b3747ce9 037cdfe0

doc/src/fix_bond_react.txt

+5 −2
Original line number Diff line number Diff line
@@ -36,10 +36,12 @@ react = mandatory argument indicating new reaction specification :l 
  template-ID(post-reacted) = ID of a molecule template containing post-reaction topology :l

  map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates :l

  zero or more individual keyword/value pairs may be appended to each react argument :l

  individual_keyword = {prob} or {stabilize_steps} or {update_edges} :l

  individual_keyword = {prob} or {max_rxn} or {stabilize_steps} or {update_edges} :l

    {prob} values = fraction seed

      fraction = initiate reaction with this probability if otherwise eligible

      seed = random number seed (positive integer)

    {max_rxn} value = N

      N = maximum number of reactions allowed to occur

    {stabilize_steps} value = timesteps

      timesteps = number of timesteps to apply the internally-created "nve/limit"_fix_nve_limit.html fix to reacting atoms

    {update_edges} value = {none} or {charges} or {custom}
@@ -285,7 +287,8 @@ The {prob} keyword can affect whether an eligible reaction actually 
occurs. The fraction setting must be a value between 0.0 and 1.0. A

uniform random number between 0.0 and 1.0 is generated and the

eligible reaction only occurs if the random number is less than the

fraction.

fraction. Up to N reactions are permitted to occur, as optionally

specified by the {max_rxn} keyword.



The {stabilize_steps} keyword allows for the specification of how many

timesteps a reaction site is stabilized before being returned to the

src/KOKKOS/pair_reaxc_kokkos.cpp

+8 −8
Original line number Diff line number Diff line
@@ -1434,14 +1434,6 @@ void PairReaxCKokkos<DeviceType>::allocate_array() 
  d_CdDelta = typename AT::t_ffloat_1d("reax/c/kk:CdDelta",nmax);

  d_sum_ovun = typename AT::t_ffloat_2d_dl("reax/c/kk:sum_ovun",nmax,3);



  // FixReaxCSpecies

  if (fixspecies_flag) {

    memoryKK->destroy_kokkos(k_tmpid,tmpid);

    memoryKK->destroy_kokkos(k_tmpbo,tmpbo);

    memoryKK->create_kokkos(k_tmpid,tmpid,nmax,MAXSPECBOND,"pair:tmpid");

    memoryKK->create_kokkos(k_tmpbo,tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");

  }



  // FixReaxCBonds

  d_abo = typename AT::t_ffloat_2d("reax/c/kk:abo",nmax,maxbo);

  d_neighid = typename AT::t_tagint_2d("reax/c/kk:neighid",nmax,maxbo);
@@ -4237,6 +4229,14 @@ void PairReaxCKokkos<DeviceType>::pack_bond_buffer_item(int i, int &j, const boo 
template<class DeviceType>

void PairReaxCKokkos<DeviceType>::FindBondSpecies()

{



  if (nmax > k_tmpid.extent(0)) {

    memoryKK->destroy_kokkos(k_tmpid,tmpid);

    memoryKK->destroy_kokkos(k_tmpbo,tmpbo);

    memoryKK->create_kokkos(k_tmpid,tmpid,nmax,MAXSPECBOND,"pair:tmpid");

    memoryKK->create_kokkos(k_tmpbo,tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");

  }



  copymode = 1;

  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondSpeciesZero>(0,nmax),*this);

src/MOLECULE/improper_umbrella.cpp

+22 −22
Original line number Diff line number Diff line
@@ -189,17 +189,17 @@ void ImproperUmbrella::compute(int eflag, int vflag) 
    dahy = ary-c*hry;

    dahz = arz-c*hrz;



    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;

    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;

    f2[2] = (dhax*vb1y - dhay*vb1x)*rar;

    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;

    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;

    f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;



    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;

    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;

    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;

    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;

    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;

    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;



    f4[0] = dahx*rhr;

    f4[1] = dahy*rhr;

    f4[2] = dahz*rhr;

    f4[0] = dahx*rhr*a;

    f4[1] = dahy*rhr*a;

    f4[2] = dahz*rhr*a;



    f1[0] = -(f2[0] + f3[0] + f4[0]);

    f1[1] = -(f2[1] + f3[1] + f4[1]);
@@ -208,27 +208,27 @@ void ImproperUmbrella::compute(int eflag, int vflag) 
    // apply force to each of 4 atoms



    if (newton_bond || i1 < nlocal) {

      f[i1][0] += f1[0]*a;

      f[i1][1] += f1[1]*a;

      f[i1][2] += f1[2]*a;

      f[i1][0] += f1[0];

      f[i1][1] += f1[1];

      f[i1][2] += f1[2];

    }



    if (newton_bond || i2 < nlocal) {

      f[i2][0] += f3[0]*a;

      f[i2][1] += f3[1]*a;

      f[i2][2] += f3[2]*a;

      f[i2][0] += f3[0];

      f[i2][1] += f3[1];

      f[i2][2] += f3[2];

    }



    if (newton_bond || i3 < nlocal) {

      f[i3][0] += f2[0]*a;

      f[i3][1] += f2[1]*a;

      f[i3][2] += f2[2]*a;

      f[i3][0] += f2[0];

      f[i3][1] += f2[1];

      f[i3][2] += f2[2];

    }



    if (newton_bond || i4 < nlocal) {

      f[i4][0] += f4[0]*a;

      f[i4][1] += f4[1]*a;

      f[i4][2] += f4[2]*a;

      f[i4][0] += f4[0];

      f[i4][1] += f4[1];

      f[i4][2] += f4[2];

    }



    if (evflag) {
@@ -247,7 +247,7 @@ void ImproperUmbrella::compute(int eflag, int vflag) 
      vb3y = x[i4][1] - x[i3][1];

      vb3z = x[i4][2] - x[i3][2];



      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,

      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,

               vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);

    }

  }

src/Purge.list

+3 −0
Original line number Diff line number Diff line
@@ -24,6 +24,9 @@ style_nstencil.h 
style_ntopo.h

# other auto-generated files

lmpinstalledpkgs.h

# renamed on 7 January 2019

pair_lebedeva.cpp

pair_lebedeva.h

# removed on 12 December 2018

pair_reax.cpp

pair_reax.h

src/USER-MISC/fix_bond_react.cpp

+81 −30
Original line number Diff line number Diff line
@@ -161,6 +161,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 
  memory->create(unreacted_mol,nreacts,"bond/react:unreacted_mol");

  memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");

  memory->create(fraction,nreacts,"bond/react:fraction");

  memory->create(max_rxn,nreacts,"bond/react:max_rxn");

  memory->create(nlocalskips,nreacts,"bond/react:nlocalskips");

  memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips");

  memory->create(seed,nreacts,"bond/react:seed");

  memory->create(limit_duration,nreacts,"bond/react:limit_duration");

  memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
@@ -179,6 +182,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 
  for (int i = 0; i < nreacts; i++) {

    fraction[i] = 1;

    seed[i] = 12345;

    max_rxn[i] = INT_MAX;

    stabilize_steps_flag[i] = 0;

    update_edges_flag[i] = 0;

    // set default limit duration to 60 timesteps
@@ -244,6 +248,13 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 
        if (seed[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command: "

                                       "probability seed must be positive");

        iarg += 3;

      } else if (strcmp(arg[iarg],"max_rxn") == 0) {

	      if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "

	      			                        "'max_rxn' has too few arguments");

	      max_rxn[rxn] = force->inumeric(FLERR,arg[iarg+1]);

	      if (max_rxn[rxn] < 0) error->all(FLERR,"Illegal fix bond/react command: "

	      				                         "'max_rxn' cannot be negative");

	      iarg += 2;

      } else if (strcmp(arg[iarg],"stabilize_steps") == 0) {

        if (stabilization_flag == 0) error->all(FLERR,"Stabilize_steps keyword "

                                                "used without stabilization keyword");
@@ -379,9 +390,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 


FixBondReact::~FixBondReact()

{

  // unregister callbacks to this fix from Atom class

  atom->delete_callback(id,0);



  for (int i = 0; i < nreacts; i++) {

    delete random[i];

  }
@@ -396,6 +404,7 @@ FixBondReact::~FixBondReact() 
  memory->destroy(edge);

  memory->destroy(equivalences);

  memory->destroy(reverse_equiv);

  memory->destroy(landlocked_atoms);

  memory->destroy(custom_edges);

  memory->destroy(delete_atoms);
@@ -405,6 +414,9 @@ FixBondReact::~FixBondReact() 
  memory->destroy(reacted_mol);

  memory->destroy(fraction);

  memory->destroy(seed);

  memory->destroy(max_rxn);

  memory->destroy(nlocalskips);

  memory->destroy(nghostlyskips);

  memory->destroy(limit_duration);

  memory->destroy(stabilize_steps_flag);

  memory->destroy(update_edges_flag);
@@ -434,7 +446,6 @@ FixBondReact::~FixBondReact() 
    memory->destroy(restore);

    memory->destroy(glove);

    memory->destroy(pioneers);

    memory->destroy(landlocked_atoms);

    memory->destroy(local_mega_glove);

    memory->destroy(ghostly_mega_glove);

  }
@@ -445,14 +456,14 @@ FixBondReact::~FixBondReact() 
    delete [] exclude_group;



    // check nfix in case all fixes have already been deleted

    if (id_fix1 == NULL && modify->nfix) modify->delete_fix(id_fix1);

    if (id_fix1 && modify->nfix) modify->delete_fix(id_fix1);

    delete [] id_fix1;



    if (id_fix3 == NULL && modify->nfix) modify->delete_fix(id_fix3);

    if (id_fix3 && modify->nfix) modify->delete_fix(id_fix3);

    delete [] id_fix3;

  }



  if (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2);

  if (id_fix2 && modify->nfix) modify->delete_fix(id_fix2);

  delete [] id_fix2;



  delete [] statted_id;
@@ -677,6 +688,8 @@ void FixBondReact::post_integrate() 
    reaction_count[i] = 0;

    local_rxn_count[i] = 0;

    ghostly_rxn_count[i] = 0;

    nlocalskips[i] = 0;

    nghostlyskips[i] = 0;

  }



  if (nevery_check) {
@@ -748,6 +761,7 @@ void FixBondReact::post_integrate() 


  int j;

  for (rxnID = 0; rxnID < nreacts; rxnID++) {

    if (max_rxn[rxnID] <= reaction_count_total[rxnID]) continue;

    for (int ii = 0; ii < nall; ii++) {

      partner[ii] = 0;

      finalpartner[ii] = 0;
@@ -1145,14 +1159,43 @@ void FixBondReact::superimpose_algorithm() 


  MPI_Allreduce(&local_rxn_count[0],&reaction_count[0],nreacts,MPI_INT,MPI_SUM,world);



  for (int i = 0; i < nreacts; i++)

    reaction_count_total[i] += reaction_count[i];



  // this assumes compute_vector is called from process 0

  // ...so doesn't bother to bcast ghostly_rxn_count

  if (me == 0)

    for (int i = 0; i < nreacts; i++)

      reaction_count_total[i] += ghostly_rxn_count[i];

      reaction_count_total[i] += reaction_count[i] + ghostly_rxn_count[i];



  MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);



  // check if we overstepped our reaction limit

  for (int i = 0; i < nreacts; i++) {

    if (reaction_count_total[i] > max_rxn[i]) {

      // let's randomly choose rxns to skip, unbiasedly from local and ghostly

      int local_rxncounts[nprocs];

      int all_localskips[nprocs];

      MPI_Gather(&local_rxn_count[i],1,MPI_INT,local_rxncounts,1,MPI_INT,0,world);

      if (me == 0) {

        int overstep = reaction_count_total[i] - max_rxn[i];

        int delta_rxn = reaction_count[i] + ghostly_rxn_count[i];

        int rxn_by_proc[delta_rxn];

        for (int j = 0; j < delta_rxn; j++)

          rxn_by_proc[j] = -1; // corresponds to ghostly

        int itemp = 0;

        for (int j = 0; j < nprocs; j++)

          for (int k = 0; k < local_rxn_count[j]; k++)

            rxn_by_proc[itemp++] = j;

        std::random_shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn]);

        for (int j = 0; j < nprocs; j++)

          all_localskips[j] = 0;

        nghostlyskips[i] = 0;

        for (int j = 0; j < overstep; j++) {

          if (rxn_by_proc[j] == -1) nghostlyskips[i]++;

          else all_localskips[rxn_by_proc[j]]++;

        }

      }

      reaction_count_total[i] = max_rxn[i];

      MPI_Scatter(&all_localskips[0],1,MPI_INT,&nlocalskips[i],1,MPI_INT,0,world);

      MPI_Bcast(&nghostlyskips[i],1,MPI_INT,0,world);

    }

  }



  // this updates topology next step

  next_reneighbor = update->ntimestep;
@@ -1948,6 +1991,7 @@ void FixBondReact::glove_ghostcheck() 
  // 'ghosts of another' indication taken from comm->sendlist



  int ghostly = 0;

  #if !defined(MPI_STUBS)

    if (comm->style == 0) {

      for (int i = 0; i < onemol->natoms; i++) {

        int ilocal = atom->map(glove[i][1]);
@@ -1957,10 +2001,9 @@ void FixBondReact::glove_ghostcheck() 
        }

      }

    } else {

    #if !defined(MPI_STUBS)

      ghostly = 1;

    #endif

    }

  #endif



  if (ghostly == 1) {

    ghostly_mega_glove[0][ghostly_num_mega] = rxnID;
@@ -2095,18 +2138,26 @@ void FixBondReact::update_everything() 
  memory->create(update_mega_glove,max_natoms+1,MAX(local_num_mega,global_megasize),"bond/react:update_mega_glove");



  for (int pass = 0; pass < 2; pass++) {



    update_num_mega = 0;

    int iskip[nreacts];

    for (int i = 0; i < nreacts; i++) iskip[i] = 0;

    if (pass == 0) {

      update_num_mega = local_num_mega;

      for (int i = 0; i < update_num_mega; i++) {

      for (int i = 0; i < local_num_mega; i++) {

        rxnID = local_mega_glove[0][i];

        // reactions already shuffled from dedup procedure, so can skip first N

        if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;

        for (int j = 0; j < max_natoms+1; j++)

          update_mega_glove[j][i] = local_mega_glove[j][i];

          update_mega_glove[j][update_num_mega] = local_mega_glove[j][i];

        update_num_mega++;

      }

    } else if (pass == 1) {

      update_num_mega = global_megasize;

      for (int i = 0; i < global_megasize; i++) {

        rxnID = global_mega_glove[0][i];

        // reactions already shuffled from dedup procedure, so can skip first N

        if (iskip[rxnID]++ < nghostlyskips[rxnID]) continue;

        for (int j = 0; j < max_natoms+1; j++)

          update_mega_glove[j][i] = global_mega_glove[j][i];

          update_mega_glove[j][update_num_mega] = global_mega_glove[j][i];

        update_num_mega++;

      }

    }

CRA Git | Maintained and supported by SUSTech CRA and CCSE