From 0288bb4b6b2d0eec681ee428d57cd81b18e1d081 Mon Sep 17 00:00:00 2001 From: wverestek Date: Wed, 22 Apr 2020 13:27:10 +0200 Subject: [PATCH] small modification to fix bond/react to allow equal style variables as probability o On branch master Your branch is up to date with 'origin/master'. Changes to be committed: modified: doc/src/fix_bond_react.rst new file: examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability modified: src/USER-REACTION/fix_bond_react.cpp modified: src/USER-REACTION/fix_bond_react.h --- doc/src/fix_bond_react.rst | 6 +- .../in.large_nylon_melt_variable_probability | 55 +++++++++++++++++++ src/USER-REACTION/fix_bond_react.cpp | 35 +++++++++++- src/USER-REACTION/fix_bond_react.h | 1 + 4 files changed, 93 insertions(+), 4 deletions(-) create mode 100644 examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 8daa52a73d..3467932500 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -419,8 +419,10 @@ it occurs: The *prob* keyword can affect whether or not 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 +1.0 or can be an equal style variable. In the later case the variable +is evaluated during runtime and adjusted to be between 0.0 and 1.0 if +necessary. 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. Up to N reactions are permitted to occur, as optionally specified by the *max_rxn* keyword. diff --git a/examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability b/examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability new file mode 100644 index 0000000000..f169f19504 --- /dev/null +++ b/examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability @@ -0,0 +1,55 @@ +# 35,000 atom nylon melt example + +units real + +boundary p p p + +atom_style full + +kspace_style pppm 1.0e-4 + +pair_style lj/class2/coul/long 8.5 + +angle_style class2 + +bond_style class2 + +dihedral_style class2 + +improper_style class2 + +read_data large_nylon_melt.data.gz + +variable runsteps equal 200 +varaible prob equal step/v_runsteps + +velocity all create 800.0 4928459 dist gaussian + +molecule mol1 rxn1_stp1_unreacted.data_template +molecule mol2 rxn1_stp1_reacted.data_template +molecule mol3 rxn1_stp2_unreacted.data_template +molecule mol4 rxn1_stp2_reacted.data_template + +thermo 50 + +# dump 1 all xyz 100 test_vis.xyz + +fix myrxns all bond/react stabilization yes statted_grp .03 & + react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map & + react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map + +# stable at 800K +fix 1 statted_grp_REACT nvt temp 800 800 100 + +# in order to customize behavior of reacting atoms, +# you can use the internally created 'bond_react_MASTER_group', like so: +# fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1 + +thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts + +# restart 100 restart1 restart2 + +run 200 + +# write_restart restart_longrun +# write_data restart_longrun.data diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 44e2bd172e..ae4c977f4d 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -39,6 +39,8 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) #include "math_extra.h" #include "memory.h" #include "error.h" +#include "input.h" +#include "variable.h" #include @@ -175,6 +177,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 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(var_fraction_flag,nreacts,"bond/react:var_fraction_flag"); + memory->create(var_fraction_id,nreacts,"bond/react:var_fraction_id"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag"); memory->create(constraints,1,MAXCONARGS,"bond/react:constraints"); @@ -193,6 +197,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fraction[i] = 1; seed[i] = 12345; max_rxn[i] = INT_MAX; + var_fraction_flag[i] = 0; + var_fraction_id[i] = 0; stabilize_steps_flag[i] = 0; update_edges_flag[i] = 0; // set default limit duration to 60 timesteps @@ -251,7 +257,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg],"prob") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/react command: " "'prob' keyword has too few arguments"); - fraction[rxn] = force->numeric(FLERR,arg[iarg+1]); + // check if probability is a variable + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + int n = strlen(&arg[iarg+1][2]) + 1; + char *fracstr = new char[n]; + strcpy(fracstr,&arg[iarg+1][2]); + var_fraction_id[rxn] = input->variable->find(fracstr); + if (var_fraction_id[rxn] < 0) + error->all(FLERR,"variable name for fix bond/react does not exist"); + if (! input->variable->equalstyle(var_fraction_id[rxn])) + error->all(FLERR,"variable in bond/react is not equal style"); + fraction[rxn] = input->variable->compute_equal(var_fraction_id[rxnID]); + var_fraction_flag[rxn] = 1.0; + delete [] fracstr; + } else { + // otherwise probability should be a number + fraction[rxn] = force->numeric(FLERR,arg[iarg+1]); + } seed[rxn] = force->inumeric(FLERR,arg[iarg+2]); if (fraction[rxn] < 0.0 || fraction[rxn] > 1.0) error->all(FLERR,"Illegal fix bond/react command: " @@ -447,6 +469,8 @@ FixBondReact::~FixBondReact() memory->destroy(nlocalskips); memory->destroy(nghostlyskips); memory->destroy(limit_duration); + memory->destroy(var_fraction_flag); + memory->destroy(var_fraction_id); memory->destroy(stabilize_steps_flag); memory->destroy(update_edges_flag); @@ -824,10 +848,17 @@ void FixBondReact::post_integrate() comm->reverse_comm_fix(this); } + // update reaction probability + if (var_fraction_flag[rxnID]) { + fraction[rxnID] = input->variable->compute_equal(var_fraction_id[rxnID]); + if (fraction[rxnID] < 0.0) fraction[rxnID] = 0.0; + if (fraction[rxnID] > 1.0) fraction[rxnID] = 1.0; + } + // each atom now knows its winning partner // for prob check, generate random value for each atom with a bond partner // forward comm of partner and random value, so ghosts have it - + if (fraction[rxnID] < 1.0) { for (int i = 0; i < nlocal; i++) if (partner[i]) probability[i] = random[rxnID]->uniform(); diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index 84cb9401f8..6f5ec9c7d7 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -64,6 +64,7 @@ class FixBondReact : public Fix { int custom_exclude_flag; int *stabilize_steps_flag; int *update_edges_flag; + int *var_fraction_flag, *var_fraction_id; int nconstraints; int narrhenius; double **constraints;