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

Fix issues in ReaxFF QEq and ACKS2

Browse Source
This commit is contained in:
Stan Gerald Moore
2022-04-04 10:24:54 -06:00
parent db00b49a50
commit b4fc86e467
9 changed files with 51 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
src/KOKKOS/fix_acks2_reaxff_kokkos.cpp
View File

@ -219,7 +219,7 @@ void FixACKS2ReaxFFKokkos<DeviceType>::pre_force(int vflag)
d_ilist = k_list->d_ilist;
nn = list->inum;
NN = list->inum + list->gnum;
NN = atom->nlocal + atom->nghost;
copymode = 1;
@ -526,7 +526,7 @@ void FixACKS2ReaxFFKokkos<DeviceType>::allocate_array()
if (efield) get_chi_field();
// init_storage
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2Zero>(0,NN),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2Zero>(0,nn),*this);
}
@ -1378,18 +1378,15 @@ template<class DeviceType>
void FixACKS2ReaxFFKokkos<DeviceType>::calculate_Q()
{
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2CalculateQ1>(0,nn),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2CalculateQ>(0,nn),*this);
pack_flag = 2;
//comm->forward_comm( this ); //Dist_vector( s );
k_s.modify<DeviceType>();
k_s.sync<LMPHostType>();
pack_flag = 4;
//comm->forward_comm( this ); //Dist_vector( q );
atomKK->k_q.modify<DeviceType>();
atomKK->k_q.sync<LMPHostType>();
comm->forward_comm(this);
k_s.modify<LMPHostType>();
k_s.sync<DeviceType>();
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2CalculateQ2>(0,NN),*this);
atomKK->k_q.modify<LMPHostType>();
atomKK->k_q.sync<DeviceType>();
}
/* ---------------------------------------------------------------------- */
@ -1822,11 +1819,13 @@ void FixACKS2ReaxFFKokkos<DeviceType>::operator() (TagACKS2Norm3, const int &ii,
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixACKS2ReaxFFKokkos<DeviceType>::operator() (TagACKS2CalculateQ1, const int &ii) const
void FixACKS2ReaxFFKokkos<DeviceType>::operator() (TagACKS2CalculateQ, const int &ii) const
{
const int i = d_ilist[ii];
if (mask[i] & groupbit) {
q(i) = d_s(i);
/* backup s */
for (int k = nprev-1; k > 0; --k) {
d_s_hist(i,k) = d_s_hist(i,k-1);
@ -1848,17 +1847,6 @@ void FixACKS2ReaxFFKokkos<DeviceType>::operator() (TagACKS2CalculateQ1, const in
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixACKS2ReaxFFKokkos<DeviceType>::operator() (TagACKS2CalculateQ2, const int &ii) const
{
const int i = d_ilist[ii];
if (mask[i] & groupbit)
q(i) = d_s(i);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixACKS2ReaxFFKokkos<DeviceType>::cleanup_copy()
{

8
src/KOKKOS/fix_acks2_reaxff_kokkos.h
View File

@ -54,8 +54,7 @@ struct TagACKS2Precon1B{};
struct TagACKS2Precon2{};
struct TagACKS2Add{};
struct TagACKS2ZeroQGhosts{};
struct TagACKS2CalculateQ1{};
struct TagACKS2CalculateQ2{};
struct TagACKS2CalculateQ{};
template<class DeviceType>
class FixACKS2ReaxFFKokkos : public FixACKS2ReaxFF {
@ -152,10 +151,7 @@ class FixACKS2ReaxFFKokkos : public FixACKS2ReaxFF {
void operator()(TagACKS2ZeroQGhosts, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagACKS2CalculateQ1, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagACKS2CalculateQ2, const int&) const;
void operator()(TagACKS2CalculateQ, const int&) const;
KOKKOS_INLINE_FUNCTION
double calculate_H_k(const F_FLOAT &r, const F_FLOAT &shld) const;

5
src/KOKKOS/fix_qeq_reaxff_kokkos.cpp
View File

@ -208,7 +208,6 @@ void FixQEqReaxFFKokkos<DeviceType>::pre_force(int /*vflag*/)
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
nn = list->inum;
NN = list->inum + list->gnum;
copymode = 1;
@ -375,7 +374,7 @@ void FixQEqReaxFFKokkos<DeviceType>::allocate_array()
if (efield) get_chi_field();
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagQEqZero>(0,NN),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagQEqZero>(0,nn),*this);
}
/* ---------------------------------------------------------------------- */
@ -872,7 +871,7 @@ void FixQEqReaxFFKokkos<DeviceType>::sparse_matvec_kokkos(typename AT::t_ffloat2
}
if (neighflag != FULL) {
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagQEqZeroQGhosts>(nn,NN),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagQEqZeroQGhosts>(atom->nlocal,atom->nghost),*this);
if (need_dup)
dup_o.reset_except(d_o);

2
src/KOKKOS/pair_reaxff_kokkos.cpp
View File

@ -715,7 +715,7 @@ void PairReaxFFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
newton_pair = force->newton_pair;
nn = list->inum;
NN = list->inum + list->gnum;
NN = atom->nlocal + atom->nghost;
const int inum = list->inum;
const int ignum = inum + list->gnum;

25
src/OPENMP/fix_qeq_reaxff_omp.cpp
View File

@ -237,7 +237,7 @@ void FixQEqReaxFFOMP::init_storage()
#if defined(_OPENMP)
#pragma omp parallel for schedule(static)
#endif
for (int i = 0; i < NN; i++) {
for (int i = 0; i < nn; i++) {
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
if (efield) b_s[i] -= chi_field[i];
@ -258,13 +258,11 @@ void FixQEqReaxFFOMP::pre_force(int /* vflag */)
if (reaxff) {
nn = reaxff->list->inum;
NN = reaxff->list->inum + reaxff->list->gnum;
ilist = reaxff->list->ilist;
numneigh = reaxff->list->numneigh;
firstneigh = reaxff->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -505,15 +503,14 @@ void FixQEqReaxFFOMP::sparse_matvec(sparse_matrix *A, double *x, double *b)
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (ii = nn; ii < NN; ++ii) {
i = ilist[ii];
for (i = atom->nlocal; i < atom->nghost; ++i) {
if (atom->mask[i] & groupbit) b[i] = 0;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i)
for (i = 0; i < atom->nghost; ++i)
for (int t=0; t<nthreads; t++) b_temp[t][i] = 0.0;
// Wait for b accumulated and b_temp zeroed.
@ -538,7 +535,7 @@ void FixQEqReaxFFOMP::sparse_matvec(sparse_matrix *A, double *x, double *b)
#pragma omp barrier
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i)
for (i = 0; i < atom->nghost; ++i)
for (int t = 0; t < nthreads; ++t) b[i] += b_temp[t][i];
} //end omp parallel
@ -838,8 +835,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x1, double *x
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (ii = nn; ii < NN; ++ii) {
i = ilist[ii];
for (i = atom->nlocal; i < atom->nghost; ++i) {
if (atom->mask[i] & groupbit) {
indxI = 2 * i;
b[indxI] = 0;
@ -850,7 +846,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x1, double *x
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i) {
for (i = 0; i < atom->nghost; ++i) {
indxI = 2 * i;
for (int t=0; t<nthreads; t++) {
b_temp[t][indxI] = 0.0;
@ -884,7 +880,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x1, double *x
#pragma omp barrier
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i) {
for (i = 0; i < atom->nghost; ++i) {
indxI = 2 * i;
for (int t = 0; t < nthreads; ++t) {
b[indxI] += b_temp[t][indxI];
@ -929,8 +925,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x, double *b)
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (ii = nn; ii < NN; ++ii) {
i = ilist[ii];
for (i = atom->nlocal; i < atom->nghost; ++i) {
if (atom->mask[i] & groupbit) {
indxI = 2 * i;
b[indxI] = 0;
@ -941,7 +936,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x, double *b)
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i) {
for (i = 0; i < atom->nghost; ++i) {
indxI = 2 * i;
for (int t=0; t<nthreads; t++) {
b_temp[t][indxI] = 0.0;
@ -975,7 +970,7 @@ void FixQEqReaxFFOMP::dual_sparse_matvec(sparse_matrix *A, double *x, double *b)
#pragma omp barrier
#pragma omp for schedule(dynamic,50)
#endif
for (i = 0; i < NN; ++i) {
for (i = 0; i < atom->nghost; ++i) {
indxI = 2 * i;
for (int t = 0; t < nthreads; ++t) {
b[indxI] += b_temp[t][indxI];

32
src/REAXFF/fix_acks2_reaxff.cpp
View File

@ -306,7 +306,7 @@ void FixACKS2ReaxFF::init_storage()
{
if (efield) get_chi_field();
for (int ii = 0; ii < NN; ii++) {
for (int ii = 0; ii < nn; ii++) {
int i = ilist[ii];
if (atom->mask[i] & groupbit) {
b_s[i] = -chi[atom->type[i]];
@ -329,17 +329,15 @@ void FixACKS2ReaxFF::pre_force(int /*vflag*/)
{
if (update->ntimestep % nevery) return;
int n = atom->nlocal;
NN = atom->nlocal + atom->nghost;
if (reaxff) {
nn = reaxff->list->inum;
NN = reaxff->list->inum + reaxff->list->gnum;
ilist = reaxff->list->ilist;
numneigh = reaxff->list->numneigh;
firstneigh = reaxff->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -349,7 +347,7 @@ void FixACKS2ReaxFF::pre_force(int /*vflag*/)
// need to be atom->nmax in length
if (atom->nmax > nmax) reallocate_storage();
if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
if (atom->nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
reallocate_matrix();
if (efield) get_chi_field();
@ -626,8 +624,7 @@ void FixACKS2ReaxFF::sparse_matvec_acks2(sparse_matrix *H, sparse_matrix *X, dou
}
}
for (ii = nn; ii < NN; ++ii) {
i = ilist[ii];
for (i = atom->nlocal; i < atom->nghost; ++i) {
if (atom->mask[i] & groupbit) {
b[i] = 0;
b[NN + i] = 0;
@ -680,6 +677,8 @@ void FixACKS2ReaxFF::calculate_Q()
i = ilist[ii];
if (atom->mask[i] & groupbit) {
atom->q[i] = s[i];
/* backup s */
for (k = nprev-1; k > 0; --k) {
s_hist[i][k] = s_hist[i][k-1];
@ -698,14 +697,8 @@ void FixACKS2ReaxFF::calculate_Q()
}
}
pack_flag = 2;
comm->forward_comm(this); //Dist_vector(s);
for (int ii = 0; ii < NN; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
atom->q[i] = s[i];
}
pack_flag = 4;
comm->forward_comm(this); //Dist_vector(q);
}
/* ---------------------------------------------------------------------- */
@ -733,6 +726,11 @@ int FixACKS2ReaxFF::pack_forward_comm(int n, int *list, double *buf,
buf[m++] = q_hat[j];
buf[m++] = q_hat[NN+j];
}
} else if (pack_flag == 4) {
for(int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = atom->q[j];
}
}
return m;
}
@ -761,6 +759,10 @@ void FixACKS2ReaxFF::unpack_forward_comm(int n, int first, double *buf)
q_hat[i] = buf[m++];
q_hat[NN+i] = buf[m++];
}
} else if (pack_flag == 4) {
for(i = first; i < last; i++) {
atom->q[i] = buf[m++];
}
}
}

2
src/REAXFF/fix_acks2_reaxff.h
View File

@ -37,7 +37,7 @@ class FixACKS2ReaxFF : public FixQEqReaxFF {
double *get_s() { return s; }
protected:
int last_rows_rank, last_rows_flag;
int NN, last_rows_rank, last_rows_flag;
double **s_hist_X, **s_hist_last;
double *bcut_acks2, bond_softness, **bcut; // acks2 parameters

18
src/REAXFF/fix_qeq_reaxff.cpp
View File

@ -104,7 +104,7 @@ FixQEqReaxFF::FixQEqReaxFF(LAMMPS *lmp, int narg, char **arg) :
shld = nullptr;
nn = n_cap = 0;
NN = nmax = 0;
nmax = 0;
m_fill = m_cap = 0;
pack_flag = 0;
s = nullptr;
@ -319,7 +319,7 @@ void FixQEqReaxFF::reallocate_storage()
void FixQEqReaxFF::allocate_matrix()
{
int i,ii,n,m;
int i,ii,m;
int mincap;
double safezone;
@ -332,8 +332,7 @@ void FixQEqReaxFF::allocate_matrix()
safezone = REAX_SAFE_ZONE;
}
n = atom->nlocal;
n_cap = MAX((int)(n * safezone), mincap);
n_cap = MAX((int)(atom->nlocal * safezone), mincap);
// determine the total space for the H matrix
@ -493,13 +492,11 @@ void FixQEqReaxFF::setup_pre_force(int vflag)
{
if (reaxff) {
nn = reaxff->list->inum;
NN = reaxff->list->inum + reaxff->list->gnum;
ilist = reaxff->list->ilist;
numneigh = reaxff->list->numneigh;
firstneigh = reaxff->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -537,7 +534,7 @@ void FixQEqReaxFF::init_storage()
{
if (efield) get_chi_field();
for (int ii = 0; ii < NN; ii++) {
for (int ii = 0; ii < nn; ii++) {
int i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[atom->type[i]];
@ -561,13 +558,11 @@ void FixQEqReaxFF::pre_force(int /*vflag*/)
if (reaxff) {
nn = reaxff->list->inum;
NN = reaxff->list->inum + reaxff->list->gnum;
ilist = reaxff->list->ilist;
numneigh = reaxff->list->numneigh;
firstneigh = reaxff->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -791,11 +786,8 @@ void FixQEqReaxFF::sparse_matvec(sparse_matrix *A, double *x, double *b)
b[i] = eta[atom->type[i]] * x[i];
}
for (ii = nn; ii < NN; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
for (i = atom->nlocal; i < atom->nghost; ++i)
b[i] = 0;
}
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];

2
src/REAXFF/fix_qeq_reaxff.h
View File

@ -59,7 +59,7 @@ class FixQEqReaxFF : public Fix {
protected:
int nevery, reaxflag;
int matvecs;
int nn, NN, m_fill;
int nn, m_fill;
int n_cap, nmax, m_cap;
int pack_flag;
int nlevels_respa;