From 18df4177f56e2834e64f8ff651760c21f04caa7c Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 9 Jan 2019 15:52:45 -0500
Subject: [PATCH] PuReMD-old: clean-up and comment initialization
optimizations. WIP: remove redundant computation of bonds list (not working
for full far nbrs list). Improve hbond list computation for full far nbr list
case.
---
PuReMD/src/allocate.c | 2 -
PuReMD/src/bond_orders.c | 166 +++-
PuReMD/src/bond_orders.h | 22 +-
PuReMD/src/forces.c | 1614 +++++++++++++++++++-------------------
PuReMD/src/reax_types.h | 16 +-
5 files changed, 1008 insertions(+), 812 deletions(-)
diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index 65474f1e..32b68ec1 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -145,7 +145,6 @@ void DeAllocate_Workspace( control_params *control, storage *workspace )
sfree( workspace->Clp, "Clp" );
sfree( workspace->vlpex, "vlpex" );
sfree( workspace->bond_mark, "bond_mark" );
- sfree( workspace->done_after, "done_after" );
/* CM storage */
sfree( workspace->Hdia_inv, "Hdia_inv" );
@@ -304,7 +303,6 @@ int Allocate_Workspace( reax_system *system, control_params *control,
workspace->Clp = smalloc( total_real, "Clp", comm );
workspace->vlpex = smalloc( total_real, "vlpex", comm );
workspace->bond_mark = scalloc( total_cap, sizeof(int), "bond_mark", comm );
- workspace->done_after = scalloc( total_cap, sizeof(int), "done_after", comm );
/* CM storage */
workspace->Hdia_inv = scalloc( total_cap, sizeof(real), "Hdia_inv", comm );
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index f01993b8..0c5950e5 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -662,6 +662,10 @@ void Add_dBond_to_Forces( int i, int pj,
}
+/* Compute the bond order term between atoms i and j,
+ * and if this term exceeds the cutoff bo_cut, then adds
+ * BOTH atoms the bonds list (i.e., compute term once
+ * and copy to avoid redundant computation) */
int BOp( storage *workspace, reax_list *bonds, real bo_cut,
int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec,
int far_nbr_list_format, single_body_parameters *sbp_i,
@@ -683,21 +687,169 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + bo_cut) * exp( C12 );
}
- else BO_s = C12 = 0.0;
+ else
+ {
+ C12 = 0.0;
+ BO_s = 0.0;
+ }
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
{
C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 );
BO_pi = exp( C34 );
}
- else BO_pi = C34 = 0.0;
+ else
+ {
+ C34 = 0.0;
+ BO_pi = 0.0;
+ }
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
{
C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = exp( C56 );
}
- else BO_pi2 = C56 = 0.0;
+ else
+ {
+ C56 = 0.0;
+ BO_pi2 = 0.0;
+ }
+
+ /* Initially BO values are the uncorrected ones, page 1 */
+ BO = BO_s + BO_pi + BO_pi2;
+
+ if ( BO >= bo_cut )
+ {
+ /****** bonds i-j and j-i ******/
+ ibond = &bonds->bond_list[btop_i];
+ btop_j = End_Index( j, bonds );
+ jbond = &bonds->bond_list[btop_j];
+
+ ibond->nbr = j;
+ ibond->d = d;
+ rvec_Copy( ibond->dvec, *dvec );
+ ivec_Copy( ibond->rel_box, *rel_box );
+ ibond->dbond_index = btop_i;
+ ibond->sym_index = btop_j;
+ jbond->nbr = i;
+ jbond->d = d;
+ rvec_Scale( jbond->dvec, -1.0, *dvec );
+ ivec_Scale( jbond->rel_box, -1.0, *rel_box );
+ jbond->dbond_index = btop_i;
+ jbond->sym_index = btop_i;
+
+ Set_End_Index( j, btop_j + 1, bonds );
+
+ bo_ij = &ibond->bo_data;
+ bo_ij->BO = BO;
+ bo_ij->BO_s = BO_s;
+ bo_ij->BO_pi = BO_pi;
+ bo_ij->BO_pi2 = BO_pi2;
+ bo_ji = &jbond->bo_data;
+ bo_ji->BO = BO;
+ bo_ji->BO_s = BO_s;
+ bo_ji->BO_pi = BO_pi;
+ bo_ji->BO_pi2 = BO_pi2;
+
+ /* Bond Order page2-3, derivative of total bond order prime */
+ Cln_BOp_s = twbp->p_bo2 * C12 / r2;
+ Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
+ Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
+
+ /* Only dln_BOp_xx wrt. dr_i is stored here, note that
+ * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
+ rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
+ rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
+ rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
+ rvec_Scale( bo_ji->dln_BOp_s, -1.0, bo_ij->dln_BOp_s );
+ rvec_Scale( bo_ji->dln_BOp_pi, -1.0, bo_ij->dln_BOp_pi );
+ rvec_Scale( bo_ji->dln_BOp_pi2, -1.0, bo_ij->dln_BOp_pi2 );
+
+ /* Only dBOp wrt. dr_i is stored here, note that
+ * dBOp/dr_i = -dBOp/dr_j and all others are 0 */
+ rvec_Scale( bo_ij->dBOp, -1.0 * (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.0, bo_ij->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;
+ workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
+ bo_ij->Cdbo = 0.0;
+ bo_ij->Cdbopi = 0.0;
+ bo_ij->Cdbopi2 = 0.0;
+ bo_ji->BO_s -= bo_cut;
+ bo_ji->BO -= bo_cut;
+ workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+ bo_ji->Cdbo = 0.0;
+ bo_ji->Cdbopi = 0.0;
+ bo_ji->Cdbopi2 = 0.0;
+
+ return 1;
+ }
+
+ return 0;
+}
+
+
+/* Compute the bond order term between atoms i and j,
+ * and if this term exceeds the cutoff bo_cut, then adds
+ * to the bond list according to the following convention:
+ * * if the far neighbor list is store in half format,
+ * add BOTH atoms to each other's portion of the bond list
+ * * if the far neighbor list is store in full format,
+ * add atom i to atom j's bonds list ONLY */
+int BOp_redundant( storage *workspace, reax_list *bonds, real bo_cut,
+ int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec,
+ int far_nbr_list_format, single_body_parameters *sbp_i,
+ single_body_parameters *sbp_j, two_body_parameters *twbp )
+{
+ real r2, C12, C34, C56;
+ real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
+ real BO, BO_s, BO_pi, BO_pi2;
+ bond_data *ibond;
+ bond_order_data *bo_ij;
+ int btop_j;
+ bond_data *jbond;
+ bond_order_data *bo_ji;
+
+ r2 = SQR(d);
+
+ if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 )
+ {
+ C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 );
+ BO_s = (1.0 + bo_cut) * exp( C12 );
+ }
+ else
+ {
+ C12 = 0.0;
+ BO_s = 0.0;
+ }
+
+ if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
+ {
+ C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 );
+ BO_pi = exp( C34 );
+ }
+ else
+ {
+ C34 = 0.0;
+ BO_pi = 0.0;
+ }
+
+ if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
+ {
+ C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 );
+ BO_pi2 = exp( C56 );
+ }
+ else
+ {
+ C56 = 0.0;
+ BO_pi2 = 0.0;
+ }
/* Initially BO values are the uncorrected ones, page 1 */
BO = BO_s + BO_pi + BO_pi2;
@@ -730,14 +882,14 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
Set_End_Index( j, btop_j + 1, bonds );
}
- bo_ij = &( ibond->bo_data );
+ bo_ij = &ibond->bo_data;
bo_ij->BO = BO;
bo_ij->BO_s = BO_s;
bo_ij->BO_pi = BO_pi;
bo_ij->BO_pi2 = BO_pi2;
if ( far_nbr_list_format == HALF_LIST )
{
- bo_ji = &( jbond->bo_data );
+ bo_ji = &jbond->bo_data;
bo_ji->BO = BO;
bo_ji->BO_s = BO_s;
bo_ji->BO_pi = BO_pi;
@@ -750,7 +902,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
/* Only dln_BOp_xx wrt. dr_i is stored here, note that
- dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
+ * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
@@ -762,7 +914,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
}
/* Only dBOp wrt. dr_i is stored here, note that
- dBOp/dr_i = -dBOp/dr_j and all others are 0 */
+ * dBOp/dr_i = -dBOp/dr_j and all others are 0 */
rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s
+ bo_ij->BO_pi * Cln_BOp_pi
+ bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
diff --git a/PuReMD/src/bond_orders.h b/PuReMD/src/bond_orders.h
index fe142fb0..fcf4d71a 100644
--- a/PuReMD/src/bond_orders.h
+++ b/PuReMD/src/bond_orders.h
@@ -24,6 +24,7 @@
#include "reax_types.h"
+
typedef struct
{
real C1dbo, C2dbo, C3dbo;
@@ -32,29 +33,42 @@ typedef struct
real C1dDelta, C2dDelta, C3dDelta;
} dbond_coefficients;
+
#ifdef TEST_FORCES
void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+
void Get_dBOpinpi2( reax_system*, reax_list**,
- int, int, real, real, rvec*, rvec* );
+ int, int, real, real, rvec*, rvec* );
void Add_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+
void Add_dBOpinpi2( reax_system*, reax_list**,
- int, int, real, real, rvec*, rvec* );
+ int, int, real, real, rvec*, rvec* );
void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, real );
+
void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**,
- int, int, real, real );
+ int, int, real, real );
void Add_dDelta( reax_system*, reax_list**, int, real, rvec* );
+
void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
#endif
void Add_dBond_to_Forces( int, int, storage*, reax_list** );
+
void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
- storage*, reax_list** );
+ storage*, reax_list** );
+
int BOp( storage*, reax_list*, real, int, int, int, ivec*, real, rvec*,
int, single_body_parameters*, single_body_parameters*,
two_body_parameters* );
+
+int BOp_redundant( storage*, reax_list*, real, int, int, int, ivec*, real, rvec*,
+ int, single_body_parameters*, single_body_parameters*,
+ two_body_parameters* );
+
void BO( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls* );
+
#endif
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 7a5f4961..4588515d 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -20,39 +20,41 @@
----------------------------------------------------------------------*/
#include "reax_types.h"
+
#if defined(PURE_REAX)
-#include "forces.h"
-#include "bond_orders.h"
-#include "bonds.h"
-#include "basic_comm.h"
-#include "hydrogen_bonds.h"
-#include "io_tools.h"
-#include "list.h"
-#include "lookup.h"
-#include "multi_body.h"
-#include "nonbonded.h"
-#include "qEq.h"
-#include "tool_box.h"
-#include "torsion_angles.h"
-#include "valence_angles.h"
-#include "vector.h"
+ #include "forces.h"
+ #include "bond_orders.h"
+ #include "bonds.h"
+ #include "basic_comm.h"
+ #include "hydrogen_bonds.h"
+ #include "io_tools.h"
+ #include "list.h"
+ #include "lookup.h"
+ #include "multi_body.h"
+ #include "nonbonded.h"
+ #include "qEq.h"
+ #include "tool_box.h"
+ #include "torsion_angles.h"
+ #include "valence_angles.h"
+ #include "vector.h"
#elif defined(LAMMPS_REAX)
-#include "reax_forces.h"
-#include "reax_bond_orders.h"
-#include "reax_bonds.h"
-#include "reax_basic_comm.h"
-#include "reax_hydrogen_bonds.h"
-#include "reax_io_tools.h"
-#include "reax_list.h"
-#include "reax_lookup.h"
-#include "reax_multi_body.h"
-#include "reax_nonbonded.h"
-#include "reax_tool_box.h"
-#include "reax_torsion_angles.h"
-#include "reax_valence_angles.h"
-#include "reax_vector.h"
+ #include "reax_forces.h"
+ #include "reax_bond_orders.h"
+ #include "reax_bonds.h"
+ #include "reax_basic_comm.h"
+ #include "reax_hydrogen_bonds.h"
+ #include "reax_io_tools.h"
+ #include "reax_list.h"
+ #include "reax_lookup.h"
+ #include "reax_multi_body.h"
+ #include "reax_nonbonded.h"
+ #include "reax_tool_box.h"
+ #include "reax_torsion_angles.h"
+ #include "reax_valence_angles.h"
+ #include "reax_vector.h"
#endif
+
interaction_function Interaction_Functions[NUM_INTRS];
@@ -70,150 +72,13 @@ static void Dummy_Interaction( reax_system *system, control_params *control,
}
-void Init_Force_Functions( control_params *control )
-{
- Interaction_Functions[0] = &BO;
- Interaction_Functions[1] = &Bonds; //Dummy_Interaction;
- Interaction_Functions[2] = &Atom_Energy; //Dummy_Interaction;
- Interaction_Functions[3] = &Valence_Angles; //Dummy_Interaction;
- Interaction_Functions[4] = &Torsion_Angles; //Dummy_Interaction;
- if ( control->hbond_cut > 0 )
- {
- Interaction_Functions[5] = &Hydrogen_Bonds;
- }
- else
- {
- Interaction_Functions[5] = &Dummy_Interaction;
- }
- Interaction_Functions[6] = &Dummy_Interaction; //empty
- Interaction_Functions[7] = &Dummy_Interaction; //empty
- Interaction_Functions[8] = &Dummy_Interaction; //empty
- Interaction_Functions[9] = &Dummy_Interaction; //empty
-}
-
-
-void Compute_Bonded_Forces( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace,
- reax_list **lists, output_controls *out_control,
- MPI_Comm comm )
-{
- int i;
-
- /* Mark beginning of a new timestep in bonded energy files */
-#if defined(TEST_ENERGY)
- Debug_Marker_Bonded( out_control, data->step );
-#endif
-
- /* Implement all force calls as function pointers */
- for ( i = 0; i < NUM_INTRS; i++ )
- {
-#if defined(DEBUG)
- fprintf( stderr, "p%d: starting f%d\n", system->my_rank, i );
- MPI_Barrier( comm );
-#endif
- (Interaction_Functions[i])( system, control, data, workspace,
- lists, out_control );
-#if defined(DEBUG)
- fprintf( stderr, "p%d: f%d done\n", system->my_rank, i );
- MPI_Barrier( comm );
-#endif
- }
-}
-
-
-void Compute_NonBonded_Forces( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace,
- reax_list **lists, output_controls *out_control,
- MPI_Comm comm )
-{
- /* Mark beginning of a new timestep in nonbonded energy files */
-#if defined(TEST_ENERGY)
- Debug_Marker_Nonbonded( out_control, data->step );
-#endif
-
- /* van der Waals and Coulomb interactions */
- if ( control->tabulate == 0 )
- vdW_Coulomb_Energy( system, control, data, workspace,
- lists, out_control );
- else
- Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
- lists, out_control );
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d: nonbonded forces done\n", system->my_rank );
- MPI_Barrier( comm );
-#endif
-}
-
-
-/* this version of Compute_Total_Force computes forces from
- * coefficients accumulated by all interaction functions.
- * Saves enormous time & space! */
-void Compute_Total_Force( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace,
- reax_list **lists, mpi_datatypes *mpi_data )
-{
- int i, pj;
- reax_list *bonds = lists[BONDS];
-
- for ( i = 0; i < system->N; ++i )
- {
- for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
- {
- if ( i < bonds->bond_list[pj].nbr )
- {
- if ( control->virial == 0 )
- {
- Add_dBond_to_Forces( i, pj, workspace, lists );
- }
- else
- {
- Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
- }
- }
- }
- }
-
- //Print_Total_Force( system, data, workspace );
-
-#if defined(PURE_REAX)
- /* now all forces are computed to their partially-final values
- * based on the neighbors information each processor has had.
- * final values of force on each atom needs to be computed by adding up
- * all partially-final pieces */
- Coll( system, mpi_data, workspace->f, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
-
- for ( i = 0; i < system->n; ++i )
- {
- rvec_Copy( system->my_atoms[i].f, workspace->f[i] );
- }
-
-#if defined(TEST_FORCES)
- Coll( system, mpi_data, workspace->f_ele, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_vdw, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_be, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_lp, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_ov, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_un, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_ang, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_coa, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_pen, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_hb, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_tor, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
- Coll( system, mpi_data, workspace->f_con, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
-#endif
-
-#endif
-}
-
-
-static void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
- int step, int n, int N, int numH, MPI_Comm comm )
+static void Validate_Lists( reax_system *system, storage *workspace,
+ reax_list **lists, int step, int n, int N, int numH, MPI_Comm comm )
{
int i, comp, Hindex;
reax_list *bonds, *hbonds;
reallocate_data *realloc;
- realloc = &(workspace->realloc);
+ realloc = &workspace->realloc;
/* bond list */
if ( N > 0 )
@@ -229,19 +94,23 @@ static void Validate_Lists( reax_system *system, storage *workspace, reax_list *
//workspace->realloc.bonds = 1;
if ( i < N - 1 )
+ {
comp = Start_Index(i + 1, bonds);
- else comp = bonds->num_intrs;
+ }
+ else
+ {
+ comp = bonds->num_intrs;
+ }
if ( End_Index(i, bonds) > comp )
{
- fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
+ fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
step, i, End_Index(i, bonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
}
}
-
/* hbonds list */
if ( numH > 0 )
{
@@ -250,6 +119,7 @@ static void Validate_Lists( reax_system *system, storage *workspace, reax_list *
for ( i = 0; i < n; ++i )
{
Hindex = system->my_atoms[i].Hindex;
+
if ( Hindex > -1 )
{
system->my_atoms[i].num_hbonds =
@@ -260,48 +130,54 @@ static void Validate_Lists( reax_system *system, storage *workspace, reax_list *
// workspace->realloc.hbonds = 1;
if ( Hindex < numH - 1 )
- comp = Start_Index(Hindex + 1, hbonds);
- else comp = hbonds->num_intrs;
+ {
+ comp = Start_Index( Hindex + 1, hbonds );
+ }
+ else
+ {
+ comp = hbonds->num_intrs;
+ }
if ( End_Index(Hindex, hbonds) > comp )
{
- fprintf(stderr, "step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
+ fprintf(stderr, "[ERROR] step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
step, Hindex, End_Index(Hindex, hbonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
}
-/*
- if ( Hindex > -1 )
- {
- system->my_atoms[i].num_hbonds =
- MAX( Num_Entries(Hindex, hbonds) * SAFER_ZONE, MIN_HBONDS );
-*/
+
+// if ( Hindex > -1 )
+// {
+// system->my_atoms[i].num_hbonds =
+// MAX( Num_Entries(Hindex, hbonds) * SAFER_ZONE, MIN_HBONDS );
+
//if( Num_Entries(i, hbonds) >=
//(Start_Index(i+1,hbonds)-Start_Index(i,hbonds))*0.90/*DANGER_ZONE*/){
// workspace->realloc.hbonds = 1;
-/*
+
//TODO
- if ( Hindex < system->n - 1 )
- comp = Start_Index(Hindex + 1, hbonds);
- else comp = hbonds->num_intrs;
-
- if ( End_Index(Hindex, hbonds) > comp )
- {
- fprintf(stderr, "step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
- step, Hindex, End_Index(Hindex, hbonds), comp );
- MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
- }
- }
-
-*/
-
-
-
+// if ( Hindex < system->n - 1 )
+// {
+// comp = Start_Index(Hindex + 1, hbonds);
+// }
+// else
+// {
+// comp = hbonds->num_intrs;
+// }
+//
+// if ( End_Index(Hindex, hbonds) > comp )
+// {
+// fprintf(stderr, "[ERROR] step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
+// step, Hindex, End_Index(Hindex, hbonds), comp );
+// MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
+// }
+// }
}
}
}
+/* Computes a charge matrix entry using the Taper function */
static real Compute_H( real r, real gamma, real *ctap )
{
real taper, dr3gamij_1, dr3gamij_3;
@@ -315,24 +191,32 @@ static real Compute_H( real r, real gamma, real *ctap )
taper = taper * r + ctap[0];
dr3gamij_1 = ( r * r * r + gamma );
- dr3gamij_3 = pow( dr3gamij_1 , 0.33333333333333 );
+ dr3gamij_3 = pow( dr3gamij_1, 1.0 / 3.0 );
+
return taper * EV_to_KCALpMOL / dr3gamij_3;
}
+/* Computes a charge matrix entry using the force tabulation
+ * (i.e., an arithmetic-reducing optimization) */
static real Compute_tabH( real r_ij, int ti, int tj )
{
int r, tmin, tmax;
real val, dif, base;
LR_lookup_table *t;
- tmin = MIN( ti, tj );
- tmax = MAX( ti, tj );
- t = &( LR[tmin][tmax] );
+ tmin = MIN( ti, tj );
+ tmax = MAX( ti, tj );
+ t = &LR[tmin][tmax];
/* cubic spline interpolation */
r = (int)(r_ij * t->inv_dx);
- if ( r == 0 ) ++r;
+
+ if ( r == 0 )
+ {
+ ++r;
+ }
+
base = (real)(r + 1) * t->dx;
dif = r_ij - base;
val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
@@ -343,6 +227,8 @@ static real Compute_tabH( real r_ij, int ti, int tj )
}
+/* Compute the distances and displacement vectors for entries
+ * in the far neighbors list if it's a NOT re-neighboring step */
static void Init_Distance( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -364,7 +250,7 @@ static void Init_Distance( reax_system *system, control_params *control,
start_i = Start_Index( i, far_nbrs );
end_i = End_Index( i, far_nbrs );
- /* update i-j distance for non-reneighboring steps */
+ /* update distance and displacement vector between atoms i and j (i-j) */
for ( pj = start_i; pj < end_i; ++pj )
{
j = far_nbrs->far_nbr_list.nbr[pj];
@@ -378,12 +264,13 @@ static void Init_Distance( reax_system *system, control_params *control,
}
}
}
-
-
}
#if defined(NEUTRAL_TERRITORY)
+/* Compute the charge matrix entries and store the matrix in half format
+ * using the far neighbors list (stored in full format) and according to
+ * the neutral territory communication method */
static void Init_CM_Half_NT( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -597,6 +484,9 @@ static void Init_CM_Half_NT( reax_system *system, control_params *control,
}
+/* Compute the charge matrix entries and store the matrix in full format
+ * using the far neighbors list (stored in full format) and according to
+ * the neutral territory communication method */
static void Init_CM_Full_NT( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -810,6 +700,9 @@ static void Init_CM_Full_NT( reax_system *system, control_params *control,
#else
+/* Compute the charge matrix entries and store the matrix in half format
+ * using the far neighbors list (stored in half format) and according to
+ * the full shell communication method */
static void Init_CM_Half_FS( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -894,6 +787,9 @@ static void Init_CM_Half_FS( reax_system *system, control_params *control,
}
+/* Compute the charge matrix entries and store the matrix in full format
+ * using the far neighbors list (stored in full format) and according to
+ * the full shell communication method */
static void Init_CM_Full_FS( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -975,6 +871,11 @@ static void Init_CM_Full_FS( reax_system *system, control_params *control,
#endif
+/* Compute entries of the bonds/hbonds lists and store the lists in full format
+ * using the far neighbors list (stored in half format)
+ *
+ * Note: this version does NOT contain an optimization to restrict the bond_mark
+ * array to at most the 3-hop neighborhood */
static void Init_Bond_Half( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -1002,8 +903,8 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
}
for ( i = system->n; i < system->N; ++i )
{
- workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
- //workspace->done_after[i] = Start_Index( i, far_nbrs );
+ /* put ghost atoms to an infinite distance (i.e., 1000) */
+ workspace->bond_mark[i] = 1000;
}
num_bonds = 0;
@@ -1098,8 +999,7 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
}
/* uncorrected bond orders */
- if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- BOp( workspace, bonds, control->bo_cut,
+ if ( BOp( workspace, bonds, control->bo_cut,
i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
&far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
&far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
@@ -1149,6 +1049,8 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
}
+/* Compute entries of the bonds/hbonds lists and store the lists in full format
+ * using the far neighbors list (stored in full format) */
static void Init_Bond_Full( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
@@ -1158,7 +1060,6 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
int type_i, type_j;
int num_bonds, num_hbonds;
int ihb, jhb, ihb_top;
- int local;
real cutoff;
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
@@ -1172,95 +1073,52 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
hbonds = lists[HBONDS];
-
num_hbonds = 0;
- for ( i = 0; i < system->N; ++i )
+ /* hbonds */
+ if ( control->hbond_cut > 0.0 )
{
- atom_i = &system->my_atoms[i];
- type_i = atom_i->type;
- start_i = Start_Index( i, far_nbrs );
- end_i = End_Index( i, far_nbrs );
-
- sbp_i = &system->reax_param.sbp[type_i];
-
- if ( i < system->n )
- {
- local = 1;
- cutoff = control->nonb_cut;
- }
- else
- {
- local = 0;
- cutoff = control->bond_cut;
- }
-
- ihb = -1;
- ihb_top = -1;
- if ( local == 1 )
+ for ( i = 0; i < system->n; ++i )
{
- if ( control->hbond_cut > 0 )
- {
- ihb = sbp_i->p_hbond;
-
- if ( ihb == 1 )
- {
- ihb_top = Start_Index( atom_i->Hindex, hbonds );
- }
- else
- {
- ihb_top = -1;
- }
- }
- }
+ atom_i = &system->my_atoms[i];
+ type_i = atom_i->type;
+ sbp_i = &system->reax_param.sbp[type_i];
+ ihb = sbp_i->p_hbond;
- /* update i-j distance - check if j is within cutoff */
- for ( pj = start_i; pj < end_i; ++pj )
- {
- j = far_nbrs->far_nbr_list.nbr[pj];
- atom_j = &system->my_atoms[j];
-
- if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
+ if ( ihb == 1 )
{
- type_j = atom_j->type;
- sbp_j = &system->reax_param.sbp[type_j];
+ start_i = Start_Index( i, far_nbrs );
+ end_i = End_Index( i, far_nbrs );
+ ihb_top = Start_Index( atom_i->Hindex, hbonds );
- if ( local == 1 )
+ /* check if j is within cutoff */
+ for ( pj = start_i; pj < end_i; ++pj )
{
- /* hydrogen bond lists */
- if ( control->hbond_cut > 0
- && (ihb == 1 || ihb == 2)
- && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
+ j = far_nbrs->far_nbr_list.nbr[pj];
+ atom_j = &system->my_atoms[j];
+
+ if ( far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut
+ && system->reax_param.sbp[atom_j->type].p_hbond == 2 )
{
- // fprintf( stderr, "%d %d\n", atom1, atom2 );
- jhb = sbp_j->p_hbond;
-
- if ( ihb == 1 && jhb == 2 )
- {
- hbonds->hbond_list[ihb_top].nbr = j;
- hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = pj;
- ++ihb_top;
- ++num_hbonds;
- }
+ hbonds->hbond_list[ihb_top].nbr = j;
+ hbonds->hbond_list[ihb_top].scl = 1;
+ hbonds->hbond_list[ihb_top].ptr = pj;
+ ++ihb_top;
+ ++num_hbonds;
}
}
+ Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
}
}
-
- if ( local == 1 && ihb == 1 )
- {
- Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
- }
}
- //Computing bond_mark
+ /* compute bond_mark */
push = 0;
num_bonds = 0;
btop_i = 0;
bonds = lists[BONDS];
- q = smalloc( sizeof( int ) * (system->N - system->n),
+ q = smalloc( sizeof(int) * (system->N - system->n),
"Init_Distance::q", MPI_COMM_WORLD );
for ( i = 0; i < system->n; ++i )
@@ -1269,11 +1127,11 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
}
for ( i = system->n; i < system->N; ++i )
{
- workspace->bond_mark[i] = 1000;// put ghost atoms to an infinite distance
- //workspace->done_after[i] = Start_Index( i, far_nbrs );
+ /* put ghost atoms to an infinite distance (i.e., 1000) */
+ workspace->bond_mark[i] = 1000;
}
- //bonds that are directly connected to local atoms
+ /* bonds that are directly connected to local atoms */
for ( i = 0; i < system->n; ++i )
{
atom_i = &system->my_atoms[i];
@@ -1287,21 +1145,30 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
{
j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
+
if ( far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
+// if ( i <= j && far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
{
type_j = atom_j->type;
sbp_j = &system->reax_param.sbp[type_j];
twbp = &system->reax_param.tbp[type_i][type_j];
- if ( BOp( workspace, bonds, control->bo_cut,
+ if ( BOp_redundant( workspace, bonds, control->bo_cut,
i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
&far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
&far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
+// if ( BOp( workspace, bonds, control->bo_cut,
+// i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+// &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+// &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+// sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
++btop_i;
+ /* if j is a non-local atom, push it on the queue
+ * to search for it's bonded neighbors later */
if ( workspace->bond_mark[j] == 1000 )
{
workspace->bond_mark[j] = 1;
@@ -1314,7 +1181,7 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
Set_End_Index( i, btop_i, bonds );
}
- //bonds that are indirectly connected to local atoms
+ /* bonds that are indirectly connected to local atoms */
for ( k = 0; k < push; ++k )
{
i = q[k];
@@ -1328,31 +1195,41 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
j = far_nbrs->far_nbr_list.nbr[pj];
- if( workspace->bond_mark[i] == 3 && workspace->bond_mark[j] == 1000 )
+
+ if ( workspace->bond_mark[i] == 3
+ && workspace->bond_mark[j] == 1000 )
{
continue;
}
+
atom_j = &system->my_atoms[j];
+
if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
+// if ( i <= j && far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
{
type_j = atom_j->type;
sbp_j = &system->reax_param.sbp[type_j];
twbp = &system->reax_param.tbp[type_i][type_j];
-
- if ( BOp( workspace, bonds, control->bo_cut,
+ if ( BOp_redundant( workspace, bonds, control->bo_cut,
i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
&far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
&far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
+// if ( BOp( workspace, bonds, control->bo_cut,
+// i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+// &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+// &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+// sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
++btop_i;
if ( workspace->bond_mark[j] == 1000 )
{
- workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
- if( workspace->bond_mark[j] < 4 )
+ workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
+
+ if ( workspace->bond_mark[j] < 4 )
{
q[ push++ ] = j;
}
@@ -1369,46 +1246,47 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
//neighborhood relationship with any local atoms
//those atoms have their bond_mark value still
//equal to 1000
- /*for ( i = system->n; i < system->N; ++i )
- {
- if ( workspace->bond_mark[i] == 1000 )
- {
- atom_i = &system->my_atoms[i];
- type_i = atom_i->type;
- btop_i = End_Index( i, bonds );
- sbp_i = &system->reax_param.sbp[type_i];
- start_i = Start_Index( i, far_nbrs );
- end_i = End_Index( i, far_nbrs );
+// for ( i = system->n; i < system->N; ++i )
+// {
+// if ( workspace->bond_mark[i] == 1000 )
+// {
+// atom_i = &system->my_atoms[i];
+// type_i = atom_i->type;
+// btop_i = End_Index( i, bonds );
+// sbp_i = &system->reax_param.sbp[type_i];
+// start_i = Start_Index( i, far_nbrs );
+// end_i = End_Index( i, far_nbrs );
+//
+// for ( pj = start_i; pj < end_i; ++pj )
+// {
+// j = far_nbrs->far_nbr_list.nbr[pj];
+// atom_j = &system->my_atoms[j];
+// if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
+// {
+// type_j = atom_j->type;
+// sbp_j = &system->reax_param.sbp[type_j];
+// twbp = &system->reax_param.tbp[type_i][type_j];
+//
+// if ( BOp( workspace, bonds, control->bo_cut,
+// i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+// &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+// &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+// sbp_i, sbp_j, twbp ) )
+// {
+// num_bonds += 2;
+// ++btop_i;
+// }
+// }
+// }
+//
+// Set_End_Index( i, btop_i, bonds );
+// }
+// }
- for ( pj = start_i; pj < end_i; ++pj )
- {
- j = far_nbrs->far_nbr_list.nbr[pj];
- atom_j = &system->my_atoms[j];
- if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
- {
- type_j = atom_j->type;
- sbp_j = &system->reax_param.sbp[type_j];
- twbp = &system->reax_param.tbp[type_i][type_j];
-
- if ( BOp( workspace, bonds, control->bo_cut,
- i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
- &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
- &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
- sbp_i, sbp_j, twbp ) )
- {
- num_bonds += 2;
- ++btop_i;
- }
- }
- }
-
- Set_End_Index( i, btop_i, bonds );
- }
- }*/
workspace->realloc.num_bonds = num_bonds;
sfree( q, "Init_Bond_Full::q" );
- for( i = 0; i < system->N; ++i )
+ for ( i = 0; i < system->N; ++i )
{
qsort( &bonds->bond_list[Start_Index(i, bonds)],
Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
@@ -1456,490 +1334,633 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
}
-void Init_Forces( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace, reax_list **lists,
- output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
+void Init_Force_Functions( control_params *control )
{
- double t_start, t_dist, t_cm, t_bond;
- double timings[3], t_total[3];
-
- t_start = MPI_Wtime( );
-
- Init_Distance( system, control, data, workspace, lists, out_control, comm, mpi_data );
-
- t_dist = MPI_Wtime( );
-
-#if defined(NEUTRAL_TERRITORY)
- if ( workspace->H->format == SYM_HALF_MATRIX )
- {
- Init_CM_Half_NT( system, control, data, workspace, lists, out_control, comm, mpi_data );
- }
- else
- {
- Init_CM_Full_NT( system, control, data, workspace, lists, out_control, comm, mpi_data );
- }
-#else
- if ( workspace->H->format == SYM_HALF_MATRIX )
- {
- Init_CM_Half_FS( system, control, data, workspace, lists, out_control, comm, mpi_data );
- }
- else
- {
- Init_CM_Full_FS( system, control, data, workspace, lists, out_control, comm, mpi_data );
- }
-#endif
-
- t_cm = MPI_Wtime();
-
- if ( lists[FAR_NBRS]->format == HALF_LIST )
+ Interaction_Functions[0] = &BO;
+ Interaction_Functions[1] = &Bonds; //Dummy_Interaction;
+ Interaction_Functions[2] = &Atom_Energy; //Dummy_Interaction;
+ Interaction_Functions[3] = &Valence_Angles; //Dummy_Interaction;
+ Interaction_Functions[4] = &Torsion_Angles; //Dummy_Interaction;
+ if ( control->hbond_cut > 0.0 )
{
- Init_Bond_Half( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ Interaction_Functions[5] = &Hydrogen_Bonds;
}
else
{
- Init_Bond_Full( system, control, data, workspace, lists, out_control, comm, mpi_data );
- }
-
- t_bond = MPI_Wtime();
-
- timings[0] = t_dist - t_start;
- timings[1] = t_cm - t_dist;
- timings[2] = t_bond - t_cm;
-
- MPI_Reduce( timings, t_total, 3, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );
-
- if ( system->my_rank == MASTER_NODE )
- {
- data->timing.init_dist += t_total[0] / control->nprocs;
- data->timing.init_cm += t_total[1] / control->nprocs;
- data->timing.init_bond += t_total[2] / control->nprocs;
- }
-
-}
-
-
-/*void Init_Forces( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace, reax_list **lists,
- output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
-{
- int i, j, pj;
- int start_i, end_i;
- int type_i, type_j;
- int Htop, btop_i, num_bonds, num_hbonds;
- int ihb, jhb, ihb_top;
- int local, flag, renbr;
- real r_ij, cutoff;
- sparse_matrix *H;
- reax_list *far_nbrs, *bonds, *hbonds;
- single_body_parameters *sbp_i, *sbp_j;
- two_body_parameters *twbp;
- reax_atom *atom_i, *atom_j;
- int jhb_top;
- int start_j, end_j;
- int btop_j;
-#if defined(NEUTRAL_TERRITORY)
- int mark[6];
- int total_cnt[6];
- int bin[6];
- int total_sum[6];
- int nt_flag;
-#endif
-
- far_nbrs = lists[FAR_NBRS];
- bonds = lists[BONDS];
- hbonds = lists[HBONDS];
-
-
- for ( i = 0; i < system->n; ++i )
- workspace->bond_mark[i] = 0;
- for ( i = system->n; i < system->N; ++i )
- {
- workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
- //workspace->done_after[i] = Start_Index( i, far_nbrs );
- }
-
- H = workspace->H;
- H->n = system->n;
- Htop = 0;
- num_bonds = 0;
- num_hbonds = 0;
- btop_i = 0;
- renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
-
-#if defined(NEUTRAL_TERRITORY)
- nt_flag = 1;
- if( renbr )
- {
- for ( i = 0; i < 6; ++i )
- {
- total_cnt[i] = 0;
- bin[i] = 0;
- total_sum[i] = 0;
- }
-
- for ( i = system->n; i < system->N; ++i )
- {
- atom_i = &system->my_atoms[i];
-
- if( atom_i->nt_dir != -1 )
- {
- total_cnt[ atom_i->nt_dir ]++;
- }
- }
-
- total_sum[0] = system->n;
- for ( i = 1; i < 6; ++i )
- {
- total_sum[i] = total_sum[i-1] + total_cnt[i-1];
- }
-
- for ( i = system->n; i < system->N; ++i )
- {
- atom_i = &system->my_atoms[i];
-
- if( atom_i->nt_dir != -1 )
- {
- atom_i->pos = total_sum[ atom_i->nt_dir ] + bin[ atom_i->nt_dir ];
- bin[ atom_i->nt_dir ]++;
- }
- }
- H->NT = total_sum[5] + total_cnt[5];
+ Interaction_Functions[5] = &Dummy_Interaction;
}
+ Interaction_Functions[6] = &Dummy_Interaction; //empty
+ Interaction_Functions[7] = &Dummy_Interaction; //empty
+ Interaction_Functions[8] = &Dummy_Interaction; //empty
+ Interaction_Functions[9] = &Dummy_Interaction; //empty
+}
- mark[0] = mark[1] = 1;
- mark[2] = mark[3] = mark[4] = mark[5] = 2;
-#endif
-
- for ( i = 0; i < system->N; ++i )
- {
- atom_i = &system->my_atoms[i];
- type_i = atom_i->type;
- start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
-
- if ( far_nbrs->format == HALF_LIST )
- {
- // start at end because other atoms
- // can add to this atom's list (half-list)
- btop_i = End_Index( i, bonds );
- }
- else if ( far_nbrs->format == FULL_LIST )
- {
- btop_i = Start_Index( i, bonds );
- }
- sbp_i = &system->reax_param.sbp[type_i];
-
- if ( i < system->n )
- {
- local = 1;
- cutoff = control->nonb_cut;
- }
-#if defined(NEUTRAL_TERRITORY)
- else if ( atom_i->nt_dir != -1 )
- {
- local = 2;
- cutoff = control->nonb_cut;
- nt_flag = 0;
- }
-#endif
- else
- {
- local = 0;
- cutoff = control->bond_cut;
- }
-
- ihb = -1;
- ihb_top = -1;
- if ( local == 1 )
- {
- H->start[i] = Htop;
- H->entries[Htop].j = i;
- H->entries[Htop].val = sbp_i->eta;
- ++Htop;
-
- if ( control->hbond_cut > 0 )
- {
- ihb = sbp_i->p_hbond;
- if ( ihb == 1 )
- {
- if ( far_nbrs->format == HALF_LIST )
- {
- // start at end because other atoms
- // can add to this atom's list (half-list)
- ihb_top = End_Index( atom_i->Hindex, hbonds );
- }
- else if ( far_nbrs->format == FULL_LIST )
- {
- ihb_top = Start_Index( atom_i->Hindex, hbonds );
- }
- }
- else
- {
- ihb_top = -1;
- }
- }
- }
-
- // update i-j distance - check if j is within cutoff
- for ( pj = start_i; pj < end_i; ++pj )
- {
- j = far_nbrs->far_nbr_list.nbr[pj];
- atom_j = &system->my_atoms[j];
-
- if ( renbr )
- {
- if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
- flag = 1;
- else
- flag = 0;
- }
- else
- {
- far_nbrs->far_nbr_list.dvec[pj][0] = atom_j->x[0] - atom_i->x[0];
- far_nbrs->far_nbr_list.dvec[pj][1] = atom_j->x[1] - atom_i->x[1];
- far_nbrs->far_nbr_list.dvec[pj][2] = atom_j->x[2] - atom_i->x[2];
- far_nbrs->far_nbr_list.d[pj] = rvec_Norm_Sqr( far_nbrs->far_nbr_list.dvec[pj] );
-
- if ( far_nbrs->far_nbr_list.d[pj] <= SQR(cutoff) )
- {
- far_nbrs->far_nbr_list.d[pj] = sqrt( far_nbrs->far_nbr_list.d[pj] );
- flag = 1;
- }
- else
- {
- flag = 0;
- }
- }
-
- if ( flag )
- {
- type_j = atom_j->type;
- r_ij = far_nbrs->far_nbr_list.d[pj];
- sbp_j = &system->reax_param.sbp[type_j];
- twbp = &system->reax_param.tbp[type_i][type_j];
-
- if ( local == 1 )
- {
- // H matrix entry
-#if defined(NEUTRAL_TERRITORY)
- if ( atom_j->nt_dir > 0 || (j < system->n
- && (H->format == SYM_FULL_MATRIX
- || (H->format == SYM_HALF_MATRIX && i < j))) )
- {
- if( j < system->n )
- {
- H->entries[Htop].j = j;
- }
- else
- {
- H->entries[Htop].j = atom_j->pos;
- }
-
- if ( control->tabulate == 0 )
- {
- H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
- }
- else
- {
- H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
- }
-
- ++Htop;
- }
-#else
- if ( (far_nbrs->format == HALF_LIST
- && (j < system->n || atom_i->orig_id < atom_j->orig_id))
- || far_nbrs->format == FULL_LIST )
- {
- H->entries[Htop].j = j;
-
- if ( control->tabulate == 0 )
- {
- H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
- }
- else
- {
- H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
- }
-
- ++Htop;
- }
-#endif
-
- // hydrogen bond lists
- if ( control->hbond_cut > 0.0
- && (ihb == 1 || ihb == 2)
- && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
- {
- // fprintf( stderr, "%d %d\n", atom1, atom2 );
- jhb = sbp_j->p_hbond;
- if ( ihb == 1 && jhb == 2 )
- {
- hbonds->hbond_list[ihb_top].nbr = j;
- hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = pj;
- ++ihb_top;
- ++num_hbonds;
- }
- // only add to list for local j (far nbrs is half-list)
- else if ( far_nbrs->format == HALF_LIST
- && (j < system->n && ihb == 2 && jhb == 1) )
- {
- jhb_top = End_Index( atom_j->Hindex, hbonds );
- hbonds->hbond_list[jhb_top].nbr = i;
- hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = pj;
- Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
- ++num_hbonds;
- }
- }
- }
-#if defined(NEUTRAL_TERRITORY)
- else if ( local == 2 )
- {
- // H matrix entry
- if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
- && ( H->format == SYM_FULL_MATRIX
- || (H->format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
- || ( j < system->n && atom_i->nt_dir != 0 && H->format == SYM_FULL_MATRIX ))
- {
- if( !nt_flag )
- {
- nt_flag = 1;
- H->start[atom_i->pos] = Htop;
- }
-
- if( j < system->n )
- {
- H->entries[Htop].j = j;
- }
- else
- {
- H->entries[Htop].j = atom_j->pos;
- }
-
- if ( control->tabulate == 0 )
- {
- H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
- }
- else
- {
- H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
- }
- ++Htop;
- }
- }
+void Compute_Bonded_Forces( reax_system *system, control_params *control,
+ simulation_data *data, storage *workspace, reax_list **lists,
+ output_controls *out_control, MPI_Comm comm )
+{
+ int i;
+
+ /* Mark beginning of a new timestep in bonded energy files */
+#if defined(TEST_ENERGY)
+ Debug_Marker_Bonded( out_control, data->step );
#endif
- // uncorrected bond orders
- if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
- && BOp( workspace, bonds, control->bo_cut,
- i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
- &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
- &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
- sbp_i, sbp_j, twbp ) )
- {
- num_bonds += 2;
- ++btop_i;
+ /* Implement all force calls as function pointers */
+ for ( i = 0; i < NUM_INTRS; i++ )
+ {
+#if defined(DEBUG)
+ fprintf( stderr, "p%d: starting f%d\n", system->my_rank, i );
+ MPI_Barrier( comm );
+#endif
- if ( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
- workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
- else if ( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 )
- {
- workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
- }
- }
- }
- }
+ (Interaction_Functions[i])( system, control, data, workspace,
+ lists, out_control );
- Set_End_Index( i, btop_i, bonds );
- if ( local == 1 )
- {
- H->end[i] = Htop;
- if ( ihb == 1 )
- Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
- }
-#if defined(NEUTRAL_TERRITORY)
- else if ( local == 2 )
- {
- if( nt_flag )
- {
- H->end[atom_i->pos] = Htop;
- }
- else
- {
- H->start[atom_i->pos] = 0;
- H->end[atom_i->pos] = 0;
- }
- }
+#if defined(DEBUG)
+ fprintf( stderr, "p%d: f%d done\n", system->my_rank, i );
+ MPI_Barrier( comm );
#endif
}
+}
- if ( far_nbrs->format == FULL_LIST )
+
+void Compute_NonBonded_Forces( reax_system *system, control_params *control,
+ simulation_data *data, storage *workspace, reax_list **lists,
+ output_controls *out_control, MPI_Comm comm )
+{
+ /* Mark beginning of a new timestep in nonbonded energy files */
+#if defined(TEST_ENERGY)
+ Debug_Marker_Nonbonded( out_control, data->step );
+#endif
+
+ /* van der Waals and Coulomb interactions */
+ if ( control->tabulate == 0 )
+ {
+ vdW_Coulomb_Energy( system, control, data, workspace,
+ lists, out_control );
+ }
+ else
{
+ Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
+ lists, out_control );
+ }
- for( i = 0; i < system->N; ++i )
- qsort( &bonds->bond_list[Start_Index(i, bonds)],
- Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
+#if defined(DEBUG)
+ fprintf( stderr, "p%d: nonbonded forces done\n", system->my_rank );
+ MPI_Barrier( comm );
+#endif
+}
- // set sym_index for bonds list (far_nbrs full list)
- for ( i = 0; i < system->N; ++i )
- {
- start_i = Start_Index( i, bonds );
- end_i = End_Index( i, bonds );
- for ( btop_i = start_i; btop_i < end_i; ++btop_i )
- {
- j = bonds->bond_list[btop_i].nbr;
- start_j = Start_Index( j, bonds );
- end_j = End_Index( j, bonds );
+/* this version of Compute_Total_Force computes forces from
+ * coefficients accumulated by all interaction functions.
+ * Saves enormous time & space! */
+void Compute_Total_Force( reax_system *system, control_params *control,
+ simulation_data *data, storage *workspace,
+ reax_list **lists, mpi_datatypes *mpi_data )
+{
+ int i, pj;
+ reax_list *bonds;
- for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+ bonds = lists[BONDS];
+
+ for ( i = 0; i < system->N; ++i )
+ {
+ for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
+ {
+ if ( i < bonds->bond_list[pj].nbr )
+ {
+ if ( control->virial == 0 )
{
- if ( bonds->bond_list[btop_j].nbr == i )
- {
- bonds->bond_list[btop_i].sym_index = btop_j;
- break;
- }
+ Add_dBond_to_Forces( i, pj, workspace, lists );
+ }
+ else
+ {
+ Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
}
}
}
}
-#if defined(DEBUG)
- Print_Sparse_Matrix2( system, H, NULL );
-#endif
+ //Print_Total_Force( system, data, workspace );
- workspace->realloc.Htop = Htop;
- workspace->realloc.num_bonds = num_bonds;
- workspace->realloc.num_hbonds = num_hbonds;
+#if defined(PURE_REAX)
+ /* now all forces are computed to their partially-final values
+ * based on the neighbors information each processor has had.
+ * final values of force on each atom needs to be computed by adding up
+ * all partially-final pieces */
+ Coll( system, mpi_data, workspace->f, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "p%d @ step%d: Htop = %d num_bonds = %d num_hbonds = %d\n",
- system->my_rank, data->step, Htop, num_bonds, num_hbonds );
- MPI_Barrier( comm );
+ for ( i = 0; i < system->n; ++i )
+ {
+ rvec_Copy( system->my_atoms[i].f, workspace->f[i] );
+ }
+
+#if defined(TEST_FORCES)
+ Coll( system, mpi_data, workspace->f_ele, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_vdw, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_be, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_lp, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_ov, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_un, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_ang, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_coa, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_pen, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_hb, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_tor, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
+ Coll( system, mpi_data, workspace->f_con, RVEC_PTR_TYPE, mpi_data->mpi_rvec );
#endif
-#if defined( DEBUG )
- Print_Bonds( system, bonds, "debugbonds.out" );
- Print_Bond_List2( system, bonds, "pbonds.out" );
- Print_Sparse_Matrix( system, H );
- for ( i = 0; i < H->n; ++i )
- for ( j = H->start[i]; j < H->end[i]; ++j )
- fprintf( stderr, "%d %d %.15e\n",
- MIN(system->my_atoms[i].orig_id,
- system->my_atoms[H->entries[j].j].orig_id),
- MAX(system->my_atoms[i].orig_id,
- system->my_atoms[H->entries[j].j].orig_id),
- H->entries[j].val );
#endif
+}
- Validate_Lists( system, workspace, lists, data->step,
- system->n, system->N, system->numH, comm );
-}*/
+void Init_Forces( reax_system *system, control_params *control,
+ simulation_data *data, storage *workspace, reax_list **lists,
+ output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
+{
+ double t_start, t_dist, t_cm, t_bond;
+ double timings[3], t_total[3];
+
+ t_start = MPI_Wtime( );
+
+ Init_Distance( system, control, data, workspace, lists, out_control, comm, mpi_data );
+
+ t_dist = MPI_Wtime( );
+
+#if defined(NEUTRAL_TERRITORY)
+ if ( workspace->H->format == SYM_HALF_MATRIX )
+ {
+ Init_CM_Half_NT( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+ else
+ {
+ Init_CM_Full_NT( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+#else
+ if ( workspace->H->format == SYM_HALF_MATRIX )
+ {
+ Init_CM_Half_FS( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+ else
+ {
+ Init_CM_Full_FS( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+#endif
+
+ t_cm = MPI_Wtime();
+
+ if ( lists[FAR_NBRS]->format == HALF_LIST )
+ {
+ Init_Bond_Half( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+ else
+ {
+ Init_Bond_Full( system, control, data, workspace, lists, out_control, comm, mpi_data );
+ }
+
+ t_bond = MPI_Wtime();
+
+ timings[0] = t_dist - t_start;
+ timings[1] = t_cm - t_dist;
+ timings[2] = t_bond - t_cm;
+
+ MPI_Reduce( timings, t_total, 3, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );
+
+ if ( system->my_rank == MASTER_NODE )
+ {
+ data->timing.init_dist += t_total[0] / control->nprocs;
+ data->timing.init_cm += t_total[1] / control->nprocs;
+ data->timing.init_bond += t_total[2] / control->nprocs;
+ }
+
+}
+
+
+//void Init_Forces( reax_system *system, control_params *control,
+// simulation_data *data, storage *workspace, reax_list **lists,
+// output_controls *out_control, MPI_Comm comm, mpi_datatypes *mpi_data )
+//{
+// int i, j, pj;
+// int start_i, end_i;
+// int type_i, type_j;
+// int Htop, btop_i, num_bonds, num_hbonds;
+// int ihb, jhb, ihb_top;
+// int local, flag, renbr;
+// real r_ij, cutoff;
+// sparse_matrix *H;
+// reax_list *far_nbrs, *bonds, *hbonds;
+// single_body_parameters *sbp_i, *sbp_j;
+// two_body_parameters *twbp;
+// reax_atom *atom_i, *atom_j;
+// int jhb_top;
+// int start_j, end_j;
+// int btop_j;
+//#if defined(NEUTRAL_TERRITORY)
+// int mark[6];
+// int total_cnt[6];
+// int bin[6];
+// int total_sum[6];
+// int nt_flag;
+//#endif
+//
+// far_nbrs = lists[FAR_NBRS];
+// bonds = lists[BONDS];
+// hbonds = lists[HBONDS];
+//
+//
+// for ( i = 0; i < system->n; ++i )
+// workspace->bond_mark[i] = 0;
+// for ( i = system->n; i < system->N; ++i )
+// {
+// /* put ghost atoms to an infinite distance (i.e., 1000) */
+// workspace->bond_mark[i] = 1000;
+// }
+//
+// H = workspace->H;
+// H->n = system->n;
+// Htop = 0;
+// num_bonds = 0;
+// num_hbonds = 0;
+// btop_i = 0;
+// renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
+//
+//#if defined(NEUTRAL_TERRITORY)
+// nt_flag = 1;
+// if( renbr )
+// {
+// for ( i = 0; i < 6; ++i )
+// {
+// total_cnt[i] = 0;
+// bin[i] = 0;
+// total_sum[i] = 0;
+// }
+//
+// for ( i = system->n; i < system->N; ++i )
+// {
+// atom_i = &system->my_atoms[i];
+//
+// if( atom_i->nt_dir != -1 )
+// {
+// total_cnt[ atom_i->nt_dir ]++;
+// }
+// }
+//
+// total_sum[0] = system->n;
+// for ( i = 1; i < 6; ++i )
+// {
+// total_sum[i] = total_sum[i-1] + total_cnt[i-1];
+// }
+//
+// for ( i = system->n; i < system->N; ++i )
+// {
+// atom_i = &system->my_atoms[i];
+//
+// if( atom_i->nt_dir != -1 )
+// {
+// atom_i->pos = total_sum[ atom_i->nt_dir ] + bin[ atom_i->nt_dir ];
+// bin[ atom_i->nt_dir ]++;
+// }
+// }
+// H->NT = total_sum[5] + total_cnt[5];
+// }
+//
+// mark[0] = mark[1] = 1;
+// mark[2] = mark[3] = mark[4] = mark[5] = 2;
+//#endif
+//
+// for ( i = 0; i < system->N; ++i )
+// {
+// atom_i = &system->my_atoms[i];
+// type_i = atom_i->type;
+// start_i = Start_Index(i, far_nbrs);
+// end_i = End_Index(i, far_nbrs);
+//
+// if ( far_nbrs->format == HALF_LIST )
+// {
+// // start at end because other atoms
+// // can add to this atom's list (half-list)
+// btop_i = End_Index( i, bonds );
+// }
+// else if ( far_nbrs->format == FULL_LIST )
+// {
+// btop_i = Start_Index( i, bonds );
+// }
+// sbp_i = &system->reax_param.sbp[type_i];
+//
+// if ( i < system->n )
+// {
+// local = 1;
+// cutoff = control->nonb_cut;
+// }
+//#if defined(NEUTRAL_TERRITORY)
+// else if ( atom_i->nt_dir != -1 )
+// {
+// local = 2;
+// cutoff = control->nonb_cut;
+// nt_flag = 0;
+// }
+//#endif
+// else
+// {
+// local = 0;
+// cutoff = control->bond_cut;
+// }
+//
+// ihb = -1;
+// ihb_top = -1;
+// if ( local == 1 )
+// {
+// H->start[i] = Htop;
+// H->entries[Htop].j = i;
+// H->entries[Htop].val = sbp_i->eta;
+// ++Htop;
+//
+// if ( control->hbond_cut > 0 )
+// {
+// ihb = sbp_i->p_hbond;
+// if ( ihb == 1 )
+// {
+// if ( far_nbrs->format == HALF_LIST )
+// {
+// // start at end because other atoms
+// // can add to this atom's list (half-list)
+// ihb_top = End_Index( atom_i->Hindex, hbonds );
+// }
+// else if ( far_nbrs->format == FULL_LIST )
+// {
+// ihb_top = Start_Index( atom_i->Hindex, hbonds );
+// }
+// }
+// else
+// {
+// ihb_top = -1;
+// }
+// }
+// }
+//
+// // update i-j distance - check if j is within cutoff
+// for ( pj = start_i; pj < end_i; ++pj )
+// {
+// j = far_nbrs->far_nbr_list.nbr[pj];
+// atom_j = &system->my_atoms[j];
+//
+// if ( renbr )
+// {
+// if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
+// flag = 1;
+// else
+// flag = 0;
+// }
+// else
+// {
+// far_nbrs->far_nbr_list.dvec[pj][0] = atom_j->x[0] - atom_i->x[0];
+// far_nbrs->far_nbr_list.dvec[pj][1] = atom_j->x[1] - atom_i->x[1];
+// far_nbrs->far_nbr_list.dvec[pj][2] = atom_j->x[2] - atom_i->x[2];
+// far_nbrs->far_nbr_list.d[pj] = rvec_Norm_Sqr( far_nbrs->far_nbr_list.dvec[pj] );
+//
+// if ( far_nbrs->far_nbr_list.d[pj] <= SQR(cutoff) )
+// {
+// far_nbrs->far_nbr_list.d[pj] = sqrt( far_nbrs->far_nbr_list.d[pj] );
+// flag = 1;
+// }
+// else
+// {
+// flag = 0;
+// }
+// }
+//
+// if ( flag )
+// {
+// type_j = atom_j->type;
+// r_ij = far_nbrs->far_nbr_list.d[pj];
+// sbp_j = &system->reax_param.sbp[type_j];
+// twbp = &system->reax_param.tbp[type_i][type_j];
+//
+// if ( local == 1 )
+// {
+// // H matrix entry
+//#if defined(NEUTRAL_TERRITORY)
+// if ( atom_j->nt_dir > 0 || (j < system->n
+// && (H->format == SYM_FULL_MATRIX
+// || (H->format == SYM_HALF_MATRIX && i < j))) )
+// {
+// if( j < system->n )
+// {
+// H->entries[Htop].j = j;
+// }
+// else
+// {
+// H->entries[Htop].j = atom_j->pos;
+// }
+//
+// if ( control->tabulate == 0 )
+// {
+// H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
+// }
+// else
+// {
+// H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
+// }
+//
+// ++Htop;
+// }
+//#else
+// if ( (far_nbrs->format == HALF_LIST
+// && (j < system->n || atom_i->orig_id < atom_j->orig_id))
+// || far_nbrs->format == FULL_LIST )
+// {
+// H->entries[Htop].j = j;
+//
+// if ( control->tabulate == 0 )
+// {
+// H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
+// }
+// else
+// {
+// H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
+// }
+//
+// ++Htop;
+// }
+//#endif
+//
+// // hydrogen bond lists
+// if ( control->hbond_cut > 0.0
+// && (ihb == 1 || ihb == 2)
+// && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
+// {
+// // fprintf( stderr, "%d %d\n", atom1, atom2 );
+// jhb = sbp_j->p_hbond;
+// if ( ihb == 1 && jhb == 2 )
+// {
+// hbonds->hbond_list[ihb_top].nbr = j;
+// hbonds->hbond_list[ihb_top].scl = 1;
+// hbonds->hbond_list[ihb_top].ptr = pj;
+// ++ihb_top;
+// ++num_hbonds;
+// }
+// // only add to list for local j (far nbrs is half-list)
+// else if ( far_nbrs->format == HALF_LIST
+// && (j < system->n && ihb == 2 && jhb == 1) )
+// {
+// jhb_top = End_Index( atom_j->Hindex, hbonds );
+// hbonds->hbond_list[jhb_top].nbr = i;
+// hbonds->hbond_list[jhb_top].scl = -1;
+// hbonds->hbond_list[jhb_top].ptr = pj;
+// Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
+// ++num_hbonds;
+// }
+// }
+// }
+//#if defined(NEUTRAL_TERRITORY)
+// else if ( local == 2 )
+// {
+// // H matrix entry
+// if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
+// && ( H->format == SYM_FULL_MATRIX
+// || (H->format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
+// || ( j < system->n && atom_i->nt_dir != 0 && H->format == SYM_FULL_MATRIX ))
+// {
+// if( !nt_flag )
+// {
+// nt_flag = 1;
+// H->start[atom_i->pos] = Htop;
+// }
+//
+// if( j < system->n )
+// {
+// H->entries[Htop].j = j;
+// }
+// else
+// {
+// H->entries[Htop].j = atom_j->pos;
+// }
+//
+// if ( control->tabulate == 0 )
+// {
+// H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
+// }
+// else
+// {
+// H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
+// }
+//
+// ++Htop;
+// }
+// }
+//#endif
+//
+// // uncorrected bond orders
+// if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
+// far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
+// && BOp( workspace, bonds, control->bo_cut,
+// i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+// &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+// &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+// sbp_i, sbp_j, twbp ) )
+// {
+// num_bonds += 2;
+// ++btop_i;
+//
+// if ( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
+// workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
+// else if ( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 )
+// {
+// workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
+// }
+// }
+// }
+// }
+//
+// Set_End_Index( i, btop_i, bonds );
+// if ( local == 1 )
+// {
+// H->end[i] = Htop;
+// if ( ihb == 1 )
+// Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
+// }
+//#if defined(NEUTRAL_TERRITORY)
+// else if ( local == 2 )
+// {
+// if( nt_flag )
+// {
+// H->end[atom_i->pos] = Htop;
+// }
+// else
+// {
+// H->start[atom_i->pos] = 0;
+// H->end[atom_i->pos] = 0;
+// }
+// }
+//#endif
+// }
+//
+// 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 )
+// {
+// start_i = Start_Index( i, bonds );
+// end_i = End_Index( i, bonds );
+//
+// for ( btop_i = start_i; btop_i < end_i; ++btop_i )
+// {
+// j = bonds->bond_list[btop_i].nbr;
+// start_j = Start_Index( j, bonds );
+// end_j = End_Index( j, bonds );
+//
+// for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+// {
+// if ( bonds->bond_list[btop_j].nbr == i )
+// {
+// bonds->bond_list[btop_i].sym_index = btop_j;
+// break;
+// }
+// }
+// }
+// }
+// }
+//
+//#if defined(DEBUG)
+// Print_Sparse_Matrix2( system, H, NULL );
+//#endif
+//
+// workspace->realloc.Htop = Htop;
+// workspace->realloc.num_bonds = num_bonds;
+// workspace->realloc.num_hbonds = num_hbonds;
+//
+//#if defined(DEBUG_FOCUS)
+// fprintf( stderr, "p%d @ step%d: Htop = %d num_bonds = %d num_hbonds = %d\n",
+// system->my_rank, data->step, Htop, num_bonds, num_hbonds );
+// MPI_Barrier( comm );
+//#endif
+//
+//#if defined( DEBUG )
+// Print_Bonds( system, bonds, "debugbonds.out" );
+// Print_Bond_List2( system, bonds, "pbonds.out" );
+// Print_Sparse_Matrix( system, H );
+// for ( i = 0; i < H->n; ++i )
+// for ( j = H->start[i]; j < H->end[i]; ++j )
+// fprintf( stderr, "%d %d %.15e\n",
+// MIN(system->my_atoms[i].orig_id,
+// system->my_atoms[H->entries[j].j].orig_id),
+// MAX(system->my_atoms[i].orig_id,
+// system->my_atoms[H->entries[j].j].orig_id),
+// H->entries[j].val );
+//#endif
+//
+// Validate_Lists( system, workspace, lists, data->step,
+// system->n, system->N, system->numH, comm );
+//
+//}
void Init_Forces_noQEq( reax_system *system, control_params *control,
@@ -1966,11 +1987,13 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
hbonds = lists[HBONDS];
for ( i = 0; i < system->n; ++i )
+ {
workspace->bond_mark[i] = 0;
+ }
for ( i = system->n; i < system->N; ++i )
{
- workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
- //workspace->done_after[i] = Start_Index( i, far_nbrs );
+ /* put ghost atoms to an infinite distance (i.e., 1000) */
+ workspace->bond_mark[i] = 1000;
}
num_bonds = 0;
@@ -2317,21 +2340,33 @@ void Estimate_Storages( reax_system *system, control_params *control,
C12 = twbp->p_bo1 * pow( r_ij / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + control->bo_cut) * exp( C12 );
}
- else BO_s = C12 = 0.0;
+ else
+ {
+ C12 = 0.0;
+ BO_s = 0.0;
+ }
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0)
{
C34 = twbp->p_bo3 * pow( r_ij / twbp->r_p, twbp->p_bo4 );
BO_pi = exp( C34 );
}
- else BO_pi = C34 = 0.0;
+ else
+ {
+ C34 = 0.0;
+ BO_pi = 0.0;
+ }
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0)
{
C56 = twbp->p_bo5 * pow( r_ij / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = exp( C56 );
}
- else BO_pi2 = C56 = 0.0;
+ else
+ {
+ C56 = 0.0;
+ BO_pi2 = 0.0;
+ }
/* Initially BO values are the uncorrected ones, page 1 */
BO = BO_s + BO_pi + BO_pi2;
@@ -2354,7 +2389,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
* Therefore, we assume it is full and divide 2 if necessary. */
if ( cm_format == SYM_HALF_MATRIX )
{
- *Htop = (*Htop + system->n) / 2;
+ *Htop = (*Htop + system->n + 1) / 2;
}
#endif
@@ -2373,10 +2408,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
- *num_3body += SQR(bond_top[i]);
- //if( i < system->n )
+ *num_3body += SQR( bond_top[i] );
+ //TODO: why x2?
bond_top[i] = MAX( bond_top[i] * 2, MIN_BONDS );
- //else bond_top[i] = MAX_BONDS;
}
#if defined(DEBUG_FOCUS)
@@ -2397,12 +2431,12 @@ void Compute_Forces( reax_system *system, control_params *control,
#if defined(LOG_PERFORMANCE)
real t_start = 0.0, t_end;
- //MPI_Barrier( mpi_data->world );
if ( system->my_rank == MASTER_NODE )
{
t_start = MPI_Wtime();
}
#endif
+
comm = mpi_data->world;
/********* init forces ************/
@@ -2467,7 +2501,6 @@ void Compute_Forces( reax_system *system, control_params *control,
}
#if defined(LOG_PERFORMANCE)
- //MPI_Barrier( mpi_data->world );
if ( system->my_rank == MASTER_NODE )
{
t_end = MPI_Wtime( );
@@ -2487,7 +2520,6 @@ void Compute_Forces( reax_system *system, control_params *control,
lists, out_control, mpi_data->world );
#if defined(LOG_PERFORMANCE)
- //MPI_Barrier( mpi_data->world );
if ( system->my_rank == MASTER_NODE )
{
t_end = MPI_Wtime( );
@@ -2506,7 +2538,6 @@ void Compute_Forces( reax_system *system, control_params *control,
Compute_Total_Force( system, control, data, workspace, lists, mpi_data );
#if defined(LOG_PERFORMANCE)
- //MPI_Barrier( mpi_data->world );
if ( system->my_rank == MASTER_NODE )
{
t_end = MPI_Wtime( );
@@ -2517,12 +2548,13 @@ void Compute_Forces( reax_system *system, control_params *control,
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: total forces computed\n",
system->my_rank, data->step );
+
//Print_Total_Force( system, data, workspace );
MPI_Barrier( mpi_data->world );
#endif
#if defined(TEST_FORCES)
Print_Force_Files( system, control, data, workspace,
- lists, out_control, mpi_data );
+ lists, out_control, mpi_data );
#endif
}
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 8a94c0c9..d93ec5c2 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -1033,9 +1033,9 @@ typedef struct
real *dDelta_lp, *dDelta_lp_temp;
real *nlp, *nlp_temp, *Clp, *vlpex;
rvec *dDeltap_self;
- int *bond_mark, *done_after;
+ int *bond_mark;
- /* QEq storage */
+ /* charge matrix storage */
sparse_matrix *H;
sparse_matrix *L;
sparse_matrix *U;
@@ -1051,23 +1051,23 @@ typedef struct
real *y, *g;
real *hc, *hs;
real **h, **v;
- /* GMRES, PIPECG storage */
+ /* GMRES, PIPECG, PIPECR storage */
real *z;
- /* CG, PIPECG storage */
+ /* CG, PIPECG, PIPECR storage */
real *d, *p, *q, *r;
- /* PIPECG storage */
+ /* PIPECG, PIPECR storage */
real *m, *n, *u, *w;
/* dual-CG storage */
rvec2 *r2, *d2, *q2, *p2;
/* Taper */
- real Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
+ real Tap[8];
/* storage for analysis */
- int *mark, *old_mark;
+ int *mark, *old_mark;
rvec *x_old;
/* storage space for bond restrictions */
- int *restricted;
+ int *restricted;
int **restricted_list;
/* integrator */
--
GitLab