From 7a138b74eb754e752ac70f2dd220bbb718bb1309 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 8 Nov 2018 10:15:13 -0800
Subject: [PATCH] Fix issue with bonds list ordering leading to different
torsion and conjugation energies due to changing the far neighbors list to a
full list.
---
PuReMD/src/bond_orders.c | 2 +-
PuReMD/src/comm_tools.c | 2 +-
PuReMD/src/forces.c | 12 ++++++++++++
PuReMD/src/torsion_angles.c | 8 ++++----
4 files changed, 18 insertions(+), 6 deletions(-)
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index b1e9b228..da07f839 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -802,7 +802,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
}
-int compare_bonds( const void *p1, const void *p2 )
+static int compare_bonds( const void *p1, const void *p2 )
{
return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
}
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 0b399825..9c2728f8 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -53,7 +53,7 @@ void Setup_NT_Comm( reax_system* system, control_params* control,
bndry_cut = system->bndry_cuts.ghost_cutoff;
/* identify my neighbors */
- system->num_nt_nbrs = MAX_NT_NBRS;
+ system->num_nt_nbrs = REAX_MAX_NT_NBRS;
for ( i = 0; i < system->num_nt_nbrs; ++i )
{
ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 30a169bc..dadb988e 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -55,6 +55,13 @@
interaction_function Interaction_Functions[NUM_INTRS];
+
+static int compare_bonds( const void *p1, const void *p2 )
+{
+ return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
+}
+
+
void Dummy_Interaction( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control )
@@ -531,6 +538,11 @@ void Init_Forces( reax_system *system, control_params *control,
if ( far_nbrs->format == FULL_LIST )
{
+
+ for( i = 0; i < system->N; ++i )
+ qsort( &bonds->bond_list[Start_Index(i, bonds)],
+ Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
+
/* set sym_index for bonds list (far_nbrs full list) */
for ( i = 0; i < system->N; ++i )
{
diff --git a/PuReMD/src/torsion_angles.c b/PuReMD/src/torsion_angles.c
index 8d9949da..4719aa33 100644
--- a/PuReMD/src/torsion_angles.c
+++ b/PuReMD/src/torsion_angles.c
@@ -342,7 +342,8 @@ void Torsion_Angles( reax_system *system, control_params *control,
fbp->V2 * exp_tor1 * (1.0 - cos2omega) +
fbp->V3 * (1.0 + cos3omega) );
- data->my_en.e_tor += e_tor = fn10 * sin_ijk * sin_jkl * CV;
+ e_tor = fn10 * sin_ijk * sin_jkl * CV;
+ data->my_en.e_tor += e_tor;
dfn11 = (-p_tor3 * exp_tor3_DjDk +
(p_tor3 * exp_tor3_DjDk - p_tor4 * exp_tor4_DjDk) *
@@ -375,9 +376,8 @@ void Torsion_Angles( reax_system *system, control_params *control,
/* 4-body conjugation energy */
fn12 = exp_cot2_ij * exp_cot2_jk * exp_cot2_kl;
- data->my_en.e_con += e_con =
- fbp->p_cot1 * fn12 *
- (1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);
+ e_con = fbp->p_cot1 * fn12 * (1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);
+ data->my_en.e_con += e_con;
Cconj = -2.0 * fn12 * fbp->p_cot1 * p_cot2 *
(1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);
--
GitLab