diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index 4a328a7b18..f30a82873f 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -223,15 +223,12 @@ void PairReaxC::coeff( int nargs, char **args ) // read args that map atom types to elements in potential file // map[i] = which element the Ith atom type is, -1 if NULL - // NOTE: for now throw an error if NULL is used to disallow use with hybrid - // REAX/C lib will have to be modified to allow this int itmp; int nreax_types = system->reax_param.num_atom_types; for (int i = 3; i < nargs; i++) { if (strcmp(args[i],"NULL") == 0) { map[i-2] = -1; - error->all("Cannot currently use pair reax/c with pair hybrid"); continue; } @@ -242,6 +239,7 @@ void PairReaxC::coeff( int nargs, char **args ) if (itmp < 0 || itmp >= nreax_types) error->all("Non-existent ReaxFF type"); + } int n = atom->ntypes; @@ -292,7 +290,7 @@ void PairReaxC::init_style( ) cutmax = MAX3(control->nonb_cut, control->hbond_cut, 2*control->bond_cut); - for(int i = 0; i < LIST_N; ++i ) + for( int i = 0; i < LIST_N; ++i ) lists[i].allocated = 0; if (fix_reax == NULL) { diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp index de2f117bbc..e31c86d491 100644 --- a/src/USER-REAXC/reaxc_allocate.cpp +++ b/src/USER-REAXC/reaxc_allocate.cpp @@ -512,7 +512,7 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds, *total_bonds = 0; *est_3body = 0; for( i = 0; i < system->N; ++i ){ - *est_3body += SQR( Num_Entries( i, bonds ) ); + *est_3body += SQR(system->my_atoms[i].num_bonds); // commented out - already updated in validate_lists in forces.c // system->my_atoms[i].num_bonds = MAX( Num_Entries(i,bonds)*2, MIN_BONDS ); *total_bonds += system->my_atoms[i].num_bonds; diff --git a/src/USER-REAXC/reaxc_nonbonded.cpp b/src/USER-REAXC/reaxc_nonbonded.cpp index e03c9ef102..ed475606fa 100644 --- a/src/USER-REAXC/reaxc_nonbonded.cpp +++ b/src/USER-REAXC/reaxc_nonbonded.cpp @@ -388,7 +388,7 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace, lr->e_vdW = Tap * twbp->D * (exp1 - 2.0 * exp2); lr->CEvd = dTap * twbp->D * (exp1 - 2.0 * exp2) - - Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2); + Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) / r_ij; } if(system->reax_param.gp.vdw_type==2 || system->reax_param.gp.vdw_type==3) @@ -397,7 +397,7 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace, lr->e_vdW += Tap * e_core; de_core = -(twbp->acore/twbp->rcore) * e_core; - lr->CEvd += dTap * e_core + Tap * de_core; + lr->CEvd += dTap * e_core + Tap * de_core / r_ij; }