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

more whitespace

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-04-20 20:14:30 -04:00
parent daf984835c
commit 2a7e662d71
28 changed files with 532 additions and 532 deletions
Download Patch File Download Diff File Expand all files Collapse all files

90
src/KOKKOS/fix_qeq_reax_kokkos.cpp
View File

@ -140,7 +140,7 @@ void FixQEqReaxKokkos<DeviceType>::init_shielding_k()
for (i = 1; i <= ntypes; ++i)
for (j = 1; j <= ntypes; ++j)
k_shield.h_view(i,j) = pow( gamma[i] * gamma[j], -1.5 );
k_shield.h_view(i,j) = pow(gamma[i] * gamma[j], -1.5);
k_shield.template modify<LMPHostType>();
k_shield.template sync<DeviceType>();
@ -260,13 +260,13 @@ void FixQEqReaxKokkos<DeviceType>::pre_force(int /*vflag*/)
FixQEqReaxKokkosMatVecFunctor<DeviceType> matvec_functor(this);
Kokkos::parallel_for(inum,matvec_functor);
// comm->forward_comm_fix(this); //Dist_vector( s );
// comm->forward_comm_fix(this); //Dist_vector(s);
pack_flag = 2;
k_s.template modify<DeviceType>();
comm->forward_comm_fix(this);
k_s.template sync<DeviceType>();
// comm->forward_comm_fix(this); //Dist_vector( t );
// comm->forward_comm_fix(this); //Dist_vector(t);
pack_flag = 3;
k_t.template modify<DeviceType>();
comm->forward_comm_fix(this);
@ -731,7 +731,7 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve1()
if (execution_space == Host) teamsize = 1;
else teamsize = 128;
// sparse_matvec( &H, x, q );
// sparse_matvec(&H, x, q);
FixQEqReaxKokkosSparse12Functor<DeviceType> sparse12_functor(this);
Kokkos::parallel_for(inum,sparse12_functor);
if (neighflag != FULL) {
@ -751,38 +751,38 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve1()
if (neighflag != FULL) {
k_o.template modify<DeviceType>();
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->reverse_comm_fix(this); //Coll_vector(q);
k_o.template sync<DeviceType>();
}
// vector_sum( r , 1., b, -1., q, nn );
// vector_sum(r , 1., b, -1., q, nn);
// preconditioning: d[j] = r[j] * Hdia_inv[j];
// b_norm = parallel_norm( b, nn );
// b_norm = parallel_norm(b, nn);
F_FLOAT my_norm = 0.0;
FixQEqReaxKokkosNorm1Functor<DeviceType> norm1_functor(this);
Kokkos::parallel_reduce(inum,norm1_functor,my_norm);
F_FLOAT norm_sqr = 0.0;
MPI_Allreduce( &my_norm, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_norm, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
b_norm = sqrt(norm_sqr);
// sig_new = parallel_dot( r, d, nn);
// sig_new = parallel_dot(r, d, nn);
F_FLOAT my_dot = 0.0;
FixQEqReaxKokkosDot1Functor<DeviceType> dot1_functor(this);
Kokkos::parallel_reduce(inum,dot1_functor,my_dot);
F_FLOAT dot_sqr = 0.0;
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
F_FLOAT sig_new = dot_sqr;
int loop;
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
// comm->forward_comm_fix(this); //Dist_vector(d);
pack_flag = 1;
k_d.template modify<DeviceType>();
comm->forward_comm_fix(this);
k_d.template sync<DeviceType>();
// sparse_matvec( &H, d, q );
// sparse_matvec(&H, d, q);
FixQEqReaxKokkosSparse22Functor<DeviceType> sparse22_functor(this);
Kokkos::parallel_for(inum,sparse22_functor);
if (neighflag != FULL) {
@ -804,36 +804,36 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve1()
if (neighflag != FULL) {
k_o.template modify<DeviceType>();
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->reverse_comm_fix(this); //Coll_vector(q);
k_o.template sync<DeviceType>();
}
// tmp = parallel_dot( d, q, nn);
// tmp = parallel_dot(d, q, nn);
my_dot = dot_sqr = 0.0;
FixQEqReaxKokkosDot2Functor<DeviceType> dot2_functor(this);
Kokkos::parallel_reduce(inum,dot2_functor,my_dot);
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
tmp = dot_sqr;
alpha = sig_new / tmp;
sig_old = sig_new;
// vector_add( s, alpha, d, nn );
// vector_add( r, -alpha, q, nn );
// vector_add(s, alpha, d, nn);
// vector_add(r, -alpha, q, nn);
my_dot = dot_sqr = 0.0;
FixQEqReaxKokkosPrecon1Functor<DeviceType> precon1_functor(this);
Kokkos::parallel_for(inum,precon1_functor);
// preconditioning: p[j] = r[j] * Hdia_inv[j];
// sig_new = parallel_dot( r, p, nn);
// sig_new = parallel_dot(r, p, nn);
FixQEqReaxKokkosPreconFunctor<DeviceType> precon_functor(this);
Kokkos::parallel_reduce(inum,precon_functor,my_dot);
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
sig_new = dot_sqr;
beta = sig_new / sig_old;
// vector_sum( d, 1., p, beta, d, nn );
// vector_sum(d, 1., p, beta, d, nn);
FixQEqReaxKokkosVecSum2Functor<DeviceType> vecsum2_functor(this);
Kokkos::parallel_for(inum,vecsum2_functor);
}
@ -859,7 +859,7 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve2()
if (execution_space == Host) teamsize = 1;
else teamsize = 128;
// sparse_matvec( &H, x, q );
// sparse_matvec(&H, x, q);
FixQEqReaxKokkosSparse32Functor<DeviceType> sparse32_functor(this);
Kokkos::parallel_for(inum,sparse32_functor);
if (neighflag != FULL) {
@ -881,38 +881,38 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve2()
if (neighflag != FULL) {
k_o.template modify<DeviceType>();
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->reverse_comm_fix(this); //Coll_vector(q);
k_o.template sync<DeviceType>();
}
// vector_sum( r , 1., b, -1., q, nn );
// vector_sum(r , 1., b, -1., q, nn);
// preconditioning: d[j] = r[j] * Hdia_inv[j];
// b_norm = parallel_norm( b, nn );
// b_norm = parallel_norm(b, nn);
F_FLOAT my_norm = 0.0;
FixQEqReaxKokkosNorm2Functor<DeviceType> norm2_functor(this);
Kokkos::parallel_reduce(inum,norm2_functor,my_norm);
F_FLOAT norm_sqr = 0.0;
MPI_Allreduce( &my_norm, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_norm, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
b_norm = sqrt(norm_sqr);
// sig_new = parallel_dot( r, d, nn);
// sig_new = parallel_dot(r, d, nn);
F_FLOAT my_dot = 0.0;
FixQEqReaxKokkosDot1Functor<DeviceType> dot1_functor(this);
Kokkos::parallel_reduce(inum,dot1_functor,my_dot);
F_FLOAT dot_sqr = 0.0;
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
F_FLOAT sig_new = dot_sqr;
int loop;
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
// comm->forward_comm_fix(this); //Dist_vector(d);
pack_flag = 1;
k_d.template modify<DeviceType>();
comm->forward_comm_fix(this);
k_d.template sync<DeviceType>();
// sparse_matvec( &H, d, q );
// sparse_matvec(&H, d, q);
FixQEqReaxKokkosSparse22Functor<DeviceType> sparse22_functor(this);
Kokkos::parallel_for(inum,sparse22_functor);
if (neighflag != FULL) {
@ -934,36 +934,36 @@ int FixQEqReaxKokkos<DeviceType>::cg_solve2()
if (neighflag != FULL) {
k_o.template modify<DeviceType>();
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->reverse_comm_fix(this); //Coll_vector(q);
k_o.template sync<DeviceType>();
}
// tmp = parallel_dot( d, q, nn);
// tmp = parallel_dot(d, q, nn);
my_dot = dot_sqr = 0.0;
FixQEqReaxKokkosDot2Functor<DeviceType> dot2_functor(this);
Kokkos::parallel_reduce(inum,dot2_functor,my_dot);
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
tmp = dot_sqr;
alpha = sig_new / tmp;
sig_old = sig_new;
// vector_add( t, alpha, d, nn );
// vector_add( r, -alpha, q, nn );
// vector_add(t, alpha, d, nn);
// vector_add(r, -alpha, q, nn);
my_dot = dot_sqr = 0.0;
FixQEqReaxKokkosPrecon2Functor<DeviceType> precon2_functor(this);
Kokkos::parallel_for(inum,precon2_functor);
// preconditioning: p[j] = r[j] * Hdia_inv[j];
// sig_new = parallel_dot( r, p, nn);
// sig_new = parallel_dot(r, p, nn);
FixQEqReaxKokkosPreconFunctor<DeviceType> precon_functor(this);
Kokkos::parallel_reduce(inum,precon_functor,my_dot);
MPI_Allreduce( &my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&my_dot, &dot_sqr, 1, MPI_DOUBLE, MPI_SUM, world);
sig_new = dot_sqr;
beta = sig_new / sig_old;
// vector_sum( d, 1., p, beta, d, nn );
// vector_sum(d, 1., p, beta, d, nn);
FixQEqReaxKokkosVecSum2Functor<DeviceType> vecsum2_functor(this);
Kokkos::parallel_for(inum,vecsum2_functor);
}
@ -984,18 +984,18 @@ void FixQEqReaxKokkos<DeviceType>::calculate_q()
F_FLOAT sum, sum_all;
const int inum = list->inum;
// s_sum = parallel_vector_acc( s, nn );
// s_sum = parallel_vector_acc(s, nn);
sum = sum_all = 0.0;
FixQEqReaxKokkosVecAcc1Functor<DeviceType> vecacc1_functor(this);
Kokkos::parallel_reduce(inum,vecacc1_functor,sum);
MPI_Allreduce(&sum, &sum_all, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&sum, &sum_all, 1, MPI_DOUBLE, MPI_SUM, world);
const F_FLOAT s_sum = sum_all;
// t_sum = parallel_vector_acc( t, nn);
// t_sum = parallel_vector_acc(t, nn);
sum = sum_all = 0.0;
FixQEqReaxKokkosVecAcc2Functor<DeviceType> vecacc2_functor(this);
Kokkos::parallel_reduce(inum,vecacc2_functor,sum);
MPI_Allreduce(&sum, &sum_all, 1, MPI_DOUBLE, MPI_SUM, world );
MPI_Allreduce(&sum, &sum_all, 1, MPI_DOUBLE, MPI_SUM, world);
const F_FLOAT t_sum = sum_all;
// u = s_sum / t_sum;
@ -1007,7 +1007,7 @@ void FixQEqReaxKokkos<DeviceType>::calculate_q()
atomKK->modified(execution_space,Q_MASK);
pack_flag = 4;
//comm->forward_comm_fix( this ); //Dist_vector( atom->q );
//comm->forward_comm_fix(this); //Dist_vector(atom->q);
comm->forward_comm_fix(this);
}
@ -1382,11 +1382,11 @@ KOKKOS_INLINE_FUNCTION
void FixQEqReaxKokkos<DeviceType>::operator()(TagFixQEqReaxUnpackForwardComm, const int &i) const {
if (pack_flag == 1)
d_d[i + first] = d_buf[i];
else if ( pack_flag == 2)
else if (pack_flag == 2)
d_s[i + first] = d_buf[i];
else if ( pack_flag == 3)
else if (pack_flag == 3)
d_t[i + first] = d_buf[i];
else if ( pack_flag == 4)
else if (pack_flag == 4)
q[i + first] = d_buf[i];
}

2
src/KOKKOS/pair_reaxc_kokkos.cpp
View File

@ -2389,7 +2389,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
exp_pen2ij = exp(-p_pen2 * (BOA_ij - 2.0)*(BOA_ij - 2.0));
exp_pen2jk = exp(-p_pen2 * (BOA_ik - 2.0)*(BOA_ik - 2.0));
exp_pen3 = exp(-p_pen3 * d_Delta[i]);
exp_pen4 = exp( p_pen4 * d_Delta[i]);
exp_pen4 = exp(p_pen4 * d_Delta[i]);
trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
f9_Dj = (2.0 + exp_pen3) / trm_pen34;
Cf9j = (-p_pen3 * exp_pen3 * trm_pen34 - (2.0 + exp_pen3) *

2
src/KOKKOS/pair_reaxc_kokkos.h
View File

@ -339,7 +339,7 @@ class PairReaxCKokkos : public PairReaxC {
void init_md();
int Init_Lookup_Tables();
void Deallocate_Lookup_Tables();
void LR_vdW_Coulomb( int i, int j, double r_ij, ReaxFF::LR_data *lr );
void LR_vdW_Coulomb(int i, int j, double r_ij, ReaxFF::LR_data *lr);
typedef Kokkos::DualView<int*,DeviceType> tdual_int_1d;
Kokkos::DualView<params_sing*,typename DeviceType::array_layout,DeviceType> k_params_sing;

2
src/QEQ/fix_qeq.cpp
View File

@ -434,7 +434,7 @@ void FixQEq::sparse_matvec(sparse_matrix *A, double *x, double *b)
for (i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit)
b[i] = eta[ atom->type[i] ] * x[i];
b[i] = eta[atom->type[i]] * x[i];
}
for (i = nlocal; i < nall; ++i) {

2
src/QEQ/fix_qeq_dynamic.cpp
View File

@ -266,7 +266,7 @@ void FixQEqDynamic::unpack_forward_comm(int n, int first, double *buf)
if (pack_flag == 1)
for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
else if ( pack_flag == 2)
else if (pack_flag == 2)
for (m = 0, i = first; m < n; m++, i++) qf[i] = buf[m];
}

2
src/QEQ/fix_qeq_fire.cpp
View File

@ -326,7 +326,7 @@ void FixQEqFire::unpack_forward_comm(int n, int first, double *buf)
if (pack_flag == 1)
for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
else if ( pack_flag == 2)
else if (pack_flag == 2)
for (m = 0, i = first; m < n; m++, i++) qf[i] = buf[m];
}

10
src/QEQ/fix_qeq_point.cpp
View File

@ -109,18 +109,18 @@ void FixQEqPoint::init_matvec()
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -(chi[atom->type[i]] + chizj[i]);
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
comm->forward_comm_fix(this); //Dist_vector(s);
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
comm->forward_comm_fix(this); //Dist_vector(t);
}
/* ---------------------------------------------------------------------- */

24
src/QEQ/fix_qeq_shielded.cpp
View File

@ -108,7 +108,7 @@ void FixQEqShielded::init_shielding()
int ntypes = atom->ntypes;
for (i = 1; i <= ntypes; ++i)
for (j = 1; j <= ntypes; ++j)
shld[i][j] = pow( gamma[i] * gamma[j], -1.5 );
shld[i][j] = pow(gamma[i] * gamma[j], -1.5);
if (fabs(swa) > 0.01 && comm->me == 0)
error->warning(FLERR,"Fix qeq has non-zero lower Taper radius cutoff");
@ -117,7 +117,7 @@ void FixQEqShielded::init_shielding()
else if (swb < 5 && comm->me == 0)
error->warning(FLERR,"Fix qeq has very low Taper radius cutoff");
d7 = pow( swb - swa, 7 );
d7 = pow(swb - swa, 7);
swa2 = swa*swa;
swa3 = swa2*swa;
swb2 = swb*swb;
@ -126,12 +126,12 @@ void FixQEqShielded::init_shielding()
Tap[7] = 20.0 / d7;
Tap[6] = -70.0 * (swa + swb) / d7;
Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7;
Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3 ) / d7;
Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3 ) / d7;
Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3) / d7;
Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3) / d7;
Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7;
Tap[1] = 140.0 * swa3 * swb3 / d7;
Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 -
7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7;
7.0*swa*swb3*swb3 + swb3*swb3*swb) / d7;
}
/* ---------------------------------------------------------------------- */
@ -171,18 +171,18 @@ void FixQEqShielded::init_matvec()
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -(chi[atom->type[i]] + chizj[i]);
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
comm->forward_comm_fix(this); //Dist_vector(s);
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
comm->forward_comm_fix(this); //Dist_vector(t);
}
/* ---------------------------------------------------------------------- */
@ -225,7 +225,7 @@ void FixQEqShielded::compute_H()
if (r_sqr <= cutoff_sq) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);
H.val[m_fill] = 0.5 * calculate_H( r, shld[type[i]][type[j]] );
H.val[m_fill] = 0.5 * calculate_H(r, shld[type[i]][type[j]]);
m_fill++;
}
}
@ -240,7 +240,7 @@ void FixQEqShielded::compute_H()
/* ---------------------------------------------------------------------- */
double FixQEqShielded::calculate_H( double r, double gamma )
double FixQEqShielded::calculate_H(double r, double gamma)
{
double Taper, denom;

14
src/QEQ/fix_qeq_slater.cpp
View File

@ -147,18 +147,18 @@ void FixQEqSlater::init_matvec()
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -(chi[atom->type[i]] + chizj[i]);
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
comm->forward_comm_fix(this); //Dist_vector(s);
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
comm->forward_comm_fix(this); //Dist_vector(t);
}
/* ---------------------------------------------------------------------- */
@ -346,7 +346,7 @@ double FixQEqSlater::calculate_H_wolf(double zei, double zej, double zj,
/* ---------------------------------------------------------------------- */
void FixQEqSlater::sparse_matvec( sparse_matrix *A, double *x, double *b )
void FixQEqSlater::sparse_matvec(sparse_matrix *A, double *x, double *b)
{
int i, j, itr_j;
@ -368,7 +368,7 @@ void FixQEqSlater::sparse_matvec( sparse_matrix *A, double *x, double *b )
for (i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
for( itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
for(itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];

98
src/USER-OMP/fix_qeq_reax_omp.cpp
View File

@ -212,7 +212,7 @@ void FixQEqReaxOMP::compute_H()
if (flag) {
H.jlist[mfill] = j;
H.val[mfill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] );
H.val[mfill] = calculate_H(sqrt(r_sqr), shld[type[i]][type[j]]);
mfill++;
}
}
@ -307,8 +307,8 @@ void FixQEqReaxOMP::init_matvec()
if (atom->mask[i] & groupbit) {
/* init pre-conditioner for H and init solution vectors */
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -chi[ atom->type[i] ];
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
// Predictor Step
@ -335,8 +335,8 @@ void FixQEqReaxOMP::init_matvec()
if (atom->mask[i] & groupbit) {
/* init pre-conditioner for H and init solution vectors */
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -chi[ atom->type[i] ];
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
/* linear extrapolation for s & t from previous solutions */
@ -344,8 +344,8 @@ void FixQEqReaxOMP::init_matvec()
//t[i] = 2 * t_hist[i][0] - t_hist[i][1];
/* quadratic extrapolation for s & t from previous solutions */
//s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] );
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
//s[i] = s_hist[i][2] + 3 * (s_hist[i][0] - s_hist[i][1]);
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
/* cubic extrapolation for s & t from previous solutions */
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
@ -355,14 +355,14 @@ void FixQEqReaxOMP::init_matvec()
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
comm->forward_comm_fix(this); //Dist_vector(s);
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
comm->forward_comm_fix(this); //Dist_vector(t);
}
/* ---------------------------------------------------------------------- */
int FixQEqReaxOMP::CG( double *b, double *x)
int FixQEqReaxOMP::CG(double *b, double *x)
{
int i;
double alpha, beta, b_norm;
@ -371,8 +371,8 @@ int FixQEqReaxOMP::CG( double *b, double *x)
double my_buf[2], buf[2];
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this); //Coll_Vector( q );
sparse_matvec(&H, x, q);
comm->reverse_comm_fix(this); //Coll_Vector(q);
double tmp1, tmp2;
tmp1 = tmp2 = 0.0;
@ -400,9 +400,9 @@ int FixQEqReaxOMP::CG( double *b, double *x)
sig_new = buf[1];
for (i = 1; i < imax && sqrt(sig_new) / b_norm > tolerance; ++i) {
comm->forward_comm_fix(this); //Dist_vector( d );
sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->forward_comm_fix(this); //Dist_vector(d);
sparse_matvec(&H, d, q);
comm->reverse_comm_fix(this); //Coll_vector(q);
tmp1 = 0.0;
#if defined(_OPENMP)
@ -471,7 +471,7 @@ int FixQEqReaxOMP::CG( double *b, double *x)
/* ---------------------------------------------------------------------- */
void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
void FixQEqReaxOMP::sparse_matvec(sparse_matrix *A, double *x, double *b)
{
#if defined(_OPENMP)
#pragma omp parallel default(shared)
@ -491,7 +491,7 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
#endif
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit) b[i] = eta[ atom->type[i] ] * x[i];
if (atom->mask[i] & groupbit) b[i] = eta[atom->type[i]] * x[i];
}
#if defined(_OPENMP)
@ -586,12 +586,12 @@ void FixQEqReaxOMP::calculate_Q()
}
pack_flag = 4;
comm->forward_comm_fix( this); //Dist_vector( atom->q );
comm->forward_comm_fix(this); //Dist_vector(atom->q);
}
/* ---------------------------------------------------------------------- */
void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
void FixQEqReaxOMP::vector_sum(double* dest, double c, double* v,
double d, double* y, int k)
{
int i;
@ -607,7 +607,7 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
/* ---------------------------------------------------------------------- */
void FixQEqReaxOMP::vector_add( double* dest, double c, double* v, int k)
void FixQEqReaxOMP::vector_add(double* dest, double c, double* v, int k)
{
int i;
@ -626,7 +626,7 @@ void FixQEqReaxOMP::vector_add( double* dest, double c, double* v, int k)
/* dual CG support */
/* ---------------------------------------------------------------------- */
int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
int FixQEqReaxOMP::dual_CG(double *b1, double *b2, double *x1, double *x2)
{
int i;
double alpha_s, alpha_t, beta_s, beta_t, b_norm_s, b_norm_t;
@ -635,8 +635,8 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
double my_buf[4], buf[4];
pack_flag = 5; // forward 2x d and reverse 2x q
dual_sparse_matvec( &H, x1, x2, q );
comm->reverse_comm_fix(this); //Coll_Vector( q );
dual_sparse_matvec(&H, x1, x2, q);
comm->reverse_comm_fix(this); //Coll_Vector(q);
double tmp1, tmp2, tmp3, tmp4;
tmp1 = tmp2 = tmp3 = tmp4 = 0.0;
@ -648,16 +648,16 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
int ii = ilist[jj];
if (atom->mask[ii] & groupbit) {
int indxI = 2 * ii;
r[indxI ] = b1[ii] - q[indxI ];
r[indxI] = b1[ii] - q[indxI];
r[indxI+1] = b2[ii] - q[indxI+1];
d[indxI ] = r[indxI ] * Hdia_inv[ii]; //pre-condition
d[indxI] = r[indxI] * Hdia_inv[ii]; //pre-condition
d[indxI+1] = r[indxI+1] * Hdia_inv[ii];
tmp1 += b1[ii] * b1[ii];
tmp2 += b2[ii] * b2[ii];
tmp3 += r[indxI ] * d[indxI ];
tmp3 += r[indxI] * d[indxI];
tmp4 += r[indxI+1] * d[indxI+1];
}
}
@ -676,9 +676,9 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
sig_new_t = buf[3];
for (i = 1; i < imax; ++i) {
comm->forward_comm_fix(this); //Dist_vector( d );
dual_sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this); //Coll_vector( q );
comm->forward_comm_fix(this); //Dist_vector(d);
dual_sparse_matvec(&H, d, q);
comm->reverse_comm_fix(this); //Coll_vector(q);
tmp1 = tmp2 = 0.0;
#if defined(_OPENMP)
@ -693,7 +693,7 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
int ii = ilist[jj];
if (atom->mask[ii] & groupbit) {
int indxI = 2 * ii;
tmp1 += d[indxI ] * q[indxI ];
tmp1 += d[indxI] * q[indxI];
tmp2 += d[indxI+1] * q[indxI+1];
}
}
@ -722,17 +722,17 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
int ii = ilist[jj];
if (atom->mask[ii] & groupbit) {
int indxI = 2 * ii;
x1[ii] += alpha_s * d[indxI ];
x1[ii] += alpha_s * d[indxI];
x2[ii] += alpha_t * d[indxI+1];
r[indxI ] -= alpha_s * q[indxI ];
r[indxI] -= alpha_s * q[indxI];
r[indxI+1] -= alpha_t * q[indxI+1];
// pre-conditioning
p[indxI ] = r[indxI ] * Hdia_inv[ii];
p[indxI] = r[indxI] * Hdia_inv[ii];
p[indxI+1] = r[indxI+1] * Hdia_inv[ii];
tmp1 += r[indxI ] * p[indxI ];
tmp1 += r[indxI] * p[indxI];
tmp2 += r[indxI+1] * p[indxI+1];
}
}
@ -763,7 +763,7 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
if (atom->mask[ii] & groupbit) {
int indxI = 2 * ii;
d[indxI ] = p[indxI ] + beta_s * d[indxI ];
d[indxI] = p[indxI] + beta_s * d[indxI];
d[indxI+1] = p[indxI+1] + beta_t * d[indxI+1];
}
}
@ -798,7 +798,7 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
/* ---------------------------------------------------------------------- */
void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2, double *b)
void FixQEqReaxOMP::dual_sparse_matvec(sparse_matrix *A, double *x1, double *x2, double *b)
{
#if defined(_OPENMP)
#pragma omp parallel default(shared)
@ -822,8 +822,8 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
i = ilist[ii];
if (atom->mask[i] & groupbit) {
indxI = 2 * i;
b[indxI ] = eta[ atom->type[i] ] * x1[i];
b[indxI+1] = eta[ atom->type[i] ] * x2[i];
b[indxI] = eta[atom->type[i]] * x1[i];
b[indxI+1] = eta[atom->type[i]] * x2[i];
}
}
@ -845,7 +845,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
for (i = 0; i < NN; ++i) {
indxI = 2 * i;
for (int t=0; t<nthreads; t++) {
b_temp[t][indxI ] = 0.0;
b_temp[t][indxI] = 0.0;
b_temp[t][indxI+1] = 0.0;
}
}
@ -862,10 +862,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
for (itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
indxJ = 2 * j;
b[indxI ] += A->val[itr_j] * x1[j];
b[indxI] += A->val[itr_j] * x1[j];
b[indxI+1] += A->val[itr_j] * x2[j];
b_temp[tid][indxJ ] += A->val[itr_j] * x1[i];
b_temp[tid][indxJ] += A->val[itr_j] * x1[i];
b_temp[tid][indxJ+1] += A->val[itr_j] * x2[i];
}
}
@ -879,7 +879,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
for (i = 0; i < NN; ++i) {
indxI = 2 * i;
for (int t = 0; t < nthreads; ++t) {
b[indxI ] += b_temp[t][indxI ];
b[indxI] += b_temp[t][indxI];
b[indxI+1] += b_temp[t][indxI+1];
}
}
@ -889,7 +889,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
/* ---------------------------------------------------------------------- */
void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
void FixQEqReaxOMP::dual_sparse_matvec(sparse_matrix *A, double *x, double *b)
{
#if defined(_OPENMP)
#pragma omp parallel default(shared)
@ -913,8 +913,8 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
i = ilist[ii];
if (atom->mask[i] & groupbit) {
indxI = 2 * i;
b[indxI ] = eta[ atom->type[i] ] * x[indxI ];
b[indxI+1] = eta[ atom->type[i] ] * x[indxI+1];
b[indxI] = eta[atom->type[i]] * x[indxI];
b[indxI+1] = eta[atom->type[i]] * x[indxI+1];
}
}
@ -936,7 +936,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
for (i = 0; i < NN; ++i) {
indxI = 2 * i;
for (int t=0; t<nthreads; t++) {
b_temp[t][indxI ] = 0.0;
b_temp[t][indxI] = 0.0;
b_temp[t][indxI+1] = 0.0;
}
}
@ -953,10 +953,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
for (itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
indxJ = 2 * j;
b[indxI ] += A->val[itr_j] * x[indxJ ];
b[indxI] += A->val[itr_j] * x[indxJ];
b[indxI+1] += A->val[itr_j] * x[indxJ+1];
b_temp[tid][indxJ ] += A->val[itr_j] * x[indxI ];
b_temp[tid][indxJ] += A->val[itr_j] * x[indxI];
b_temp[tid][indxJ+1] += A->val[itr_j] * x[indxI+1];
}
}
@ -970,7 +970,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
for (i = 0; i < NN; ++i) {
indxI = 2 * i;
for (int t = 0; t < nthreads; ++t) {
b[indxI ] += b_temp[t][indxI ];
b[indxI] += b_temp[t][indxI];
b[indxI+1] += b_temp[t][indxI+1];
}
}

184
src/USER-OMP/reaxc_bond_orders_omp.cpp
View File

@ -58,7 +58,7 @@ namespace ReaxFF {
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);
bo_ji = &(bonds->select.bond_list[nbr_j->sym_index].bo_data);
double c = bo_ij->Cdbo + bo_ji->Cdbo;
coef.C1dbo = bo_ij->C1dbo * c;
@ -83,19 +83,19 @@ namespace ReaxFF {
coef.C3dDelta = bo_ij->C3dbo * c;
c = (coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2);
rvec_Scale( temp, c, bo_ij->dBOp );
rvec_Scale( temp, c, bo_ij->dBOp);
c = (coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2);
rvec_ScaledAdd( temp, c, workspace->dDeltap_self[i] );
rvec_ScaledAdd(temp, c, workspace->dDeltap_self[i]);
rvec_ScaledAdd( temp, coef.C1dbopi, bo_ij->dln_BOp_pi );
rvec_ScaledAdd( temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2 );
rvec_ScaledAdd(temp, coef.C1dbopi, bo_ij->dln_BOp_pi);
rvec_ScaledAdd(temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2);
rvec_Add(workspace->forceReduction[reductionOffset+i],temp );
rvec_Add(workspace->forceReduction[reductionOffset+i],temp);
if (system->pair_ptr->vflag_atom) {
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,
fi_tmp[0],fi_tmp[1],fi_tmp[2],
@ -103,20 +103,20 @@ namespace ReaxFF {
}
c = -(coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2);
rvec_Scale( temp, c, bo_ij->dBOp );
rvec_Scale( temp, c, bo_ij->dBOp);
c = (coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2);
rvec_ScaledAdd( temp, c, workspace->dDeltap_self[j] );
rvec_ScaledAdd(temp, c, workspace->dDeltap_self[j]);
rvec_ScaledAdd( temp, -coef.C1dbopi, bo_ij->dln_BOp_pi );
rvec_ScaledAdd( temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2 );
rvec_ScaledAdd(temp, -coef.C1dbopi, bo_ij->dln_BOp_pi);
rvec_ScaledAdd(temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2);
rvec_Add(workspace->forceReduction[reductionOffset+j],temp );
rvec_Add(workspace->forceReduction[reductionOffset+j],temp);
if (system->pair_ptr->vflag_atom) {
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,
fj_tmp[0],fj_tmp[1],fj_tmp[2],
@ -128,15 +128,15 @@ namespace ReaxFF {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
// rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd( temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd( temp, -coef.C3dbopi, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd( temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);
// rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C3dbopi, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);
const double c = -(coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2);
rvec_Scale(temp, c, nbr_k->bo_data.dBOp);
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
rvec_Add(workspace->forceReduction[reductionOffset+k],temp);
if (system->pair_ptr->vflag_atom) {
rvec_Scale(fk_tmp, -1.0, temp);
@ -158,15 +158,15 @@ namespace ReaxFF {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
// rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp );
// rvec_ScaledAdd( temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd( temp, -coef.C4dbopi, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd( temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);
// rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C4dbopi, nbr_k->bo_data.dBOp);
// rvec_ScaledAdd(temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);
const double c = -(coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2);
rvec_Scale(temp, c, nbr_k->bo_data.dBOp);
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
rvec_Add(workspace->forceReduction[reductionOffset+k],temp);
if (system->pair_ptr->vflag_atom) {
rvec_Scale(fk_tmp, -1.0, temp);
@ -203,7 +203,7 @@ namespace ReaxFF {
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);
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);
@ -238,59 +238,59 @@ namespace ReaxFF {
rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);/*3rd, dBOpi2*/
/* force */
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
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_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*/
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 );
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*/
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 );
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_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.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*/
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 );
rvec_Add(workspace->forceReduction[reductionOffset+j],temp);
}
/* ---------------------------------------------------------------------- */
int BOp_OMP( storage * /* workspace */, reax_list *bonds, double bo_cut,
int BOp_OMP(storage * /* workspace */, reax_list *bonds, double bo_cut,
int i, int btop_i, far_neighbor_data *nbr_pj,
single_body_parameters * /* sbp_i */, single_body_parameters * /* sbp_j */,
two_body_parameters *twbp,
@ -309,24 +309,24 @@ namespace ReaxFF {
/* Initially BO values are the uncorrected ones, page 1 */
/****** bonds i-j and j-i ******/
ibond = &( bonds->select.bond_list[btop_i] );
jbond = &( bonds->select.bond_list[btop_j] );
ibond = &(bonds->select.bond_list[btop_i]);
jbond = &(bonds->select.bond_list[btop_j]);
ibond->nbr = j;
jbond->nbr = i;
ibond->d = nbr_pj->d;
jbond->d = nbr_pj->d;
rvec_Copy( ibond->dvec, nbr_pj->dvec );
rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
rvec_Copy(ibond->dvec, nbr_pj->dvec);
rvec_Scale(jbond->dvec, -1, nbr_pj->dvec);
ivec_Copy(ibond->rel_box, nbr_pj->rel_box);
ivec_Scale(jbond->rel_box, -1, nbr_pj->rel_box);
ibond->dbond_index = btop_i;
jbond->dbond_index = btop_i;
ibond->sym_index = btop_j;
jbond->sym_index = btop_i;
bo_ij = &( ibond->bo_data );
bo_ji = &( jbond->bo_data );
bo_ij = &(ibond->bo_data);
bo_ji = &(jbond->bo_data);
bo_ji->BO = bo_ij->BO = BO;
bo_ji->BO_s = bo_ij->BO_s = BO_s;
bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
@ -344,14 +344,14 @@ namespace ReaxFF {
rvec_Scale(bo_ij->dln_BOp_pi2,
-bo_ij->BO_pi2*Cln_BOp_pi2,ibond->dvec);
rvec_Scale(bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s);
rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi);
rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2);
rvec_Scale( bo_ij->dBOp,
rvec_Scale(bo_ij->dBOp,
-(bo_ij->BO_s * Cln_BOp_s +
bo_ij->BO_pi * Cln_BOp_pi +
bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec);
rvec_Scale(bo_ji->dBOp, -1., bo_ij->dBOp);
bo_ij->BO_s -= bo_cut;
bo_ij->BO -= bo_cut;
@ -366,7 +366,7 @@ namespace ReaxFF {
/* ---------------------------------------------------------------------- */
void BOOMP( reax_system *system, storage *workspace, reax_list **lists)
void BOOMP(reax_system *system, storage *workspace, reax_list **lists)
{
double p_lp1 = system->reax_param.gp.l[15];
double p_boc1 = system->reax_param.gp.l[0];
@ -428,10 +428,10 @@ namespace ReaxFF {
j = bonds->select.bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
bo_ij = &( bonds->select.bond_list[pj].bo_data );
bo_ij = &(bonds->select.bond_list[pj].bo_data);
if (i < j || workspace->bond_mark[j] > 3) {
twbp = &( system->reax_param.tbp[type_i][type_j] );
twbp = &(system->reax_param.tbp[type_i][type_j]);
if (twbp->ovc < 0.001 && twbp->v13cor < 0.001) {
bo_ij->C1dbo = 1.000000;
@ -456,38 +456,38 @@ namespace ReaxFF {
/* on page 1 */
if (twbp->ovc >= 0.001) {
/* Correction for overcoordination */
exp_p1i = exp( -p_boc1 * Deltap_i );
exp_p2i = exp( -p_boc2 * Deltap_i );
exp_p1j = exp( -p_boc1 * Deltap_j );
exp_p2j = exp( -p_boc2 * Deltap_j );
exp_p1i = exp(-p_boc1 * Deltap_i);
exp_p2i = exp(-p_boc2 * Deltap_i);
exp_p1j = exp(-p_boc1 * Deltap_j);
exp_p2j = exp(-p_boc2 * Deltap_j);
f2 = exp_p1i + exp_p1j;
f3 = -1.0 / p_boc2 * log( 0.5 * ( exp_p2i + exp_p2j ) );
f1 = 0.5 * ( ( val_i + f2 )/( val_i + f2 + f3 ) +
( val_j + f2 )/( val_j + f2 + f3 ) );
f3 = -1.0 / p_boc2 * log(0.5 * (exp_p2i + exp_p2j));
f1 = 0.5 * ((val_i + f2)/(val_i + f2 + f3) +
(val_j + f2)/(val_j + f2 + f3));
/* Now come the derivates */
/* Bond Order pages 5-7, derivative of f1 */
temp = f2 + f3;
u1_ij = val_i + temp;
u1_ji = val_j + temp;
Cf1A_ij = 0.5 * f3 * (1.0 / SQR( u1_ij ) +
1.0 / SQR( u1_ji ));
Cf1B_ij = -0.5 * (( u1_ij - f3 ) / SQR( u1_ij ) +
( u1_ji - f3 ) / SQR( u1_ji ));
Cf1A_ij = 0.5 * f3 * (1.0 / SQR(u1_ij) +
1.0 / SQR(u1_ji));
Cf1B_ij = -0.5 * ((u1_ij - f3) / SQR(u1_ij) +
(u1_ji - f3) / SQR(u1_ji));
Cf1_ij = 0.50 * ( -p_boc1 * exp_p1i / u1_ij -
Cf1_ij = 0.50 * (-p_boc1 * exp_p1i / u1_ij -
((val_i+f2) / SQR(u1_ij)) *
( -p_boc1 * exp_p1i +
exp_p2i / ( exp_p2i + exp_p2j ) ) +
(-p_boc1 * exp_p1i +
exp_p2i / (exp_p2i + exp_p2j)) +
-p_boc1 * exp_p1i / u1_ji -
((val_j+f2) / SQR(u1_ji)) *
( -p_boc1 * exp_p1i +
exp_p2i / ( exp_p2i + exp_p2j ) ));
(-p_boc1 * exp_p1i +
exp_p2i / (exp_p2i + exp_p2j)));
Cf1_ji = -Cf1A_ij * p_boc1 * exp_p1j +
Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j );
Cf1B_ij * exp_p2j / (exp_p2i + exp_p2j);
} else {
/* No overcoordination correction! */
f1 = 1.0;
@ -496,9 +496,9 @@ namespace ReaxFF {
if (twbp->v13cor >= 0.001) {
/* Correction for 1-3 bond orders */
exp_f4 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
exp_f4 =exp(-(twbp->p_boc4 * SQR(bo_ij->BO) -
Deltap_boc_i) * twbp->p_boc3 + twbp->p_boc5);
exp_f5 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
exp_f5 =exp(-(twbp->p_boc4 * SQR(bo_ij->BO) -
Deltap_boc_j) * twbp->p_boc3 + twbp->p_boc5);
f4 = 1. / (1. + exp_f4);
@ -526,7 +526,7 @@ namespace ReaxFF {
bo_ij->BO = bo_ij->BO * A0_ij;
bo_ij->BO_pi = bo_ij->BO_pi * A0_ij *f1;
bo_ij->BO_pi2= bo_ij->BO_pi2* A0_ij *f1;
bo_ij->BO_s = bo_ij->BO - ( bo_ij->BO_pi + bo_ij->BO_pi2 );
bo_ij->BO_s = bo_ij->BO - (bo_ij->BO_pi + bo_ij->BO_pi2);
bo_ij->C1dbo = A0_ij + bo_ij->BO * A1_ij;
bo_ij->C2dbo = bo_ij->BO * A2_ij;
@ -585,8 +585,8 @@ namespace ReaxFF {
everything else is set in uncorrected_bo calculations */
sym_index = bonds->select.bond_list[pj].sym_index;
bo_ij = &( bonds->select.bond_list[pj].bo_data );
bo_ji = &(bonds->select.bond_list[ sym_index ].bo_data);
bo_ij = &(bonds->select.bond_list[pj].bo_data);
bo_ji = &(bonds->select.bond_list[sym_index].bo_data);
bo_ij->BO = bo_ji->BO;
bo_ij->BO_s = bo_ji->BO_s;
bo_ij->BO_pi = bo_ji->BO_pi;
@ -612,7 +612,7 @@ namespace ReaxFF {
for (j = 0; j < system->N; ++j) {
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
sbp_j = &(system->reax_param.sbp[ type_j ]);
sbp_j = &(system->reax_param.sbp[type_j]);
workspace->Delta[j] = workspace->total_bond_order[j] - sbp_j->valency;
workspace->Delta_e[j] = workspace->total_bond_order[j] - sbp_j->valency_e;

72
src/USER-OMP/reaxc_hydrogen_bonds_omp.cpp
View File

@ -41,7 +41,7 @@ namespace ReaxFF {
/* ---------------------------------------------------------------------- */
void Hydrogen_BondsOMP( reax_system *system, control_params *control,
void Hydrogen_BondsOMP(reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists)
{
@ -105,20 +105,20 @@ namespace ReaxFF {
type_j = system->my_atoms[j].type;
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
hb_start_j = Start_Index( system->my_atoms[j].Hindex, hbonds );
hb_end_j = End_Index( system->my_atoms[j].Hindex, hbonds );
hb_start_j = Start_Index(system->my_atoms[j].Hindex, hbonds);
hb_end_j = End_Index(system->my_atoms[j].Hindex, hbonds);
if (type_j < 0) continue;
top = 0;
for (pi = start_j; pi < end_j; ++pi) {
pbond_ij = &( bond_list[pi] );
pbond_ij = &(bond_list[pi]);
i = pbond_ij->nbr;
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
bo_ij = &(pbond_ij->bo_data);
if ( system->reax_param.sbp[type_i].p_hbond == 2 &&
bo_ij->BO >= HB_THRESHOLD )
if (system->reax_param.sbp[type_i].p_hbond == 2 &&
bo_ij->BO >= HB_THRESHOLD)
hblist[top++] = pi;
}
@ -129,35 +129,35 @@ namespace ReaxFF {
if (type_k < 0) continue;
nbr_jk = hbond_list[pk].ptr;
r_jk = nbr_jk->d;
rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
rvec_Scale(dvec_jk, hbond_list[pk].scl, nbr_jk->dvec);
for (itr = 0; itr < top; ++itr) {
pi = hblist[itr];
pbond_ij = &( bonds->select.bond_list[pi] );
pbond_ij = &(bonds->select.bond_list[pi]);
i = pbond_ij->nbr;
if (system->my_atoms[i].orig_id != system->my_atoms[k].orig_id) {
bo_ij = &(pbond_ij->bo_data);
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]);
hbp = &(system->reax_param.hbp[type_i][type_j][type_k]);
++num_hb_intrs;
Calculate_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&theta, &cos_theta );
Calculate_Theta(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&theta, &cos_theta);
/* the derivative of cos(theta) */
Calculate_dCos_ThetaOMP( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
Calculate_dCos_ThetaOMP(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&dcos_theta_di, &dcos_theta_dj,
&dcos_theta_dk );
&dcos_theta_dk);
/* hydrogen bond energy*/
sin_theta2 = sin( theta/2.0 );
sin_theta2 = sin(theta/2.0);
sin_xhz4 = SQR(sin_theta2);
sin_xhz4 *= sin_xhz4;
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 ) );
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));
e_hb_thr += e_hb = hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
@ -171,38 +171,38 @@ namespace ReaxFF {
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 );
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 );
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_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 );
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 );
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_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 );
rvec_Scale(force, CEhb3/r_jk, dvec_jk);
rvec_Add(workspace->forceReduction[reductionOffset+k], force);
}
/* 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 );
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);

8
src/USER-OMP/reaxc_multi_body_omp.cpp
View File

@ -93,7 +93,7 @@ namespace ReaxFF {
for (i = 0; i < system->n; ++i) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[ type_i ]);
sbp_i = &(system->reax_param.sbp[type_i]);
/* lone-pair Energy */
p_lp2 = sbp_i->p_lp2;
@ -158,7 +158,7 @@ namespace ReaxFF {
for (i = 0; i < system->n; ++i) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[ type_i ]);
sbp_i = &(system->reax_param.sbp[type_i]);
/* over-coordination energy */
if (sbp_i->mass > 21.0)
@ -172,7 +172,7 @@ namespace ReaxFF {
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
bo_ij = &(bonds->select.bond_list[pj].bo_data);
twbp = &(system->reax_param.tbp[ type_i ][ type_j ]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
sum_ovun1 += twbp->p_ovun1 * twbp->De_s * bo_ij->BO;
sum_ovun2 += (workspace->Delta[j] - dfvl*workspace->Delta_lp_temp[j])*
@ -243,7 +243,7 @@ namespace ReaxFF {
pbond = &(bonds->select.bond_list[pj]);
j = pbond->nbr;
bo_ij = &(pbond->bo_data);
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
twbp = &(system->reax_param.tbp[system->my_atoms[i].type]
[system->my_atoms[pbond->nbr].type]);
bo_ij->Cdbo += CEover1 * twbp->p_ovun1 * twbp->De_s; // OvCoor-1st

4
src/USER-OMP/reaxc_nonbonded_omp.cpp
View File

@ -112,8 +112,8 @@ namespace ReaxFF {
if (flag) {
r_ij = nbr_pj->d;
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
[ system->my_atoms[j].type ]);
twbp = &(system->reax_param.tbp[system->my_atoms[i].type]
[system->my_atoms[j].type]);
/* Calculate Taper and its derivative */
// Tap = nbr_pj->Tap; -- precomputed during compte_H

8
src/USER-OMP/reaxc_valence_angles_omp.cpp
View File

@ -276,8 +276,8 @@ namespace ReaxFF {
my_offset = _my_offset[j];
p_val3 = system->reax_param.sbp[ type_j ].p_val3;
p_val5 = system->reax_param.sbp[ type_j ].p_val5;
p_val3 = system->reax_param.sbp[type_j].p_val3;
p_val5 = system->reax_param.sbp[type_j].p_val5;
SBOp = 0, prod_SBO = 1;
for (t = start_j; t < end_j; ++t) {
@ -390,7 +390,7 @@ namespace ReaxFF {
(bo_ij->BO > control->thb_cut) &&
(bo_jk->BO > control->thb_cut) &&
(bo_ij->BO * bo_jk->BO > control->thb_cutsq)) {
thbh = &(system->reax_param.thbp[ type_i ][ type_j ][ type_k ]);
thbh = &(system->reax_param.thbp[type_i][type_j][type_k]);
for (cnt = 0; cnt < thbh->cnt; ++cnt) {
@ -457,7 +457,7 @@ namespace ReaxFF {
exp_pen2ij = exp(-p_pen2 * SQR(BOA_ij - 2.0));
exp_pen2jk = exp(-p_pen2 * SQR(BOA_jk - 2.0));
exp_pen3 = exp(-p_pen3 * workspace->Delta[j]);
exp_pen4 = exp( p_pen4 * workspace->Delta[j]);
exp_pen4 = exp(p_pen4 * workspace->Delta[j]);
trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
f9_Dj = (2.0 + exp_pen3) / trm_pen34;
Cf9j = (-p_pen3 * exp_pen3 * trm_pen34 -

14
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -593,8 +593,8 @@ void FixQEqReax::init_matvec()
if (atom->mask[i] & groupbit) {
/* init pre-conditioner for H and init solution vectors */
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -chi[ atom->type[i] ];
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
/* quadratic extrapolation for s & t from previous solutions */
@ -761,7 +761,7 @@ void FixQEqReax::sparse_matvec(sparse_matrix *A, double *x, double *b)
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
b[i] = eta[ atom->type[i] ] * x[i];
b[i] = eta[atom->type[i]] * x[i];
}
for (ii = nn; ii < NN; ++ii) {
@ -835,7 +835,7 @@ int FixQEqReax::pack_forward_comm(int n, int *list, double *buf,
m = 0;
for (int i = 0; i < n; i++) {
int j = 2 * list[i];
buf[m++] = d[j ];
buf[m++] = d[j];
buf[m++] = d[j+1];
}
return m;
@ -862,7 +862,7 @@ void FixQEqReax::unpack_forward_comm(int n, int first, double *buf)
m = 0;
for (i = first; i < last; i++) {
int j = 2 * i;
d[j ] = buf[m++];
d[j] = buf[m++];
d[j+1] = buf[m++];
}
}
@ -878,7 +878,7 @@ int FixQEqReax::pack_reverse_comm(int n, int first, double *buf)
int last = first + n;
for (i = first; i < last; i++) {
int indxI = 2 * i;
buf[m++] = q[indxI ];
buf[m++] = q[indxI];
buf[m++] = q[indxI+1];
}
return m;
@ -896,7 +896,7 @@ void FixQEqReax::unpack_reverse_comm(int n, int *list, double *buf)
int m = 0;
for (int i = 0; i < n; i++) {
int indxI = 2 * list[i];
q[indxI ] += buf[m++];
q[indxI] += buf[m++];
q[indxI+1] += buf[m++];
}
} else {

6
src/USER-REAXC/fix_qeq_reax.h
View File

@ -127,9 +127,9 @@ class FixQEqReax : public Fix {
virtual int pack_exchange(int, double *);
virtual int unpack_exchange(int, double *);
virtual double parallel_norm( double*, int );
virtual double parallel_dot( double*, double*, int );
virtual double parallel_vector_acc( double*, int );
virtual double parallel_norm(double*, int);
virtual double parallel_dot(double*, double*, int);
virtual double parallel_vector_acc(double*, int);
virtual void vector_sum(double*,double,double*,double,double*,int);
virtual void vector_add(double*, double, double*,int);

6
src/USER-REAXC/fix_reaxc_bonds.cpp
View File

@ -47,7 +47,7 @@ FixReaxCBonds::FixReaxCBonds(LAMMPS *lmp, int narg, char **arg) :
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0 )
if (nevery <= 0)
error->all(FLERR,"Illegal fix reax/c/bonds command");
if (me == 0) {
@ -194,7 +194,7 @@ void FixReaxCBonds::FindBond(struct _reax_list * /*lists*/, int &numbonds)
nj = 0;
for (pj = Start_Index(i, reaxc->api->lists); pj < End_Index(i, reaxc->api->lists); ++pj) {
bo_ij = &( reaxc->api->lists->select.bond_list[pj] );
bo_ij = &(reaxc->api->lists->select.bond_list[pj]);
j = bo_ij->nbr;
jtag = tag[j];
bo_tmp = bo_ij->bo_data.BO;
@ -233,7 +233,7 @@ void FixReaxCBonds::PassBuffer(double *buf, int &nbuf_local)
}
j += (5+numbonds);
if (atom->molecule == nullptr ) buf[j] = 0.0;
if (atom->molecule == nullptr) buf[j] = 0.0;
else buf[j] = atom->molecule[i];
j ++;

228
src/USER-REAXC/reaxc_bond_orders.cpp
View File

@ -31,7 +31,7 @@
#include <cmath>
namespace ReaxFF {
void Add_dBond_to_Forces_NPT( int i, int pj, storage *workspace, reax_list **lists )
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;
@ -44,7 +44,7 @@ namespace ReaxFF {
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);
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);
@ -74,58 +74,58 @@ namespace ReaxFF {
rvec_ScaledAdd(temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);/*3rd, dBOpi2*/
/* force */
rvec_Add( workspace->f[k], temp );
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_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*/
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 );
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*/
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 );
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_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.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*/
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 );
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;
@ -140,7 +140,7 @@ namespace ReaxFF {
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);
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);
@ -161,40 +161,40 @@ namespace ReaxFF {
coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]);
// forces on i
rvec_Scale( temp, coef.C1dbo, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C2dbo, workspace->dDeltap_self[i] );
rvec_ScaledAdd( temp, coef.C1dDelta, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C2dDelta, workspace->dDeltap_self[i] );
rvec_ScaledAdd( temp, coef.C1dbopi, bo_ij->dln_BOp_pi );
rvec_ScaledAdd( temp, coef.C2dbopi, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C3dbopi, workspace->dDeltap_self[i]);
rvec_ScaledAdd( temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2 );
rvec_ScaledAdd( temp, coef.C2dbopi2, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C3dbopi2, workspace->dDeltap_self[i] );
rvec_Add( workspace->f[i], temp );
rvec_Scale( temp, coef.C1dbo, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C2dbo, workspace->dDeltap_self[i]);
rvec_ScaledAdd(temp, coef.C1dDelta, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C2dDelta, workspace->dDeltap_self[i]);
rvec_ScaledAdd(temp, coef.C1dbopi, bo_ij->dln_BOp_pi);
rvec_ScaledAdd(temp, coef.C2dbopi, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C3dbopi, workspace->dDeltap_self[i]);
rvec_ScaledAdd(temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2);
rvec_ScaledAdd(temp, coef.C2dbopi2, bo_ij->dBOp);
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);
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);
}
// forces on j
rvec_Scale( temp, -coef.C1dbo, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C3dbo, workspace->dDeltap_self[j] );
rvec_ScaledAdd( temp, -coef.C1dDelta, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C3dDelta, workspace->dDeltap_self[j]);
rvec_ScaledAdd( temp, -coef.C1dbopi, bo_ij->dln_BOp_pi );
rvec_ScaledAdd( temp, -coef.C2dbopi, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C4dbopi, workspace->dDeltap_self[j]);
rvec_ScaledAdd( temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2 );
rvec_ScaledAdd( temp, -coef.C2dbopi2, bo_ij->dBOp );
rvec_ScaledAdd( temp, coef.C4dbopi2, workspace->dDeltap_self[j]);
rvec_Add( workspace->f[j], temp );
rvec_Scale( temp, -coef.C1dbo, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C3dbo, workspace->dDeltap_self[j]);
rvec_ScaledAdd(temp, -coef.C1dDelta, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C3dDelta, workspace->dDeltap_self[j]);
rvec_ScaledAdd(temp, -coef.C1dbopi, bo_ij->dln_BOp_pi);
rvec_ScaledAdd(temp, -coef.C2dbopi, bo_ij->dBOp);
rvec_ScaledAdd(temp, coef.C4dbopi, workspace->dDeltap_self[j]);
rvec_ScaledAdd(temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2);
rvec_ScaledAdd(temp, -coef.C2dbopi2, bo_ij->dBOp);
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);
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);
}
@ -203,11 +203,11 @@ namespace ReaxFF {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp);
rvec_ScaledAdd( temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);
rvec_ScaledAdd( temp, -coef.C3dbopi, nbr_k->bo_data.dBOp);
rvec_ScaledAdd( temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);
rvec_Add( workspace->f[k], temp );
rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp);
rvec_ScaledAdd(temp, -coef.C2dDelta, nbr_k->bo_data.dBOp);
rvec_ScaledAdd(temp, -coef.C3dbopi, nbr_k->bo_data.dBOp);
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);
@ -223,11 +223,11 @@ namespace ReaxFF {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp );
rvec_ScaledAdd( temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);
rvec_ScaledAdd( temp, -coef.C4dbopi, nbr_k->bo_data.dBOp);
rvec_ScaledAdd( temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);
rvec_Add( workspace->f[k], temp );
rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp);
rvec_ScaledAdd(temp, -coef.C3dDelta, nbr_k->bo_data.dBOp);
rvec_ScaledAdd(temp, -coef.C4dbopi, nbr_k->bo_data.dBOp);
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);
@ -254,18 +254,18 @@ namespace ReaxFF {
r2 = SQR(nbr_pj->d);
if (sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) {
C12 = twbp->p_bo1 * pow( nbr_pj->d / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + bo_cut) * exp( C12 );
C12 = twbp->p_bo1 * pow(nbr_pj->d / twbp->r_s, twbp->p_bo2);
BO_s = (1.0 + bo_cut) * exp(C12);
} else BO_s = C12 = 0.0;
if (sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0) {
C34 = twbp->p_bo3 * pow( nbr_pj->d / twbp->r_p, twbp->p_bo4 );
BO_pi = exp( C34 );
C34 = twbp->p_bo3 * pow(nbr_pj->d / twbp->r_p, twbp->p_bo4);
BO_pi = exp(C34);
} else BO_pi = C34 = 0.0;
if (sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0) {
C56 = twbp->p_bo5 * pow( nbr_pj->d / twbp->r_pp, twbp->p_bo6 );
BO_pi2= exp( C56 );
C56 = twbp->p_bo5 * pow(nbr_pj->d / twbp->r_pp, twbp->p_bo6);
BO_pi2= exp(C56);
} else BO_pi2 = C56 = 0.0;
/* Initially BO values are the uncorrected ones, page 1 */
@ -273,26 +273,26 @@ namespace ReaxFF {
if (BO >= bo_cut) {
/****** bonds i-j and j-i ******/
ibond = &( bonds->select.bond_list[btop_i] );
btop_j = End_Index( j, bonds );
ibond = &(bonds->select.bond_list[btop_i]);
btop_j = End_Index(j, bonds);
jbond = &(bonds->select.bond_list[btop_j]);
ibond->nbr = j;
jbond->nbr = i;
ibond->d = nbr_pj->d;
jbond->d = nbr_pj->d;
rvec_Copy( ibond->dvec, nbr_pj->dvec );
rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
rvec_Copy(ibond->dvec, nbr_pj->dvec);
rvec_Scale(jbond->dvec, -1, nbr_pj->dvec);
ivec_Copy(ibond->rel_box, nbr_pj->rel_box);
ivec_Scale(jbond->rel_box, -1, nbr_pj->rel_box);
ibond->dbond_index = btop_i;
jbond->dbond_index = btop_i;
ibond->sym_index = btop_j;
jbond->sym_index = btop_i;
Set_End_Index( j, btop_j+1, bonds );
Set_End_Index(j, btop_j+1, bonds);
bo_ij = &( ibond->bo_data );
bo_ji = &( jbond->bo_data );
bo_ij = &(ibond->bo_data);
bo_ji = &(jbond->bo_data);
bo_ji->BO = bo_ij->BO = BO;
bo_ji->BO_s = bo_ij->BO_s = BO_s;
bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
@ -310,17 +310,17 @@ namespace ReaxFF {
rvec_Scale(bo_ij->dln_BOp_pi2,
-bo_ij->BO_pi2*Cln_BOp_pi2,ibond->dvec);
rvec_Scale(bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s);
rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi);
rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2);
rvec_Scale( bo_ij->dBOp,
rvec_Scale(bo_ij->dBOp,
-(bo_ij->BO_s * Cln_BOp_s +
bo_ij->BO_pi * Cln_BOp_pi +
bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec);
rvec_Scale(bo_ji->dBOp, -1., bo_ij->dBOp);
rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
rvec_Add(workspace->dDeltap_self[i], bo_ij->dBOp);
rvec_Add(workspace->dDeltap_self[j], bo_ji->dBOp);
bo_ij->BO_s -= bo_cut;
bo_ij->BO -= bo_cut;
@ -383,11 +383,11 @@ namespace ReaxFF {
j = bonds->select.bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
bo_ij = &( bonds->select.bond_list[pj].bo_data );
// fprintf( stderr, "\tj:%d - ubo: %8.3f\n", j+1, bo_ij->BO );
bo_ij = &(bonds->select.bond_list[pj].bo_data);
// fprintf(stderr, "\tj:%d - ubo: %8.3f\n", j+1, bo_ij->BO);
if (i < j || workspace->bond_mark[j] > 3) {
twbp = &( system->reax_param.tbp[type_i][type_j] );
twbp = &(system->reax_param.tbp[type_i][type_j]);
if (twbp->ovc < 0.001 && twbp->v13cor < 0.001) {
bo_ij->C1dbo = 1.000000;
@ -412,36 +412,36 @@ namespace ReaxFF {
/* on page 1 */
if (twbp->ovc >= 0.001) {
/* Correction for overcoordination */
exp_p1i = exp( -p_boc1 * Deltap_i );
exp_p2i = exp( -p_boc2 * Deltap_i );
exp_p1j = exp( -p_boc1 * Deltap_j );
exp_p2j = exp( -p_boc2 * Deltap_j );
exp_p1i = exp(-p_boc1 * Deltap_i);
exp_p2i = exp(-p_boc2 * Deltap_i);
exp_p1j = exp(-p_boc1 * Deltap_j);
exp_p2j = exp(-p_boc2 * Deltap_j);
f2 = exp_p1i + exp_p1j;
f3 = -1.0 / p_boc2 * log( 0.5 * ( exp_p2i + exp_p2j ) );
f1 = 0.5 * ( ( val_i + f2 )/( val_i + f2 + f3 ) +
( val_j + f2 )/( val_j + f2 + f3 ) );
f3 = -1.0 / p_boc2 * log(0.5 * (exp_p2i + exp_p2j));
f1 = 0.5 * ((val_i + f2)/(val_i + f2 + f3) +
(val_j + f2)/(val_j + f2 + f3));
temp = f2 + f3;
u1_ij = val_i + temp;
u1_ji = val_j + temp;
Cf1A_ij = 0.5 * f3 * (1.0 / SQR( u1_ij ) +
1.0 / SQR( u1_ji ));
Cf1B_ij = -0.5 * (( u1_ij - f3 ) / SQR( u1_ij ) +
( u1_ji - f3 ) / SQR( u1_ji ));
Cf1A_ij = 0.5 * f3 * (1.0 / SQR(u1_ij) +
1.0 / SQR(u1_ji));
Cf1B_ij = -0.5 * ((u1_ij - f3) / SQR(u1_ij) +
(u1_ji - f3) / SQR(u1_ji));
Cf1_ij = 0.50 * ( -p_boc1 * exp_p1i / u1_ij -
Cf1_ij = 0.50 * (-p_boc1 * exp_p1i / u1_ij -
((val_i+f2) / SQR(u1_ij)) *
( -p_boc1 * exp_p1i +
exp_p2i / ( exp_p2i + exp_p2j ) ) +
(-p_boc1 * exp_p1i +
exp_p2i / (exp_p2i + exp_p2j)) +
-p_boc1 * exp_p1i / u1_ji -
((val_j+f2) / SQR(u1_ji)) *
( -p_boc1 * exp_p1i +
exp_p2i / ( exp_p2i + exp_p2j ) ));
(-p_boc1 * exp_p1i +
exp_p2i / (exp_p2i + exp_p2j)));
Cf1_ji = -Cf1A_ij * p_boc1 * exp_p1j +
Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j );
Cf1B_ij * exp_p2j / (exp_p2i + exp_p2j);
} else {
/* No overcoordination correction! */
@ -451,9 +451,9 @@ namespace ReaxFF {
if (twbp->v13cor >= 0.001) {
/* Correction for 1-3 bond orders */
exp_f4 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
exp_f4 =exp(-(twbp->p_boc4 * SQR(bo_ij->BO) -
Deltap_boc_i) * twbp->p_boc3 + twbp->p_boc5);
exp_f5 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
exp_f5 =exp(-(twbp->p_boc4 * SQR(bo_ij->BO) -
Deltap_boc_j) * twbp->p_boc3 + twbp->p_boc5);
f4 = 1. / (1. + exp_f4);
@ -481,7 +481,7 @@ namespace ReaxFF {
bo_ij->BO = bo_ij->BO * A0_ij;
bo_ij->BO_pi = bo_ij->BO_pi * A0_ij *f1;
bo_ij->BO_pi2= bo_ij->BO_pi2* A0_ij *f1;
bo_ij->BO_s = bo_ij->BO - ( bo_ij->BO_pi + bo_ij->BO_pi2 );
bo_ij->BO_s = bo_ij->BO - (bo_ij->BO_pi + bo_ij->BO_pi2);
bo_ij->C1dbo = A0_ij + bo_ij->BO * A1_ij;
bo_ij->C2dbo = bo_ij->BO * A2_ij;
@ -515,7 +515,7 @@ namespace ReaxFF {
/* We only need to update bond orders from bo_ji
everything else is set in uncorrected_bo calculations */
sym_index = bonds->select.bond_list[pj].sym_index;
bo_ji = &(bonds->select.bond_list[ sym_index ].bo_data);
bo_ji = &(bonds->select.bond_list[sym_index].bo_data);
bo_ij->BO = bo_ji->BO;
bo_ij->BO_s = bo_ji->BO_s;
bo_ij->BO_pi = bo_ji->BO_pi;
@ -531,7 +531,7 @@ namespace ReaxFF {
for (j = 0; j < system->N; ++j) {
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
sbp_j = &(system->reax_param.sbp[ type_j ]);
sbp_j = &(system->reax_param.sbp[type_j]);
workspace->Delta[j] = workspace->total_bond_order[j] - sbp_j->valency;
workspace->Delta_e[j] = workspace->total_bond_order[j] - sbp_j->valency_e;

20
src/USER-REAXC/reaxc_bonds.cpp
View File

@ -75,17 +75,17 @@ namespace ReaxFF {
/* set the pointers */
type_i = system->my_atoms[i].type;
type_j = system->my_atoms[j].type;
sbp_i = &( system->reax_param.sbp[type_i] );
sbp_j = &( system->reax_param.sbp[type_j] );
twbp = &( system->reax_param.tbp[type_i][type_j] );
bo_ij = &( bonds->select.bond_list[pj].bo_data );
sbp_i = &(system->reax_param.sbp[type_i]);
sbp_j = &(system->reax_param.sbp[type_j]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
bo_ij = &(bonds->select.bond_list[pj].bo_data);
/* calculate the constants */
if (bo_ij->BO_s == 0.0) pow_BOs_be2 = 0.0;
else pow_BOs_be2 = pow( bo_ij->BO_s, twbp->p_be2 );
exp_be12 = exp( twbp->p_be1 * ( 1.0 - pow_BOs_be2 ) );
else pow_BOs_be2 = pow(bo_ij->BO_s, twbp->p_be2);
exp_be12 = exp(twbp->p_be1 * (1.0 - pow_BOs_be2));
CEbo = -twbp->De_s * exp_be12 *
( 1.0 - twbp->p_be1 * twbp->p_be2 * pow_BOs_be2 );
(1.0 - twbp->p_be1 * twbp->p_be2 * pow_BOs_be2);
/* calculate the Bond Energy */
data->my_en.e_bond += ebond =
@ -104,10 +104,10 @@ namespace ReaxFF {
/* Stabilisation terminal triple bond */
if (bo_ij->BO >= 1.00) {
if ( gp37 == 2 ||
if (gp37 == 2 ||
(sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990) ||
(sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990)) {
exphu = exp( -gp7 * SQR(bo_ij->BO - 2.50) );
exphu = exp(-gp7 * SQR(bo_ij->BO - 2.50));
exphua1 = exp(-gp3 * (workspace->total_bond_order[i]-bo_ij->BO));
exphub1 = exp(-gp3 * (workspace->total_bond_order[j]-bo_ij->BO));
exphuov = exp(gp4 * (workspace->Delta[i] + workspace->Delta[j]));
@ -117,7 +117,7 @@ namespace ReaxFF {
data->my_en.e_bond += estriph;
decobdbo = gp10 * exphu * hulpov * (exphua1 + exphub1) *
( gp3 - 2.0 * gp7 * (bo_ij->BO-2.50) );
(gp3 - 2.0 * gp7 * (bo_ij->BO-2.50));
decobdboua = -gp10 * exphu * hulpov *
(gp3*exphua1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
decobdboub = -gp10 * exphu * hulpov *

70
src/USER-REAXC/reaxc_hydrogen_bonds.cpp
View File

@ -67,20 +67,20 @@ namespace ReaxFF {
type_j = system->my_atoms[j].type;
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
hb_start_j = Start_Index( system->my_atoms[j].Hindex, hbonds );
hb_end_j = End_Index( system->my_atoms[j].Hindex, hbonds );
hb_start_j = Start_Index(system->my_atoms[j].Hindex, hbonds);
hb_end_j = End_Index(system->my_atoms[j].Hindex, hbonds);
if (type_j < 0) continue;
top = 0;
for (pi = start_j; pi < end_j; ++pi) {
pbond_ij = &( bond_list[pi] );
pbond_ij = &(bond_list[pi]);
i = pbond_ij->nbr;
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
bo_ij = &(pbond_ij->bo_data);
if ( system->reax_param.sbp[type_i].p_hbond == 2 &&
bo_ij->BO >= HB_THRESHOLD )
if (system->reax_param.sbp[type_i].p_hbond == 2 &&
bo_ij->BO >= HB_THRESHOLD)
hblist[top++] = pi;
}
@ -91,36 +91,36 @@ namespace ReaxFF {
if (type_k < 0) continue;
nbr_jk = hbond_list[pk].ptr;
r_jk = nbr_jk->d;
rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
rvec_Scale(dvec_jk, hbond_list[pk].scl, nbr_jk->dvec);
for (itr = 0; itr < top; ++itr) {
pi = hblist[itr];
pbond_ij = &( bonds->select.bond_list[pi] );
pbond_ij = &(bonds->select.bond_list[pi]);
i = pbond_ij->nbr;
if (system->my_atoms[i].orig_id != system->my_atoms[k].orig_id) {
bo_ij = &(pbond_ij->bo_data);
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]);
hbp = &(system->reax_param.hbp[type_i][type_j][type_k]);
if (hbp->r0_hb <= 0.0) continue;
++num_hb_intrs;
Calculate_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&theta, &cos_theta );
Calculate_Theta(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&theta, &cos_theta);
/* the derivative of cos(theta) */
Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
Calculate_dCos_Theta(pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
&dcos_theta_di, &dcos_theta_dj,
&dcos_theta_dk );
&dcos_theta_dk);
/* hyrogen bond energy*/
sin_theta2 = sin( theta/2.0 );
sin_theta2 = sin(theta/2.0);
sin_xhz4 = SQR(sin_theta2);
sin_xhz4 *= sin_xhz4;
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 ) );
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));
data->my_en.e_hb += e_hb =
hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
@ -135,36 +135,36 @@ namespace ReaxFF {
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 );
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 );
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_Scale(force, +CEhb2, dcos_theta_di); // dcos terms
rvec_Add(workspace->f[i], force);
rvec_ScaledAdd( workspace->f[j], +CEhb2, dcos_theta_dj );
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 );
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_ScaledAdd(workspace->f[j], -CEhb3/r_jk, dvec_jk);
rvec_Scale( force, CEhb3/r_jk, dvec_jk );
rvec_Add( workspace->f[k], force );
rvec_Scale(force, CEhb3/r_jk, dvec_jk);
rvec_Add(workspace->f[k], force);
}
/* 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 );
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);

2
src/USER-REAXC/reaxc_init_md.cpp
View File

@ -54,7 +54,7 @@ namespace ReaxFF {
if (control->hbond_cut > 0)
for (i = 0; i < system->n; ++i) {
atom = &(system->my_atoms[i]);
if (system->reax_param.sbp[ atom->type ].p_hbond == 1 && atom->type >= 0)
if (system->reax_param.sbp[atom->type].p_hbond == 1 && atom->type >= 0)
atom->Hindex = system->numH++;
else atom->Hindex = -1;
}

26
src/USER-REAXC/reaxc_multi_body.cpp
View File

@ -72,8 +72,8 @@ namespace ReaxFF {
/* lone-pair Energy */
p_lp2 = sbp_i->p_lp2;
expvd2 = exp( -75 * workspace->Delta_lp[i] );
inv_expvd2 = 1. / (1. + expvd2 );
expvd2 = exp(-75 * workspace->Delta_lp[i]);
inv_expvd2 = 1. / (1. + expvd2);
numbonds = 0;
e_lp = 0.0;
@ -103,9 +103,9 @@ namespace ReaxFF {
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
if (!strcmp( system->reax_param.sbp[type_j].name, "C" )) {
twbp = &( system->reax_param.tbp[type_i][type_j]);
bo_ij = &( bonds->select.bond_list[pj].bo_data );
if (!strcmp(system->reax_param.sbp[type_j].name, "C")) {
twbp = &(system->reax_param.tbp[type_i][type_j]);
bo_ij = &(bonds->select.bond_list[pj].bo_data);
Di = workspace->Delta[i];
vov3 = bo_ij->BO - Di - 0.040*pow(Di, 4.);
@ -148,16 +148,16 @@ namespace ReaxFF {
sum_ovun1 += twbp->p_ovun1 * twbp->De_s * bo_ij->BO;
sum_ovun2 += (workspace->Delta[j] - dfvl*workspace->Delta_lp_temp[j])*
( bo_ij->BO_pi + bo_ij->BO_pi2 );
(bo_ij->BO_pi + bo_ij->BO_pi2);
}
exp_ovun1 = p_ovun3 * exp( p_ovun4 * sum_ovun2 );
exp_ovun1 = p_ovun3 * exp(p_ovun4 * sum_ovun2);
inv_exp_ovun1 = 1.0 / (1 + exp_ovun1);
Delta_lpcorr = workspace->Delta[i] -
(dfvl * workspace->Delta_lp_temp[i]) * inv_exp_ovun1;
exp_ovun2 = exp( p_ovun2 * Delta_lpcorr );
exp_ovun2 = exp(p_ovun2 * Delta_lpcorr);
inv_exp_ovun2 = 1.0 / (1.0 + exp_ovun2);
DlpVi = 1.0 / (Delta_lpcorr + sbp_i->valency + 1e-8);
@ -166,9 +166,9 @@ namespace ReaxFF {
data->my_en.e_ov += e_ov = sum_ovun1 * CEover1;
CEover2 = sum_ovun1 * DlpVi * inv_exp_ovun2 *
(1.0 - Delta_lpcorr * ( DlpVi + p_ovun2 * exp_ovun2 * inv_exp_ovun2 ));
(1.0 - Delta_lpcorr * (DlpVi + p_ovun2 * exp_ovun2 * inv_exp_ovun2));
CEover3 = CEover2 * (1.0 - dfvl * workspace->dDelta_lp[i] * inv_exp_ovun1 );
CEover3 = CEover2 * (1.0 - dfvl * workspace->dDelta_lp[i] * inv_exp_ovun1);
CEover4 = CEover2 * (dfvl * workspace->Delta_lp_temp[i]) *
p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1);
@ -179,7 +179,7 @@ namespace ReaxFF {
p_ovun5 = sbp_i->p_ovun5;
exp_ovun2n = 1.0 / exp_ovun2;
exp_ovun6 = exp( p_ovun6 * Delta_lpcorr );
exp_ovun6 = exp(p_ovun6 * Delta_lpcorr);
exp_ovun8 = p_ovun7 * exp(p_ovun8 * sum_ovun2);
inv_exp_ovun2n = 1.0 / (1.0 + exp_ovun2n);
inv_exp_ovun8 = 1.0 / (1.0 + exp_ovun8);
@ -194,8 +194,8 @@ namespace ReaxFF {
-p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
CEunder1 = inv_exp_ovun2n *
( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
p_ovun2 * e_un * exp_ovun2n );
(p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
p_ovun2 * e_un * exp_ovun2n);
CEunder2 = -e_un * p_ovun8 * exp_ovun8 * inv_exp_ovun8;
CEunder3 = CEunder1 * (1.0 - dfvl*workspace->dDelta_lp[i]*inv_exp_ovun1);
CEunder4 = CEunder1 * (dfvl*workspace->Delta_lp_temp[i]) *

4
src/USER-REAXC/reaxc_nonbonded.cpp
View File

@ -111,8 +111,8 @@ namespace ReaxFF {
if (flag) {
r_ij = nbr_pj->d;
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
[ system->my_atoms[j].type ]);
twbp = &(system->reax_param.tbp[system->my_atoms[i].type]
[system->my_atoms[j].type]);
Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
Tap = Tap * r_ij + workspace->Tap[5];

2
src/USER-REAXC/reaxc_reset_tools.cpp
View File

@ -43,7 +43,7 @@ namespace ReaxFF {
for (i = 0; i < system->n; ++i) {
atom = &(system->my_atoms[i]);
if (atom->type < 0) continue;
if (system->reax_param.sbp[ atom->type ].p_hbond == 1)
if (system->reax_param.sbp[atom->type].p_hbond == 1)
atom->Hindex = system->numH++;
else atom->Hindex = -1;
}

158
src/USER-REAXC/reaxc_valence_angles.cpp
View File

@ -32,20 +32,20 @@
#include "error.h"
namespace ReaxFF {
void Calculate_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk,
double *theta, double *cos_theta )
void Calculate_Theta(rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk,
double *theta, double *cos_theta)
{
(*cos_theta) = rvec_Dot(dvec_ji,dvec_jk) / ( d_ji * d_jk );
(*cos_theta) = rvec_Dot(dvec_ji,dvec_jk) / (d_ji * d_jk);
if (*cos_theta > 1.) *cos_theta = 1.0;
if (*cos_theta < -1.) *cos_theta = -1.0;
(*theta) = acos( *cos_theta );
(*theta) = acos(*cos_theta);
}
void Calculate_dCos_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk,
void Calculate_dCos_Theta(rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk,
rvec* dcos_theta_di,
rvec* dcos_theta_dj,
rvec* dcos_theta_dk )
rvec* dcos_theta_dk)
{
int t;
double sqr_d_ji = SQR(d_ji);
@ -59,14 +59,14 @@ namespace ReaxFF {
(*dcos_theta_di)[t] = dvec_jk[t] * inv_dists -
Cdot_inv3 * sqr_d_jk * dvec_ji[t];
(*dcos_theta_dj)[t] = -(dvec_jk[t] + dvec_ji[t]) * inv_dists +
Cdot_inv3 * ( sqr_d_jk * dvec_ji[t] + sqr_d_ji * dvec_jk[t] );
Cdot_inv3 * (sqr_d_jk * dvec_ji[t] + sqr_d_ji * dvec_jk[t]);
(*dcos_theta_dk)[t] = dvec_ji[t] * inv_dists -
Cdot_inv3 * sqr_d_ji * dvec_jk[t];
}
}
void Valence_Angles( reax_system *system, control_params *control,
void Valence_Angles(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists)
{
@ -119,17 +119,17 @@ namespace ReaxFF {
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
p_val3 = system->reax_param.sbp[ type_j ].p_val3;
p_val5 = system->reax_param.sbp[ type_j ].p_val5;
p_val3 = system->reax_param.sbp[type_j].p_val3;
p_val5 = system->reax_param.sbp[type_j].p_val5;
SBOp = 0, prod_SBO = 1;
for (t = start_j; t < end_j; ++t) {
bo_jt = &(bonds->select.bond_list[t].bo_data);
SBOp += (bo_jt->BO_pi + bo_jt->BO_pi2);
temp = SQR( bo_jt->BO );
temp = SQR(bo_jt->BO);
temp *= temp;
temp *= temp;
prod_SBO *= exp( -temp );
prod_SBO *= exp(-temp);
}
if (workspace->vlpex[j] >= 0) {
@ -141,50 +141,50 @@ namespace ReaxFF {
}
SBO = SBOp + (1 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
dSBO1 = -8 * prod_SBO * ( workspace->Delta_boc[j] + p_val8 * vlpadj );
dSBO1 = -8 * prod_SBO * (workspace->Delta_boc[j] + p_val8 * vlpadj);
if (SBO <= 0)
SBO2 = 0, CSBO2 = 0;
else if (SBO > 0 && SBO <= 1) {
SBO2 = pow( SBO, p_val9 );
CSBO2 = p_val9 * pow( SBO, p_val9 - 1 );
SBO2 = pow(SBO, p_val9);
CSBO2 = p_val9 * pow(SBO, p_val9 - 1);
}
else if (SBO > 1 && SBO < 2) {
SBO2 = 2 - pow( 2-SBO, p_val9 );
CSBO2 = p_val9 * pow( 2 - SBO, p_val9 - 1 );
SBO2 = 2 - pow(2-SBO, p_val9);
CSBO2 = p_val9 * pow(2 - SBO, p_val9 - 1);
}
else
SBO2 = 2, CSBO2 = 0;
expval6 = exp( p_val6 * workspace->Delta_boc[j] );
expval6 = exp(p_val6 * workspace->Delta_boc[j]);
for (pi = start_j; pi < end_j; ++pi) {
Set_Start_Index( pi, num_thb_intrs, thb_intrs );
Set_Start_Index(pi, num_thb_intrs, thb_intrs);
pbond_ij = &(bonds->select.bond_list[pi]);
bo_ij = &(pbond_ij->bo_data);
BOA_ij = bo_ij->BO - control->thb_cut;
if ( BOA_ij/*bo_ij->BO*/ > 0.0 &&
( j < system->n || pbond_ij->nbr < system->n )) {
if (BOA_ij/*bo_ij->BO*/ > 0.0 &&
(j < system->n || pbond_ij->nbr < system->n)) {
i = pbond_ij->nbr;
type_i = system->my_atoms[i].type;
for (pk = start_j; pk < pi; ++pk) {
start_pk = Start_Index( pk, thb_intrs );
end_pk = End_Index( pk, thb_intrs );
start_pk = Start_Index(pk, thb_intrs);
end_pk = End_Index(pk, thb_intrs);
for (t = start_pk; t < end_pk; ++t)
if (thb_intrs->select.three_body_list[t].thb == i) {
p_ijk = &(thb_intrs->select.three_body_list[num_thb_intrs] );
p_ijk = &(thb_intrs->select.three_body_list[num_thb_intrs]);
p_kji = &(thb_intrs->select.three_body_list[t]);
p_ijk->thb = bonds->select.bond_list[pk].nbr;
p_ijk->pthb = pk;
p_ijk->theta = p_kji->theta;
rvec_Copy( p_ijk->dcos_di, p_kji->dcos_dk );
rvec_Copy( p_ijk->dcos_dj, p_kji->dcos_dj );
rvec_Copy( p_ijk->dcos_dk, p_kji->dcos_di );
rvec_Copy(p_ijk->dcos_di, p_kji->dcos_dk);
rvec_Copy(p_ijk->dcos_dj, p_kji->dcos_dj);
rvec_Copy(p_ijk->dcos_dk, p_kji->dcos_di);
++num_thb_intrs;
break;
@ -197,21 +197,21 @@ namespace ReaxFF {
BOA_jk = bo_jk->BO - control->thb_cut;
k = pbond_jk->nbr;
type_k = system->my_atoms[k].type;
p_ijk = &( thb_intrs->select.three_body_list[num_thb_intrs] );
p_ijk = &(thb_intrs->select.three_body_list[num_thb_intrs]);
Calculate_Theta( pbond_ij->dvec, pbond_ij->d,
Calculate_Theta(pbond_ij->dvec, pbond_ij->d,
pbond_jk->dvec, pbond_jk->d,
&theta, &cos_theta );
&theta, &cos_theta);
Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d,
Calculate_dCos_Theta(pbond_ij->dvec, pbond_ij->d,
pbond_jk->dvec, pbond_jk->d,
&(p_ijk->dcos_di), &(p_ijk->dcos_dj),
&(p_ijk->dcos_dk) );
&(p_ijk->dcos_dk));
p_ijk->thb = k;
p_ijk->pthb = pk;
p_ijk->theta = theta;
sin_theta = sin( theta );
sin_theta = sin(theta);
if (sin_theta < 1.0e-5)
sin_theta = 1.0e-5;
@ -222,11 +222,11 @@ namespace ReaxFF {
(bo_ij->BO > control->thb_cut) &&
(bo_jk->BO > control->thb_cut) &&
(bo_ij->BO * bo_jk->BO > control->thb_cutsq)) {
thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
thbh = &(system->reax_param.thbp[type_i][type_j][type_k]);
for (cnt = 0; cnt < thbh->cnt; ++cnt) {
if (fabs(thbh->prm[cnt].p_val1) > 0.001) {
thbp = &( thbh->prm[cnt] );
thbp = &(thbh->prm[cnt]);
/* ANGLE ENERGY */
p_val1 = thbp->p_val1;
@ -235,26 +235,26 @@ namespace ReaxFF {
p_val7 = thbp->p_val7;
theta_00 = thbp->theta_00;
exp3ij = exp( -p_val3 * pow( BOA_ij, p_val4 ) );
exp3ij = exp(-p_val3 * pow(BOA_ij, p_val4));
f7_ij = 1.0 - exp3ij;
Cf7ij = p_val3 * p_val4 * pow( BOA_ij, p_val4 - 1.0 ) * exp3ij;
Cf7ij = p_val3 * p_val4 * pow(BOA_ij, p_val4 - 1.0) * exp3ij;
exp3jk = exp( -p_val3 * pow( BOA_jk, p_val4 ) );
exp3jk = exp(-p_val3 * pow(BOA_jk, p_val4));
f7_jk = 1.0 - exp3jk;
Cf7jk = p_val3 * p_val4 * pow( BOA_jk, p_val4 - 1.0 ) * exp3jk;
Cf7jk = p_val3 * p_val4 * pow(BOA_jk, p_val4 - 1.0) * exp3jk;
expval7 = exp( -p_val7 * workspace->Delta_boc[j] );
expval7 = exp(-p_val7 * workspace->Delta_boc[j]);
trm8 = 1.0 + expval6 + expval7;
f8_Dj = p_val5 - ( (p_val5 - 1.0) * (2.0 + expval6) / trm8 );
Cf8j = ( (1.0 - p_val5) / SQR(trm8) ) *
( p_val6 * expval6 * trm8 -
(2.0 + expval6) * ( p_val6*expval6 - p_val7*expval7 ) );
f8_Dj = p_val5 - ((p_val5 - 1.0) * (2.0 + expval6) / trm8);
Cf8j = ((1.0 - p_val5) / SQR(trm8)) *
(p_val6 * expval6 * trm8 -
(2.0 + expval6) * (p_val6*expval6 - p_val7*expval7));
theta_0 = 180.0 - theta_00 * (1.0 -
exp(-p_val10 * (2.0 - SBO2)));
theta_0 = DEG2RAD( theta_0 );
theta_0 = DEG2RAD(theta_0);
expval2theta = exp( -p_val2 * SQR(theta_0 - theta) );
expval2theta = exp(-p_val2 * SQR(theta_0 - theta));
if (p_val1 >= 0)
expval12theta = p_val1 * (1.0 - expval2theta);
else // To avoid linear Me-H-Me angles (6/6/06)
@ -267,7 +267,7 @@ namespace ReaxFF {
expval2theta * (theta_0 - theta);
Ctheta_0 = p_val10 * DEG2RAD(theta_00) *
exp( -p_val10 * (2.0 - SBO2) );
exp(-p_val10 * (2.0 - SBO2));
CEval5 = -CEval4 * Ctheta_0 * CSBO2;
CEval6 = CEval5 * dSBO1;
@ -284,16 +284,16 @@ namespace ReaxFF {
p_pen3 = system->reax_param.gp.l[20];
p_pen4 = system->reax_param.gp.l[21];
exp_pen2ij = exp( -p_pen2 * SQR( BOA_ij - 2.0 ) );
exp_pen2jk = exp( -p_pen2 * SQR( BOA_jk - 2.0 ) );
exp_pen3 = exp( -p_pen3 * workspace->Delta[j] );
exp_pen4 = exp( p_pen4 * workspace->Delta[j] );
exp_pen2ij = exp(-p_pen2 * SQR(BOA_ij - 2.0));
exp_pen2jk = exp(-p_pen2 * SQR(BOA_jk - 2.0));
exp_pen3 = exp(-p_pen3 * workspace->Delta[j]);
exp_pen4 = exp( p_pen4 * workspace->Delta[j]);
trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
f9_Dj = ( 2.0 + exp_pen3 ) / trm_pen34;
Cf9j = ( -p_pen3 * exp_pen3 * trm_pen34 -
(2.0 + exp_pen3) * ( -p_pen3 * exp_pen3 +
p_pen4 * exp_pen4 ) ) /
SQR( trm_pen34 );
f9_Dj = (2.0 + exp_pen3) / trm_pen34;
Cf9j = (-p_pen3 * exp_pen3 * trm_pen34 -
(2.0 + exp_pen3) * (-p_pen3 * exp_pen3 +
p_pen4 * exp_pen4)) /
SQR(trm_pen34);
data->my_en.e_pen += e_pen =
p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
@ -310,13 +310,13 @@ namespace ReaxFF {
p_coa3 = system->reax_param.gp.l[38];
p_coa4 = system->reax_param.gp.l[30];
exp_coa2 = exp( p_coa2 * workspace->Delta_val[j] );
exp_coa2 = exp(p_coa2 * workspace->Delta_val[j]);
data->my_en.e_coa += e_coa =
p_coa1 / (1. + exp_coa2) *
exp( -p_coa3 * SQR(workspace->total_bond_order[i]-BOA_ij) ) *
exp( -p_coa3 * SQR(workspace->total_bond_order[k]-BOA_jk) ) *
exp( -p_coa4 * SQR(BOA_ij - 1.5) ) *
exp( -p_coa4 * SQR(BOA_jk - 1.5) );
exp(-p_coa3 * SQR(workspace->total_bond_order[i]-BOA_ij)) *
exp(-p_coa3 * SQR(workspace->total_bond_order[k]-BOA_jk)) *
exp(-p_coa4 * SQR(BOA_ij - 1.5)) *
exp(-p_coa4 * SQR(BOA_jk - 1.5));
CEcoa1 = -2 * p_coa4 * (BOA_ij - 1.5) * e_coa;
CEcoa2 = -2 * p_coa4 * (BOA_jk - 1.5) * e_coa;
@ -335,10 +335,10 @@ namespace ReaxFF {
workspace->CdDelta[k] += CEcoa5;
for (t = start_j; t < end_j; ++t) {
pbond_jt = &( bonds->select.bond_list[t] );
pbond_jt = &(bonds->select.bond_list[t]);
bo_jt = &(pbond_jt->bo_data);
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
temp = CUBE(temp_bo_jt);
pBOjt7 = temp * temp * temp_bo_jt;
bo_jt->Cdbo += (CEval6 * pBOjt7);
@ -347,31 +347,31 @@ namespace ReaxFF {
}
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 );
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_Scale(force, CEval8, p_ijk->dcos_di);
rvec_Add(workspace->f[i], force);
rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
rvec_ScaledAdd(workspace->f[j], CEval8, p_ijk->dcos_dj);
rvec_Scale( force, CEval8, p_ijk->dcos_dk );
rvec_Add( workspace->f[k], force );
rvec_Scale(force, CEval8, p_ijk->dcos_dk);
rvec_Add(workspace->f[k], force);
}
/* tally into per-atom virials */
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
/* 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 );
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;
@ -386,7 +386,7 @@ namespace ReaxFF {
}
}
Set_End_Index(pi, num_thb_intrs, thb_intrs );
Set_End_Index(pi, num_thb_intrs, thb_intrs);
}
}

4
src/USER-REAXC/reaxff_api.h
View File

@ -184,8 +184,8 @@ namespace ReaxFF
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
inline double rvec_Norm( rvec v ) {
return sqrt( SQR(v[0]) + SQR(v[1]) + SQR(v[2]) );
inline double rvec_Norm(rvec v) {
return sqrt(SQR(v[0]) + SQR(v[1]) + SQR(v[2]));
}
inline void rvec_Scale(rvec ret, double c, rvec v) {

2
src/USER-REAXC/reaxff_defs.h
View File

@ -29,7 +29,7 @@
#define CUBE(x) ((x)*(x)*(x))
#define DEG2RAD(a) ((a)*constPI/180.0)
#define RAD2DEG(a) ((a)*180.0/constPI)
#define MAX3(x,y,z) MAX( MAX(x,y), z)
#define MAX3(x,y,z) MAX(MAX(x,y), z)
#define constPI 3.14159265
#define C_ele 332.06371