ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

recover virial and nofdotr related changes by @athomps and @akohlmey

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-07-12 15:58:27 -04:00
parent 757e2f8cff
commit f4239530bd
26 changed files with 244 additions and 675 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/ML-IAP/pair_mliap.cpp
View File

@ -264,8 +264,7 @@ void PairMLIAP::e_tally(MLIAPData* data)
add virial contribution into global and per-atom accumulators add virial contribution into global and per-atom accumulators
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMLIAP::v_tally(int i, int j, void PairMLIAP::v_tally(int i, int j, double *fij, double *rij)
double *fij, double *rij)
{ {
double v[6]; double v[6];

3
src/OPENMP/pair_reaxc_omp.cpp
View File

@ -233,9 +233,6 @@ void PairReaxCOMP::compute(int eflag, int vflag)
evdwl = ecoul = 0.0; evdwl = ecoul = 0.0;
ev_init(eflag,vflag); ev_init(eflag,vflag);
if (vflag_global) api->control->virial = 1;
else api->control->virial = 0;
if (vflag_atom) if (vflag_atom)
error->all(FLERR,"Pair style reax/c/omp does not support " error->all(FLERR,"Pair style reax/c/omp does not support "
"computing per-atom stress"); "computing per-atom stress");

28
src/OPENMP/pair_reaxc_omp.h
View File

@ -47,33 +47,49 @@ class PairReaxCOMP : public PairReaxC, public ThrOMP {
reduce_thr(styleparm, eflagparm, vflagparm, thrparm); 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 int nlocalparm, const int newton_pairparm,
const double evdwlparm, const double ecoulparm, const double evdwlparm, const double ecoulparm,
const double fpairparm, const double delxparm, const double fpairparm, const double delxparm,
const double delyparm, const double delzparm, const double delyparm, const double delzparm,
ThrData *const thrparm) 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); 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 int nlocalparm, const int newton_pairparm,
const double evdwlparm, const double ecoulparm, const double evdwlparm, const double ecoulparm,
const double fxparm, const double fyparm, const double fzparm, const double fxparm, const double fyparm, const double fzparm,
const double delxparm, const double delyparm, const double delxparm, const double delyparm,
const double delzparm, ThrData *const thrparm) 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); 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, double ecoul, double *fj, double *fk, double *drji, double *drki,
ThrData *const thrparm) 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: protected:

123
src/OPENMP/reaxc_bond_orders_omp.cpp
View File

@ -37,109 +37,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
namespace ReaxFF { 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, void Add_dBond_to_ForcesOMP(reax_system *system, int i, int pj,
storage *workspace, reax_list **lists) { storage *workspace, reax_list **lists) {
reax_list *bonds = (*lists) + BONDS; reax_list *bonds = (*lists) + BONDS;
@ -196,11 +93,11 @@ namespace ReaxFF {
rvec_Add(workspace->forceReduction[reductionOffset+i],temp); 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_Scale(fi_tmp, -1.0, temp);
rvec_ScaledSum(delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x); 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], fi_tmp[0],fi_tmp[1],fi_tmp[2],
delij[0],delij[1],delij[2],thr); delij[0],delij[1],delij[2],thr);
} }
@ -217,11 +114,11 @@ namespace ReaxFF {
rvec_Add(workspace->forceReduction[reductionOffset+j],temp); 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_Scale(fj_tmp, -1.0, temp);
rvec_ScaledSum(delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x); 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], fj_tmp[0],fj_tmp[1],fj_tmp[2],
delji[0],delji[1],delji[2],thr); delji[0],delji[1],delji[2],thr);
} }
@ -236,16 +133,16 @@ namespace ReaxFF {
rvec_Add(workspace->forceReduction[reductionOffset+k],temp); 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_Scale(fk_tmp, -1.0, temp);
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); 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], fk_tmp[0],fk_tmp[1],fk_tmp[2],
delki[0],delki[1],delki[2],thr); delki[0],delki[1],delki[2],thr);
rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); 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], fk_tmp[0],fk_tmp[1],fk_tmp[2],
delkj[0],delkj[1],delkj[2],thr); delkj[0],delkj[1],delkj[2],thr);
} }
@ -261,17 +158,17 @@ namespace ReaxFF {
rvec_Add(workspace->forceReduction[reductionOffset+k],temp); 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_Scale(fk_tmp, -1.0, temp);
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); 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], fk_tmp[0],fk_tmp[1],fk_tmp[2],
delki[0],delki[1],delki[2],thr); delki[0],delki[1],delki[2],thr);
rvec_ScaledSum(delkj,1.,system->my_atoms[k].x,-1.,system->my_atoms[j].x); 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], fk_tmp[0],fk_tmp[1],fk_tmp[2],
delkj[0],delkj[1],delkj[2],thr); delkj[0],delkj[1],delkj[2],thr);
} }

8
src/OPENMP/reaxc_bonds_omp.cpp
View File

@ -121,8 +121,8 @@ namespace ReaxFF {
-twbp->De_pp * bo_ij->BO_pi2; -twbp->De_pp * bo_ij->BO_pi2;
/* tally into per-atom energy */ /* tally into per-atom energy */
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, 1, pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1,
ebond, 0.0, 0.0, 0.0, 0.0, 0.0, thr); ebond, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
/* calculate derivatives of Bond Orders */ /* calculate derivatives of Bond Orders */
@ -152,8 +152,8 @@ namespace ReaxFF {
(gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1)); (gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
/* tally into per-atom energy */ /* tally into per-atom energy */
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, natoms, 1, pair_reax_ptr->ev_tally_thr_proxy( i, j, natoms, 1,
estriph, 0.0, 0.0, 0.0, 0.0, 0.0, thr); estriph, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
bo_ij->Cdbo += decobdbo; bo_ij->Cdbo += decobdbo;

17
src/OPENMP/reaxc_forces_omp.cpp
View File

@ -119,8 +119,6 @@ namespace ReaxFF {
} }
} }
if (control->virial == 0) {
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp for schedule(dynamic,50) #pragma omp for schedule(dynamic,50)
#endif #endif
@ -132,21 +130,6 @@ namespace ReaxFF {
Add_dBond_to_ForcesOMP(system, i, pj, workspace, lists); 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)
pair_reax_ptr->reduce_thr_proxy(system->pair_ptr, 0, 1, thr); pair_reax_ptr->reduce_thr_proxy(system->pair_ptr, 0, 1, thr);
#if defined(_OPENMP) #if defined(_OPENMP)

34
src/OPENMP/reaxc_hydrogen_bonds_omp.cpp
View File

@ -58,11 +58,10 @@ namespace ReaxFF {
int hblist[MAX_BONDS]; int hblist[MAX_BONDS];
int itr, top; int itr, top;
int num_hb_intrs = 0; int num_hb_intrs = 0;
ivec rel_jk;
double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; 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; 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 dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
rvec dvec_jk, force; rvec dvec_jk;
hbond_parameters *hbp; hbond_parameters *hbp;
bond_order_data *bo_ij; bond_order_data *bo_ij;
bond_data *pbond_ij; bond_data *pbond_ij;
@ -170,7 +169,6 @@ namespace ReaxFF {
/* hydrogen bond forces */ /* hydrogen bond forces */
bo_ij->Cdbo += CEhb1; // dbo term bo_ij->Cdbo += CEhb1; // dbo term
if (control->virial == 0) {
// dcos terms // dcos terms
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], +CEhb2, dcos_theta_di); 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+j], +CEhb2, dcos_theta_dj);
@ -178,43 +176,21 @@ namespace ReaxFF {
// dr terms // dr terms
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], -CEhb3/r_jk, dvec_jk); rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], -CEhb3/r_jk, dvec_jk);
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], +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);
}
/* tally into per-atom virials */ /* tally into per-atom virials */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { if (system->pair_ptr->vflag_either || system->pair_ptr->eflag_either) {
rvec_ScaledSum(delij, 1., system->my_atoms[j].x, rvec_ScaledSum(delij, 1., system->my_atoms[j].x, -1., system->my_atoms[i].x);
-1., system->my_atoms[i].x); rvec_ScaledSum(delkj, 1., system->my_atoms[j].x, -1., system->my_atoms[k].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(fi_tmp, CEhb2, dcos_theta_di);
rvec_Scale(fk_tmp, CEhb2, dcos_theta_dk); rvec_Scale(fk_tmp, CEhb2, dcos_theta_dk);
rvec_ScaledAdd(fk_tmp, CEhb3/r_jk, dvec_jk); 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) #if defined(_OPENMP)

12
src/OPENMP/reaxc_multi_body_omp.cpp
View File

@ -118,8 +118,8 @@ namespace ReaxFF {
if (numbonds > 0) workspace->CdDelta[i] += CElp; // lp - 1st term if (numbonds > 0) workspace->CdDelta[i] += CElp; // lp - 1st term
/* tally into per-atom energy */ /* tally into per-atom energy */
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
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,
e_lp, 0.0, 0.0, 0.0, 0.0, 0.0, thr); e_lp, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
/* correction for C2 */ /* correction for C2 */
@ -145,8 +145,8 @@ namespace ReaxFF {
workspace->CdDelta[i] += deahu2dsbo; workspace->CdDelta[i] += deahu2dsbo;
/* tally into per-atom energy */ /* tally into per-atom energy */
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, system->n, 1, 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); 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; p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
/* tally into per-atom energy */ /* tally into per-atom energy */
if (system->pair_ptr->evflag) { if (system->pair_ptr->eflag_either) {
eng_tmp = e_ov; eng_tmp = e_ov;
if (numbonds > 0) eng_tmp+= e_un; 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); eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
} }

26
src/OPENMP/reaxc_nonbonded_omp.cpp
View File

@ -64,7 +64,6 @@ namespace ReaxFF {
double e_ele, e_vdW, e_core; double e_ele, e_vdW, e_core;
const double SMALL = 0.0001; const double SMALL = 0.0001;
double e_lg, de_lg, r_ij5, r_ij6, re6; double e_lg, de_lg, r_ij5, r_ij6, re6;
rvec temp;
two_body_parameters *twbp; two_body_parameters *twbp;
far_neighbor_data *nbr_pj; far_neighbor_data *nbr_pj;
@ -203,23 +202,14 @@ namespace ReaxFF {
rvec_ScaledSum(delij, 1., system->my_atoms[i].x, rvec_ScaledSum(delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x); -1., system->my_atoms[j].x);
f_tmp = -(CEvd + CEclmb); 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, 1, pe_vdw, e_ele, f_tmp,
delij[0], delij[1], delij[2], thr); delij[0], delij[1], delij[2], thr);
} }
if (control->virial == 0) {
rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec); rvec_ScaledAdd(workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec);
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j],
+(CEvd + CEclmb), nbr_pj->dvec); +(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);
}
} }
} }
} }
@ -257,7 +247,6 @@ namespace ReaxFF {
double e_vdW, e_ele; double e_vdW, e_ele;
double CEvd, CEclmb; double CEvd, CEclmb;
double f_tmp, delij[3]; double f_tmp, delij[3];
rvec temp;
far_neighbor_data *nbr_pj; far_neighbor_data *nbr_pj;
LR_lookup_table *t; LR_lookup_table *t;
@ -332,26 +321,17 @@ namespace ReaxFF {
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q; CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
/* tally into per-atom energy */ /* 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, rvec_ScaledSum(delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x); -1., system->my_atoms[j].x);
f_tmp = -(CEvd + CEclmb); 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); 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->f[i], -(CEvd + CEclmb), nbr_pj->dvec);
rvec_ScaledAdd(workspace->forceReduction[froffset+j], rvec_ScaledAdd(workspace->forceReduction[froffset+j],
+(CEvd + CEclmb), nbr_pj->dvec); +(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);
}
} }
} }
} }

53
src/OPENMP/reaxc_torsion_angles_omp.cpp
View File

@ -87,8 +87,6 @@ namespace ReaxFF {
double CEconj4, CEconj5, CEconj6; double CEconj4, CEconj5, CEconj6;
double e_tor, e_con; double e_tor, e_con;
rvec dvec_li; rvec dvec_li;
rvec force;
ivec rel_box_jl;
four_body_header *fbh; four_body_header *fbh;
four_body_parameters *fbp; four_body_parameters *fbp;
bond_data *pbond_ij, *pbond_jk, *pbond_kl; bond_data *pbond_ij, *pbond_jk, *pbond_kl;
@ -319,7 +317,6 @@ namespace ReaxFF {
//bo_kl->Cdbo += (CEtors6 + CEconj3); //bo_kl->Cdbo += (CEtors6 + CEconj3);
bo_kl->CdboReduction[tid] += (CEtors6 + CEconj3); bo_kl->CdboReduction[tid] += (CEtors6 + CEconj3);
if (control->virial == 0) {
/* dcos_theta_ijk */ /* dcos_theta_ijk */
rvec_ScaledAdd(workspace->f[j], rvec_ScaledAdd(workspace->f[j],
CEtors7 + CEconj4, p_ijk->dcos_dj); CEtors7 + CEconj4, p_ijk->dcos_dj);
@ -345,46 +342,9 @@ namespace ReaxFF {
CEtors9 + CEconj6, dcos_omega_dk); CEtors9 + CEconj6, dcos_omega_dk);
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l], rvec_ScaledAdd(workspace->forceReduction[reductionOffset+l],
CEtors9 + CEconj6, dcos_omega_dl); 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);
}
/* tally into per-atom virials */ /* tally into per-atom virials */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { if (system->pair_ptr->evflag) {
// acquire vectors // acquire vectors
rvec_ScaledSum(delil, 1., system->my_atoms[l].x, rvec_ScaledSum(delil, 1., system->my_atoms[l].x,
@ -410,14 +370,13 @@ namespace ReaxFF {
// tally // tally
eng_tmp = e_tor + e_con; eng_tmp = e_tor + e_con;
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, j, k, system->n, 1, 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); 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_either)
if (system->pair_ptr->vflag_atom) pair_reax_ptr->v_tally4_thr_proxy(i, j, k, l, fi_tmp, fj_tmp, fk_tmp,
system->pair_ptr->v_tally4(i, j, k, l, fi_tmp, fj_tmp, fk_tmp, delil, deljl, delkl, thr);
delil, deljl, delkl);
} }
} // pl check ends } // pl check ends

30
src/OPENMP/reaxc_valence_angles_omp.cpp
View File

@ -140,7 +140,6 @@ namespace ReaxFF {
double f7_ij, f7_jk, f8_Dj, f9_Dj; double f7_ij, f7_jk, f8_Dj, f9_Dj;
double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
double BOA_ij, BOA_jk; double BOA_ij, BOA_jk;
rvec force;
// Tallying variables // Tallying variables
double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
@ -519,26 +518,12 @@ namespace ReaxFF {
bo_jt->Cdbopi2 += CEval5; bo_jt->Cdbopi2 += CEval5;
} }
if (control->virial == 0) {
rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj);
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], CEval8, p_ijk->dcos_di);
CEval8, p_ijk->dcos_di); rvec_ScaledAdd(workspace->forceReduction[reductionOffset+k], CEval8, p_ijk->dcos_dk);
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);
}
/* tally into per-atom virials */ /* tally into per-atom virials */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { if (system->pair_ptr->evflag) {
/* Acquire vectors */ /* Acquire vectors */
rvec_ScaledSum(delij, 1., system->my_atoms[i].x, rvec_ScaledSum(delij, 1., system->my_atoms[i].x,
@ -552,12 +537,11 @@ namespace ReaxFF {
eng_tmp = e_ang + e_pen + e_coa; eng_tmp = e_ang + e_pen + e_coa;
if (system->pair_ptr->evflag) if (system->pair_ptr->eflag_either)
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, j, j, system->N, 1, 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); eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
if (system->pair_ptr->vflag_atom) if (system->pair_ptr->vflag_either)
// NEED TO MAKE AN OMP VERSION OF THIS CALL! pair_reax_ptr->v_tally3_thr_proxy(i, j, k, fi_tmp, fk_tmp, delij, delkj, thr);
system->pair_ptr->v_tally3(i, j, k, fi_tmp, fk_tmp, delij, delkj);
} }
} // if (p_val1>0.001) } // if (p_val1>0.001)

3
src/REAXFF/pair_reaxc.cpp
View File

@ -464,9 +464,6 @@ void PairReaxC::compute(int eflag, int vflag)
evdwl = ecoul = 0.0; evdwl = ecoul = 0.0;
ev_init(eflag,vflag); 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; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system

129
src/REAXFF/reaxc_bond_orders.cpp
View File

@ -33,106 +33,7 @@
#include <cmath> #include <cmath>
namespace ReaxFF { namespace ReaxFF {
void Add_dBond_to_Forces_NPT(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;
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)
{ {
reax_list *bonds = (*lists) + BONDS; reax_list *bonds = (*lists) + BONDS;
bond_data *nbr_j, *nbr_k; bond_data *nbr_j, *nbr_k;
@ -180,10 +81,10 @@ namespace ReaxFF {
rvec_ScaledAdd(temp, coef.C3dbopi2, workspace->dDeltap_self[i]); rvec_ScaledAdd(temp, coef.C3dbopi2, workspace->dDeltap_self[i]);
rvec_Add(workspace->f[i], temp); rvec_Add(workspace->f[i], temp);
if (system->pair_ptr->vflag_atom) { if (system->pair_ptr->vflag_either) {
rvec_Scale(fi_tmp, -1.0, temp); rvec_Scale(fi_tmp, -0.5, temp);
rvec_ScaledSum(delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x); 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 // forces on j
@ -199,10 +100,10 @@ namespace ReaxFF {
rvec_ScaledAdd(temp, coef.C4dbopi2, workspace->dDeltap_self[j]); rvec_ScaledAdd(temp, coef.C4dbopi2, workspace->dDeltap_self[j]);
rvec_Add(workspace->f[j], temp); rvec_Add(workspace->f[j], temp);
if (system->pair_ptr->vflag_atom) { if (system->pair_ptr->vflag_either) {
rvec_Scale(fj_tmp, -1.0, temp); rvec_Scale(fj_tmp, -0.5, temp);
rvec_ScaledSum(delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x); 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 // forces on k: i neighbor
@ -216,12 +117,12 @@ namespace ReaxFF {
rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp); rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);
rvec_Add(workspace->f[k], temp); rvec_Add(workspace->f[k], temp);
if (system->pair_ptr->vflag_atom) { if (system->pair_ptr->vflag_either) {
rvec_Scale(fk_tmp, -1.0, temp); rvec_Scale(fk_tmp, -0.5, temp);
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); 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); 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_ScaledAdd(temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);
rvec_Add(workspace->f[k], temp); rvec_Add(workspace->f[k], temp);
if (system->pair_ptr->vflag_atom) { if (system->pair_ptr->vflag_either) {
rvec_Scale(fk_tmp, -1.0, temp); rvec_Scale(fk_tmp, -0.5, temp);
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x); 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); 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);
} }
} }
} }

11
src/REAXFF/reaxc_bonds.cpp
View File

@ -32,8 +32,7 @@
#include <cmath> #include <cmath>
namespace ReaxFF { namespace ReaxFF {
void Bonds(reax_system *system, simulation_data *data, void Bonds(reax_system *system, simulation_data *data, storage *workspace, reax_list **lists)
storage *workspace, reax_list **lists)
{ {
int i, j, pj, natoms; int i, j, pj, natoms;
int start_i, end_i; int start_i, end_i;
@ -94,8 +93,8 @@ namespace ReaxFF {
-twbp->De_p * bo_ij->BO_pi -twbp->De_p * bo_ij->BO_pi
-twbp->De_pp * bo_ij->BO_pi2; -twbp->De_pp * bo_ij->BO_pi2;
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) 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); 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 */ /* calculate derivatives of Bond Orders */
@ -124,8 +123,8 @@ namespace ReaxFF {
decobdboub = -gp10 * exphu * hulpov * decobdboub = -gp10 * exphu * hulpov *
(gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1)); (gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) 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); system->pair_ptr->ev_tally(i,j,natoms,1,estriph,0.0,0.0,0.0,0.0,0.0);
bo_ij->Cdbo += decobdbo; bo_ij->Cdbo += decobdbo;

1
src/REAXFF/reaxc_control.cpp
View File

@ -71,7 +71,6 @@ namespace ReaxFF {
/* assign default values */ /* assign default values */
control->nthreads = 1; control->nthreads = 1;
control->tabulate = 0; control->tabulate = 0;
control->virial = 0;
control->bond_cut = 5.0; control->bond_cut = 5.0;
control->bg_cut = 0.3; control->bg_cut = 0.3;
control->thb_cut = 0.001; control->thb_cut = 0.001;

14
src/REAXFF/reaxc_forces.cpp
View File

@ -47,7 +47,7 @@ namespace ReaxFF {
Valence_Angles(system, control, data, workspace, lists); Valence_Angles(system, control, data, workspace, lists);
Torsion_Angles(system, control, data, workspace, lists); Torsion_Angles(system, control, data, workspace, lists);
if (control->hbond_cut > 0) 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, static void Compute_NonBonded_Forces(reax_system *system,
@ -64,21 +64,15 @@ namespace ReaxFF {
Tabulated_vdW_Coulomb_Energy(system, control, data, workspace, lists); Tabulated_vdW_Coulomb_Energy(system, control, data, workspace, lists);
} }
static void Compute_Total_Force(reax_system *system, control_params *control, static void Compute_Total_Force(reax_system *system, storage *workspace, reax_list **lists)
storage *workspace, reax_list **lists)
{ {
int i, pj; int i, pj;
reax_list *bonds = (*lists) + BONDS; reax_list *bonds = (*lists) + BONDS;
for (i = 0; i < system->N; ++i) for (i = 0; i < system->N; ++i)
for (pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj) for (pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj)
if (i < bonds->select.bond_list[pj].nbr) { if (i < bonds->select.bond_list[pj].nbr)
if (control->virial == 0)
Add_dBond_to_Forces(system, i, pj, workspace, lists); 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, static void Validate_Lists(reax_system *system, reax_list **lists,
@ -387,6 +381,6 @@ namespace ReaxFF {
Compute_NonBonded_Forces(system, control, data, workspace, lists); Compute_NonBonded_Forces(system, control, data, workspace, lists);
/*********** total force ***************/ /*********** total force ***************/
Compute_Total_Force(system, control, workspace, lists); Compute_Total_Force(system, workspace, lists);
} }
} }

27
src/REAXFF/reaxc_hydrogen_bonds.cpp
View File

@ -32,9 +32,7 @@
#include "pair.h" #include "pair.h"
namespace ReaxFF { namespace ReaxFF {
void Hydrogen_Bonds(reax_system *system, control_params *control, void Hydrogen_Bonds(reax_system *system, simulation_data *data, storage *workspace, reax_list **lists)
simulation_data *data, storage *workspace,
reax_list **lists)
{ {
int i, j, k, pi, pk; int i, j, k, pi, pk;
int type_i, type_j, type_k; int type_i, type_j, type_k;
@ -42,11 +40,10 @@ namespace ReaxFF {
int hblist[MAX_BONDS]; int hblist[MAX_BONDS];
int itr, top; int itr, top;
int num_hb_intrs = 0; int num_hb_intrs = 0;
ivec rel_jk;
double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk; rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
rvec dvec_jk, force; rvec dvec_jk;
hbond_parameters *hbp; hbond_parameters *hbp;
bond_order_data *bo_ij; bond_order_data *bo_ij;
bond_data *pbond_ij; bond_data *pbond_ij;
@ -134,7 +131,6 @@ namespace ReaxFF {
/* hydrogen bond forces */ /* hydrogen bond forces */
bo_ij->Cdbo += CEhb1; // dbo term bo_ij->Cdbo += CEhb1; // dbo term
if (control->virial == 0) {
// dcos terms // dcos terms
rvec_ScaledAdd(workspace->f[i], +CEhb2, dcos_theta_di); rvec_ScaledAdd(workspace->f[i], +CEhb2, dcos_theta_di);
rvec_ScaledAdd(workspace->f[j], +CEhb2, dcos_theta_dj); rvec_ScaledAdd(workspace->f[j], +CEhb2, dcos_theta_dj);
@ -142,26 +138,9 @@ namespace ReaxFF {
// dr terms // dr terms
rvec_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk); rvec_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk);
rvec_ScaledAdd(workspace->f[k], +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);
}
/* tally into per-atom virials */ /* 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, rvec_ScaledSum(delij, 1., system->my_atoms[j].x,
-1., system->my_atoms[i].x); -1., system->my_atoms[i].x);
rvec_ScaledSum(delkj, 1., system->my_atoms[j].x, rvec_ScaledSum(delkj, 1., system->my_atoms[j].x,

12
src/REAXFF/reaxc_multi_body.cpp
View File

@ -93,8 +93,8 @@ namespace ReaxFF {
if (numbonds > 0 || control->enobondsflag) if (numbonds > 0 || control->enobondsflag)
workspace->CdDelta[i] += CElp; // lp - 1st term workspace->CdDelta[i] += CElp; // lp - 1st term
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) 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); 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 */ /* correction for C2 */
@ -119,8 +119,8 @@ namespace ReaxFF {
bo_ij->Cdbo += deahu2dbo; bo_ij->Cdbo += deahu2dbo;
workspace->CdDelta[i] += deahu2dsbo; workspace->CdDelta[i] += deahu2dsbo;
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) 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); 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]) * CEunder4 = CEunder1 * (dfvl*workspace->Delta_lp_temp[i]) *
p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2; p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) { if (system->pair_ptr->eflag_either) {
eng_tmp = e_ov; eng_tmp = e_ov;
if (numbonds > 0 || control->enobondsflag) if (numbonds > 0 || control->enobondsflag)
eng_tmp += e_un; eng_tmp += e_un;

24
src/REAXFF/reaxc_nonbonded.cpp
View File

@ -47,8 +47,8 @@ namespace ReaxFF {
(system->reax_param.sbp[type_i].eta / 2.) * SQR(q)); (system->reax_param.sbp[type_i].eta / 2.) * SQR(q));
data->my_en.e_pol += en_tmp; data->my_en.e_pol += en_tmp;
/* tally into per-atom energy */ /* tally energy into global or per-atom energy accumulators */
if (system->pair_ptr->evflag) 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); 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 dr3gamij_1, dr3gamij_3;
double e_ele, e_vdW, e_core, SMALL = 0.0001; double e_ele, e_vdW, e_core, SMALL = 0.0001;
double e_lg, de_lg, r_ij5, r_ij6, re6; double e_lg, de_lg, r_ij5, r_ij6, re6;
rvec temp;
two_body_parameters *twbp; two_body_parameters *twbp;
far_neighbor_data *nbr_pj; far_neighbor_data *nbr_pj;
reax_list *far_nbrs; reax_list *far_nbrs;
@ -193,7 +192,7 @@ namespace ReaxFF {
(dTap - Tap * r_ij / dr3gamij_1) / dr3gamij_3; (dTap - Tap * r_ij / dr3gamij_1) / dr3gamij_3;
/* tally into per-atom energy */ /* 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); pe_vdw = Tap * (e_vdW + e_core + e_lg);
rvec_ScaledSum(delij, 1., system->my_atoms[i].x, rvec_ScaledSum(delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x); -1., system->my_atoms[j].x);
@ -202,15 +201,8 @@ namespace ReaxFF {
f_tmp,delij[0],delij[1],delij[2]); 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[i], -(CEvd + CEclmb), nbr_pj->dvec);
rvec_ScaledAdd(workspace->f[j], +(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);
}
} }
} }
} }
@ -231,7 +223,6 @@ namespace ReaxFF {
double CEvd, CEclmb, SMALL = 0.0001; double CEvd, CEclmb, SMALL = 0.0001;
double f_tmp, delij[3]; double f_tmp, delij[3];
rvec temp;
far_neighbor_data *nbr_pj; far_neighbor_data *nbr_pj;
reax_list *far_nbrs; reax_list *far_nbrs;
LR_lookup_table *t; LR_lookup_table *t;
@ -301,7 +292,7 @@ namespace ReaxFF {
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q; CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
/* tally into per-atom energy */ /* 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, rvec_ScaledSum(delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x); -1., system->my_atoms[j].x);
f_tmp = -(CEvd + CEclmb); f_tmp = -(CEvd + CEclmb);
@ -309,15 +300,8 @@ namespace ReaxFF {
f_tmp,delij[0],delij[1],delij[2]); 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[i], -(CEvd + CEclmb), nbr_pj->dvec);
rvec_ScaledAdd(workspace->f[j], +(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);
}
} }
} }
} }

86
src/REAXFF/reaxc_torsion_angles.cpp
View File

@ -144,8 +144,6 @@ namespace ReaxFF {
double CEconj4, CEconj5, CEconj6; double CEconj4, CEconj5, CEconj6;
double e_tor, e_con; double e_tor, e_con;
rvec dvec_li; rvec dvec_li;
rvec force;
ivec rel_box_jl;
four_body_header *fbh; four_body_header *fbh;
four_body_parameters *fbp; four_body_parameters *fbp;
bond_data *pbond_ij, *pbond_jk, *pbond_kl; bond_data *pbond_ij, *pbond_jk, *pbond_kl;
@ -351,81 +349,29 @@ namespace ReaxFF {
bo_jk->Cdbo += (CEtors5 + CEconj2); bo_jk->Cdbo += (CEtors5 + CEconj2);
bo_kl->Cdbo += (CEtors6 + CEconj3); bo_kl->Cdbo += (CEtors6 + CEconj3);
if (control->virial == 0) {
/* dcos_theta_ijk */ /* dcos_theta_ijk */
rvec_ScaledAdd(workspace->f[i], rvec_ScaledAdd(workspace->f[i], CEtors7 + CEconj4, p_ijk->dcos_dk);
CEtors7 + CEconj4, p_ijk->dcos_dk); rvec_ScaledAdd(workspace->f[j], CEtors7 + CEconj4, p_ijk->dcos_dj);
rvec_ScaledAdd(workspace->f[j], rvec_ScaledAdd(workspace->f[k], CEtors7 + CEconj4, p_ijk->dcos_di);
CEtors7 + CEconj4, p_ijk->dcos_dj);
rvec_ScaledAdd(workspace->f[k],
CEtors7 + CEconj4, p_ijk->dcos_di);
/* dcos_theta_jkl */ /* dcos_theta_jkl */
rvec_ScaledAdd(workspace->f[j], rvec_ScaledAdd(workspace->f[j], CEtors8 + CEconj5, p_jkl->dcos_di);
CEtors8 + CEconj5, p_jkl->dcos_di); rvec_ScaledAdd(workspace->f[k], CEtors8 + CEconj5, p_jkl->dcos_dj);
rvec_ScaledAdd(workspace->f[k], rvec_ScaledAdd(workspace->f[l], CEtors8 + CEconj5, p_jkl->dcos_dk);
CEtors8 + CEconj5, p_jkl->dcos_dj);
rvec_ScaledAdd(workspace->f[l],
CEtors8 + CEconj5, p_jkl->dcos_dk);
/* dcos_omega */ /* dcos_omega */
rvec_ScaledAdd(workspace->f[i], rvec_ScaledAdd(workspace->f[i], CEtors9 + CEconj6, dcos_omega_di);
CEtors9 + CEconj6, dcos_omega_di); rvec_ScaledAdd(workspace->f[j], CEtors9 + CEconj6, dcos_omega_dj);
rvec_ScaledAdd(workspace->f[j], rvec_ScaledAdd(workspace->f[k], CEtors9 + CEconj6, dcos_omega_dk);
CEtors9 + CEconj6, dcos_omega_dj); rvec_ScaledAdd(workspace->f[l], CEtors9 + CEconj6, dcos_omega_dl);
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);
}
/* tally into per-atom virials */ /* tally into per-atom virials */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { if (system->pair_ptr->evflag) {
// acquire vectors // acquire vectors
rvec_ScaledSum(delil, 1., system->my_atoms[l].x, rvec_ScaledSum(delil, 1., system->my_atoms[l].x, -1., system->my_atoms[i].x);
-1., system->my_atoms[i].x); rvec_ScaledSum(deljl, 1., system->my_atoms[l].x, -1., system->my_atoms[j].x);
rvec_ScaledSum(deljl, 1., system->my_atoms[l].x, rvec_ScaledSum(delkl, 1., system->my_atoms[l].x, -1., system->my_atoms[k].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 // dcos_theta_ijk
rvec_Scale(fi_tmp, CEtors7 + CEconj4, p_ijk->dcos_dk); rvec_Scale(fi_tmp, CEtors7 + CEconj4, p_ijk->dcos_dk);
rvec_Scale(fj_tmp, CEtors7 + CEconj4, p_ijk->dcos_dj); rvec_Scale(fj_tmp, CEtors7 + CEconj4, p_ijk->dcos_dj);
@ -442,9 +388,9 @@ namespace ReaxFF {
// tally // tally
eng_tmp = e_tor + e_con; 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); 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); system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl);
} }
} // pl check ends } // pl check ends

38
src/REAXFF/reaxc_valence_angles.cpp
View File

@ -67,9 +67,8 @@ namespace ReaxFF {
} }
void Valence_Angles(reax_system *system, control_params *control, void Valence_Angles(reax_system *system, control_params *control, simulation_data *data,
simulation_data *data, storage *workspace, storage *workspace, reax_list **lists)
reax_list **lists)
{ {
int i, j, pi, k, pk, t; int i, j, pi, k, pk, t;
int type_i, type_j, type_k; int type_i, type_j, type_k;
@ -92,7 +91,6 @@ namespace ReaxFF {
double f7_ij, f7_jk, f8_Dj, f9_Dj; double f7_ij, f7_jk, f8_Dj, f9_Dj;
double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
double BOA_ij, BOA_jk; double BOA_ij, BOA_jk;
rvec force;
// Tallying variables // Tallying variables
double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
@ -114,7 +112,7 @@ namespace ReaxFF {
num_thb_intrs = 0; 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; type_j = system->my_atoms[j].type;
if (type_j < 0) continue; if (type_j < 0) continue;
start_j = Start_Index(j, bonds); start_j = Start_Index(j, bonds);
@ -218,7 +216,6 @@ namespace ReaxFF {
++num_thb_intrs; ++num_thb_intrs;
if ((j < system->n) && (BOA_jk > 0.0) && if ((j < system->n) && (BOA_jk > 0.0) &&
(bo_ij->BO > control->thb_cut) && (bo_ij->BO > control->thb_cut) &&
(bo_jk->BO > control->thb_cut) && (bo_jk->BO > control->thb_cut) &&
@ -347,38 +344,25 @@ namespace ReaxFF {
bo_jt->Cdbopi2 += CEval5; bo_jt->Cdbopi2 += CEval5;
} }
if (control->virial == 0) {
rvec_ScaledAdd(workspace->f[i], CEval8, p_ijk->dcos_di); rvec_ScaledAdd(workspace->f[i], CEval8, p_ijk->dcos_di);
rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj); rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj);
rvec_ScaledAdd(workspace->f[k], CEval8, p_ijk->dcos_dk); 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[j], CEval8, p_ijk->dcos_dj); /* tally energy */
if (system->pair_ptr->eflag_either) {
rvec_Scale(force, CEval8, p_ijk->dcos_dk); eng_tmp = e_ang + e_pen + e_coa;
rvec_Add(workspace->f[k], force); 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 */ /* tally virial */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) { if (system->pair_ptr->vflag_either) {
/* Acquire vectors */ /* Acquire vectors */
rvec_ScaledSum(delij, 1., system->my_atoms[i].x, rvec_ScaledSum(delij, 1., system->my_atoms[i].x, -1., system->my_atoms[j].x);
-1., system->my_atoms[j].x); rvec_ScaledSum(delkj, 1., system->my_atoms[k].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(fi_tmp, -CEval8, p_ijk->dcos_di);
rvec_Scale(fj_tmp, -CEval8, p_ijk->dcos_dj); rvec_Scale(fj_tmp, -CEval8, p_ijk->dcos_dj);
rvec_Scale(fk_tmp, -CEval8, p_ijk->dcos_dk); 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);
} }
} }

5
src/REAXFF/reaxff_api.h
View File

@ -54,7 +54,6 @@ namespace ReaxFF
single_body_parameters *, single_body_parameters *, single_body_parameters *, single_body_parameters *,
two_body_parameters *); two_body_parameters *);
extern void Add_dBond_to_Forces(reax_system*, int, int, storage*, reax_list**); 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 // bonds
@ -78,8 +77,8 @@ namespace ReaxFF
// hydrogen bonds // hydrogen bonds
extern void Hydrogen_Bonds(reax_system *, control_params *, extern void Hydrogen_Bonds(reax_system *, simulation_data *, storage *, reax_list **);
simulation_data *, storage *, reax_list **);
// init md // init md
extern void Init_System(reax_system *, control_params *); extern void Init_System(reax_system *, control_params *);

1
src/REAXFF/reaxff_types.h
View File

@ -240,7 +240,6 @@ namespace ReaxFF
double thb_cutsq; double thb_cutsq;
int tabulate; int tabulate;
int virial;
int lgflag; int lgflag;
int enobondsflag; int enobondsflag;

33
src/USER-MISC/pair_drip.cpp
View File

@ -366,10 +366,7 @@ void PairDRIP::compute(int eflag, int vflag)
f[j][0] += fj[0]; f[j][0] += fj[0];
f[j][1] += fj[1]; f[j][1] += fj[1];
f[j][2] += fj[2]; f[j][2] += fj[2];
if (vflag_either) v_tally2_newton(j, fj, x[j]);
// 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]);
} }
} //loop over jj } //loop over jj
@ -377,10 +374,7 @@ void PairDRIP::compute(int eflag, int vflag)
f[i][0] += fi[0]; f[i][0] += fi[0];
f[i][1] += fi[1]; f[i][1] += fi[1];
f[i][2] += fi[2]; f[i][2] += fi[2];
if (vflag_either) v_tally2_newton(i, fi, x[i]);
// 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]);
} // loop over ii } // loop over ii
@ -529,22 +523,13 @@ double PairDRIP::calc_repulsive(int const i, int const j, Param& p,
f[nbj3][k] += fnbj3[k]; f[nbj3][k] += fnbj3[k];
} }
if (vflag_atom) { if (vflag_either) {
// multiply since v_tally has a 0.5 coeff v_tally2_newton(nbi1, fnbi1, x[nbi1]);
for (int k = 0; k < DIM; k++) { v_tally2_newton(nbi2, fnbi2, x[nbi2]);
fnbi1[k] *= 2; v_tally2_newton(nbi3, fnbi3, x[nbi3]);
fnbi2[k] *= 2; v_tally2_newton(nbj1, fnbj1, x[nbj1]);
fnbi3[k] *= 2; v_tally2_newton(nbj2, fnbj2, x[nbj2]);
fnbj1[k] *= 2; v_tally2_newton(nbj3, fnbj3, x[nbj3]);
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]);
} }
return phi; return phi;

54
src/pair.cpp
View File

@ -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 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 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 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 for virial, have delx,dely,delz and fx,fy,fz
called when using full neighbor lists 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 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 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 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 tally virial into global or per-atom accumulators
called by REAX/C potential, newton_pair is always on called by ReaxFF potential, newton_pair is always on
fi is magnitude of force on atom i 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]; double v[6];
v[0] = 0.5*deli[0]*fi[0]; v[0] = deli[0]*fi[0];
v[1] = 0.5*deli[1]*fi[1]; v[1] = deli[1]*fi[1];
v[2] = 0.5*deli[2]*fi[2]; v[2] = deli[2]*fi[2];
v[3] = 0.5*deli[0]*fi[1]; v[3] = deli[0]*fi[1];
v[4] = 0.5*deli[0]*fi[2]; v[4] = deli[0]*fi[2];
v[5] = 0.5*deli[1]*fi[2]; v[5] = deli[1]*fi[2];
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][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]; 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 called by AIREBO potential, newton_pair is always on
fpair is magnitude of force on atom I 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 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) 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 tally virial into global or per-atom accumulators
called by AIREBO potential, newton_pair is always on called by AIREBO potential, Tersoff, ReaxFF potentials, newton_pair is always on
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Pair::v_tally4(int i, int j, int k, int m, 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 called by pair lubricate potential with 6 tensor components
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

2
src/pair.h
View File

@ -137,9 +137,9 @@ class Pair : protected Pointers {
// need to be public, so can be called by pair_style reaxc // 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_tally(int, int, int, int, double, double, double, double, double, double);
void ev_tally3(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_tally3(int, int, int, double *, double *, double *, double *);
void v_tally4(int, int, int, int, double *, double *, 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, void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double,