diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 0a6ed3e34a..b70092296f 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -1475,15 +1475,14 @@ void FixBondReact::superimpose_algorithm() int rxnflag = 0; int *delta_rxn; - memory->create(delta_rxn,nreacts,"bond/react:delta_rxn"); + memory->create(delta_rxn, nreacts, "bond/react:delta_rxn"); if (comm->me == 0) for (int i = 0; i < nreacts; i++) { delta_rxn[i] = reaction_count[i] + ghostly_rxn_count[i]; - reaction_count_total[i] += delta_rxn[i]; rxnflag += delta_rxn[i]; } - MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world); + MPI_Bcast(&delta_rxn[0], nreacts, MPI_INT, 0, world); MPI_Bcast(&rxnflag, 1, MPI_INT, 0, world); if (!rxnflag) return; @@ -1495,12 +1494,12 @@ void FixBondReact::superimpose_algorithm() // check if we overstepped our reaction limit, via either max_rxn or rate_limit for (int i = 0; i < nreacts; i++) { int overstep = 0; - int max_rxn_overstep = reaction_count_total[i] - max_rxn[i]; + int max_rxn_overstep = reaction_count_total[i] + delta_rxn[i] - max_rxn[i]; overstep = MAX(overstep,max_rxn_overstep); if (rate_limit[0][i] == 1) { int myrxn_count = store_rxn_count[rate_limit[2][i]-1][i]; if (myrxn_count != -1) { - int nrxn_delta = reaction_count_total[i] - myrxn_count; + int nrxn_delta = reaction_count_total[i] + delta_rxn[i] - myrxn_count; int my_nrate; if (var_flag[NRATE][i] == 1) { my_nrate = input->variable->compute_equal(var_id[NRATE][i]); @@ -1540,7 +1539,6 @@ void FixBondReact::superimpose_algorithm() else all_localkeep[rxn_by_proc[j]]++; } memory->destroy(rxn_by_proc); - reaction_count_total[i] -= overstep; } MPI_Scatter(&all_localkeep[0],1,MPI_INT,&nlocalkeep[i],1,MPI_INT,0,world); MPI_Bcast(&nghostlykeep[i],1,MPI_INT,0,world); @@ -1548,7 +1546,6 @@ void FixBondReact::superimpose_algorithm() memory->destroy(all_localkeep); } } - MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world); memory->destroy(delta_rxn); // this updates topology next step @@ -3103,18 +3100,16 @@ void FixBondReact::update_everything() if (create_atoms_flag[rxnID] == 1) { onemol = atom->molecules[unreacted_mol[rxnID]]; twomol = atom->molecules[reacted_mol[rxnID]]; - if (insert_atoms_setup(update_mega_glove,update_num_mega)) { - inserted_atoms_flag = 1; - } else { // create aborted - reaction_count_total[rxnID]--; - continue; - } + if (insert_atoms_setup(update_mega_glove,update_num_mega)) inserted_atoms_flag = 1; + else continue; } noccur[rxnID]++; if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = local_mega_glove[1][i]; update_num_mega++; } + MPI_Allreduce(MPI_IN_PLACE, &noccur[0], nreacts, MPI_INT, MPI_SUM, world); + reaction_count_total[rxnID] += noccur[rxnID]; } else if (pass == 1) { for (int i = 0; i < global_megasize; i++) { rxnID = (int) global_mega_glove[0][i]; @@ -3129,17 +3124,15 @@ void FixBondReact::update_everything() // we can insert atoms here, now that reactions are finalized // can't do it any earlier, due to skipped reactions (max_rxn) // for MPI build, reactions that create atoms are always treated as 'global' + if (create_atoms_flag[rxnID] == 1) { onemol = atom->molecules[unreacted_mol[rxnID]]; twomol = atom->molecules[reacted_mol[rxnID]]; - if (insert_atoms_setup(update_mega_glove,update_num_mega)) { - inserted_atoms_flag = 1; - } else { // create aborted - reaction_count_total[rxnID]--; - continue; - } + if (insert_atoms_setup(update_mega_glove,update_num_mega)) inserted_atoms_flag = 1; + else continue; } noccur[rxnID]++; + reaction_count_total[rxnID]++; if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = global_mega_glove[1][i]; update_num_mega++;