diff --git a/src/ML-IAP/pair_mliap.cpp b/src/ML-IAP/pair_mliap.cpp index 0d2a80635a..64c5a547be 100644 --- a/src/ML-IAP/pair_mliap.cpp +++ b/src/ML-IAP/pair_mliap.cpp @@ -264,8 +264,7 @@ void PairMLIAP::e_tally(MLIAPData* data) add virial contribution into global and per-atom accumulators ------------------------------------------------------------------------- */ -void PairMLIAP::v_tally(int i, int j, - double *fij, double *rij) +void PairMLIAP::v_tally(int i, int j, double *fij, double *rij) { double v[6]; diff --git a/src/OPENMP/pair_reaxc_omp.cpp b/src/OPENMP/pair_reaxc_omp.cpp index 87190f5285..2596f8390b 100644 --- a/src/OPENMP/pair_reaxc_omp.cpp +++ b/src/OPENMP/pair_reaxc_omp.cpp @@ -233,9 +233,6 @@ void PairReaxCOMP::compute(int eflag, int vflag) evdwl = ecoul = 0.0; ev_init(eflag,vflag); - if (vflag_global) api->control->virial = 1; - else api->control->virial = 0; - if (vflag_atom) error->all(FLERR,"Pair style reax/c/omp does not support " "computing per-atom stress"); diff --git a/src/OPENMP/pair_reaxc_omp.h b/src/OPENMP/pair_reaxc_omp.h index 895dc646af..d8152656ff 100644 --- a/src/OPENMP/pair_reaxc_omp.h +++ b/src/OPENMP/pair_reaxc_omp.h @@ -47,33 +47,49 @@ class PairReaxCOMP : public PairReaxC, public ThrOMP { reduce_thr(styleparm, eflagparm, vflagparm, thrparm); } - inline void ev_tally_thr_proxy(Pair *const pairparm, const int iparm, const int jparm, + inline void ev_tally_thr_proxy(const int iparm, const int jparm, const int nlocalparm, const int newton_pairparm, const double evdwlparm, const double ecoulparm, const double fpairparm, const double delxparm, const double delyparm, const double delzparm, ThrData *const thrparm) { - ev_tally_thr(pairparm, iparm, jparm, nlocalparm, newton_pairparm, evdwlparm, ecoulparm, + ev_tally_thr(this, iparm, jparm, nlocalparm, newton_pairparm, evdwlparm, ecoulparm, fpairparm, delxparm, delyparm, delzparm, thrparm); } - inline void ev_tally_xyz_thr_proxy(Pair *const pairparm, const int iparm, const int jparm, + inline void ev_tally_xyz_thr_proxy(const int iparm, const int jparm, const int nlocalparm, const int newton_pairparm, const double evdwlparm, const double ecoulparm, const double fxparm, const double fyparm, const double fzparm, const double delxparm, const double delyparm, const double delzparm, ThrData *const thrparm) { - ev_tally_xyz_thr(pairparm, iparm, jparm, nlocalparm, newton_pairparm, evdwlparm, ecoulparm, + ev_tally_xyz_thr(this, iparm, jparm, nlocalparm, newton_pairparm, evdwlparm, ecoulparm, fxparm, fyparm, fzparm, delxparm, delyparm, delzparm, thrparm); } - inline void ev_tally3_thr_proxy(Pair *const pairparm, int i, int j, int k, double evdwl, + inline void ev_tally3_thr_proxy(int i, int j, int k, double evdwl, double ecoul, double *fj, double *fk, double *drji, double *drki, ThrData *const thrparm) { - ev_tally3_thr(pairparm, i, j, k, evdwl, ecoul, fj, fk, drji, drki, thrparm); + ev_tally3_thr(this, i, j, k, evdwl, ecoul, fj, fk, drji, drki, thrparm); + } + + inline void v_tally3_thr_proxy(const int i, const int j, const int k, const double *const fi, + const double *const fk, const double *const drij, + const double *const drkj, ThrData *const thrparm) + { + v_tally3_thr(this, i, j, k, fi, fk, drij, drkj, thrparm); + } + + inline void v_tally4_thr_proxy(const int i, const int j, const int k, const int l, + const double *const fi, const double *const fj, + const double *const fk, const double *const dril, + const double *const drjl, const double *const drkl, + ThrData *const thrparm) + { + v_tally4_thr(this, i, j, k, l, fi, fj, fk, dril, drjl, drkl, thrparm); } protected: diff --git a/src/OPENMP/reaxc_bond_orders_omp.cpp b/src/OPENMP/reaxc_bond_orders_omp.cpp index 22d7fe1e6e..9c026ee95c 100644 --- a/src/OPENMP/reaxc_bond_orders_omp.cpp +++ b/src/OPENMP/reaxc_bond_orders_omp.cpp @@ -37,109 +37,6 @@ using namespace LAMMPS_NS; namespace ReaxFF { - void Add_dBond_to_Forces_NPTOMP(reax_system *system, int i, int pj, - storage *workspace, reax_list **lists) { - reax_list *bonds = (*lists) + BONDS; - bond_data *nbr_j, *nbr_k; - bond_order_data *bo_ij, *bo_ji; - dbond_coefficients coef; - rvec temp; - int pk, k, j; - - int tid = get_tid(); - long reductionOffset = (system->N * tid); - - /* Initializations */ - nbr_j = &(bonds->select.bond_list[pj]); - j = nbr_j->nbr; - bo_ij = &(nbr_j->bo_data); - bo_ji = &(bonds->select.bond_list[nbr_j->sym_index].bo_data); - - coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - coef.C2dbo = bo_ij->C2dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - coef.C3dbo = bo_ij->C3dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - - coef.C1dbopi = bo_ij->C1dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C2dbopi = bo_ij->C2dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C3dbopi = bo_ij->C3dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C4dbopi = bo_ij->C4dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - - coef.C1dbopi2 = bo_ij->C1dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C2dbopi2 = bo_ij->C2dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C3dbopi2 = bo_ij->C3dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C4dbopi2 = bo_ij->C4dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - - coef.C1dDelta = bo_ij->C1dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - coef.C2dDelta = bo_ij->C2dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - - - /************************************ - * forces related to atom i * - * first neighbors of atom i * - ************************************/ - for (pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk) { - nbr_k = &(bonds->select.bond_list[pk]); - k = nbr_k->nbr; - - rvec_Scale(temp, -coef.C2dbo, nbr_k->bo_data.dBOp); /*2nd, dBO*/ - rvec_ScaledAdd(temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);/*dDelta*/ - rvec_ScaledAdd(temp, -coef.C3dbopi, nbr_k->bo_data.dBOp); /*3rd, dBOpi*/ - rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);/*3rd, dBOpi2*/ - - /* force */ - rvec_Add(workspace->forceReduction[reductionOffset+k],temp); - } - - /* then atom i itself */ - rvec_Scale(temp, coef.C1dbo, bo_ij->dBOp); /*1st,dBO*/ - rvec_ScaledAdd(temp, coef.C2dbo, workspace->dDeltap_self[i]); /*2nd,dBO*/ - rvec_ScaledAdd(temp, coef.C1dDelta, bo_ij->dBOp); /*1st,dBO*/ - rvec_ScaledAdd(temp, coef.C2dDelta, workspace->dDeltap_self[i]);/*2nd,dBO*/ - rvec_ScaledAdd(temp, coef.C1dbopi, bo_ij->dln_BOp_pi); /*1st,dBOpi*/ - rvec_ScaledAdd(temp, coef.C2dbopi, bo_ij->dBOp); /*2nd,dBOpi*/ - rvec_ScaledAdd(temp, coef.C3dbopi, workspace->dDeltap_self[i]);/*3rd,dBOpi*/ - - rvec_ScaledAdd(temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2); /*1st,dBO_pi2*/ - rvec_ScaledAdd(temp, coef.C2dbopi2, bo_ij->dBOp); /*2nd,dBO_pi2*/ - rvec_ScaledAdd(temp, coef.C3dbopi2, workspace->dDeltap_self[i]);/*3rd*/ - - /* force */ - rvec_Add(workspace->forceReduction[reductionOffset+i],temp); - - for (pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk) { - nbr_k = &(bonds->select.bond_list[pk]); - k = nbr_k->nbr; - - rvec_Scale(temp, -coef.C3dbo, nbr_k->bo_data.dBOp); /*3rd,dBO*/ - rvec_ScaledAdd(temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);/*dDelta*/ - rvec_ScaledAdd(temp, -coef.C4dbopi, nbr_k->bo_data.dBOp); /*4th,dBOpi*/ - rvec_ScaledAdd(temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);/*4th,dBOpi2*/ - - /* force */ - rvec_Add(workspace->forceReduction[reductionOffset+k],temp); - } - - /* then atom j itself */ - rvec_Scale(temp, -coef.C1dbo, bo_ij->dBOp); /*1st, dBO*/ - rvec_ScaledAdd(temp, coef.C3dbo, workspace->dDeltap_self[j]); /*2nd, dBO*/ - rvec_ScaledAdd(temp, -coef.C1dDelta, bo_ij->dBOp); /*1st, dBO*/ - rvec_ScaledAdd(temp, coef.C3dDelta, workspace->dDeltap_self[j]);/*2nd, dBO*/ - - rvec_ScaledAdd(temp, -coef.C1dbopi, bo_ij->dln_BOp_pi); /*1st,dBOpi*/ - rvec_ScaledAdd(temp, -coef.C2dbopi, bo_ij->dBOp); /*2nd,dBOpi*/ - rvec_ScaledAdd(temp, coef.C4dbopi, workspace->dDeltap_self[j]);/*3rd,dBOpi*/ - - rvec_ScaledAdd(temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2); /*1st,dBOpi2*/ - rvec_ScaledAdd(temp, -coef.C2dbopi2, bo_ij->dBOp); /*2nd,dBOpi2*/ - rvec_ScaledAdd(temp,coef.C4dbopi2,workspace->dDeltap_self[j]);/*3rd,dBOpi2*/ - - /* force */ - rvec_Add(workspace->forceReduction[reductionOffset+j],temp); - } - -/* ---------------------------------------------------------------------- */ - void Add_dBond_to_ForcesOMP(reax_system *system, int i, int pj, storage *workspace, reax_list **lists) { reax_list *bonds = (*lists) + BONDS; @@ -196,11 +93,11 @@ namespace ReaxFF { rvec_Add(workspace->forceReduction[reductionOffset+i],temp); - if (system->pair_ptr->vflag_atom) { + if (system->pair_ptr->vflag_either) { rvec_Scale(fi_tmp, -1.0, temp); rvec_ScaledSum(delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,i,j,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(i,j,system->N,0,0,0, fi_tmp[0],fi_tmp[1],fi_tmp[2], delij[0],delij[1],delij[2],thr); } @@ -217,11 +114,11 @@ namespace ReaxFF { rvec_Add(workspace->forceReduction[reductionOffset+j],temp); - if (system->pair_ptr->vflag_atom) { + if (system->pair_ptr->vflag_either) { rvec_Scale(fj_tmp, -1.0, temp); rvec_ScaledSum(delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,j,i,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(j,i,system->N,0,0,0, fj_tmp[0],fj_tmp[1],fj_tmp[2], delji[0],delji[1],delji[2],thr); } @@ -236,16 +133,16 @@ namespace ReaxFF { rvec_Add(workspace->forceReduction[reductionOffset+k],temp); - if (system->pair_ptr->vflag_atom) { + if (system->pair_ptr->vflag_either) { rvec_Scale(fk_tmp, -1.0, temp); rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,k,i,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(k,i,system->N,0,0,0, fk_tmp[0],fk_tmp[1],fk_tmp[2], delki[0],delki[1],delki[2],thr); rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,k,j,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(k,j,system->N,0,0,0, fk_tmp[0],fk_tmp[1],fk_tmp[2], delkj[0],delkj[1],delkj[2],thr); } @@ -261,17 +158,17 @@ namespace ReaxFF { rvec_Add(workspace->forceReduction[reductionOffset+k],temp); - if (system->pair_ptr->vflag_atom) { + if (system->pair_ptr->vflag_either) { rvec_Scale(fk_tmp, -1.0, temp); rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,k,i,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(k,i,system->N,0,0,0, fk_tmp[0],fk_tmp[1],fk_tmp[2], delki[0],delki[1],delki[2],thr); rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); - pair_reax_ptr->ev_tally_xyz_thr_proxy(system->pair_ptr,k,j,system->N,0,0,0, + pair_reax_ptr->ev_tally_xyz_thr_proxy(k,j,system->N,0,0,0, fk_tmp[0],fk_tmp[1],fk_tmp[2], delkj[0],delkj[1],delkj[2],thr); } diff --git a/src/OPENMP/reaxc_bonds_omp.cpp b/src/OPENMP/reaxc_bonds_omp.cpp index 9380aad5af..2357d20690 100644 --- a/src/OPENMP/reaxc_bonds_omp.cpp +++ b/src/OPENMP/reaxc_bonds_omp.cpp @@ -121,8 +121,8 @@ namespace ReaxFF { -twbp->De_pp * bo_ij->BO_pi2; /* tally into per-atom energy */ - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1, ebond, 0.0, 0.0, 0.0, 0.0, 0.0, thr); /* calculate derivatives of Bond Orders */ @@ -152,8 +152,8 @@ namespace ReaxFF { (gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1)); /* tally into per-atom energy */ - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1, estriph, 0.0, 0.0, 0.0, 0.0, 0.0, thr); bo_ij->Cdbo += decobdbo; diff --git a/src/OPENMP/reaxc_forces_omp.cpp b/src/OPENMP/reaxc_forces_omp.cpp index 18c9419723..26f36f1357 100644 --- a/src/OPENMP/reaxc_forces_omp.cpp +++ b/src/OPENMP/reaxc_forces_omp.cpp @@ -119,33 +119,16 @@ namespace ReaxFF { } } - if (control->virial == 0) { - #if defined(_OPENMP) #pragma omp for schedule(dynamic,50) #endif - for (i = 0; i < system->N; ++i) { - const int startj = Start_Index(i, bonds); - const int endj = End_Index(i, bonds); - for (pj = startj; pj < endj; ++pj) - if (i < bonds->select.bond_list[pj].nbr) - Add_dBond_to_ForcesOMP(system, i, pj, workspace, lists); - } - - } else { - -#if defined(_OPENMP) -#pragma omp for schedule(dynamic,50) -#endif - for (i = 0; i < system->N; ++i) { - const int startj = Start_Index(i, bonds); - const int endj = End_Index(i, bonds); - for (pj = startj; pj < endj; ++pj) - if (i < bonds->select.bond_list[pj].nbr) - Add_dBond_to_Forces_NPTOMP(system, i, pj, workspace, lists); - } - - } // if (virial == 0) + for (i = 0; i < system->N; ++i) { + const int startj = Start_Index(i, bonds); + const int endj = End_Index(i, bonds); + for (pj = startj; pj < endj; ++pj) + if (i < bonds->select.bond_list[pj].nbr) + Add_dBond_to_ForcesOMP(system, i, pj, workspace, lists); + } pair_reax_ptr->reduce_thr_proxy(system->pair_ptr, 0, 1, thr); diff --git a/src/OPENMP/reaxc_hydrogen_bonds_omp.cpp b/src/OPENMP/reaxc_hydrogen_bonds_omp.cpp index d3adedd5b1..1b174aef81 100644 --- a/src/OPENMP/reaxc_hydrogen_bonds_omp.cpp +++ b/src/OPENMP/reaxc_hydrogen_bonds_omp.cpp @@ -43,7 +43,7 @@ namespace ReaxFF { /* ---------------------------------------------------------------------- */ void Hydrogen_BondsOMP(reax_system *system, control_params *control, - simulation_data *data, storage *workspace, reax_list **lists) + simulation_data *data, storage *workspace, reax_list **lists) { const int nthreads = control->nthreads; @@ -58,11 +58,10 @@ namespace ReaxFF { int hblist[MAX_BONDS]; int itr, top; int num_hb_intrs = 0; - ivec rel_jk; double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; double e_hb, e_hb_thr = 0.0, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk; - rvec dvec_jk, force; + rvec dvec_jk; hbond_parameters *hbp; bond_order_data *bo_ij; bond_data *pbond_ij; @@ -119,7 +118,7 @@ namespace ReaxFF { bo_ij = &(pbond_ij->bo_data); if (system->reax_param.sbp[type_i].p_hbond == 2 && - bo_ij->BO >= HB_THRESHOLD) + bo_ij->BO >= HB_THRESHOLD) hblist[top++] = pi; } @@ -145,11 +144,11 @@ namespace ReaxFF { ++num_hb_intrs; Calculate_Theta(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk, - &theta, &cos_theta); + &theta, &cos_theta); /* the derivative of cos(theta) */ Calculate_dCos_ThetaOMP(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk, - &dcos_theta_di, &dcos_theta_dj, - &dcos_theta_dk); + &dcos_theta_di, &dcos_theta_dj, + &dcos_theta_dk); /* hydrogen bond energy*/ sin_theta2 = sin(theta/2.0); @@ -158,7 +157,7 @@ namespace ReaxFF { cos_xhz1 = (1.0 - cos_theta); exp_hb2 = exp(-hbp->p_hb2 * bo_ij->BO); exp_hb3 = exp(-hbp->p_hb3 * (hbp->r0_hb / r_jk + - r_jk / hbp->r0_hb - 2.0)); + r_jk / hbp->r0_hb - 2.0)); e_hb_thr += e_hb = hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4; @@ -170,51 +169,28 @@ namespace ReaxFF { /* hydrogen bond forces */ bo_ij->Cdbo += CEhb1; // dbo term - if (control->virial == 0) { - // dcos terms - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], +CEhb2, dcos_theta_di); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], +CEhb2, dcos_theta_dj); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], +CEhb2, dcos_theta_dk); - // dr terms - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], -CEhb3/r_jk, dvec_jk); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], +CEhb3/r_jk, dvec_jk); - } - else { - /* for pressure coupling, terms that are not related to bond order - derivatives are added directly into pressure vector/tensor */ - rvec_Scale(force, +CEhb2, dcos_theta_di); // dcos terms - rvec_Add(workspace->forceReduction[reductionOffset+i], force); - - - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], +CEhb2, dcos_theta_dj); - - ivec_Scale(rel_jk, hbond_list[pk].scl, nbr_jk->rel_box); - rvec_Scale(force, +CEhb2, dcos_theta_dk); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - // dr terms - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j],-CEhb3/r_jk, dvec_jk); - - rvec_Scale(force, CEhb3/r_jk, dvec_jk); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - } + // dcos terms + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], +CEhb2, dcos_theta_di); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], +CEhb2, dcos_theta_dj); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], +CEhb2, dcos_theta_dk); + // dr terms + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], -CEhb3/r_jk, dvec_jk); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], +CEhb3/r_jk, dvec_jk); /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { - rvec_ScaledSum(delij, 1., system->my_atoms[j].x, - -1., system->my_atoms[i].x); - rvec_ScaledSum(delkj, 1., system->my_atoms[j].x, - -1., system->my_atoms[k].x); + if (system->pair_ptr->vflag_either || system->pair_ptr->eflag_either) { + rvec_ScaledSum(delij, 1., system->my_atoms[j].x, -1., system->my_atoms[i].x); + rvec_ScaledSum(delkj, 1., system->my_atoms[j].x, -1., system->my_atoms[k].x); rvec_Scale(fi_tmp, CEhb2, dcos_theta_di); rvec_Scale(fk_tmp, CEhb2, dcos_theta_dk); rvec_ScaledAdd(fk_tmp, CEhb3/r_jk, dvec_jk); - pair_reax_ptr->ev_tally3_thr_proxy(system->pair_ptr,i,j,k,e_hb,0.0,fi_tmp,fk_tmp,delij,delkj,thr); + pair_reax_ptr->ev_tally3_thr_proxy(i,j,k,e_hb,0.0,fi_tmp,fk_tmp,delij,delkj,thr); } } } } - } } #if defined(_OPENMP) diff --git a/src/OPENMP/reaxc_multi_body_omp.cpp b/src/OPENMP/reaxc_multi_body_omp.cpp index 7d338511d6..6b7a7c76f4 100644 --- a/src/OPENMP/reaxc_multi_body_omp.cpp +++ b/src/OPENMP/reaxc_multi_body_omp.cpp @@ -118,8 +118,8 @@ namespace ReaxFF { if (numbonds > 0) workspace->CdDelta[i] += CElp; // lp - 1st term /* tally into per-atom energy */ - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, i, system->n, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( i, i, system->n, 1, e_lp, 0.0, 0.0, 0.0, 0.0, 0.0, thr); /* correction for C2 */ @@ -145,8 +145,8 @@ namespace ReaxFF { workspace->CdDelta[i] += deahu2dsbo; /* tally into per-atom energy */ - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, system->n, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( i, j, system->n, 1, e_lph, 0.0, 0.0, 0.0, 0.0, 0.0, thr); } } @@ -229,10 +229,10 @@ namespace ReaxFF { p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2; /* tally into per-atom energy */ - if (system->pair_ptr->evflag) { + if (system->pair_ptr->eflag_either) { eng_tmp = e_ov; if (numbonds > 0) eng_tmp+= e_un; - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, i, system->n, 1, + pair_reax_ptr->ev_tally_thr_proxy( i, i, system->n, 1, eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr); } diff --git a/src/OPENMP/reaxc_nonbonded_omp.cpp b/src/OPENMP/reaxc_nonbonded_omp.cpp index e6feb0b497..595731b005 100644 --- a/src/OPENMP/reaxc_nonbonded_omp.cpp +++ b/src/OPENMP/reaxc_nonbonded_omp.cpp @@ -64,7 +64,6 @@ namespace ReaxFF { double e_ele, e_vdW, e_core; const double SMALL = 0.0001; double e_lg, de_lg, r_ij5, r_ij6, re6; - rvec temp; two_body_parameters *twbp; far_neighbor_data *nbr_pj; @@ -203,23 +202,14 @@ namespace ReaxFF { rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x); f_tmp = -(CEvd + CEclmb); - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, + pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1, pe_vdw, e_ele, f_tmp, delij[0], delij[1], delij[2], thr); } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], - +(CEvd + CEclmb), nbr_pj->dvec); - } else { /* NPT, iNPT or sNPT */ - /* for pressure coupling, terms not related to bond order - derivatives are added directly into pressure vector/tensor */ - - rvec_Scale(temp, CEvd + CEclmb, nbr_pj->dvec); - rvec_ScaledAdd(workspace->f[reductionOffset+i], -1., temp); - rvec_Add(workspace->forceReduction[reductionOffset+j], temp); - } + rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], + +(CEvd + CEclmb), nbr_pj->dvec); } } } @@ -257,7 +247,6 @@ namespace ReaxFF { double e_vdW, e_ele; double CEvd, CEclmb; double f_tmp, delij[3]; - rvec temp; far_neighbor_data *nbr_pj; LR_lookup_table *t; @@ -332,26 +321,17 @@ namespace ReaxFF { CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q; /* tally into per-atom energy */ - if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) { + if (system->pair_ptr->evflag) { rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x); f_tmp = -(CEvd + CEclmb); - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, 1, e_vdW, e_ele, + pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1, e_vdW, e_ele, f_tmp, delij[0], delij[1], delij[2], thr); } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); - rvec_ScaledAdd(workspace->forceReduction[froffset+j], - +(CEvd + CEclmb), nbr_pj->dvec); - } else { // NPT, iNPT or sNPT - /* for pressure coupling, terms not related to bond order derivatives - are added directly into pressure vector/tensor */ - rvec_Scale(temp, CEvd + CEclmb, nbr_pj->dvec); - - rvec_ScaledAdd(workspace->f[i], -1., temp); - rvec_Add(workspace->forceReduction[froffset+j], temp); - } + rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); + rvec_ScaledAdd(workspace->forceReduction[froffset+j], + +(CEvd + CEclmb), nbr_pj->dvec); } } } diff --git a/src/OPENMP/reaxc_torsion_angles_omp.cpp b/src/OPENMP/reaxc_torsion_angles_omp.cpp index f3283ef93c..83f49e96de 100644 --- a/src/OPENMP/reaxc_torsion_angles_omp.cpp +++ b/src/OPENMP/reaxc_torsion_angles_omp.cpp @@ -87,8 +87,6 @@ namespace ReaxFF { double CEconj4, CEconj5, CEconj6; double e_tor, e_con; rvec dvec_li; - rvec force; - ivec rel_box_jl; four_body_header *fbh; four_body_parameters *fbp; bond_data *pbond_ij, *pbond_jk, *pbond_kl; @@ -319,72 +317,34 @@ namespace ReaxFF { //bo_kl->Cdbo += (CEtors6 + CEconj3); bo_kl->CdboReduction[tid] += (CEtors6 + CEconj3); - if (control->virial == 0) { - /* dcos_theta_ijk */ - rvec_ScaledAdd(workspace->f[j], - CEtors7 + CEconj4, p_ijk->dcos_dj); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], - CEtors7 + CEconj4, p_ijk->dcos_dk); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], - CEtors7 + CEconj4, p_ijk->dcos_di); + /* dcos_theta_ijk */ + rvec_ScaledAdd(workspace->f[j], + CEtors7 + CEconj4, p_ijk->dcos_dj); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], + CEtors7 + CEconj4, p_ijk->dcos_dk); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], + CEtors7 + CEconj4, p_ijk->dcos_di); - /* dcos_theta_jkl */ - rvec_ScaledAdd(workspace->f[j], - CEtors8 + CEconj5, p_jkl->dcos_di); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], - CEtors8 + CEconj5, p_jkl->dcos_dj); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l], - CEtors8 + CEconj5, p_jkl->dcos_dk); + /* dcos_theta_jkl */ + rvec_ScaledAdd(workspace->f[j], + CEtors8 + CEconj5, p_jkl->dcos_di); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], + CEtors8 + CEconj5, p_jkl->dcos_dj); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l], + CEtors8 + CEconj5, p_jkl->dcos_dk); - /* dcos_omega */ - rvec_ScaledAdd(workspace->f[j], - CEtors9 + CEconj6, dcos_omega_dj); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], - CEtors9 + CEconj6, dcos_omega_di); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], - CEtors9 + CEconj6, dcos_omega_dk); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l], - CEtors9 + CEconj6, dcos_omega_dl); - } - else { - ivec_Sum(rel_box_jl, pbond_jk->rel_box, pbond_kl->rel_box); - - /* dcos_theta_ijk */ - rvec_Scale(force, CEtors7 + CEconj4, p_ijk->dcos_dk); - rvec_Add(workspace->forceReduction[reductionOffset+i], force); - - rvec_ScaledAdd(workspace->f[j], - CEtors7 + CEconj4, p_ijk->dcos_dj); - - rvec_Scale(force, CEtors7 + CEconj4, p_ijk->dcos_di); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - - /* dcos_theta_jkl */ - rvec_ScaledAdd(workspace->f[j], - CEtors8 + CEconj5, p_jkl->dcos_di); - - rvec_Scale(force, CEtors8 + CEconj5, p_jkl->dcos_dj); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - - rvec_Scale(force, CEtors8 + CEconj5, p_jkl->dcos_dk); - rvec_Add(workspace->forceReduction[reductionOffset+l], force); - - /* dcos_omega */ - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_di); - rvec_Add(workspace->forceReduction[reductionOffset+i], force); - - rvec_ScaledAdd(workspace->f[j], - CEtors9 + CEconj6, dcos_omega_dj); - - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_dk); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_dl); - rvec_Add(workspace->forceReduction[reductionOffset+i], force); - } + /* dcos_omega */ + rvec_ScaledAdd(workspace->f[j], + CEtors9 + CEconj6, dcos_omega_dj); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], + CEtors9 + CEconj6, dcos_omega_di); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], + CEtors9 + CEconj6, dcos_omega_dk); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l], + CEtors9 + CEconj6, dcos_omega_dl); /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { + if (system->pair_ptr->evflag) { // acquire vectors rvec_ScaledSum(delil, 1., system->my_atoms[l].x, @@ -410,14 +370,13 @@ namespace ReaxFF { // tally eng_tmp = e_tor + e_con; - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, j, k, system->n, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( j, k, system->n, 1, eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr); - // NEED TO MAKE AN OMP VERSION OF THIS CALL! - if (system->pair_ptr->vflag_atom) - system->pair_ptr->v_tally4(i, j, k, l, fi_tmp, fj_tmp, fk_tmp, - delil, deljl, delkl); + if (system->pair_ptr->vflag_either) + pair_reax_ptr->v_tally4_thr_proxy(i, j, k, l, fi_tmp, fj_tmp, fk_tmp, + delil, deljl, delkl, thr); } } // pl check ends diff --git a/src/OPENMP/reaxc_valence_angles_omp.cpp b/src/OPENMP/reaxc_valence_angles_omp.cpp index 5785e3c966..ed457c96fe 100644 --- a/src/OPENMP/reaxc_valence_angles_omp.cpp +++ b/src/OPENMP/reaxc_valence_angles_omp.cpp @@ -140,7 +140,6 @@ namespace ReaxFF { double f7_ij, f7_jk, f8_Dj, f9_Dj; double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; double BOA_ij, BOA_jk; - rvec force; // Tallying variables double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; @@ -519,26 +518,12 @@ namespace ReaxFF { bo_jt->Cdbopi2 += CEval5; } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], - CEval8, p_ijk->dcos_di); - rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], - CEval8, p_ijk->dcos_dk); - } else { - /* terms not related to bond order derivatives are - added directly into forces and pressure vector/tensor */ - rvec_Scale(force, CEval8, p_ijk->dcos_di); - rvec_Add(workspace->forceReduction[reductionOffset+i], force); - - rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); - - rvec_Scale(force, CEval8, p_ijk->dcos_dk); - rvec_Add(workspace->forceReduction[reductionOffset+k], force); - } + rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], CEval8, p_ijk->dcos_di); + rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], CEval8, p_ijk->dcos_dk); /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { + if (system->pair_ptr->evflag) { /* Acquire vectors */ rvec_ScaledSum(delij, 1., system->my_atoms[i].x, @@ -552,12 +537,11 @@ namespace ReaxFF { eng_tmp = e_ang + e_pen + e_coa; - if (system->pair_ptr->evflag) - pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, j, j, system->N, 1, + if (system->pair_ptr->eflag_either) + pair_reax_ptr->ev_tally_thr_proxy( j, j, system->N, 1, eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr); - if (system->pair_ptr->vflag_atom) - // NEED TO MAKE AN OMP VERSION OF THIS CALL! - system->pair_ptr->v_tally3(i, j, k, fi_tmp, fk_tmp, delij, delkj); + if (system->pair_ptr->vflag_either) + pair_reax_ptr->v_tally3_thr_proxy(i, j, k, fi_tmp, fk_tmp, delij, delkj, thr); } } // if (p_val1>0.001) diff --git a/src/REAXFF/pair_reaxc.cpp b/src/REAXFF/pair_reaxc.cpp index 3b1d09a002..389b7d9e13 100644 --- a/src/REAXFF/pair_reaxc.cpp +++ b/src/REAXFF/pair_reaxc.cpp @@ -464,9 +464,6 @@ void PairReaxC::compute(int eflag, int vflag) evdwl = ecoul = 0.0; ev_init(eflag,vflag); - if (vflag_global) api->control->virial = 1; - else api->control->virial = 0; - api->system->n = atom->nlocal; // my atoms api->system->N = atom->nlocal + atom->nghost; // mine + ghosts api->system->bigN = static_cast (atom->natoms); // all atoms in the system diff --git a/src/REAXFF/reaxc_bond_orders.cpp b/src/REAXFF/reaxc_bond_orders.cpp index 3866c61110..0e4227922b 100644 --- a/src/REAXFF/reaxc_bond_orders.cpp +++ b/src/REAXFF/reaxc_bond_orders.cpp @@ -33,106 +33,7 @@ #include namespace ReaxFF { - void Add_dBond_to_Forces_NPT(int i, int pj, storage *workspace, reax_list **lists) - { - reax_list *bonds = (*lists) + BONDS; - bond_data *nbr_j, *nbr_k; - bond_order_data *bo_ij, *bo_ji; - dbond_coefficients coef; - rvec temp; - int pk, k, j; - - /* Initializations */ - nbr_j = &(bonds->select.bond_list[pj]); - j = nbr_j->nbr; - bo_ij = &(nbr_j->bo_data); - bo_ji = &(bonds->select.bond_list[nbr_j->sym_index].bo_data); - - coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - coef.C2dbo = bo_ij->C2dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - coef.C3dbo = bo_ij->C3dbo * (bo_ij->Cdbo + bo_ji->Cdbo); - - coef.C1dbopi = bo_ij->C1dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C2dbopi = bo_ij->C2dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C3dbopi = bo_ij->C3dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - coef.C4dbopi = bo_ij->C4dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); - - coef.C1dbopi2 = bo_ij->C1dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C2dbopi2 = bo_ij->C2dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C3dbopi2 = bo_ij->C3dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - coef.C4dbopi2 = bo_ij->C4dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); - - coef.C1dDelta = bo_ij->C1dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - coef.C2dDelta = bo_ij->C2dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]); - - - /************************************ - * forces related to atom i * - * first neighbors of atom i * - ************************************/ - for (pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk) { - nbr_k = &(bonds->select.bond_list[pk]); - k = nbr_k->nbr; - - rvec_Scale(temp, -coef.C2dbo, nbr_k->bo_data.dBOp); /*2nd, dBO*/ - rvec_ScaledAdd(temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);/*dDelta*/ - rvec_ScaledAdd(temp, -coef.C3dbopi, nbr_k->bo_data.dBOp); /*3rd, dBOpi*/ - rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);/*3rd, dBOpi2*/ - - /* force */ - rvec_Add(workspace->f[k], temp); - } - - /* then atom i itself */ - rvec_Scale(temp, coef.C1dbo, bo_ij->dBOp); /*1st,dBO*/ - rvec_ScaledAdd(temp, coef.C2dbo, workspace->dDeltap_self[i]); /*2nd,dBO*/ - rvec_ScaledAdd(temp, coef.C1dDelta, bo_ij->dBOp); /*1st,dBO*/ - rvec_ScaledAdd(temp, coef.C2dDelta, workspace->dDeltap_self[i]);/*2nd,dBO*/ - rvec_ScaledAdd(temp, coef.C1dbopi, bo_ij->dln_BOp_pi); /*1st,dBOpi*/ - rvec_ScaledAdd(temp, coef.C2dbopi, bo_ij->dBOp); /*2nd,dBOpi*/ - rvec_ScaledAdd(temp, coef.C3dbopi, workspace->dDeltap_self[i]);/*3rd,dBOpi*/ - - rvec_ScaledAdd(temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2); /*1st,dBO_pi2*/ - rvec_ScaledAdd(temp, coef.C2dbopi2, bo_ij->dBOp); /*2nd,dBO_pi2*/ - rvec_ScaledAdd(temp, coef.C3dbopi2, workspace->dDeltap_self[i]);/*3rd*/ - - /* force */ - rvec_Add(workspace->f[i], temp); - - for (pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk) { - nbr_k = &(bonds->select.bond_list[pk]); - k = nbr_k->nbr; - - rvec_Scale(temp, -coef.C3dbo, nbr_k->bo_data.dBOp); /*3rd,dBO*/ - rvec_ScaledAdd(temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);/*dDelta*/ - rvec_ScaledAdd(temp, -coef.C4dbopi, nbr_k->bo_data.dBOp); /*4th,dBOpi*/ - rvec_ScaledAdd(temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);/*4th,dBOpi2*/ - - /* force */ - rvec_Add(workspace->f[k], temp); - } - - /* then atom j itself */ - rvec_Scale(temp, -coef.C1dbo, bo_ij->dBOp); /*1st, dBO*/ - rvec_ScaledAdd(temp, coef.C3dbo, workspace->dDeltap_self[j]); /*2nd, dBO*/ - rvec_ScaledAdd(temp, -coef.C1dDelta, bo_ij->dBOp); /*1st, dBO*/ - rvec_ScaledAdd(temp, coef.C3dDelta, workspace->dDeltap_self[j]);/*2nd, dBO*/ - - rvec_ScaledAdd(temp, -coef.C1dbopi, bo_ij->dln_BOp_pi); /*1st,dBOpi*/ - rvec_ScaledAdd(temp, -coef.C2dbopi, bo_ij->dBOp); /*2nd,dBOpi*/ - rvec_ScaledAdd(temp, coef.C4dbopi, workspace->dDeltap_self[j]);/*3rd,dBOpi*/ - - rvec_ScaledAdd(temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2); /*1st,dBOpi2*/ - rvec_ScaledAdd(temp, -coef.C2dbopi2, bo_ij->dBOp); /*2nd,dBOpi2*/ - rvec_ScaledAdd(temp,coef.C4dbopi2,workspace->dDeltap_self[j]);/*3rd,dBOpi2*/ - - /* force */ - rvec_Add(workspace->f[j], temp); - } - - void Add_dBond_to_Forces(reax_system *system, int i, int pj, - storage *workspace, reax_list **lists) + void Add_dBond_to_Forces(reax_system *system, int i, int pj, storage *workspace, reax_list **lists) { reax_list *bonds = (*lists) + BONDS; bond_data *nbr_j, *nbr_k; @@ -180,10 +81,10 @@ namespace ReaxFF { rvec_ScaledAdd(temp, coef.C3dbopi2, workspace->dDeltap_self[i]); rvec_Add(workspace->f[i], temp); - if (system->pair_ptr->vflag_atom) { - rvec_Scale(fi_tmp, -1.0, temp); + if (system->pair_ptr->vflag_either) { + rvec_Scale(fi_tmp, -0.5, temp); rvec_ScaledSum(delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x); - system->pair_ptr->v_tally(i,fi_tmp,delij); + system->pair_ptr->v_tally2_newton(i,fi_tmp,delij); } // forces on j @@ -199,10 +100,10 @@ namespace ReaxFF { rvec_ScaledAdd(temp, coef.C4dbopi2, workspace->dDeltap_self[j]); rvec_Add(workspace->f[j], temp); - if (system->pair_ptr->vflag_atom) { - rvec_Scale(fj_tmp, -1.0, temp); + if (system->pair_ptr->vflag_either) { + rvec_Scale(fj_tmp, -0.5, temp); rvec_ScaledSum(delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x); - system->pair_ptr->v_tally(j,fj_tmp,delji); + system->pair_ptr->v_tally2_newton(j,fj_tmp,delji); } // forces on k: i neighbor @@ -216,12 +117,12 @@ namespace ReaxFF { rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp); rvec_Add(workspace->f[k], temp); - if (system->pair_ptr->vflag_atom) { - rvec_Scale(fk_tmp, -1.0, temp); + if (system->pair_ptr->vflag_either) { + rvec_Scale(fk_tmp, -0.5, temp); rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); - system->pair_ptr->v_tally(k,fk_tmp,delki); + system->pair_ptr->v_tally2_newton(k,fk_tmp,delki); rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); - system->pair_ptr->v_tally(k,fk_tmp,delkj); + system->pair_ptr->v_tally2_newton(k,fk_tmp,delkj); } } @@ -236,12 +137,12 @@ namespace ReaxFF { rvec_ScaledAdd(temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp); rvec_Add(workspace->f[k], temp); - if (system->pair_ptr->vflag_atom) { - rvec_Scale(fk_tmp, -1.0, temp); + if (system->pair_ptr->vflag_either) { + rvec_Scale(fk_tmp, -0.5, temp); rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); - system->pair_ptr->v_tally(k,fk_tmp,delki); + system->pair_ptr->v_tally2_newton(k,fk_tmp,delki); rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); - system->pair_ptr->v_tally(k,fk_tmp,delkj); + system->pair_ptr->v_tally2_newton(k,fk_tmp,delkj); } } } diff --git a/src/REAXFF/reaxc_bonds.cpp b/src/REAXFF/reaxc_bonds.cpp index 2c83ce8144..7feea043fa 100644 --- a/src/REAXFF/reaxc_bonds.cpp +++ b/src/REAXFF/reaxc_bonds.cpp @@ -32,8 +32,7 @@ #include namespace ReaxFF { - void Bonds(reax_system *system, simulation_data *data, - storage *workspace, reax_list **lists) + void Bonds(reax_system *system, simulation_data *data, storage *workspace, reax_list **lists) { int i, j, pj, natoms; int start_i, end_i; @@ -94,8 +93,8 @@ namespace ReaxFF { -twbp->De_p * bo_ij->BO_pi -twbp->De_pp * bo_ij->BO_pi2; - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(i,j,natoms,1,ebond,0.0,0.0,0.0,0.0,0.0); /* calculate derivatives of Bond Orders */ @@ -124,8 +123,8 @@ namespace ReaxFF { decobdboub = -gp10 * exphu * hulpov * (gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1)); - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(i,j,natoms,1,estriph,0.0,0.0,0.0,0.0,0.0); bo_ij->Cdbo += decobdbo; diff --git a/src/REAXFF/reaxc_control.cpp b/src/REAXFF/reaxc_control.cpp index 97fd65df9d..c84c3afe9d 100644 --- a/src/REAXFF/reaxc_control.cpp +++ b/src/REAXFF/reaxc_control.cpp @@ -71,7 +71,6 @@ namespace ReaxFF { /* assign default values */ control->nthreads = 1; control->tabulate = 0; - control->virial = 0; control->bond_cut = 5.0; control->bg_cut = 0.3; control->thb_cut = 0.001; diff --git a/src/REAXFF/reaxc_forces.cpp b/src/REAXFF/reaxc_forces.cpp index 5bb0dcca88..ec783102b4 100644 --- a/src/REAXFF/reaxc_forces.cpp +++ b/src/REAXFF/reaxc_forces.cpp @@ -47,7 +47,7 @@ namespace ReaxFF { Valence_Angles(system, control, data, workspace, lists); Torsion_Angles(system, control, data, workspace, lists); if (control->hbond_cut > 0) - Hydrogen_Bonds(system, control, data, workspace, lists); + Hydrogen_Bonds(system, data, workspace, lists); } static void Compute_NonBonded_Forces(reax_system *system, @@ -64,21 +64,15 @@ namespace ReaxFF { Tabulated_vdW_Coulomb_Energy(system, control, data, workspace, lists); } - static void Compute_Total_Force(reax_system *system, control_params *control, - storage *workspace, reax_list **lists) + static void Compute_Total_Force(reax_system *system, storage *workspace, reax_list **lists) { int i, pj; reax_list *bonds = (*lists) + BONDS; for (i = 0; i < system->N; ++i) for (pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj) - if (i < bonds->select.bond_list[pj].nbr) { - if (control->virial == 0) + if (i < bonds->select.bond_list[pj].nbr) Add_dBond_to_Forces(system, i, pj, workspace, lists); - else - Add_dBond_to_Forces_NPT(i, pj, workspace, lists); - } - } static void Validate_Lists(reax_system *system, reax_list **lists, @@ -387,6 +381,6 @@ namespace ReaxFF { Compute_NonBonded_Forces(system, control, data, workspace, lists); /*********** total force ***************/ - Compute_Total_Force(system, control, workspace, lists); + Compute_Total_Force(system, workspace, lists); } } diff --git a/src/REAXFF/reaxc_hydrogen_bonds.cpp b/src/REAXFF/reaxc_hydrogen_bonds.cpp index 3072968d58..5c7629ffb7 100644 --- a/src/REAXFF/reaxc_hydrogen_bonds.cpp +++ b/src/REAXFF/reaxc_hydrogen_bonds.cpp @@ -32,9 +32,7 @@ #include "pair.h" namespace ReaxFF { - void Hydrogen_Bonds(reax_system *system, control_params *control, - simulation_data *data, storage *workspace, - reax_list **lists) + void Hydrogen_Bonds(reax_system *system, simulation_data *data, storage *workspace, reax_list **lists) { int i, j, k, pi, pk; int type_i, type_j, type_k; @@ -42,11 +40,10 @@ namespace ReaxFF { int hblist[MAX_BONDS]; int itr, top; int num_hb_intrs = 0; - ivec rel_jk; double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk; - rvec dvec_jk, force; + rvec dvec_jk; hbond_parameters *hbp; bond_order_data *bo_ij; bond_data *pbond_ij; @@ -134,34 +131,16 @@ namespace ReaxFF { /* hydrogen bond forces */ bo_ij->Cdbo += CEhb1; // dbo term - if (control->virial == 0) { - // dcos terms - rvec_ScaledAdd(workspace->f[i], +CEhb2, dcos_theta_di); - rvec_ScaledAdd(workspace->f[j], +CEhb2, dcos_theta_dj); - rvec_ScaledAdd(workspace->f[k], +CEhb2, dcos_theta_dk); - // dr terms - rvec_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk); - rvec_ScaledAdd(workspace->f[k], +CEhb3/r_jk, dvec_jk); - } - else { - rvec_Scale(force, +CEhb2, dcos_theta_di); // dcos terms - rvec_Add(workspace->f[i], force); - - rvec_ScaledAdd(workspace->f[j], +CEhb2, dcos_theta_dj); - - ivec_Scale(rel_jk, hbond_list[pk].scl, nbr_jk->rel_box); - rvec_Scale(force, +CEhb2, dcos_theta_dk); - rvec_Add(workspace->f[k], force); - - // dr terms - rvec_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk); - - rvec_Scale(force, CEhb3/r_jk, dvec_jk); - rvec_Add(workspace->f[k], force); - } + // dcos terms + rvec_ScaledAdd(workspace->f[i], +CEhb2, dcos_theta_di); + rvec_ScaledAdd(workspace->f[j], +CEhb2, dcos_theta_dj); + rvec_ScaledAdd(workspace->f[k], +CEhb2, dcos_theta_dk); + // dr terms + rvec_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk); + rvec_ScaledAdd(workspace->f[k], +CEhb3/r_jk, dvec_jk); /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { + if (system->pair_ptr->evflag) { rvec_ScaledSum(delij, 1., system->my_atoms[j].x, -1., system->my_atoms[i].x); rvec_ScaledSum(delkj, 1., system->my_atoms[j].x, diff --git a/src/REAXFF/reaxc_multi_body.cpp b/src/REAXFF/reaxc_multi_body.cpp index bd1a5ccaa6..2390b54474 100644 --- a/src/REAXFF/reaxc_multi_body.cpp +++ b/src/REAXFF/reaxc_multi_body.cpp @@ -93,8 +93,8 @@ namespace ReaxFF { if (numbonds > 0 || control->enobondsflag) workspace->CdDelta[i] += CElp; // lp - 1st term - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0); /* correction for C2 */ @@ -119,8 +119,8 @@ namespace ReaxFF { bo_ij->Cdbo += deahu2dbo; workspace->CdDelta[i] += deahu2dsbo; - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0); } @@ -202,8 +202,8 @@ namespace ReaxFF { CEunder4 = CEunder1 * (dfvl*workspace->Delta_lp_temp[i]) * p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2; - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) { + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) { eng_tmp = e_ov; if (numbonds > 0 || control->enobondsflag) eng_tmp += e_un; diff --git a/src/REAXFF/reaxc_nonbonded.cpp b/src/REAXFF/reaxc_nonbonded.cpp index e86ffe3a8b..8a1e41da29 100644 --- a/src/REAXFF/reaxc_nonbonded.cpp +++ b/src/REAXFF/reaxc_nonbonded.cpp @@ -47,8 +47,8 @@ namespace ReaxFF { (system->reax_param.sbp[type_i].eta / 2.) * SQR(q)); data->my_en.e_pol += en_tmp; - /* tally into per-atom energy */ - if (system->pair_ptr->evflag) + /* tally energy into global or per-atom energy accumulators */ + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(i,i,system->n,1,0.0,en_tmp,0.0,0.0,0.0,0.0); } } @@ -67,7 +67,6 @@ namespace ReaxFF { double dr3gamij_1, dr3gamij_3; double e_ele, e_vdW, e_core, SMALL = 0.0001; double e_lg, de_lg, r_ij5, r_ij6, re6; - rvec temp; two_body_parameters *twbp; far_neighbor_data *nbr_pj; reax_list *far_nbrs; @@ -193,7 +192,7 @@ namespace ReaxFF { (dTap - Tap * r_ij / dr3gamij_1) / dr3gamij_3; /* tally into per-atom energy */ - if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) { + if (system->pair_ptr->evflag) { pe_vdw = Tap * (e_vdW + e_core + e_lg); rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x); @@ -202,15 +201,8 @@ namespace ReaxFF { f_tmp,delij[0],delij[1],delij[2]); } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); - rvec_ScaledAdd(workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec); - } else { /* NPT, iNPT or sNPT */ - rvec_Scale(temp, CEvd + CEclmb, nbr_pj->dvec); - - rvec_ScaledAdd(workspace->f[i], -1., temp); - rvec_Add(workspace->f[j], temp); - } + rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); + rvec_ScaledAdd(workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec); } } } @@ -231,7 +223,6 @@ namespace ReaxFF { double CEvd, CEclmb, SMALL = 0.0001; double f_tmp, delij[3]; - rvec temp; far_neighbor_data *nbr_pj; reax_list *far_nbrs; LR_lookup_table *t; @@ -301,7 +292,7 @@ namespace ReaxFF { CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q; /* tally into per-atom energy */ - if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) { + if (system->pair_ptr->evflag) { rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x); f_tmp = -(CEvd + CEclmb); @@ -309,15 +300,8 @@ namespace ReaxFF { f_tmp,delij[0],delij[1],delij[2]); } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); - rvec_ScaledAdd(workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec); - } else { // NPT, iNPT or sNPT - rvec_Scale(temp, CEvd + CEclmb, nbr_pj->dvec); - - rvec_ScaledAdd(workspace->f[i], -1., temp); - rvec_Add(workspace->f[j], temp); - } + rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); + rvec_ScaledAdd(workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec); } } } diff --git a/src/REAXFF/reaxc_torsion_angles.cpp b/src/REAXFF/reaxc_torsion_angles.cpp index adbb6ffc06..329f7b8a7a 100644 --- a/src/REAXFF/reaxc_torsion_angles.cpp +++ b/src/REAXFF/reaxc_torsion_angles.cpp @@ -144,8 +144,6 @@ namespace ReaxFF { double CEconj4, CEconj5, CEconj6; double e_tor, e_con; rvec dvec_li; - rvec force; - ivec rel_box_jl; four_body_header *fbh; four_body_parameters *fbp; bond_data *pbond_ij, *pbond_jk, *pbond_kl; @@ -351,81 +349,29 @@ namespace ReaxFF { bo_jk->Cdbo += (CEtors5 + CEconj2); bo_kl->Cdbo += (CEtors6 + CEconj3); - if (control->virial == 0) { - /* dcos_theta_ijk */ - rvec_ScaledAdd(workspace->f[i], - CEtors7 + CEconj4, p_ijk->dcos_dk); - rvec_ScaledAdd(workspace->f[j], - CEtors7 + CEconj4, p_ijk->dcos_dj); - rvec_ScaledAdd(workspace->f[k], - CEtors7 + CEconj4, p_ijk->dcos_di); + /* dcos_theta_ijk */ + rvec_ScaledAdd(workspace->f[i], CEtors7 + CEconj4, p_ijk->dcos_dk); + rvec_ScaledAdd(workspace->f[j], CEtors7 + CEconj4, p_ijk->dcos_dj); + rvec_ScaledAdd(workspace->f[k], CEtors7 + CEconj4, p_ijk->dcos_di); - /* dcos_theta_jkl */ - rvec_ScaledAdd(workspace->f[j], - CEtors8 + CEconj5, p_jkl->dcos_di); - rvec_ScaledAdd(workspace->f[k], - CEtors8 + CEconj5, p_jkl->dcos_dj); - rvec_ScaledAdd(workspace->f[l], - CEtors8 + CEconj5, p_jkl->dcos_dk); + /* dcos_theta_jkl */ + rvec_ScaledAdd(workspace->f[j], CEtors8 + CEconj5, p_jkl->dcos_di); + rvec_ScaledAdd(workspace->f[k], CEtors8 + CEconj5, p_jkl->dcos_dj); + rvec_ScaledAdd(workspace->f[l], CEtors8 + CEconj5, p_jkl->dcos_dk); - /* dcos_omega */ - rvec_ScaledAdd(workspace->f[i], - CEtors9 + CEconj6, dcos_omega_di); - rvec_ScaledAdd(workspace->f[j], - CEtors9 + CEconj6, dcos_omega_dj); - rvec_ScaledAdd(workspace->f[k], - CEtors9 + CEconj6, dcos_omega_dk); - rvec_ScaledAdd(workspace->f[l], - CEtors9 + CEconj6, dcos_omega_dl); - } - else { - ivec_Sum(rel_box_jl, pbond_jk->rel_box, pbond_kl->rel_box); - - /* dcos_theta_ijk */ - rvec_Scale(force, CEtors7 + CEconj4, p_ijk->dcos_dk); - rvec_Add(workspace->f[i], force); - - rvec_ScaledAdd(workspace->f[j], - CEtors7 + CEconj4, p_ijk->dcos_dj); - - rvec_Scale(force, CEtors7 + CEconj4, p_ijk->dcos_di); - rvec_Add(workspace->f[k], force); - - /* dcos_theta_jkl */ - rvec_ScaledAdd(workspace->f[j], - CEtors8 + CEconj5, p_jkl->dcos_di); - - rvec_Scale(force, CEtors8 + CEconj5, p_jkl->dcos_dj); - rvec_Add(workspace->f[k], force); - - rvec_Scale(force, CEtors8 + CEconj5, p_jkl->dcos_dk); - rvec_Add(workspace->f[l], force); - - - /* dcos_omega */ - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_di); - rvec_Add(workspace->f[i], force); - - rvec_ScaledAdd(workspace->f[j], - CEtors9 + CEconj6, dcos_omega_dj); - - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_dk); - rvec_Add(workspace->f[k], force); - - rvec_Scale(force, CEtors9 + CEconj6, dcos_omega_dl); - rvec_Add(workspace->f[l], force); - } + /* dcos_omega */ + rvec_ScaledAdd(workspace->f[i], CEtors9 + CEconj6, dcos_omega_di); + rvec_ScaledAdd(workspace->f[j], CEtors9 + CEconj6, dcos_omega_dj); + rvec_ScaledAdd(workspace->f[k], CEtors9 + CEconj6, dcos_omega_dk); + rvec_ScaledAdd(workspace->f[l], CEtors9 + CEconj6, dcos_omega_dl); /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { + if (system->pair_ptr->evflag) { // acquire vectors - rvec_ScaledSum(delil, 1., system->my_atoms[l].x, - -1., system->my_atoms[i].x); - rvec_ScaledSum(deljl, 1., system->my_atoms[l].x, - -1., system->my_atoms[j].x); - rvec_ScaledSum(delkl, 1., system->my_atoms[l].x, - -1., system->my_atoms[k].x); + rvec_ScaledSum(delil, 1., system->my_atoms[l].x, -1., system->my_atoms[i].x); + rvec_ScaledSum(deljl, 1., system->my_atoms[l].x, -1., system->my_atoms[j].x); + rvec_ScaledSum(delkl, 1., system->my_atoms[l].x, -1., system->my_atoms[k].x); // dcos_theta_ijk rvec_Scale(fi_tmp, CEtors7 + CEconj4, p_ijk->dcos_dk); rvec_Scale(fj_tmp, CEtors7 + CEconj4, p_ijk->dcos_dj); @@ -442,9 +388,9 @@ namespace ReaxFF { // tally eng_tmp = e_tor + e_con; - if (system->pair_ptr->evflag) + if (system->pair_ptr->eflag_either) system->pair_ptr->ev_tally(j,k,natoms,1,eng_tmp,0.0,0.0,0.0,0.0,0.0); - if (system->pair_ptr->vflag_atom) + if (system->pair_ptr->vflag_either) system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl); } } // pl check ends diff --git a/src/REAXFF/reaxc_valence_angles.cpp b/src/REAXFF/reaxc_valence_angles.cpp index 6f609646e3..5f75b6637d 100644 --- a/src/REAXFF/reaxc_valence_angles.cpp +++ b/src/REAXFF/reaxc_valence_angles.cpp @@ -67,9 +67,8 @@ namespace ReaxFF { } - void Valence_Angles(reax_system *system, control_params *control, - simulation_data *data, storage *workspace, - reax_list **lists) + void Valence_Angles(reax_system *system, control_params *control, simulation_data *data, + storage *workspace, reax_list **lists) { int i, j, pi, k, pk, t; int type_i, type_j, type_k; @@ -92,7 +91,6 @@ namespace ReaxFF { double f7_ij, f7_jk, f8_Dj, f9_Dj; double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; double BOA_ij, BOA_jk; - rvec force; // Tallying variables double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; @@ -114,7 +112,7 @@ namespace ReaxFF { num_thb_intrs = 0; - for (j = 0; j < system->N; ++j) { // Ray: the first one with system->N + for (j = 0; j < system->N; ++j) { type_j = system->my_atoms[j].type; if (type_j < 0) continue; start_j = Start_Index(j, bonds); @@ -218,7 +216,6 @@ namespace ReaxFF { ++num_thb_intrs; - if ((j < system->n) && (BOA_jk > 0.0) && (bo_ij->BO > control->thb_cut) && (bo_jk->BO > control->thb_cut) && @@ -347,39 +344,26 @@ namespace ReaxFF { bo_jt->Cdbopi2 += CEval5; } - if (control->virial == 0) { - rvec_ScaledAdd(workspace->f[i], CEval8, p_ijk->dcos_di); - rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); - rvec_ScaledAdd(workspace->f[k], CEval8, p_ijk->dcos_dk); - } else { - rvec_Scale(force, CEval8, p_ijk->dcos_di); - rvec_Add(workspace->f[i], force); + rvec_ScaledAdd(workspace->f[i], CEval8, p_ijk->dcos_di); + rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); + rvec_ScaledAdd(workspace->f[k], CEval8, p_ijk->dcos_dk); - rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); - - rvec_Scale(force, CEval8, p_ijk->dcos_dk); - rvec_Add(workspace->f[k], force); + /* tally energy */ + if (system->pair_ptr->eflag_either) { + eng_tmp = e_ang + e_pen + e_coa; + system->pair_ptr->ev_tally(j,j,system->N,1,eng_tmp,0.0,0.0,0.0,0.0,0.0); } - /* tally into per-atom virials */ - if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { - + /* tally virial */ + if (system->pair_ptr->vflag_either) { + /* Acquire vectors */ - rvec_ScaledSum(delij, 1., system->my_atoms[i].x, - -1., system->my_atoms[j].x); - rvec_ScaledSum(delkj, 1., system->my_atoms[k].x, - -1., system->my_atoms[j].x); - + rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x); + rvec_ScaledSum(delkj, 1., system->my_atoms[k].x, -1., system->my_atoms[j].x); rvec_Scale(fi_tmp, -CEval8, p_ijk->dcos_di); rvec_Scale(fj_tmp, -CEval8, p_ijk->dcos_dj); rvec_Scale(fk_tmp, -CEval8, p_ijk->dcos_dk); - - eng_tmp = e_ang + e_pen + e_coa; - - if (system->pair_ptr->evflag) - system->pair_ptr->ev_tally(j,j,system->N,1,eng_tmp,0.0,0.0,0.0,0.0,0.0); - if (system->pair_ptr->vflag_atom) - system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj); + system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj); } } } diff --git a/src/REAXFF/reaxff_api.h b/src/REAXFF/reaxff_api.h index 96febe70f9..6179da9c72 100644 --- a/src/REAXFF/reaxff_api.h +++ b/src/REAXFF/reaxff_api.h @@ -54,7 +54,6 @@ namespace ReaxFF single_body_parameters *, single_body_parameters *, two_body_parameters *); extern void Add_dBond_to_Forces(reax_system*, int, int, storage*, reax_list**); - extern void Add_dBond_to_Forces_NPT(int, int, storage*, reax_list**); // bonds @@ -78,8 +77,8 @@ namespace ReaxFF // hydrogen bonds - extern void Hydrogen_Bonds(reax_system *, control_params *, - simulation_data *, storage *, reax_list **); + extern void Hydrogen_Bonds(reax_system *, simulation_data *, storage *, reax_list **); + // init md extern void Init_System(reax_system *, control_params *); diff --git a/src/REAXFF/reaxff_types.h b/src/REAXFF/reaxff_types.h index c88e0f5a6a..aec48072df 100644 --- a/src/REAXFF/reaxff_types.h +++ b/src/REAXFF/reaxff_types.h @@ -240,7 +240,6 @@ namespace ReaxFF double thb_cutsq; int tabulate; - int virial; int lgflag; int enobondsflag; diff --git a/src/USER-MISC/pair_drip.cpp b/src/USER-MISC/pair_drip.cpp index 97f5ca439e..a59d74bab9 100644 --- a/src/USER-MISC/pair_drip.cpp +++ b/src/USER-MISC/pair_drip.cpp @@ -366,10 +366,7 @@ void PairDRIP::compute(int eflag, int vflag) f[j][0] += fj[0]; f[j][1] += fj[1]; f[j][2] += fj[2]; - - // multiply 2 since v_tally has a 0.5 coeff - fj[0] *= 2; fj[1] *= 2; fj[2] *= 2; - if (vflag_atom) v_tally(j, fj, x[j]); + if (vflag_either) v_tally2_newton(j, fj, x[j]); } } //loop over jj @@ -377,10 +374,7 @@ void PairDRIP::compute(int eflag, int vflag) f[i][0] += fi[0]; f[i][1] += fi[1]; f[i][2] += fi[2]; - - // multiply 2 since v_tally has a 0.5 coeff - fi[0] *= 2; fi[1] *= 2; fi[2] *= 2; - if (vflag_atom) v_tally(i, fi, x[i]); + if (vflag_either) v_tally2_newton(i, fi, x[i]); } // loop over ii @@ -529,22 +523,13 @@ double PairDRIP::calc_repulsive(int const i, int const j, Param& p, f[nbj3][k] += fnbj3[k]; } - if (vflag_atom) { - // multiply since v_tally has a 0.5 coeff - for (int k = 0; k < DIM; k++) { - fnbi1[k] *= 2; - fnbi2[k] *= 2; - fnbi3[k] *= 2; - fnbj1[k] *= 2; - fnbj2[k] *= 2; - fnbj3[k] *= 2; - } - v_tally(nbi1, fnbi1, x[nbi1]); - v_tally(nbi2, fnbi2, x[nbi2]); - v_tally(nbi3, fnbi3, x[nbi3]); - v_tally(nbj1, fnbj1, x[nbj1]); - v_tally(nbj2, fnbj2, x[nbj2]); - v_tally(nbj3, fnbj3, x[nbj3]); + if (vflag_either) { + v_tally2_newton(nbi1, fnbi1, x[nbi1]); + v_tally2_newton(nbi2, fnbi2, x[nbi2]); + v_tally2_newton(nbi3, fnbi3, x[nbi3]); + v_tally2_newton(nbj1, fnbj1, x[nbj1]); + v_tally2_newton(nbj2, fnbj2, x[nbj2]); + v_tally2_newton(nbj3, fnbj3, x[nbj3]); } return phi; diff --git a/src/pair.cpp b/src/pair.cpp index 5e8c666743..bc8cb2ecd1 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -991,7 +991,7 @@ void Pair::ev_unset() } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators need i < nlocal test since called by bond_quartic and dihedral_charmm ------------------------------------------------------------------------- */ @@ -1092,7 +1092,7 @@ void Pair::ev_tally(int i, int j, int nlocal, int newton_pair, } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators can use this version with full neighbor lists ------------------------------------------------------------------------- */ @@ -1138,7 +1138,7 @@ void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz ------------------------------------------------------------------------- */ @@ -1232,7 +1232,7 @@ void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair, } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz called when using full neighbor lists ------------------------------------------------------------------------- */ @@ -1285,7 +1285,7 @@ void Pair::ev_tally_xyz_full(int i, double evdwl, double ecoul, } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators called by SW and hbond potentials, newton_pair is always on virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ @@ -1342,7 +1342,7 @@ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, } /* ---------------------------------------------------------------------- - tally eng_vdwl and virial into global and per-atom accumulators + tally eng_vdwl and virial into global or per-atom accumulators called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ @@ -1487,28 +1487,40 @@ void Pair::ev_tally_tip4p(int key, int *list, double *v, } /* ---------------------------------------------------------------------- - tally virial into per-atom accumulators - called by REAX/C potential, newton_pair is always on - fi is magnitude of force on atom i + tally virial into global or per-atom accumulators + called by ReaxFF potential, newton_pair is always on + fi is magnitude of force on atom i, deli is the direction + note that the other atom (j) is not updated, due to newton on ------------------------------------------------------------------------- */ -void Pair::v_tally(int i, double *fi, double *deli) +void Pair::v_tally2_newton(int i, double *fi, double *deli) { double v[6]; - v[0] = 0.5*deli[0]*fi[0]; - v[1] = 0.5*deli[1]*fi[1]; - v[2] = 0.5*deli[2]*fi[2]; - v[3] = 0.5*deli[0]*fi[1]; - v[4] = 0.5*deli[0]*fi[2]; - v[5] = 0.5*deli[1]*fi[2]; + v[0] = deli[0]*fi[0]; + v[1] = deli[1]*fi[1]; + v[2] = deli[2]*fi[2]; + v[3] = deli[0]*fi[1]; + v[4] = deli[0]*fi[2]; + v[5] = deli[1]*fi[2]; - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + } } /* ---------------------------------------------------------------------- - tally virial into per-atom accumulators + tally virial into global or per-atom accumulators called by AIREBO potential, newton_pair is always on fpair is magnitude of force on atom I ------------------------------------------------------------------------- */ @@ -1549,7 +1561,7 @@ void Pair::v_tally2(int i, int j, double fpair, double *drij) /* ---------------------------------------------------------------------- tally virial into per-atom accumulators - called by AIREBO and Tersoff potential, newton_pair is always on + called by AIREBO and Tersoff potentials, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, double *drjk) @@ -1589,8 +1601,8 @@ void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, d } /* ---------------------------------------------------------------------- - tally virial into per-atom accumulators - called by AIREBO potential, newton_pair is always on + tally virial into global or per-atom accumulators + called by AIREBO potential, Tersoff, ReaxFF potentials, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally4(int i, int j, int k, int m, @@ -1634,7 +1646,7 @@ void Pair::v_tally4(int i, int j, int k, int m, } /* ---------------------------------------------------------------------- - tally virial into global and per-atom accumulators + tally virial into global or per-atom accumulators called by pair lubricate potential with 6 tensor components ------------------------------------------------------------------------- */ diff --git a/src/pair.h b/src/pair.h index ad158bba9e..c4a5d3678c 100644 --- a/src/pair.h +++ b/src/pair.h @@ -137,9 +137,9 @@ class Pair : protected Pointers { // need to be public, so can be called by pair_style reaxc - void v_tally(int, double *, double *); void ev_tally(int, int, int, int, double, double, double, double, double, double); void ev_tally3(int, int, int, double, double, double *, double *, double *, double *); + void v_tally2_newton(int, double *, double *); void v_tally3(int, int, int, double *, double *, double *, double *); void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double,