From ac3b81a42c99a959a409f00c8f798eff50fd5a54 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 30 May 2018 00:15:15 -0400
Subject: [PATCH] PG-PuReMD: fix function pointers for bonded forces. Declare
constants appropriately in bonded force computations. Fix sym_index
computation for bond list. Other code clean-up.
---
PG-PuReMD/src/bond_orders.c | 62 ++++++++++++++++-----------------
PG-PuReMD/src/bonds.c | 43 +++++++++++------------
PG-PuReMD/src/ffield.c | 27 +++++++++------
PG-PuReMD/src/forces.c | 38 ++++++++++++--------
PG-PuReMD/src/multi_body.c | 63 +++++++++++++++++-----------------
PG-PuReMD/src/puremd.c | 5 ---
PG-PuReMD/src/tool_box.c | 15 ++++++--
PG-PuReMD/src/torsion_angles.c | 16 ++++-----
PG-PuReMD/src/valence_angles.c | 48 +++++++++++++++-----------
PG-PuReMD/src/valence_angles.h | 4 +--
10 files changed, 172 insertions(+), 149 deletions(-)
diff --git a/PG-PuReMD/src/bond_orders.c b/PG-PuReMD/src/bond_orders.c
index 281ea2c8..22336cb4 100644
--- a/PG-PuReMD/src/bond_orders.c
+++ b/PG-PuReMD/src/bond_orders.c
@@ -460,8 +460,9 @@ int compare_bonds( const void *p1, const void *p2 )
#endif
-void BO( reax_system * const system, control_params * const control, simulation_data * const data,
- storage * const workspace, reax_list ** const lists, output_controls * const out_control )
+void BO( reax_system * const system, control_params * const control,
+ simulation_data * const data, storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control )
{
int i, j, pj, type_i, type_j;
int start_i, end_i, sym_index;
@@ -470,9 +471,12 @@ void BO( reax_system * const system, control_params * const control, simulation_
real f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5;
real exp_p1i, exp_p2i, exp_p1j, exp_p2j;
real temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji;
- real Cf45_ij, Cf45_ji, p_lp1; //u_ij, u_ji
+ real Cf45_ij, Cf45_ji; //u_ij, u_ji
+ const real p_lp1 = system->reax_param.gp.l[15];
real A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji;
- real explp1, p_boc1, p_boc2;
+ real explp1;
+ const real p_boc1 = system->reax_param.gp.l[0];
+ const real p_boc2 = system->reax_param.gp.l[1];
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
bond_order_data *bo_ij, *bo_ji;
@@ -489,10 +493,6 @@ void BO( reax_system * const system, control_params * const control, simulation_
top_dDelta = 0;
#endif
- p_boc1 = system->reax_param.gp.l[0];
- p_boc2 = system->reax_param.gp.l[1];
- p_lp1 = system->reax_param.gp.l[15];
-
#ifdef TEST_FORCES
for( i = 0; i < bond_list->n; ++i )
{
@@ -666,8 +666,8 @@ void BO( reax_system * const system, control_params * const control, simulation_
/* Bond Order page 10, derivative of total bond order */
A0_ij = f1 * f4f5;
- A1_ij = -2.0 * twbp->p_boc3 * twbp->p_boc4 * bo_ij->BO *
- (Cf45_ij + Cf45_ji);
+ A1_ij = -2.0 * twbp->p_boc3 * twbp->p_boc4 * bo_ij->BO
+ * (Cf45_ij + Cf45_ji);
A2_ij = Cf1_ij / f1 + twbp->p_boc3 * Cf45_ij;
A2_ji = Cf1_ji / f1 + twbp->p_boc3 * Cf45_ji;
A3_ij = A2_ij + Cf1_ij / f1;
@@ -784,39 +784,39 @@ void BO( reax_system * const system, control_params * const control, simulation_
/* Calculate some helper variables that are used at many places
* throughout force calculations */
- for ( j = 0; j < system->N; ++j )
+ for ( i = 0; i < system->N; ++i )
{
- type_j = system->my_atoms[j].type;
+ type_j = system->my_atoms[i].type;
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;
- workspace->Delta_boc[j] = workspace->total_bond_order[j] -
+ workspace->Delta[i] = workspace->total_bond_order[i] - sbp_j->valency;
+ workspace->Delta_e[i] = workspace->total_bond_order[i] - sbp_j->valency_e;
+ workspace->Delta_boc[i] = workspace->total_bond_order[i] -
sbp_j->valency_boc;
- workspace->vlpex[j] = workspace->Delta_e[j] -
- 2.0 * (int)(workspace->Delta_e[j] / 2.0);
- explp1 = EXP( -1.0 * p_lp1 * SQR(2.0 + workspace->vlpex[j]) );
- workspace->nlp[j] = explp1 - (int)(workspace->Delta_e[j] / 2.0);
- workspace->Delta_lp[j] = sbp_j->nlp_opt - workspace->nlp[j];
- workspace->Clp[j] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[j]);
+ workspace->vlpex[i] = workspace->Delta_e[i]
+ - 2.0 * (int)(workspace->Delta_e[i] / 2.0);
+ explp1 = EXP( -1.0 * p_lp1 * SQR(2.0 + workspace->vlpex[i]) );
+ workspace->nlp[i] = explp1 - (int)(workspace->Delta_e[i] / 2.0);
+ workspace->Delta_lp[i] = sbp_j->nlp_opt - workspace->nlp[i];
+ workspace->Clp[i] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[i]);
/* Adri uses different dDelta_lp values than the ones in notes... */
- workspace->dDelta_lp[j] = workspace->Clp[j];
- //workspace->dDelta_lp[j] = workspace->Clp[j] + (0.5-workspace->Clp[j]) *
- //((FABS(workspace->Delta_e[j]/2.0 -
- // (int)(workspace->Delta_e[j]/2.0)) < 0.1) ? 1 : 0 );
+ workspace->dDelta_lp[i] = workspace->Clp[i];
+ //workspace->dDelta_lp[i] = workspace->Clp[i] + (0.5-workspace->Clp[i]) *
+ //((FABS(workspace->Delta_e[i]/2.0 -
+ // (int)(workspace->Delta_e[i]/2.0)) < 0.1) ? 1 : 0 );
if ( sbp_j->mass > 21.0 )
{
- workspace->nlp_temp[j] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
- workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
- workspace->dDelta_lp_temp[j] = 0.0;
+ workspace->nlp_temp[i] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
+ workspace->Delta_lp_temp[i] = sbp_j->nlp_opt - workspace->nlp_temp[i];
+ workspace->dDelta_lp_temp[i] = 0.0;
}
else
{
- workspace->nlp_temp[j] = workspace->nlp[j];
- workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
- workspace->dDelta_lp_temp[j] = workspace->Clp[j];
+ workspace->nlp_temp[i] = workspace->nlp[i];
+ workspace->Delta_lp_temp[i] = sbp_j->nlp_opt - workspace->nlp_temp[i];
+ workspace->dDelta_lp_temp[i] = workspace->Clp[i];
}
}
}
diff --git a/PG-PuReMD/src/bonds.c b/PG-PuReMD/src/bonds.c
index c6a39b8a..acad8963 100644
--- a/PG-PuReMD/src/bonds.c
+++ b/PG-PuReMD/src/bonds.c
@@ -42,11 +42,15 @@ void Bonds( reax_system * const system, control_params * const control,
simulation_data * const data, storage * const workspace, reax_list ** const lists,
output_controls * const out_control )
{
- int i, j, pj, natoms;
+ int i, j, pj;
int start_i, end_i;
int type_i, type_j;
real ebond, pow_BOs_be2, exp_be12, CEbo;
- real gp3, gp4, gp7, gp10, gp37;
+ const real gp3 = system->reax_param.gp.l[3];
+ const real gp4 = system->reax_param.gp.l[4];
+ const real gp7 = system->reax_param.gp.l[7];
+ const real gp10 = system->reax_param.gp.l[10];
+ const real gp37 = (int) system->reax_param.gp.l[37];
real exphu, exphua1, exphub1, exphuov, hulpov, estriph;
real decobdbo, decobdboua, decobdboub;
single_body_parameters *sbp_i, *sbp_j;
@@ -54,14 +58,7 @@ void Bonds( reax_system * const system, control_params * const control,
bond_order_data *bo_ij;
reax_list * const bond_list = lists[BONDS];
- gp3 = system->reax_param.gp.l[3];
- gp4 = system->reax_param.gp.l[4];
- gp7 = system->reax_param.gp.l[7];
- gp10 = system->reax_param.gp.l[10];
- gp37 = (int) system->reax_param.gp.l[37];
- natoms = system->n;
-
- for ( i = 0; i < natoms; ++i )
+ for ( i = 0; i < system->n; ++i )
{
start_i = Start_Index( i, bond_list );
end_i = End_Index( i, bond_list );
@@ -72,7 +69,6 @@ void Bonds( reax_system * const system, control_params * const control,
if ( system->my_atoms[i].orig_id <= system->my_atoms[j].orig_id )
{
- /* 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];
@@ -116,25 +112,26 @@ void Bonds( reax_system * const system, control_params * const control,
/* Stabilisation terminal triple bond */
if ( bo_ij->BO >= 1.00 )
{
- if ( gp37 == 2 ||
- (sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990) ||
- (sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) )
+ /* changed to avoid equality checks of doubles, but really need better approach... */
+ if ( gp37 == 2
+ || (sbp_i->mass - 12.0000 < 1e-6 && sbp_j->mass - 15.9990 < 1e-6)
+ || (sbp_j->mass - 12.0000 < 1e-6 && sbp_i->mass - 15.9990 < 1e-6) )
{
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]));
+ 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]) );
hulpov = 1.0 / (1.0 + 25.0 * exphuov);
estriph = gp10 * exphu * hulpov * (exphua1 + exphub1);
data->my_en.e_bond += estriph;
- decobdbo = gp10 * exphu * hulpov * (exphua1 + exphub1) *
- ( 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 *
- (gp3 * exphub1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
+ decobdbo = gp10 * exphu * hulpov * (exphua1 + exphub1)
+ * ( 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 * (gp3 * exphub1
+ + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
bo_ij->Cdbo += decobdbo;
workspace->CdDelta[i] += decobdboua;
diff --git a/PG-PuReMD/src/ffield.c b/PG-PuReMD/src/ffield.c
index e83142e8..1daadd8f 100644
--- a/PG-PuReMD/src/ffield.c
+++ b/PG-PuReMD/src/ffield.c
@@ -64,6 +64,10 @@ void Read_Force_Field_File( const char * const ffield_file, reax_interaction * c
{
n = atoi(tmp[0]);
}
+ else
+ {
+ n = 0;
+ }
if ( n < 1 )
{
@@ -346,7 +350,7 @@ void Read_Force_Field_File( const char * const ffield_file, reax_interaction * c
/* a line of comments */
fgets( s, MAX_LINE, fp );
- for (i = 0; i < l; i++)
+ for ( i = 0; i < l; i++ )
{
/* line 1 */
fgets(s, MAX_LINE, fp);
@@ -415,9 +419,9 @@ void Read_Force_Field_File( const char * const ffield_file, reax_interaction * c
}
/* calculating combination rules and filling up remaining fields. */
- for (i = 0; i < reax->num_atom_types; i++)
+ for ( i = 0; i < reax->num_atom_types; i++ )
{
- for (j = i; j < reax->num_atom_types; j++)
+ for ( j = i; j < reax->num_atom_types; j++ )
{
index1 = i * __N + j;
index2 = j * __N + i;
@@ -454,43 +458,44 @@ void Read_Force_Field_File( const char * const ffield_file, reax_interaction * c
reax->tbp[index1].D =
SQRT(reax->sbp[i].epsilon * reax->sbp[j].epsilon);
-
reax->tbp[index2].D =
SQRT(reax->sbp[j].epsilon * reax->sbp[i].epsilon);
reax->tbp[index1].alpha =
SQRT(reax->sbp[i].alpha * reax->sbp[j].alpha);
-
reax->tbp[index2].alpha =
SQRT(reax->sbp[j].alpha * reax->sbp[i].alpha);
reax->tbp[index1].r_vdW =
2.0 * SQRT(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);
-
reax->tbp[index2].r_vdW =
2.0 * SQRT(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);
reax->tbp[index1].gamma_w =
SQRT(reax->sbp[i].gamma_w * reax->sbp[j].gamma_w);
-
reax->tbp[index2].gamma_w =
SQRT(reax->sbp[j].gamma_w * reax->sbp[i].gamma_w);
reax->tbp[index1].gamma =
POW(reax->sbp[i].gamma * reax->sbp[j].gamma, -1.5);
-
reax->tbp[index2].gamma =
POW(reax->sbp[j].gamma * reax->sbp[i].gamma, -1.5);
/* additions for additional vdWaals interaction types - inner core */
- reax->tbp[index1].rcore = reax->tbp[index2].rcore =
+ reax->tbp[index1].rcore =
SQRT( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );
+ reax->tbp[index2].rcore =
+ SQRT( reax->sbp[j].rcore2 * reax->sbp[i].rcore2 );
- reax->tbp[index1].ecore = reax->tbp[index2].ecore =
+ reax->tbp[index1].ecore =
SQRT( reax->sbp[i].ecore2 * reax->sbp[j].ecore2 );
+ reax->tbp[index2].ecore =
+ SQRT( reax->sbp[j].ecore2 * reax->sbp[i].ecore2 );
- reax->tbp[index1].acore = reax->tbp[index2].acore =
+ reax->tbp[index1].acore =
SQRT( reax->sbp[i].acore2 * reax->sbp[j].acore2 );
+ reax->tbp[index2].acore =
+ SQRT( reax->sbp[j].acore2 * reax->sbp[i].acore2 );
}
}
diff --git a/PG-PuReMD/src/forces.c b/PG-PuReMD/src/forces.c
index 346e6167..0e18c533 100644
--- a/PG-PuReMD/src/forces.c
+++ b/PG-PuReMD/src/forces.c
@@ -64,6 +64,12 @@ typedef enum
} MATRIX_ENTRY_POSITION;
+int compare_bonds( const void *p1, const void *p2 )
+{
+ return ((bond_data *) p1)->nbr - ((bond_data *) p2)->nbr;
+}
+
+
/* placeholder for unused interactions in interaction list
* Interaction_Functions, which is initialized in Init_Force_Functions */
void Dummy_Interaction( reax_system *system, control_params *control,
@@ -75,23 +81,23 @@ void Dummy_Interaction( reax_system *system, control_params *control,
void Init_Force_Functions( control_params * const control )
{
- control->intr_funcs[0] = BO;
- control->intr_funcs[1] = Bonds;
- control->intr_funcs[2] = Atom_Energy;
- control->intr_funcs[3] = Valence_Angles;
- control->intr_funcs[4] = Torsion_Angles;
+ control->intr_funcs[0] = &BO;
+ control->intr_funcs[1] = &Bonds;
+ control->intr_funcs[2] = &Atom_Energy;
+ control->intr_funcs[3] = &Valence_Angles;
+ control->intr_funcs[4] = &Torsion_Angles;
if ( control->hbond_cut > 0.0 )
{
- control->intr_funcs[5] = Hydrogen_Bonds;
+ control->intr_funcs[5] = &Hydrogen_Bonds;
}
else
{
- control->intr_funcs[5] = Dummy_Interaction;
+ control->intr_funcs[5] = &Dummy_Interaction;
}
- control->intr_funcs[6] = Dummy_Interaction;
- control->intr_funcs[7] = Dummy_Interaction;
- control->intr_funcs[8] = Dummy_Interaction;
- control->intr_funcs[9] = Dummy_Interaction;
+ control->intr_funcs[6] = &Dummy_Interaction;
+ control->intr_funcs[7] = &Dummy_Interaction;
+ control->intr_funcs[8] = &Dummy_Interaction;
+ control->intr_funcs[9] = &Dummy_Interaction;
}
@@ -503,7 +509,7 @@ int Init_Forces( reax_system * const system, control_params * const control,
#if !defined(HALF_LIST)
/* set sym_index for bonds list (far_nbrs full list) */
- for ( i = 0; i < system->n; ++i )
+ for ( i = 0; i < system->N; ++i )
{
start_i = Start_Index( i, bond_list );
end_i = End_Index( i, bond_list );
@@ -534,7 +540,8 @@ int Init_Forces( reax_system * const system, control_params * const control,
workspace->realloc.bonds = TRUE;
}
- if ( system->reax_param.sbp[ system->my_atoms[i].type ].p_hbond == H_ATOM
+ if ( i < system->n
+ && system->reax_param.sbp[ system->my_atoms[i].type ].p_hbond == H_ATOM
&& Num_Entries( system->my_atoms[i].Hindex, hbond_list )
> system->max_hbonds[system->my_atoms[i].Hindex] )
{
@@ -740,7 +747,7 @@ int Init_Forces_No_Charges( reax_system * const system, control_params * const c
#if !defined(HALF_LIST)
/* set sym_index for bonds list (far_nbrs full list) */
- for ( i = 0; i < system->n; ++i )
+ for ( i = 0; i < system->N; ++i )
{
start_i = Start_Index( i, bond_list );
end_i = End_Index( i, bond_list );
@@ -771,7 +778,8 @@ int Init_Forces_No_Charges( reax_system * const system, control_params * const c
workspace->realloc.bonds = TRUE;
}
- if ( system->reax_param.sbp[ system->my_atoms[i].type ].p_hbond == H_ATOM
+ if ( i < system->n
+ && system->reax_param.sbp[ system->my_atoms[i].type ].p_hbond == H_ATOM
&& Num_Entries( system->my_atoms[i].Hindex, hbond_list )
> system->max_hbonds[system->my_atoms[i].Hindex] )
{
diff --git a/PG-PuReMD/src/multi_body.c b/PG-PuReMD/src/multi_body.c
index 6c3e3886..7214d6b4 100644
--- a/PG-PuReMD/src/multi_body.c
+++ b/PG-PuReMD/src/multi_body.c
@@ -100,7 +100,7 @@ void Atom_Energy( reax_system * const system, control_params * const control,
/* correction for C2 */
if ( system->reax_param.gp.l[5] > 0.001 &&
- !strncmp( system->reax_param.sbp[type_i].name, "C", 15 ) )
+ strncmp( system->reax_param.sbp[type_i].name, "C", 15 ) != 0 )
{
for ( pj = Start_Index(i, bond_list); pj < End_Index(i, bond_list); ++pj )
{
@@ -110,15 +110,15 @@ void Atom_Energy( reax_system * const system, control_params * const control,
j = bond_list->bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
- if ( !strncmp( system->reax_param.sbp[type_j].name, "C", 15 ) )
+ if ( strncmp( system->reax_param.sbp[type_j].name, "C", 15 ) != 0 )
{
twbp = &system->reax_param.tbp[
index_tbp(type_i, type_j, system->reax_param.num_atom_types) ];
bo_ij = &bond_list->bond_list[pj].bo_data;
Di = workspace->Delta[i];
- vov3 = bo_ij->BO - Di - 0.040 * POW(Di, 4.);
+ vov3 = bo_ij->BO - Di - 0.040 * POW( Di, 4.0 );
- if ( vov3 > 3. )
+ if ( vov3 > 3.0 )
{
e_lph = p_lp3 * SQR( vov3 - 3.0 );
data->my_en.e_lp += e_lph;
@@ -131,7 +131,7 @@ void Atom_Energy( reax_system * const system, control_params * const control,
workspace->CdDelta[i] += deahu2dsbo;
#ifdef TEST_ENERGY
- fprintf(out_control->elp, "C2cor%6d%6d%12.6f%12.6f%12.6f\n",
+ fprintf( out_control->elp, "C2cor%6d%6d%12.6f%12.6f%12.6f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
e_lph, deahu2dbo, deahu2dsbo );
#endif
@@ -164,7 +164,9 @@ void Atom_Energy( reax_system * const system, control_params * const control,
}
p_ovun2 = sbp_i->p_ovun2;
- sum_ovun1 = sum_ovun2 = 0;
+ sum_ovun1 = 0.0;
+ sum_ovun2 = 0.0;
+
for ( pj = Start_Index(i, bond_list); pj < End_Index(i, bond_list); ++pj )
{
j = bond_list->bond_list[pj].nbr;
@@ -180,8 +182,8 @@ void Atom_Energy( reax_system * const system, control_params * const control,
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;
+ Delta_lpcorr = workspace->Delta[i]
+ - (dfvl * workspace->Delta_lp_temp[i]) * inv_exp_ovun1;
exp_ovun2 = EXP( p_ovun2 * Delta_lpcorr );
inv_exp_ovun2 = 1.0 / (1.0 + exp_ovun2);
@@ -192,13 +194,13 @@ void Atom_Energy( reax_system * const system, control_params * const control,
e_ov = sum_ovun1 * CEover1;
data->my_en.e_ov += e_ov;
- CEover2 = sum_ovun1 * DlpVi * inv_exp_ovun2 *
- (1.0 - Delta_lpcorr * ( DlpVi + p_ovun2 * exp_ovun2 * inv_exp_ovun2 ));
+ CEover2 = sum_ovun1 * DlpVi * 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 );
- CEover4 = CEover2 * (dfvl * workspace->Delta_lp_temp[i]) *
- p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1);
+ CEover4 = CEover2 * (dfvl * workspace->Delta_lp_temp[i])
+ * p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1);
/* under-coordination potential */
p_ovun2 = sbp_i->p_ovun2;
@@ -213,13 +215,12 @@ void Atom_Energy( reax_system * const system, control_params * const control,
e_un = -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
data->my_en.e_un += e_un;
- CEunder1 = inv_exp_ovun2n *
- ( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
- p_ovun2 * e_un * exp_ovun2n );
+ CEunder1 = inv_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]) *
- p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
+ CEunder4 = CEunder1 * (dfvl * workspace->Delta_lp_temp[i])
+ * p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
/* forces */
workspace->CdDelta[i] += CEover3; // OvCoor - 2nd term
@@ -240,20 +241,20 @@ void Atom_Energy( reax_system * const system, control_params * const control,
system->reax_param.num_atom_types) ];
bo_ij->Cdbo += CEover1 * twbp->p_ovun1 * twbp->De_s;// OvCoor-1st
- workspace->CdDelta[j] += CEover4 * (1.0 - dfvl * workspace->dDelta_lp[j]) *
- (bo_ij->BO_pi + bo_ij->BO_pi2); // OvCoor-3a
- bo_ij->Cdbopi += CEover4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); // OvCoor-3b
- bo_ij->Cdbopi2 += CEover4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); // OvCoor-3b
-
-
- workspace->CdDelta[j] += CEunder4 * (1.0 - dfvl * workspace->dDelta_lp[j]) *
- (bo_ij->BO_pi + bo_ij->BO_pi2); // UnCoor - 2a
- bo_ij->Cdbopi += CEunder4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); // UnCoor-2b
- bo_ij->Cdbopi2 += CEunder4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); // UnCoor-2b
+ workspace->CdDelta[j] += CEover4 * (1.0 - dfvl * workspace->dDelta_lp[j])
+ * (bo_ij->BO_pi + bo_ij->BO_pi2); // OvCoor-3a
+ bo_ij->Cdbopi += CEover4 * (workspace->Delta[j] - dfvl
+ * workspace->Delta_lp_temp[j]); // OvCoor-3b
+ bo_ij->Cdbopi2 += CEover4 * (workspace->Delta[j] - dfvl
+ * workspace->Delta_lp_temp[j]); // OvCoor-3b
+
+
+ workspace->CdDelta[j] += CEunder4 * (1.0 - dfvl * workspace->dDelta_lp[j])
+ * (bo_ij->BO_pi + bo_ij->BO_pi2); // UnCoor - 2a
+ bo_ij->Cdbopi += CEunder4 * (workspace->Delta[j] - dfvl
+ * workspace->Delta_lp_temp[j]); // UnCoor-2b
+ bo_ij->Cdbopi2 += CEunder4 * (workspace->Delta[j] - dfvl
+ * workspace->Delta_lp_temp[j]); // UnCoor-2b
#ifdef TEST_ENERGY
/* fprintf( out_control->eov, "%6d%12.6f\n",
diff --git a/PG-PuReMD/src/puremd.c b/PG-PuReMD/src/puremd.c
index 326d4296..2b74e001 100644
--- a/PG-PuReMD/src/puremd.c
+++ b/PG-PuReMD/src/puremd.c
@@ -462,11 +462,6 @@ int simulate( const void * const handle )
Reset( system, control, data, workspace, lists );
- if ( ret == FAILURE )
- {
- ReAllocate( system, control, data, workspace, lists, mpi_data );
- }
-
ret = Generate_Neighbor_Lists( system, data, workspace, lists );
if ( ret != SUCCESS )
diff --git a/PG-PuReMD/src/tool_box.c b/PG-PuReMD/src/tool_box.c
index 96e57b80..36e6b43b 100644
--- a/PG-PuReMD/src/tool_box.c
+++ b/PG-PuReMD/src/tool_box.c
@@ -255,15 +255,24 @@ int Tokenize( const char* s, char*** tok )
{
char test[MAX_LINE];
char *sep = "\t \n!=";
- char *word, *saveptr;
+ char *word;
+ char *saveptr = NULL;
int count = 0;
+ if ( s == NULL )
+ {
+ fprintf( stderr, "[WARNING] passed null string to tokenizer. Returning...\n" );
+ return count;
+ }
+
strncpy( test, s, MAX_LINE );
- for ( word = strtok_r(test, sep, &saveptr); word; word = strtok_r(NULL, sep, &saveptr) )
+ for ( word = strtok_r( test, sep, &saveptr );
+ word != NULL;
+ word = strtok_r( NULL, sep, &saveptr ) )
{
strncpy( (*tok)[count], word, MAX_LINE );
- count++;
+ ++count;
}
return count;
diff --git a/PG-PuReMD/src/torsion_angles.c b/PG-PuReMD/src/torsion_angles.c
index a6722b21..7c1c428b 100644
--- a/PG-PuReMD/src/torsion_angles.c
+++ b/PG-PuReMD/src/torsion_angles.c
@@ -40,8 +40,8 @@
#define MIN_SINE 1e-10
-static real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
- rvec dvec_kl, real r_kl, rvec dvec_li, real r_li,
+static real Calculate_Omega( const rvec dvec_ij, real r_ij, const rvec dvec_jk, real r_jk,
+ const rvec dvec_kl, real r_kl, const rvec dvec_li, real r_li,
const three_body_interaction_data * const p_ijk, const three_body_interaction_data * const p_jkl,
rvec dcos_omega_di, rvec dcos_omega_dj, rvec dcos_omega_dk, rvec dcos_omega_dl )
{
@@ -101,25 +101,25 @@ static real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
arg = -1.0;
}
- if ( sin_ijk >= 0 && sin_ijk <= MIN_SINE )
+ if ( sin_ijk >= 0.0 && sin_ijk <= MIN_SINE )
{
sin_ijk = MIN_SINE;
}
- else if ( sin_ijk <= 0 && sin_ijk >= -MIN_SINE )
+ else if ( sin_ijk <= 0.0 && sin_ijk >= -MIN_SINE )
{
sin_ijk = -MIN_SINE;
}
- if ( sin_jkl >= 0 && sin_jkl <= MIN_SINE )
+ if ( sin_jkl >= 0.0 && sin_jkl <= MIN_SINE )
{
sin_jkl = MIN_SINE;
}
- else if ( sin_jkl <= 0 && sin_jkl >= -MIN_SINE )
+ else if ( sin_jkl <= 0.0 && sin_jkl >= -MIN_SINE )
{
sin_jkl = -MIN_SINE;
}
/* dcos_omega_di */
- rvec_ScaledSum( dcos_omega_di, (htra - arg * hnra) / r_ij, dvec_ij, -1., dvec_li );
+ rvec_ScaledSum( dcos_omega_di, (htra - arg * hnra) / r_ij, dvec_ij, -1.0, dvec_li );
rvec_ScaledAdd( dcos_omega_di, -(hthd - arg * hnhd) / sin_ijk, p_ijk->dcos_dk );
rvec_Scale( dcos_omega_di, 2.0 / poem, dcos_omega_di );
@@ -138,7 +138,7 @@ static real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
rvec_Scale( dcos_omega_dk, 2.0 / poem, dcos_omega_dk );
/* dcos_omega_dl */
- rvec_ScaledSum( dcos_omega_dl, (htrc - arg * hnrc) / r_kl, dvec_kl, 1., dvec_li );
+ rvec_ScaledSum( dcos_omega_dl, (htrc - arg * hnrc) / r_kl, dvec_kl, 1.0, dvec_li );
rvec_ScaledAdd( dcos_omega_dl, -(hthe - arg * hnhe) / sin_jkl, p_jkl->dcos_dk );
rvec_Scale( dcos_omega_dl, 2.0 / poem, dcos_omega_dl );
diff --git a/PG-PuReMD/src/valence_angles.c b/PG-PuReMD/src/valence_angles.c
index 7fd5f497..63582c2c 100644
--- a/PG-PuReMD/src/valence_angles.c
+++ b/PG-PuReMD/src/valence_angles.c
@@ -35,18 +35,18 @@
/* calculates the theta angle between i-j-k */
-void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
+void Calculate_Theta( const rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
real * const theta, real * const cos_theta )
{
(*cos_theta) = Dot( dvec_ji, dvec_jk, 3 ) / ( d_ji * d_jk );
if ( *cos_theta > 1.0 )
{
- *cos_theta = 1.0;
+ *cos_theta = 1.0;
}
if ( *cos_theta < -1.0 )
{
- *cos_theta = -1.0;
+ *cos_theta = -1.0;
}
(*theta) = ACOS( *cos_theta );
@@ -54,7 +54,7 @@ void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
/* calculates the derivative of the cosine of the angle between i-j-k */
-void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
+void Calculate_dCos_Theta( const rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
rvec * const dcos_theta_di, rvec * const dcos_theta_dj, rvec * const dcos_theta_dk )
{
int t;
@@ -101,8 +101,14 @@ void Valence_Angles( reax_system * const system, control_params * const control,
const real p_val8 = system->reax_param.gp.l[33];
const real p_val9 = system->reax_param.gp.l[16];
const real p_val10 = system->reax_param.gp.l[17];
- real p_pen1, p_pen2, p_pen3, p_pen4;
- real p_coa1, p_coa2, p_coa3, p_coa4;
+ real p_pen1;
+ const real p_pen2 = system->reax_param.gp.l[19];
+ const real p_pen3 = system->reax_param.gp.l[20];
+ const real p_pen4 = system->reax_param.gp.l[21];
+ real p_coa1;
+ const real p_coa2 = system->reax_param.gp.l[2];
+ const real p_coa3 = system->reax_param.gp.l[38];
+ const real p_coa4 = system->reax_param.gp.l[30];
real trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk;
real exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2;
real dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO, vlpadj;
@@ -141,7 +147,7 @@ void Valence_Angles( reax_system * const system, control_params * const control,
for ( t = start_j; t < end_j; ++t )
{
bo_jt = &bond_list->bond_list[t].bo_data;
- SBOp += (bo_jt->BO_pi + bo_jt->BO_pi2);
+ SBOp += bo_jt->BO_pi + bo_jt->BO_pi2;
temp = SQR( bo_jt->BO );
temp *= temp;
temp *= temp;
@@ -193,8 +199,8 @@ void Valence_Angles( reax_system * const system, control_params * const control,
bo_ij = &pbond_ij->bo_data;
BOA_ij = bo_ij->BO - control->thb_cut;
- if ( BOA_ij > 0.0 &&
- ( j < system->n || pbond_ij->nbr < system->n ) )
+ if ( BOA_ij > 0.0
+ && ( j < system->n || pbond_ij->nbr < system->n ) )
{
i = pbond_ij->nbr;
type_i = system->my_atoms[i].type;
@@ -326,9 +332,6 @@ void Valence_Angles( reax_system * const system, control_params * const control,
/* PENALTY ENERGY */
p_pen1 = thbp->p_pen1;
- p_pen2 = system->reax_param.gp.l[19];
- 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 ) );
@@ -351,9 +354,6 @@ void Valence_Angles( reax_system * const system, control_params * const control,
/* COALITION ENERGY */
p_coa1 = thbp->p_coa1;
- p_coa2 = system->reax_param.gp.l[2];
- p_coa3 = system->reax_param.gp.l[38];
- p_coa4 = system->reax_param.gp.l[30];
exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
e_coa = p_coa1 / (1.0 + exp_coa2)
@@ -373,9 +373,9 @@ void Valence_Angles( reax_system * const system, control_params * const control,
/* END COALITION ENERGY */
/* FORCES */
- bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
- bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
- workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
+ 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;
@@ -387,7 +387,7 @@ void Valence_Angles( reax_system * const system, control_params * const control,
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
- bo_jt->Cdbo += (CEval6 * pBOjt7);
+ bo_jt->Cdbo += CEval6 * pBOjt7;
bo_jt->Cdbopi += CEval5;
bo_jt->Cdbopi2 += CEval5;
}
@@ -523,10 +523,18 @@ void Valence_Angles( reax_system * const system, control_params * const control,
}
}
- if ( num_thb_intrs >= thb_list->max_intrs * DANGER_ZONE )
+ if ( num_thb_intrs >= (int)(thb_list->max_intrs * DANGER_ZONE) )
{
system->total_thbodies = MAX( num_thb_intrs * SAFE_ZONE, MIN_3BODIES );
workspace->realloc.thbody = TRUE;
+
+ //TODO: need to refactor Compute_Bonded_Forces to allow retry logic
+ if ( num_thb_intrs > thb_list->max_intrs )
+ {
+ fprintf( stderr, "[ERROR] step%d-ran out of space on angle_list: top=%d, max=%d",
+ data->step, num_thb_intrs, thb_list->max_intrs );
+ MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
+ }
}
#if defined(DEBUG)
diff --git a/PG-PuReMD/src/valence_angles.h b/PG-PuReMD/src/valence_angles.h
index f973d42d..5dccbaa5 100644
--- a/PG-PuReMD/src/valence_angles.h
+++ b/PG-PuReMD/src/valence_angles.h
@@ -29,10 +29,10 @@
extern "C" {
#endif
-void Calculate_Theta( rvec, real, rvec, real,
+void Calculate_Theta( const rvec, real, rvec, real,
real * const, real * const );
-void Calculate_dCos_Theta( rvec, real, rvec, real,
+void Calculate_dCos_Theta( const rvec, real, rvec, real,
rvec * const, rvec * const, rvec * const );
void Valence_Angles( reax_system * const, control_params * const, simulation_data * const,
--
GitLab