diff --git a/PG-PuReMD/src/cuda/cuda_valence_angles.cu b/PG-PuReMD/src/cuda/cuda_valence_angles.cu
index 996ab0ef0b736900bf08689a4ba9aca03ea57b08..09c288b7a4905a022b8511f0a3fde8b456a46635 100644
--- a/PG-PuReMD/src/cuda/cuda_valence_angles.cu
+++ b/PG-PuReMD/src/cuda/cuda_valence_angles.cu
@@ -114,7 +114,7 @@ CUDA_GLOBAL void Cuda_Valence_Angles_Part1( reax_atom *my_atoms,
else
{
vlpadj = workspace.nlp[j];
- dSBO2 = (prod_SBO - 1.0) * (1.0 - p_val8 * workspace.dDelta_lp[j]);
+ dSBO2 = (prod_SBO - 1.0) * (1.0 + p_val8 * workspace.dDelta_lp[j]);
}
SBO = SBOp + (1.0 - prod_SBO) * (-workspace.Delta_boc[j] - p_val8 * vlpadj);
diff --git a/PG-PuReMD/src/nonbonded.c b/PG-PuReMD/src/nonbonded.c
index cc86a3a2c60d165abe045308b74f20d7a7767c47..55b3cbefbcaad6290ffe9e44b302c9941ae008a8 100644
--- a/PG-PuReMD/src/nonbonded.c
+++ b/PG-PuReMD/src/nonbonded.c
@@ -206,7 +206,7 @@ void vdW_Coulomb_Energy( reax_system const * const system,
if ( control->virial == 0 )
{
rvec_ScaledAdd( workspace->f[i],
- -(CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
+ -1.0 * (CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
rvec_ScaledAdd( workspace->f[j],
(CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
}
diff --git a/PG-PuReMD/src/valence_angles.c b/PG-PuReMD/src/valence_angles.c
index f227372a2136d619e07c159b9fbc69ea4cb83f1e..e4e8478ca098ae49b5a4fe1a07123ec993dd3e33 100644
--- a/PG-PuReMD/src/valence_angles.c
+++ b/PG-PuReMD/src/valence_angles.c
@@ -186,7 +186,7 @@ void Valence_Angles( reax_system * const system, control_params * const control,
else
{
vlpadj = workspace->nlp[j];
- dSBO2 = (prod_SBO - 1.0) * (1.0 - p_val8 * workspace->dDelta_lp[j]);
+ dSBO2 = (prod_SBO - 1.0) * (1.0 + p_val8 * workspace->dDelta_lp[j]);
}
SBO = SBOp + (1.0 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
@@ -300,259 +300,261 @@ void Valence_Angles( reax_system * const system, control_params * const control,
/* Fortran ReaxFF code hard-codes the constant below
* as of 2019-02-27, so use that for now */
- if ( j < system->n && BOA_jk >= 0.0 && (bo_ij->BO * bo_jk->BO) >= 0.00001 )
-// if ( j < system->n && BOA_jk >= 0.0 && (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut) )
+ if ( j >= system->n || BOA_jk < 0.0 || (bo_ij->BO * bo_jk->BO) < 0.00001 )
+// if ( j < system->n || BOA_jk < 0.0 || (bo_ij->BO * bo_jk->BO) < SQR(control->thb_cut) )
{
- thbh = &system->reax_param.thbp[
- index_thbp( type_i, type_j, type_k, system->reax_param.num_atom_types ) ];
+ continue;
+ }
+
+ thbh = &system->reax_param.thbp[
+ index_thbp( type_i, type_j, type_k, system->reax_param.num_atom_types ) ];
+
+ for ( cnt = 0; cnt < thbh->cnt; ++cnt )
+ {
+ /* valence angle does not exist in the force field */
+ if ( FABS(thbh->prm[cnt].p_val1) < 0.001 )
+ {
+ continue;
+ }
+
+ thbp = &thbh->prm[cnt];
+
+ /* calculate valence angle energy */
+ p_val1 = thbp->p_val1;
+ p_val2 = thbp->p_val2;
+ p_val4 = thbp->p_val4;
+ p_val7 = thbp->p_val7;
+ theta_00 = thbp->theta_00;
+
+ 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;
+
+ 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;
+
+ 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 ));
+
+ theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
+ theta_0 = DEG2RAD( theta_0 );
+
+ expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
+ if ( p_val1 >= 0.0 )
+ {
+ expval12theta = p_val1 - expval2theta;
+ }
+ /* To avoid linear Me-H-Me angles (6/6/06) */
+ else
+ {
+ expval12theta = -expval2theta;
+ }
- for ( cnt = 0; cnt < thbh->cnt; ++cnt )
+ CEval1 = Cf7ij * f7_jk * f8_Dj * expval12theta;
+ CEval2 = Cf7jk * f7_ij * f8_Dj * expval12theta;
+ CEval3 = Cf8j * f7_ij * f7_jk * expval12theta;
+ CEval4 = 2.0 * p_val2 * f7_ij * f7_jk * f8_Dj
+ * expval2theta * (theta_0 - theta);
+
+ Ctheta_0 = p_val10 * DEG2RAD(theta_00)
+ * EXP( -p_val10 * (2.0 - SBO2) );
+
+ CEval5 = CEval4 * Ctheta_0 * CSBO2;
+ CEval6 = CEval5 * dSBO1;
+ CEval7 = CEval5 * dSBO2;
+ CEval8 = CEval4 / sin_theta;
+
+ e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
+ data->my_en.e_ang += e_ang;
+
+ /* calculate penalty for double bonds in valency angles */
+ p_pen1 = thbp->p_pen1;
+
+ 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 );
+
+ e_pen = p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
+ data->my_en.e_pen += e_pen;
+
+ CEpen1 = e_pen * Cf9j / f9_Dj;
+ temp = -2.0 * p_pen2 * e_pen;
+ CEpen2 = temp * (BOA_ij - 2.0);
+ CEpen3 = temp * (BOA_jk - 2.0);
+
+ /* calculate valency angle conjugation energy */
+ p_coa1 = thbp->p_coa1;
+
+ exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
+ e_coa = p_coa1
+ * 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) )
+ / (1.0 + exp_coa2);
+ data->my_en.e_coa += e_coa;
+
+ CEcoa1 = -2.0 * p_coa4 * (BOA_ij - 1.5) * e_coa;
+ CEcoa2 = -2.0 * p_coa4 * (BOA_jk - 1.5) * e_coa;
+ CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1.0 + exp_coa2);
+ CEcoa4 = -2.0 * p_coa3 * (workspace->total_bond_order[i] - BOA_ij) * e_coa;
+ CEcoa5 = -2.0 * p_coa3 * (workspace->total_bond_order[k] - BOA_jk) * e_coa;
+
+ /* calculate force contributions */
+ bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
+ bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
+ workspace->CdDelta[j] += (CEval3 + CEval7) + CEpen1 + CEcoa3;
+ workspace->CdDelta[i] += CEcoa4;
+ workspace->CdDelta[k] += CEcoa5;
+
+ for ( t = start_j; t < end_j; ++t )
{
- /* valence angle does not exist in the force field */
- if ( FABS(thbh->prm[cnt].p_val1) < 0.001 )
- {
- continue;
- }
-
- thbp = &thbh->prm[cnt];
-
- /* calculate valence angle energy */
- p_val1 = thbp->p_val1;
- p_val2 = thbp->p_val2;
- p_val4 = thbp->p_val4;
- p_val7 = thbp->p_val7;
- theta_00 = thbp->theta_00;
-
- 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;
-
- 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;
-
- 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 ));
-
- theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
- theta_0 = DEG2RAD( theta_0 );
-
- expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
- if ( p_val1 >= 0.0 )
- {
- expval12theta = p_val1 - expval2theta;
- }
- /* To avoid linear Me-H-Me angles (6/6/06) */
- else
- {
- expval12theta = -expval2theta;
- }
-
- CEval1 = Cf7ij * f7_jk * f8_Dj * expval12theta;
- CEval2 = Cf7jk * f7_ij * f8_Dj * expval12theta;
- CEval3 = Cf8j * f7_ij * f7_jk * expval12theta;
- CEval4 = 2.0 * p_val2 * f7_ij * f7_jk * f8_Dj
- * expval2theta * (theta_0 - theta);
-
- Ctheta_0 = p_val10 * DEG2RAD(theta_00)
- * EXP( -p_val10 * (2.0 - SBO2) );
-
- CEval5 = CEval4 * Ctheta_0 * CSBO2;
- CEval6 = CEval5 * dSBO1;
- CEval7 = CEval5 * dSBO2;
- CEval8 = CEval4 / sin_theta;
-
- e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
- data->my_en.e_ang += e_ang;
-
- /* calculate penalty for double bonds in valency angles */
- p_pen1 = thbp->p_pen1;
-
- 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 );
-
- e_pen = p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
- data->my_en.e_pen += e_pen;
-
- CEpen1 = e_pen * Cf9j / f9_Dj;
- temp = -2.0 * p_pen2 * e_pen;
- CEpen2 = temp * (BOA_ij - 2.0);
- CEpen3 = temp * (BOA_jk - 2.0);
-
- /* calculate valency angle conjugation energy */
- p_coa1 = thbp->p_coa1;
-
- exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
- e_coa = p_coa1
- * 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) )
- / (1.0 + exp_coa2);
- data->my_en.e_coa += e_coa;
-
- CEcoa1 = -2.0 * p_coa4 * (BOA_ij - 1.5) * e_coa;
- CEcoa2 = -2.0 * p_coa4 * (BOA_jk - 1.5) * e_coa;
- CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1.0 + exp_coa2);
- CEcoa4 = -2.0 * p_coa3 * (workspace->total_bond_order[i] - BOA_ij) * e_coa;
- CEcoa5 = -2.0 * p_coa3 * (workspace->total_bond_order[k] - BOA_jk) * e_coa;
-
- /* calculate force contributions */
- bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
- bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
- workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
- workspace->CdDelta[i] += CEcoa4;
- workspace->CdDelta[k] += CEcoa5;
-
- for ( t = start_j; t < end_j; ++t )
- {
- pbond_jt = &bond_list->bond_list[t];
- bo_jt = &pbond_jt->bo_data;
- temp_bo_jt = bo_jt->BO;
- temp = CUBE( temp_bo_jt );
- pBOjt7 = temp * temp * temp_bo_jt;
-
- bo_jt->Cdbo += CEval6 * pBOjt7;
- bo_jt->Cdbopi += CEval5;
- bo_jt->Cdbopi2 += CEval5;
- }
-
- if ( control->virial == 0 )
- {
- rvec_ScaledAdd( workspace->f[i], CEval8, p_ijk->dcos_di );
- rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
- rvec_ScaledAdd( workspace->f[k], CEval8, p_ijk->dcos_dk );
- }
- else
- {
- /* terms not related to bond order derivatives are
- * added directly into forces and pressure vector/tensor */
- rvec_Scale( force, CEval8, p_ijk->dcos_di );
- rvec_Add( workspace->f[i], force );
- rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
- rvec_Add( data->my_ext_press, ext_press );
-
- 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_iMultiply( ext_press, pbond_jk->rel_box, force );
- rvec_Add( data->my_ext_press, ext_press );
- }
+ pbond_jt = &bond_list->bond_list[t];
+ bo_jt = &pbond_jt->bo_data;
+ temp_bo_jt = bo_jt->BO;
+ temp = CUBE( temp_bo_jt );
+ pBOjt7 = temp * temp * temp_bo_jt;
+
+ bo_jt->Cdbo += CEval6 * pBOjt7;
+ bo_jt->Cdbopi += CEval5;
+ bo_jt->Cdbopi2 += CEval5;
+ }
+
+ if ( control->virial == 0 )
+ {
+ rvec_ScaledAdd( workspace->f[i], CEval8, p_ijk->dcos_di );
+ rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
+ rvec_ScaledAdd( workspace->f[k], CEval8, p_ijk->dcos_dk );
+ }
+ else
+ {
+ /* terms not related to bond order derivatives are
+ * added directly into forces and pressure vector/tensor */
+ rvec_Scale( force, CEval8, p_ijk->dcos_di );
+ rvec_Add( workspace->f[i], force );
+ rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
+ rvec_Add( data->my_ext_press, ext_press );
+
+ 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_iMultiply( ext_press, pbond_jk->rel_box, force );
+ rvec_Add( data->my_ext_press, ext_press );
+ }
#if defined(TEST_ENERGY)
- /*fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
- p_val3, p_val4, BOA_ij, BOA_jk );
- fprintf(out_control->eval, "%13.8f%13.8f%13.8f%13.8f%13.8f\n",
- workspace->Delta_e[j], workspace->vlpex[j],
- dSBO1, dSBO2, vlpadj );
- fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
- f7_ij, f7_jk, f8_Dj, expval12theta );
- fprintf( out_control->eval,
- "%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
- CEval1, CEval2, CEval3, CEval4,
- CEval5, CEval6, CEval7, CEval8 );
-
- fprintf( out_control->eval,
- "%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n",
- p_ijk->dcos_di[0]/sin_theta, p_ijk->dcos_di[1]/sin_theta,
- p_ijk->dcos_di[2]/sin_theta,
- p_ijk->dcos_dj[0]/sin_theta, p_ijk->dcos_dj[1]/sin_theta,
- p_ijk->dcos_dj[2]/sin_theta,
- p_ijk->dcos_dk[0]/sin_theta, p_ijk->dcos_dk[1]/sin_theta,
- p_ijk->dcos_dk[2]/sin_theta);
-
- fprintf( out_control->eval,
- "%6d%6d%6d%15.8f%15.8f\n",
+ /*fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
+ p_val3, p_val4, BOA_ij, BOA_jk );
+ fprintf(out_control->eval, "%13.8f%13.8f%13.8f%13.8f%13.8f\n",
+ workspace->Delta_e[j], workspace->vlpex[j],
+ dSBO1, dSBO2, vlpadj );
+ fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
+ f7_ij, f7_jk, f8_Dj, expval12theta );
+ fprintf( out_control->eval,
+ "%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
+ CEval1, CEval2, CEval3, CEval4,
+ CEval5, CEval6, CEval7, CEval8 );
+
+ fprintf( out_control->eval,
+ "%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n",
+ p_ijk->dcos_di[0]/sin_theta, p_ijk->dcos_di[1]/sin_theta,
+ p_ijk->dcos_di[2]/sin_theta,
+ p_ijk->dcos_dj[0]/sin_theta, p_ijk->dcos_dj[1]/sin_theta,
+ p_ijk->dcos_dj[2]/sin_theta,
+ p_ijk->dcos_dk[0]/sin_theta, p_ijk->dcos_dk[1]/sin_theta,
+ p_ijk->dcos_dk[2]/sin_theta);
+
+ fprintf( out_control->eval,
+ "%6d%6d%6d%15.8f%15.8f\n",
+ system->my_atoms[i].orig_id,
+ system->my_atoms[j].orig_id,
+ system->my_atoms[k].orig_id,
+ RAD2DEG(theta), e_ang );*/
+
+ fprintf( out_control->eval,
+ //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
+ "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
+ system->my_atoms[i].orig_id,
+ system->my_atoms[j].orig_id,
+ system->my_atoms[k].orig_id,
+ RAD2DEG(theta), theta_0, BOA_ij, BOA_jk,
+ e_ang, data->my_en.e_ang );
+
+ fprintf( out_control->epen,
+ //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
+ "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
- RAD2DEG(theta), e_ang );*/
-
- fprintf( out_control->eval,
- //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
- "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[j].orig_id,
- system->my_atoms[k].orig_id,
- RAD2DEG(theta), theta_0, BOA_ij, BOA_jk,
- e_ang, data->my_en.e_ang );
-
- fprintf( out_control->epen,
- //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
- "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[j].orig_id,
- system->my_atoms[k].orig_id,
- RAD2DEG(theta), BOA_ij, BOA_jk, e_pen,
- data->my_en.e_pen );
-
- fprintf( out_control->ecoa,
- //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
- "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[j].orig_id,
- system->my_atoms[k].orig_id,
- RAD2DEG(theta), BOA_ij, BOA_jk,
- e_coa, data->my_en.e_coa );
+ RAD2DEG(theta), BOA_ij, BOA_jk, e_pen,
+ data->my_en.e_pen );
+
+ fprintf( out_control->ecoa,
+ //"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
+ "%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
+ system->my_atoms[i].orig_id,
+ system->my_atoms[j].orig_id,
+ system->my_atoms[k].orig_id,
+ RAD2DEG(theta), BOA_ij, BOA_jk,
+ e_coa, data->my_en.e_coa );
#endif
#if defined(TEST_FORCES)
- /* angle forces */
- Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
- Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
- Add_dDelta( system, lists, j,
- CEval3 + CEval7, workspace->f_ang );
-
- for ( t = start_j; t < end_j; ++t )
- {
- pbond_jt = &bond_list->bond_list[t];
- bo_jt = &pbond_jt->bo_data;
- temp_bo_jt = bo_jt->BO;
- temp = CUBE( temp_bo_jt );
- pBOjt7 = temp * temp * temp_bo_jt;
-
- Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
- workspace->f_ang );
- Add_dBOpinpi2( system, lists, j, t, CEval5, CEval5,
- workspace->f_ang, workspace->f_ang );
- }
-
- rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
- rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
- rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
- /* end angle forces */
-
- /* penalty forces */
- Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
- Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
- Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
- /* end penalty forces */
-
- /* coalition forces */
- Add_dBO( system, lists, j, pi, CEcoa1 - CEcoa4,
- workspace->f_coa );
- Add_dBO( system, lists, j, pk, CEcoa2 - CEcoa5,
- workspace->f_coa );
- Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
- Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
- Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
- /* end coalition forces */
-#endif
+ /* angle forces */
+ Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
+ Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
+ Add_dDelta( system, lists, j,
+ CEval3 + CEval7, workspace->f_ang );
+
+ for ( t = start_j; t < end_j; ++t )
+ {
+ pbond_jt = &bond_list->bond_list[t];
+ bo_jt = &pbond_jt->bo_data;
+ temp_bo_jt = bo_jt->BO;
+ temp = CUBE( temp_bo_jt );
+ pBOjt7 = temp * temp * temp_bo_jt;
+
+ Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
+ workspace->f_ang );
+ Add_dBOpinpi2( system, lists, j, t, CEval5, CEval5,
+ workspace->f_ang, workspace->f_ang );
}
+
+ rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
+ rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
+ rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
+ /* end angle forces */
+
+ /* penalty forces */
+ Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
+ Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
+ Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
+ /* end penalty forces */
+
+ /* coalition forces */
+ Add_dBO( system, lists, j, pi, CEcoa1 - CEcoa4,
+ workspace->f_coa );
+ Add_dBO( system, lists, j, pk, CEcoa2 - CEcoa5,
+ workspace->f_coa );
+ Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
+ Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
+ Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
+ /* end coalition forces */
+#endif
}
}
}
diff --git a/sPuReMD/src/nonbonded.c b/sPuReMD/src/nonbonded.c
index b26e600a9bbdc71ff215d922c1f5dc9f127fdb1f..2c3ae783932df0dd8a21ec6b53cab08343adb49c 100644
--- a/sPuReMD/src/nonbonded.c
+++ b/sPuReMD/src/nonbonded.c
@@ -164,7 +164,8 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
dTap = dTap * r_ij + workspace->Tap[1];
/* vdWaals Calculations */
- if ( system->reax_param.gp.vdw_type == 1 || system->reax_param.gp.vdw_type == 3 )
+ if ( system->reax_param.gp.vdw_type == 1
+ || system->reax_param.gp.vdw_type == 3 )
{
/* shielding */
powr_vdW1 = POW( r_ij, p_vdW1 );
@@ -221,8 +222,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
#endif
/* Coulomb Calculations */
- dr3gamij_1 = r_ij * r_ij * r_ij
- + POW( twbp->gamma, -3.0 );
+ dr3gamij_1 = r_ij * r_ij * r_ij + POW( twbp->gamma, -3.0 );
dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
e_clb = C_ELE * (system->atoms[i].q * system->atoms[j].q) / dr3gamij_3;
e_ele = self_coef * (e_clb * Tap);
diff --git a/sPuReMD/src/valence_angles.c b/sPuReMD/src/valence_angles.c
index 7d7520bf589994639c5a3ca278931531b876c60f..d8d7e5b13ab536afe29b2ffc1c1add78a69488fc 100644
--- a/sPuReMD/src/valence_angles.c
+++ b/sPuReMD/src/valence_angles.c
@@ -64,7 +64,7 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
sqr_d_ji = SQR( d_ji );
sqr_d_jk = SQR( d_jk );
inv_dists = 1.0 / (d_ji * d_jk);
- inv_dists3 = POW( inv_dists, 3 );
+ inv_dists3 = POW( inv_dists, 3.0 );
dot_dvecs = rvec_Dot( dvec_ji, dvec_jk );
Cdot_inv3 = dot_dvecs * inv_dists3;
@@ -73,7 +73,7 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
(*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
+ (*dcos_theta_dj)[t] = -1.0 * (dvec_jk[t] + dvec_ji[t]) * inv_dists
+ Cdot_inv3 * ( sqr_d_jk * dvec_ji[t] + sqr_d_ji * dvec_jk[t] );
(*dcos_theta_dk)[t] = dvec_ji[t] * inv_dists
@@ -377,13 +377,13 @@ void Valence_Angles( reax_system *system, control_params *control,
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 ));
+ - (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 );
expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
- if ( p_val1 >= 0 )
+ if ( p_val1 >= 0.0 )
{
expval12theta = p_val1 - expval2theta;
}
@@ -477,15 +477,15 @@ void Valence_Angles( reax_system *system, control_params *control,
#if defined(_OPENMP)
// #pragma omp atomic
#endif
- bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
+ bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
#if defined(_OPENMP)
// #pragma omp atomic
#endif
- bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
+ bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
#if defined(_OPENMP)
// #pragma omp atomic
#endif
- workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
+ workspace->CdDelta[j] += (CEval3 + CEval7) + CEpen1 + CEcoa3;
#if defined(_OPENMP)
// #pragma omp atomic
#endif