diff --git a/PG-PuReMD/src/bond_orders.c b/PG-PuReMD/src/bond_orders.c
index 2b494cac938689518e9cd67c655f5996a3ea14b5..c6928d7a3beeadd2143b7fce4c7d83aa7c0ffa9a 100644
--- a/PG-PuReMD/src/bond_orders.c
+++ b/PG-PuReMD/src/bond_orders.c
@@ -919,19 +919,19 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
of bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients such that
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
- bo_ij->C1dbo = 1.000000;
- bo_ij->C2dbo = 0.000000;
- bo_ij->C3dbo = 0.000000;
-
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
-
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbo = 1.0;
+ bo_ij->C2dbo = 0.0;
+ bo_ij->C3dbo = 0.0;
+
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
+
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
diff --git a/PG-PuReMD/src/cuda/cuda_bond_orders.cu b/PG-PuReMD/src/cuda/cuda_bond_orders.cu
index 8bcd31264c579b7bf859567db8e1f8ef3cf3dbf3..bb13428ce4b20019ae742c902ac885ed84e1c3c2 100644
--- a/PG-PuReMD/src/cuda/cuda_bond_orders.cu
+++ b/PG-PuReMD/src/cuda/cuda_bond_orders.cu
@@ -113,15 +113,15 @@ CUDA_GLOBAL void Cuda_Calculate_BO( reax_atom *my_atoms, global_parameters gp,
bo_ij->C2dbo = 0.000000;
bo_ij->C3dbo = 0.000000;
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
@@ -247,12 +247,12 @@ CUDA_GLOBAL void Cuda_Calculate_BO( reax_atom *my_atoms, global_parameters gp,
bo_ij->C2dbo = bo_ij->BO * A2_ij;
bo_ij->C3dbo = bo_ij->BO * A2_ji;
- bo_ij->C1dbopi = f1*f1*f4*f5;
+ bo_ij->C1dbopi = f1 * f1 * f4 * f5;
bo_ij->C2dbopi = bo_ij->BO_pi * A1_ij;
bo_ij->C3dbopi = bo_ij->BO_pi * A3_ij;
bo_ij->C4dbopi = bo_ij->BO_pi * A3_ji;
- bo_ij->C1dbopi2 = f1*f1*f4*f5;
+ bo_ij->C1dbopi2 = f1 * f1 * f4 * f5;
bo_ij->C2dbopi2 = bo_ij->BO_pi2 * A1_ij;
bo_ij->C3dbopi2 = bo_ij->BO_pi2 * A3_ij;
bo_ij->C4dbopi2 = bo_ij->BO_pi2 * A3_ji;
diff --git a/PuReMD-GPU/src/bond_orders.c b/PuReMD-GPU/src/bond_orders.c
index 49eaed6532449ef34ee3dac26a71462b2edbdd36..05386a87615000dd85bab82885b893f0604a6476 100644
--- a/PuReMD-GPU/src/bond_orders.c
+++ b/PuReMD-GPU/src/bond_orders.c
@@ -765,19 +765,19 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients s.t.
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
- bo_ij->C1dbo = 1.000000;
- bo_ij->C2dbo = 0.000000;
- bo_ij->C3dbo = 0.000000;
-
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
-
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbo = 1.0;
+ bo_ij->C2dbo = 0.0;
+ bo_ij->C3dbo = 0.0;
+
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
+
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
diff --git a/PuReMD-GPU/src/cuda_bond_orders.cu b/PuReMD-GPU/src/cuda_bond_orders.cu
index 81a7462033a3da94d8ae601288efb83a93e387f2..ed6b73e4c088deda7553f2248e947437230abf81 100644
--- a/PuReMD-GPU/src/cuda_bond_orders.cu
+++ b/PuReMD-GPU/src/cuda_bond_orders.cu
@@ -415,19 +415,19 @@ GLOBAL void Cuda_Calculate_Bond_Orders ( reax_atom *atoms, global_parameters g_
bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients s.t.
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
- bo_ij->C1dbo = 1.000000;
- bo_ij->C2dbo = 0.000000;
- bo_ij->C3dbo = 0.000000;
-
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
-
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbo = 1.0;
+ bo_ij->C2dbo = 0.0;
+ bo_ij->C3dbo = 0.0;
+
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
+
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs.select.dbo_list[ top_dbo ]);
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index ca54db6f77f6701bdc05b851920c9f9cde5f4471..cf6c69911f889bb528ef32c39e6542960ba0d7e0 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -894,19 +894,19 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
of bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients such that
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
- bo_ij->C1dbo = 1.000000;
- bo_ij->C2dbo = 0.000000;
- bo_ij->C3dbo = 0.000000;
-
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
-
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbo = 1.0;
+ bo_ij->C2dbo = 0.0;
+ bo_ij->C3dbo = 0.0;
+
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
+
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs->dbo_list[ top_dbo ]);
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 6c07fe4f2e8ae5a25f4fc0108057d74a7b00b98e..907fb1a0b50859381961ec56eb8da6ab6af94f74 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -24,38 +24,37 @@
#include "list.h"
#include "tool_box.h"
+
/* allocate space for atoms */
-int PreAllocate_Space( reax_system *system, control_params *control,
+void PreAllocate_Space( reax_system *system, control_params *control,
static_storage *workspace )
{
int i;
system->atoms = (reax_atom*) scalloc( system->N,
- sizeof(reax_atom), "atoms" );
+ sizeof(reax_atom), "PreAllocate_Space::system->atoms" );
workspace->orig_id = (int*) scalloc( system->N,
- sizeof(int), "orid_id" );
+ sizeof(int), "PreAllocate_Space::workspace->orid_id" );
/* space for keeping restriction info, if any */
if ( control->restrict_bonds )
{
workspace->restricted = (int*) scalloc( system->N,
- sizeof(int), "restricted_atoms" );
+ sizeof(int), "PreAllocate_Space::workspace->restricted_atoms" );
workspace->restricted_list = (int**) scalloc( system->N,
- sizeof(int*), "restricted_list" );
+ sizeof(int*), "PreAllocate_Space::workspace->restricted_list" );
for ( i = 0; i < system->N; ++i )
{
workspace->restricted_list[i] = (int*) scalloc( MAX_RESTRICT,
- sizeof(int), "restricted_list[i]" );
+ sizeof(int), "PreAllocate_Space::workspace->restricted_list[i]" );
}
}
-
- return SUCCESS;
}
-void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
+static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
{
Delete_List( TYP_FAR_NEIGHBOR, far_nbrs );
@@ -70,36 +69,36 @@ void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
}
-/* dynamic allocation of memory for matrix in CSR format */
-int Allocate_Matrix( sparse_matrix **pH, int n, int m )
+/* Dynamic allocation of memory for matrix in CSR format
+ *
+ * pH (output): pointer to sparse matrix for which to allocate
+ * n: dimension of the matrix
+ * m: number of nonzeros to allocate space for in matrix
+ * */
+void Allocate_Matrix( sparse_matrix **pH, int n, int m )
{
sparse_matrix *H;
- if ( (*pH = (sparse_matrix*) smalloc( sizeof(sparse_matrix),
- "Allocate_Matrix::pH" )) == NULL )
- {
- return FAILURE;
- }
+ *pH = (sparse_matrix*) smalloc( sizeof(sparse_matrix),
+ "Allocate_Matrix::pH" );
H = *pH;
H->n = n;
H->m = m;
- if ( (H->start = (unsigned int*) smalloc( sizeof(unsigned int) * (n + 1),
- "Allocate_Matrix::H->start" )) == NULL
- || (H->j = (unsigned int*) smalloc( sizeof(unsigned int) * m,
- "Allocate_Matrix::H->j" )) == NULL
- || (H->val = (real*) smalloc( sizeof(real) * m,
- "Allocate_Matrix::H->val" )) == NULL )
- {
- return FAILURE;
- }
-
- return SUCCESS;
+ H->start = (unsigned int*) smalloc( sizeof(unsigned int) * (n + 1),
+ "Allocate_Matrix::H->start" );
+ H->j = (unsigned int*) smalloc( sizeof(unsigned int) * m,
+ "Allocate_Matrix::H->j" );
+ H->val = (real*) smalloc( sizeof(real) * m,
+ "Allocate_Matrix::H->val" );
}
-/* deallocate memory for matrix in CSR format */
+/* Deallocate memory for matrix in CSR format
+ *
+ * H (output): pointer to sparse matrix for which to allocate
+ * */
void Deallocate_Matrix( sparse_matrix *H )
{
sfree( H->start, "Deallocate_Matrix::H->start" );
@@ -109,28 +108,15 @@ void Deallocate_Matrix( sparse_matrix *H )
}
-int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name )
+static void Reallocate_Matrix( sparse_matrix **H, int n, int m )
{
Deallocate_Matrix( *H );
- if ( Allocate_Matrix( H, n, m ) == FAILURE )
- {
- fprintf(stderr, "not enough space for %s matrix. terminating!\n", name);
- exit( INSUFFICIENT_MEMORY );
- }
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "reallocating %s matrix, n = %d, m = %d\n",
- name, n, m );
- fprintf( stderr, "memory allocated: %s = %ldMB\n",
- name, m * sizeof(sparse_matrix_entry) / (1024 * 1024) );
-#endif
-
- return SUCCESS;
+ Allocate_Matrix( H, n, m );
}
-int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
+void Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
reax_list *hbonds )
{
int i, num_hbonds;
@@ -164,19 +150,14 @@ int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
#endif
- return SUCCESS;
}
-int Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
+static void Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
{
int i;
int *hb_top;
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "reallocating hbonds\n" );
-#endif
-
hb_top = (int *) scalloc( n, sizeof(int), "Reallocate_HBonds_List::hb_top" );
for ( i = 0; i < n; ++i )
{
@@ -191,12 +172,10 @@ int Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
Allocate_HBond_List( n, num_h, h_index, hb_top, hbonds );
sfree( hb_top, "Reallocate_HBonds_List::hb_top" );
-
- return SUCCESS;
}
-int Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
+void Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
{
int i, num_bonds;
@@ -223,11 +202,10 @@ int Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
fprintf( stderr, "memory allocated: bonds = %ldMB\n",
num_bonds * sizeof(bond_data) / (1024 * 1024) );
#endif
- return SUCCESS;
}
-int Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
+static void Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
{
int i;
int *bond_top;
@@ -251,8 +229,6 @@ int Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3bo
*num_bonds = bond_top[n - 1];
sfree( bond_top, "Reallocate_Bonds_List::bond_top" );
-
- return SUCCESS;
}
@@ -276,7 +252,8 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
if ( realloc->Htop > 0 )
{
- Reallocate_Matrix(&(workspace->H), system->N_cm, realloc->Htop * SAFE_ZONE, "H");
+ Reallocate_Matrix( &(workspace->H), system->N_cm,
+ realloc->Htop * SAFE_ZONE );
realloc->Htop = -1;
Deallocate_Matrix( workspace->L );
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index b961521eccce628544e254289c896f3846544cfb..10a401de7b0e6b959df0417e1b2ac72961c28fc2 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -24,16 +24,18 @@
#include "mytypes.h"
-int PreAllocate_Space( reax_system*, control_params*, static_storage* );
+
+void PreAllocate_Space( reax_system*, control_params*, static_storage* );
void Reallocate( reax_system*, control_params*, static_storage*, reax_list**, int );
-int Allocate_Matrix( sparse_matrix**, int, int );
+void Allocate_Matrix( sparse_matrix**, int, int );
void Deallocate_Matrix( sparse_matrix* );
-int Allocate_HBond_List( int, int, int*, int*, reax_list* );
+void Allocate_HBond_List( int, int, int*, int*, reax_list* );
+
+void Allocate_Bond_List( int, int*, reax_list* );
-int Allocate_Bond_List( int, int*, reax_list* );
#endif
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index 1b260eea95d5e16ab443a8ba3d5ccae8a749e22a..3f63f09d9c151c304c57be1aca93ace4d0885147 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -60,7 +60,7 @@ typedef struct
// copy bond list into old bond list
-void Copy_Bond_List( reax_system *system, control_params *control,
+static void Copy_Bond_List( reax_system *system, control_params *control,
reax_list **lists )
{
int i, j, top_old;
@@ -91,7 +91,7 @@ void Copy_Bond_List( reax_system *system, control_params *control,
// ASSUMPTION: Bond lists are sorted
-int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
+static int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
{
int oldp, newp;
reax_list *new_bonds = &(*lists)[BONDS];
@@ -169,7 +169,7 @@ int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
}
-void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
+static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
control_params *control, reax_list *bonds, int print,
FILE *fout )
{
@@ -193,7 +193,7 @@ void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
}
-void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
+static void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
size_t size )
{
int j, atom;
@@ -220,7 +220,7 @@ void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
}
-void Analyze_Molecules( reax_system *system, control_params *control,
+static void Analyze_Molecules( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, FILE *fout )
{
@@ -318,7 +318,7 @@ void Analyze_Molecules( reax_system *system, control_params *control,
}
-void Report_Bond_Change( reax_system *system, control_params *control,
+static void Report_Bond_Change( reax_system *system, control_params *control,
static_storage *workspace, reax_list *old_bonds,
reax_list *new_bonds, int a1, int a2, int flag,
FILE *fout )
@@ -389,7 +389,7 @@ void Report_Bond_Change( reax_system *system, control_params *control,
/* ASSUMPTION: Bond lists are sorted */
-void Compare_Bonding( int atom, reax_system *system, control_params *control,
+static void Compare_Bonding( int atom, reax_system *system, control_params *control,
static_storage *workspace, reax_list *old_bonds,
reax_list *new_bonds, FILE *fout )
{
@@ -507,7 +507,7 @@ void Compare_Bonding( int atom, reax_system *system, control_params *control,
}
-void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
+static void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
control_params *control, reax_list *bonds, int ignore )
{
int i, t, start, end, nbr;
@@ -531,8 +531,7 @@ void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
}
-
-void Analyze_Fragments( reax_system *system, control_params *control,
+static void Analyze_Fragments( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, FILE *fout, int ignore )
{
@@ -601,7 +600,7 @@ void Analyze_Fragments( reax_system *system, control_params *control,
}
-void Analyze_Silica( reax_system *system, control_params *control,
+static void Analyze_Silica( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, FILE *fout )
{
@@ -724,7 +723,7 @@ void Analyze_Silica( reax_system *system, control_params *control,
}
-int Get_Type_of_Molecule( molecule *m )
+static int Get_Type_of_Molecule( molecule *m )
{
if ( m->atom_count == 3 && m->mtypes[1] == 2 && m->mtypes[2] == 1 )
{
@@ -735,7 +734,7 @@ int Get_Type_of_Molecule( molecule *m )
}
-void Calculate_Dipole_Moment( reax_system *system, control_params *control,
+static void Calculate_Dipole_Moment( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list *bonds, FILE *fout )
{
@@ -781,7 +780,7 @@ void Calculate_Dipole_Moment( reax_system *system, control_params *control,
}
-void Copy_Positions( reax_system *system, static_storage *workspace )
+static void Copy_Positions( reax_system *system, static_storage *workspace )
{
int i;
@@ -790,7 +789,7 @@ void Copy_Positions( reax_system *system, static_storage *workspace )
}
-void Calculate_Drift( reax_system *system, control_params *control,
+static void Calculate_Drift( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
FILE *fout )
{
@@ -846,7 +845,7 @@ void Calculate_Drift( reax_system *system, control_params *control,
}
-void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
+static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
FILE *fout )
{
int i, j, k;
@@ -913,7 +912,7 @@ void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
}
-void Calculate_Density_Slice( reax_system *system, simulation_data *data,
+static void Calculate_Density_Slice( reax_system *system, simulation_data *data,
FILE *fout )
{
real slice_thickness = 0.5;
@@ -957,14 +956,17 @@ void Analysis( reax_system *system, control_params *control,
int steps;
steps = data->step - data->prev_steps;
- // fprintf( stderr, "prev_steps: %d\n", data->prev_steps );
if ( steps == 1 )
{
if ( control->molec_anal == REACTIONS )
+ {
Copy_Bond_List( system, control, lists );
+ }
if ( control->diffusion_coef )
+ {
Copy_Positions( system, workspace );
+ }
}
/****** Molecular Analysis ******/
@@ -977,25 +979,29 @@ void Analysis( reax_system *system, control_params *control,
out_control->mol, 0 );
/* discover fragments without the ignored atoms */
if ( control->num_ignored )
+ {
Analyze_Fragments( system, control, data, workspace, lists,
out_control->ign, 1 );
+ }
}
else if ( control->molec_anal == REACTIONS )
+ {
/* discover molecular changes - reactions */
Analyze_Molecules( system, control, data, workspace,
lists, out_control->mol );
+ }
}
/****** Electric Dipole Moment ******/
if ( control->dipole_anal && steps % control->freq_dipole_anal == 0 )
+ {
Calculate_Dipole_Moment( system, control, data, workspace,
&(*lists)[BONDS], out_control->dpl );
+ }
/****** Drift ******/
if ( control->diffusion_coef && steps % control->freq_diffusion_coef == 0 )
+ {
Calculate_Drift( system, control, data, workspace, out_control->drft );
-
-#if defined(DEBUG)
- fprintf( stderr, "analysis... done\n" );
-#endif
+ }
}
diff --git a/sPuReMD/src/analyze.h b/sPuReMD/src/analyze.h
index af2b056f5aabd598f9660a5956dd40f027ab3b4c..c8328240315db071b996c4bf19a1bd10ab1fc097 100644
--- a/sPuReMD/src/analyze.h
+++ b/sPuReMD/src/analyze.h
@@ -24,6 +24,7 @@
#include "mytypes.h"
+
void Analysis( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
@@ -50,4 +51,5 @@ void Analysis( reax_system*, control_params*, simulation_data*,
//void Calculate_Density_Slice( reax_system*, simulation_data*, FILE* );
+
#endif
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index d8f71ee014c4d6b072476a4e5fe1138fec31e425..81311ffa9dcc3acf8d1295345e54926ab2185ffc 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "bond_orders.h"
+
#include "list.h"
#include "lookup.h"
#include "print_utils.h"
@@ -671,7 +672,7 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
/* Locate j on i's list.
This function assumes that j is there for sure!
And this is the case given our method of neighbor generation*/
-int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
+static int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
{
int start = Start_Index(i, bonds);
int end = End_Index(i, bonds);
@@ -722,7 +723,7 @@ static inline void Copy_Bond_Order_Data( bond_order_data *dest, bond_order_data
}
-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;
}
@@ -832,19 +833,19 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients s.t.
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
- bo_ij->C1dbo = 1.000000;
- bo_ij->C2dbo = 0.000000;
- bo_ij->C3dbo = 0.000000;
-
- bo_ij->C1dbopi = bo_ij->BO_pi;
- bo_ij->C2dbopi = 0.000000;
- bo_ij->C3dbopi = 0.000000;
- bo_ij->C4dbopi = 0.000000;
-
- bo_ij->C1dbopi2 = bo_ij->BO_pi2;
- bo_ij->C2dbopi2 = 0.000000;
- bo_ij->C3dbopi2 = 0.000000;
- bo_ij->C4dbopi2 = 0.000000;
+ bo_ij->C1dbo = 1.0;
+ bo_ij->C2dbo = 0.0;
+ bo_ij->C3dbo = 0.0;
+
+ bo_ij->C1dbopi = 1.0;
+ bo_ij->C2dbopi = 0.0;
+ bo_ij->C3dbopi = 0.0;
+ bo_ij->C4dbopi = 0.0;
+
+ bo_ij->C1dbopi2 = 1.0;
+ bo_ij->C2dbopi2 = 0.0;
+ bo_ij->C3dbopi2 = 0.0;
+ bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
diff --git a/sPuReMD/src/bond_orders.h b/sPuReMD/src/bond_orders.h
index 67592ed87452d69c35aeeb767ad8bf96bdd4dc22..37cd62428e92c822cb1453e122808db4d42b5ac8 100644
--- a/sPuReMD/src/bond_orders.h
+++ b/sPuReMD/src/bond_orders.h
@@ -24,6 +24,7 @@
#include "mytypes.h"
+
typedef struct
{
real C1dbo, C2dbo, C3dbo;
@@ -32,6 +33,7 @@ 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* );
@@ -52,4 +54,5 @@ void Add_dBond_to_Forces_NPT( int, int, reax_system*, simulation_data*,
static_storage*, reax_list** );
void Calculate_Bond_Orders( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
#endif
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index f3646ae370f6a838baac9846e9167c8533f805bf..fba8bd3223c7f5284ca106d69d843041bcdeb572 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "box.h"
+
#include "tool_box.h"
#include "vector.h"
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index 4902a3404bc60be6dc18a4f4185cf7da9f25d5a9..8282c14cf972725dfb1c648322ddbd374955fb03 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -24,16 +24,17 @@
#include "mytypes.h"
-void Setup_Box( real, real, real, real, real, real, simulation_box* );
-
-/* Initializes box from box rtensor */
-void Update_Box(rtensor, simulation_box* /*, int*/);
-void Update_Box_Isotropic(simulation_box*, real /*, int*/);
-void Update_Box_SemiIsotropic( simulation_box*, rvec /*, int*/ );
/* Computes all the transformations,
metric and other quantities from box rtensor */
-void Make_Consistent(simulation_box*/*, int*/ );
+void Make_Consistent( simulation_box* );
+
+void Setup_Box( real, real, real, real, real, real, simulation_box* );
+
+/* Initializes box from box rtensor */
+void Update_Box( rtensor, simulation_box* );
+void Update_Box_Isotropic( simulation_box*, real );
+void Update_Box_SemiIsotropic( simulation_box*, rvec );
int Are_Far_Neighbors( rvec, rvec, simulation_box*, real, far_neighbor_data* );
void Get_NonPeriodic_Far_Neighbors( rvec, rvec, simulation_box*,
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 5c2d3bed963b38ea951da040cefe96b8dc9b14de..5d2faa75397beec670a3dedc38c758d02b4d5a2e 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -743,7 +743,6 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
}
-
static void Setup_Preconditioner_QEq( const reax_system * const system,
const control_params * const control,
simulation_data * const data, static_storage * const workspace,
@@ -791,17 +790,16 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
-
+ Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
}
- else
+ else if ( workspace->L->m < fillin )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
}
break;
@@ -809,16 +807,17 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -827,17 +826,20 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ /* TODO: safest storage estimate is ILU(0)
+ * (same as lower triangular portion of H), could improve later */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* TODO: safest storage estimate is ILU(0)
+ * (same as lower triangular portion of H), could improve later */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -845,16 +847,17 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -865,7 +868,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
break;
default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
@@ -928,17 +931,16 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
if ( workspace->L == NULL )
{
- if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
-
+ Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm );
+ Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm );
}
- else
+ else if ( workspace->L->m < fillin + system->N_cm )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm );
+ Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm );
}
break;
@@ -946,16 +948,17 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -964,17 +967,20 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
if ( workspace->L == NULL )
{
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ /* TODO: safest storage estimate is ILU(0)
+ * (same as lower triangular portion of H), could improve later */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* TODO: safest storage estimate is ILU(0)
+ * (same as lower triangular portion of H), could improve later */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -982,17 +988,22 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
+ }
+ else if ( workspace->L->m < Hptr->m )
+ {
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
exit( INSUFFICIENT_MEMORY );
}
}
- else
- {
- //TODO: reallocate
- }
break;
case SAI_PC:
@@ -1002,7 +1013,7 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
break;
default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
@@ -1073,16 +1084,17 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
if ( workspace->L == NULL )
{
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
}
- else
+ else if ( workspace->L->m < fillin )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
}
break;
@@ -1090,16 +1102,17 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -1108,17 +1121,19 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
if ( workspace->L == NULL )
{
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ /* TODO: safest storage estimate is ILU(0)
+ * (same as lower triangular portion of H), could improve later */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -1126,16 +1141,17 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
if ( workspace->L == NULL )
{
/* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
- else
+ else if ( workspace->L->m < Hptr->m )
{
- //TODO: reallocate
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+
+ /* factors have sparsity pattern as H */
+ Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
}
break;
@@ -1146,7 +1162,7 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
break;
default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
@@ -1212,8 +1228,6 @@ static void Calculate_Charges_EE( const reax_system * const system,
/* Main driver method for QEq kernel
- *
- * Rough outline:
* 1) init / setup routines for preconditioning of linear solver
* 2) compute preconditioner
* 3) extrapolate charges
@@ -1297,8 +1311,6 @@ static void QEq( reax_system * const system, control_params * const control,
/* Main driver method for EE kernel
- *
- * Rough outline:
* 1) init / setup routines for preconditioning of linear solver
* 2) compute preconditioner
* 3) extrapolate charges
@@ -1365,8 +1377,6 @@ static void EE( reax_system * const system, control_params * const control,
/* Main driver method for ACKS2 kernel
- *
- * Rough outline:
* 1) init / setup routines for preconditioning of linear solver
* 2) compute preconditioner
* 3) extrapolate charges
@@ -1391,6 +1401,21 @@ static void ACKS2( reax_system * const system, control_params * const control,
Extrapolate_Charges_EE( system, control, data, workspace );
+#if defined(DEBUG_FOCUS)
+#define SIZE (200)
+ char fname[SIZE];
+ FILE * fp;
+
+ if ( data->step % 10 == 0 )
+ {
+ snprintf( fname, SIZE + 11, "s_%d_%s.out", data->step, control->sim_name );
+ fp = fopen( fname, "w" );
+ Vector_Print( fp, NULL, workspace->s[0], system->N_cm );
+ fclose( fp );
+ }
+#undef SIZE
+#endif
+
switch ( control->cm_solver_type )
{
case GMRES_S:
@@ -1440,18 +1465,23 @@ void Compute_Charges( reax_system * const system, control_params * const control
char fname[SIZE];
FILE * fp;
- if ( data->step >= 100 )
+ if ( data->step % 10 == 0 )
{
- snprintf( fname, SIZE + 11, "s_%d_%s.out", data->step, control->sim_name );
- fp = fopen( fname, "w" );
- Vector_Print( fp, NULL, workspace->s[0], system->N_cm );
- fclose( fp );
+ snprintf( fname, SIZE + 11, "H_%d_%s.out", data->step, control->sim_name );
+ Print_Sparse_Matrix2( workspace->H, fname, NULL );
+// Print_Sparse_Matrix_Binary( workspace->H, fname );
- snprintf( fname, SIZE + 11, "t_%d_%s.out", data->step, control->sim_name );
+ snprintf( fname, SIZE + 11, "b_s_%d_%s.out", data->step, control->sim_name );
fp = fopen( fname, "w" );
- Vector_Print( fp, NULL, workspace->t[0], system->N_cm );
+ Vector_Print( fp, NULL, workspace->b_s, system->N_cm );
fclose( fp );
+
+// snprintf( fname, SIZE + 11, "b_t_%d_%s.out", data->step, control->sim_name );
+// fp = fopen( fname, "w" );
+// Vector_Print( fp, NULL, workspace->b_t, system->N_cm );
+// fclose( fp );
}
+#undef SIZE
#endif
switch ( control->charge_method )
@@ -1473,24 +1503,4 @@ void Compute_Charges( reax_system * const system, control_params * const control
exit( INVALID_INPUT );
break;
}
-
-#if defined(DEBUG_FOCUS)
- if ( data->step >= 100 )
- {
- snprintf( fname, SIZE + 11, "H_%d_%s.out", data->step, control->sim_name );
- Print_Sparse_Matrix2( workspace->H, fname, NULL );
-// Print_Sparse_Matrix_Binary( workspace->H, fname );
-
- snprintf( fname, SIZE + 11, "b_s_%d_%s.out", data->step, control->sim_name );
- fp = fopen( fname, "w" );
- Vector_Print( fp, NULL, workspace->b_s, system->N_cm );
- fclose( fp );
-
- snprintf( fname, SIZE + 11, "b_t_%d_%s.out", data->step, control->sim_name );
- fp = fopen( fname, "w" );
- Vector_Print( fp, NULL, workspace->b_t, system->N_cm );
- fclose( fp );
- }
-#undef SIZE
-#endif
}
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index e0520fc6963f1bbf7eebe091de6c217680b6512b..48028b0c8fe9fa355c7dcd6f658d00e08d0db74b 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -24,8 +24,10 @@
#include "mytypes.h"
+
void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
static_storage* const, const reax_list* const,
const output_controls* const );
+
#endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 4f5616ab0d983290d219c57d92544a752d1dcba2..783c73852c03a4a0aa7790830d71009c881a7fff 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -19,14 +19,15 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
+#include "control.h"
+
#include <ctype.h>
-#include "control.h"
#include "traj.h"
#include "tool_box.h"
-char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
+void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
output_controls *out_control )
{
char *s, **tmp;
@@ -54,9 +55,6 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->restrict_bonds = 0;
control->periodic_boundaries = 1;
- control->periodic_images[0] = 0;
- control->periodic_images[1] = 0;
- control->periodic_images[2] = 0;
control->reneighbor = 1;
control->vlist_cut = 0;
@@ -198,15 +196,6 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
ival = atoi( tmp[1] );
control->periodic_boundaries = ival;
}
- else if ( strcmp(tmp[0], "periodic_images") == 0 )
- {
- ival = atoi(tmp[1]);
- control->periodic_images[0] = ival;
- ival = atoi(tmp[2]);
- control->periodic_images[1] = ival;
- ival = atoi(tmp[3]);
- control->periodic_images[2] = ival;
- }
else if ( strcmp(tmp[0], "geo_format") == 0 )
{
ival = atoi( tmp[1] );
@@ -599,6 +588,4 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
fprintf(stderr, "control file read\n" );
#endif
-
- return SUCCESS;
}
diff --git a/sPuReMD/src/control.h b/sPuReMD/src/control.h
index 66d0dde7b4901d7a7b42512414328a8e6b256d83..da79fe527aa670fb62269a3b8b35ef5c7d09afe9 100644
--- a/sPuReMD/src/control.h
+++ b/sPuReMD/src/control.h
@@ -24,6 +24,7 @@
#include "mytypes.h"
-char Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
+
+void Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
#endif
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 1fb05c45f755e49b7d756f0a36d243dce6e63f6e..f8b39c8e6552bd8f097a1cd1f255a40fd2ddef11 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -19,13 +19,14 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
+#include "ffield.h"
+
#include <ctype.h>
-#include "ffield.h"
#include "tool_box.h"
-char Read_Force_Field( FILE* fp, reax_interaction* reax )
+void Read_Force_Field( FILE* fp, reax_interaction* reax )
{
char *s;
char **tmp;
@@ -46,7 +47,6 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
/* reading first header comment */
fgets( s, MAX_LINE, fp );
-
/* line 2 is number of global parameters */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp );
@@ -55,8 +55,8 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
n = atoi(tmp[0]);
if (n < 1)
{
- fprintf( stderr, "WARNING: number of globals in ffield file is 0!\n" );
- return 1;
+ fprintf( stderr, "[WARNING] number of globals in ffield file is 0!\n" );
+ return;
}
reax->gp.n_global = n;
@@ -74,19 +74,16 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
reax->gp.l[i] = val;
}
-
/* next line is number of atom types and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp );
reax->num_atom_types = atoi(tmp[0]);
-
/* 3 lines of comments */
fgets(s, MAX_LINE, fp);
fgets(s, MAX_LINE, fp);
fgets(s, MAX_LINE, fp);
-
/* Allocating structures in reax_interaction */
reax->sbp = (single_body_parameters*)
scalloc( reax->num_atom_types, sizeof(single_body_parameters),
@@ -743,7 +740,7 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
k = atoi(tmp[1]) - 1;
m = atoi(tmp[2]) - 1;
- if (j < reax->num_atom_types && m < reax->num_atom_types)
+ if ( j < reax->num_atom_types && m < reax->num_atom_types )
{
val = atof(tmp[3]);
reax->hbp[j][k][m].r0_hb = val;
@@ -785,10 +782,4 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
}
sfree( tor_flag, "Read_Force_Field::tor_flag" );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "force field read\n" );
-#endif
-
- return SUCCESS;
}
diff --git a/sPuReMD/src/ffield.h b/sPuReMD/src/ffield.h
index 4aaf32a644861b069e8cf87e2eec68aadf4d3c84..e2e88a10a6e7eee68402c85a5eba12eea2dc9023 100644
--- a/sPuReMD/src/ffield.h
+++ b/sPuReMD/src/ffield.h
@@ -23,6 +23,9 @@
#define __FFIELD_H_
#include "mytypes.h"
-char Read_Force_Field( FILE*, reax_interaction* );
+
+
+void Read_Force_Field( FILE*, reax_interaction* );
+
#endif
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 24e5c2914a3e76b365a7dbf10acab4ee677f05b3..787837988b059dd97043fda691404c8bb77ae4a9 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -41,7 +41,7 @@ typedef enum
} MATRIX_ENTRY_POSITION;
-void Dummy_Interaction( reax_system *system, control_params *control,
+static void Dummy_Interaction( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
@@ -72,13 +72,11 @@ void Init_Bonded_Force_Functions( control_params *control,
}
-void Compute_Bonded_Forces( reax_system *system, control_params *control,
+static void Compute_Bonded_Forces( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control, interaction_function *interaction_functions )
{
-
int i;
- //real t_start, t_end, t_elapsed;
#ifdef TEST_ENERGY
/* Mark beginning of a new timestep in each energy file */
@@ -116,10 +114,6 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
(interaction_functions[i])( system, control, data, workspace,
lists, out_control );
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "f%d-", i );
-#endif
-
#ifdef TEST_FORCES
(Print_Interactions[i])(system, control, data, workspace,
lists, out_control);
@@ -128,7 +122,7 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
}
-void Compute_NonBonded_Forces( reax_system *system, control_params *control,
+static void Compute_NonBonded_Forces( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list** lists, output_controls *out_control )
{
@@ -146,10 +140,6 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
t_elapsed = Get_Timing_Info( t_start );
data->timing.cm += t_elapsed;
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "qeq - " );
-#endif
-
if ( control->tabulate <= 0 )
{
vdW_Coulomb_Energy( system, control, data, workspace, lists, out_control );
@@ -160,10 +150,6 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
lists, out_control );
}
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "nonb forces - " );
-#endif
-
#ifdef TEST_FORCES
Print_vdW_Coulomb_Forces( system, control, data, workspace,
lists, out_control );
@@ -173,7 +159,7 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
/* 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,
+static void Compute_Total_Force( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists )
{
int i;
@@ -212,8 +198,6 @@ void Compute_Total_Force( reax_system *system, control_params *control,
}
#ifdef _OPENMP
- #pragma omp barrier
-
#pragma omp for schedule(static)
for ( i = 0; i < system->N; ++i )
{
@@ -227,7 +211,7 @@ void Compute_Total_Force( reax_system *system, control_params *control,
}
-void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
+static void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
int Hmax, int Htop, int num_bonds, int num_hbonds )
{
int i, flag;
@@ -243,7 +227,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
if ( Htop > Hmax )
{
fprintf( stderr,
- "step%d - ran out of space on H matrix: Htop=%d, max = %d",
+ "[ERROR] step%d - ran out of space on H matrix: Htop=%d, max = %d",
step, Htop, Hmax );
exit( INSUFFICIENT_MEMORY );
}
@@ -266,7 +250,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
if ( flag > -1 )
{
- 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, flag, End_Index(flag, bonds), Start_Index(flag + 1, bonds) );
exit( INSUFFICIENT_MEMORY );
}
@@ -277,7 +261,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
if ( End_Index(i, bonds) > bonds->total_intrs )
{
- fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
+ fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
step, flag, End_Index(i, bonds), bonds->total_intrs );
exit( INSUFFICIENT_MEMORY );
}
@@ -304,7 +288,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
if ( flag > -1 )
{
- fprintf( stderr, "step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
+ fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
step, flag, End_Index(flag, hbonds), Start_Index(flag + 1, hbonds) );
exit( INSUFFICIENT_MEMORY );
}
@@ -316,7 +300,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
if ( End_Index(i, hbonds) > hbonds->total_intrs )
{
- fprintf( stderr, "step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
+ fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
step, flag, End_Index(i, hbonds), hbonds->total_intrs );
exit( INSUFFICIENT_MEMORY );
}
@@ -366,14 +350,14 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
break;
default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
break;
default:
- fprintf( stderr, "Invalid charge method. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
@@ -419,14 +403,14 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
break;
default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
break;
default:
- fprintf( stderr, "Invalid charge method. Terminating...\n" );
+ fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
exit( INVALID_INPUT );
break;
}
@@ -602,7 +586,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
}
-void Init_Forces( reax_system *system, control_params *control,
+static void Init_Forces( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
@@ -835,8 +819,8 @@ void Init_Forces( reax_system *system, control_params *control,
bo_ij->BO -= control->bo_cut;
bo_ji->BO_s -= control->bo_cut;
bo_ji->BO -= control->bo_cut;
- workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
- workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+ workspace->total_bond_order[i] += bo_ij->BO;
+ workspace->total_bond_order[j] += bo_ji->BO;
bo_ij->Cdbo = 0.0;
bo_ij->Cdbopi = 0.0;
bo_ij->Cdbopi2 = 0.0;
@@ -844,35 +828,6 @@ void Init_Forces( reax_system *system, control_params *control,
bo_ji->Cdbopi = 0.0;
bo_ji->Cdbopi2 = 0.0;
- /*fprintf( stderr, "%d %d %g %g %g\n",
- i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/
-
- /*fprintf( stderr, "Cln_BOp_s: %f, pbo2: %f, C12:%f\n",
- Cln_BOp_s, twbp->p_bo2, C12 );
- fprintf( stderr, "Cln_BOp_pi: %f, pbo4: %f, C34:%f\n",
- Cln_BOp_pi, twbp->p_bo4, C34 );
- fprintf( stderr, "Cln_BOp_pi2: %f, pbo6: %f, C56:%f\n",
- Cln_BOp_pi2, twbp->p_bo6, C56 );*/
- /*fprintf(stderr, "pbo1: %f, pbo2:%f\n", twbp->p_bo1, twbp->p_bo2);
- fprintf(stderr, "pbo3: %f, pbo4:%f\n", twbp->p_bo3, twbp->p_bo4);
- fprintf(stderr, "pbo5: %f, pbo6:%f\n", twbp->p_bo5, twbp->p_bo6);
- fprintf( stderr, "r_s: %f, r_p: %f, r_pp: %f\n",
- twbp->r_s, twbp->r_p, twbp->r_pp );
- fprintf( stderr, "C12: %g, C34:%g, C56:%g\n", C12, C34, C56 );*/
-
- /*fprintf( stderr, "\tfactors: %g %g %g\n",
- -(bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi +
- bo_ij->BO_pi2 * Cln_BOp_pp),
- -bo_ij->BO_pi * Cln_BOp_pi, -bo_ij->BO_pi2 * Cln_BOp_pi2 );*/
- /*fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
- bo_ij->dBOp[0], bo_ij->dBOp[1], bo_ij->dBOp[2] );
- fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
- bo_ij->dln_BOp_pi[0], bo_ij->dln_BOp_pi[1],
- bo_ij->dln_BOp_pi[2] );
- fprintf( stderr, "dBOpi2:\t[%g, %g, %g]\n\n",
- bo_ij->dln_BOp_pi2[0], bo_ij->dln_BOp_pi2[1],
- bo_ij->dln_BOp_pi2[2] );*/
-
Set_End_Index( j, btop_j + 1, bonds );
}
}
@@ -915,7 +870,7 @@ void Init_Forces( reax_system *system, control_params *control,
}
-void Init_Forces_Tab( reax_system *system, control_params *control,
+static void Init_Forces_Tab( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
@@ -1148,8 +1103,8 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
bo_ij->BO -= control->bo_cut;
bo_ji->BO_s -= control->bo_cut;
bo_ji->BO -= control->bo_cut;
- workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
- workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+ workspace->total_bond_order[i] += bo_ij->BO;
+ workspace->total_bond_order[j] += bo_ji->BO;
bo_ij->Cdbo = 0.0;
bo_ij->Cdbopi = 0.0;
bo_ij->Cdbopi2 = 0.0;
@@ -1322,34 +1277,21 @@ void Compute_Forces( reax_system *system, control_params *control,
t_elapsed = Get_Timing_Info( t_start );
data->timing.init_forces += t_elapsed;
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "init_forces - ");
-#endif
-
t_start = Get_Time( );
Compute_Bonded_Forces( system, control, data, workspace, lists, out_control,
interaction_functions );
t_elapsed = Get_Timing_Info( t_start );
data->timing.bonded += t_elapsed;
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "bonded_forces - ");
-#endif
-
t_start = Get_Time( );
Compute_NonBonded_Forces( system, control, data, workspace,
lists, out_control );
t_elapsed = Get_Timing_Info( t_start );
data->timing.nonb += t_elapsed;
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "nonbondeds - ");
-#endif
-
Compute_Total_Force( system, control, data, workspace, lists );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "totalforces - ");
//Print_Total_Force( system, control, data, workspace, lists, out_control );
#endif
@@ -1357,8 +1299,4 @@ void Compute_Forces( reax_system *system, control_params *control,
Print_Total_Force( system, control, data, workspace, lists, out_control );
Compare_Total_Forces( system, control, data, workspace, lists, out_control );
#endif
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "forces - ");
-#endif
}
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index eeccac22c003250eae1bad325fc5ccdad08d9d88..01c789b90ffe299c2d6f8ecde4641d3de8ad66f1 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -24,13 +24,16 @@
#include "mytypes.h"
+
void Init_Bonded_Force_Functions( control_params*,
interaction_function* );
void Compute_Forces( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ static_storage*, reax_list**, output_controls*,
+ interaction_function* );
void Estimate_Storage_Sizes( reax_system*, control_params*, reax_list**,
- int*, int*, int*, int* );
+ int*, int*, int*, int* );
+
+
#endif
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index 22d104763a7aaaca6c97dfed1619ff96e11672a3..c46e9bc6dfc3938b606cfe84f29cb84165161fd8 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -20,17 +20,17 @@
----------------------------------------------------------------------*/
#include "four_body_interactions.h"
+
#include "bond_orders.h"
#include "box.h"
#include "list.h"
#include "lookup.h"
#include "vector.h"
-#include "math.h"
#define MIN_SINE 1e-10
-real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
+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,
three_body_interaction_data *p_ijk,
three_body_interaction_data *p_jkl,
@@ -118,10 +118,22 @@ real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
-p_jkl->dcos_dk[1]/sin_jkl,
-p_jkl->dcos_dk[2]/sin_jkl );*/
- if ( sin_ijk >= 0 && sin_ijk <= MIN_SINE ) sin_ijk = MIN_SINE;
- else if ( sin_ijk <= 0 && sin_ijk >= -MIN_SINE ) sin_ijk = -MIN_SINE;
- if ( sin_jkl >= 0 && sin_jkl <= MIN_SINE ) sin_jkl = MIN_SINE;
- else if ( sin_jkl <= 0 && sin_jkl >= -MIN_SINE ) sin_jkl = -MIN_SINE;
+ if ( sin_ijk >= 0 && sin_ijk <= MIN_SINE )
+ {
+ sin_ijk = MIN_SINE;
+ }
+ else if ( sin_ijk <= 0 && sin_ijk >= -MIN_SINE )
+ {
+ sin_ijk = -MIN_SINE;
+ }
+ if ( sin_jkl >= 0 && sin_jkl <= MIN_SINE )
+ {
+ sin_jkl = MIN_SINE;
+ }
+ else if ( sin_jkl <= 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 );
diff --git a/sPuReMD/src/four_body_interactions.h b/sPuReMD/src/four_body_interactions.h
index 88261c511c334d2bfcd78a4e49f972ee16fdb34e..e7a8e64fc55616b3853f2ee3ab9a5d1fb064ac31 100644
--- a/sPuReMD/src/four_body_interactions.h
+++ b/sPuReMD/src/four_body_interactions.h
@@ -24,7 +24,9 @@
#include "mytypes.h"
+
void Four_Body_Interactions( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
+
#endif
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 549d29ffce7efe6e64fa889215c7b3f131170441..05ece8c81c4ea2b2543593b98cffb27d0d378f6a 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -27,8 +27,7 @@
#include "vector.h"
-/********************* geo format routines ******************/
-void Count_Geo_Atoms( FILE *geo, reax_system *system )
+static void Count_Geo_Atoms( FILE *geo, reax_system *system )
{
int i, serial;
rvec x;
@@ -45,7 +44,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
Fit_to_Periodic_Box( &(system->box), &x );
}
- fseek( geo, 0, SEEK_SET ); // set the pointer to the beginning of the file
+ fseek( geo, 0, SEEK_SET );
fgets( line, MAX_LINE, geo );
fgets( line, MAX_LINE, geo );
@@ -55,75 +54,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
}
-char Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
-{
-
- FILE *geo;
- char descriptor[9];
- int i, serial, top;
- real box_x, box_y, box_z, alpha, beta, gamma;
- rvec x;
- char element[3], name[9];
- reax_atom *atom;
-
- /* open the geometry file */
- if ( (geo = fopen(geo_file, "r")) == NULL )
- {
- fprintf( stderr, "Error opening the geo file! terminating...\n" );
- exit( FILE_NOT_FOUND );
- }
-
- /* read box information */
- fscanf( geo, CUSTOM_BOXGEO_FORMAT,
- descriptor, &box_x, &box_y, &box_z, &alpha, &beta, &gamma );
- /* initialize the box */
- Setup_Box( box_x, box_y, box_z, alpha, beta, gamma, &(system->box) );
-
- /* count my atoms & allocate storage */
- Count_Geo_Atoms( geo, system );
- if ( PreAllocate_Space( system, control, workspace ) == FAILURE )
- {
- fprintf( stderr, "PreAllocate_Space: not enough memory!" );
- fprintf( stderr, "terminating...\n" );
- exit( INSUFFICIENT_MEMORY );
- }
-
- /* read in my atom info */
- top = 0;
- for ( i = 0; i < system->N; ++i )
- {
- fscanf( geo, CUSTOM_ATOM_FORMAT,
- &serial, element, name, &x[0], &x[1], &x[2] );
- Fit_to_Periodic_Box( &(system->box), &x );
-#if defined(DEBUG)
- fprintf( stderr, "atom%d: %s %s %f %f %f\n",
- serial, element, name, x[0], x[1], x[2] );
-#endif
-
- atom = &(system->atoms[top]);
- workspace->orig_id[i] = serial;
- atom->type = Get_Atom_Type( &(system->reaxprm), element );
- strcpy( atom->name, name );
- rvec_Copy( atom->x, x );
- rvec_MakeZero( atom->v );
- rvec_MakeZero( atom->f );
- atom->q = 0.;
-
- top++;
- }
-
- fclose( geo );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "finished reading the geo file\n" );
-#endif
-
- return SUCCESS;
-}
-
-
-int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
+static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
{
int ret;
char *cryst;
@@ -181,7 +112,64 @@ int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
}
-void Count_PDB_Atoms( FILE *geo, reax_system *system )
+void Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
+ simulation_data *data, static_storage *workspace )
+{
+
+ FILE *geo;
+ char descriptor[9];
+ int i, serial, top;
+ real box_x, box_y, box_z, alpha, beta, gamma;
+ rvec x;
+ char element[3], name[9];
+ reax_atom *atom;
+
+ /* open the geometry file */
+ if ( (geo = fopen(geo_file, "r")) == NULL )
+ {
+ fprintf( stderr, "Error opening the geo file! terminating...\n" );
+ exit( FILE_NOT_FOUND );
+ }
+
+ /* read box information */
+ fscanf( geo, CUSTOM_BOXGEO_FORMAT,
+ descriptor, &box_x, &box_y, &box_z, &alpha, &beta, &gamma );
+ /* initialize the box */
+ Setup_Box( box_x, box_y, box_z, alpha, beta, gamma, &(system->box) );
+
+ /* count my atoms & allocate storage */
+ Count_Geo_Atoms( geo, system );
+ PreAllocate_Space( system, control, workspace );
+
+ /* read in my atom info */
+ top = 0;
+ for ( i = 0; i < system->N; ++i )
+ {
+ fscanf( geo, CUSTOM_ATOM_FORMAT,
+ &serial, element, name, &x[0], &x[1], &x[2] );
+ Fit_to_Periodic_Box( &(system->box), &x );
+#if defined(DEBUG)
+ fprintf( stderr, "atom%d: %s %s %f %f %f\n",
+ serial, element, name, x[0], x[1], x[2] );
+#endif
+
+ atom = &(system->atoms[top]);
+ workspace->orig_id[i] = serial;
+ atom->type = Get_Atom_Type( &(system->reaxprm), element );
+ strcpy( atom->name, name );
+ rvec_Copy( atom->x, x );
+ rvec_MakeZero( atom->v );
+ rvec_MakeZero( atom->f );
+ atom->q = 0.;
+
+ top++;
+ }
+
+ fclose( geo );
+}
+
+
+static void Count_PDB_Atoms( FILE *geo, reax_system *system )
{
char *endptr = NULL;
char line[MAX_LINE + 1];
@@ -221,7 +209,7 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
}
-char Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
+void Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
simulation_data *data, static_storage *workspace )
{
@@ -258,18 +246,9 @@ char Read_PDB( const char * const pdb_file, reax_system* system, control_params
}
Count_PDB_Atoms( pdb, system );
- if ( PreAllocate_Space( system, control, workspace ) == FAILURE )
- {
- fprintf( stderr, "PreAllocate_Space: not enough memory!" );
- fprintf( stderr, "terminating...\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ PreAllocate_Space( system, control, workspace );
/* start reading and processing the pdb file */
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "starting to read the pdb file\n" );
-#endif
-
fseek( pdb, 0, SEEK_SET );
c = 0;
c1 = 0;
@@ -420,26 +399,21 @@ char Read_PDB( const char * const pdb_file, reax_system* system, control_params
}
if ( ferror( pdb ) )
{
- return FAILURE;
+ fprintf( stderr, "[ERROR] Unable to read PDB file. Terminating...\n" );
+ exit( INVALID_INPUT );
}
fclose( pdb );
Deallocate_Tokenizer_Space( &s, &s1, &tmp );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "finished reading the pdb file\n" );
-#endif
-
- return SUCCESS;
}
/* PDB serials are written without regard to the order, we'll see if this
- cause trouble, if so we'll have to rethink this approach
- Also, we do not write connect lines yet.
-*/
-char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
+ * cause trouble, if so we'll have to rethink this approach
+ * Also, we do not write connect lines yet.
+ * */
+void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
control_params *control, static_storage *workspace, output_controls *out_control )
{
int i, buffer_req, buffer_len;
@@ -476,7 +450,7 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
system->box.box[2][2] * system->box.box[1][2]) /
(system->box.box_norms[2] * system->box.box_norms[1]) );
- /*open pdb and write header*/
+ /* open pdb and write header */
snprintf( fname, MAX_STR + 9, "%s-%d.pdb", control->sim_name, data->step );
pdb = fopen(fname, "w");
fprintf( pdb, PDB_CRYST1_FORMAT_O,
@@ -487,7 +461,7 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
fprintf( out_control->log, "Box written\n" );
fflush( out_control->log );
- /*write atom lines to buffer*/
+ /* write atom lines to buffer */
for ( i = 0; i < system->N; i++)
{
p_atom = &(system->atoms[i]);
@@ -508,6 +482,12 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
buffer_len = system->N * PDB_ATOM_FORMAT_O_LENGTH;
buffer[buffer_len] = 0;
+
+ if ( ferror( pdb ) )
+ {
+ fprintf( stderr, "[ERROR] Unable to write PDB file. Terminating...\n" );
+ exit( INVALID_INPUT );
+ }
fprintf( pdb, "%s", buffer );
fclose( pdb );
@@ -533,12 +513,10 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
sfree( buffer, "Write_PDB::buffer" );
sfree( line, "Write_PDB::line" );
-
- return SUCCESS;
}
-char Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
+void Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
simulation_data *data, static_storage *workspace )
{
FILE *bgf;
@@ -584,17 +562,22 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
token_cnt = Tokenize( line, &tokens );
if ( !strcmp( tokens[0], "ATOM" ) || !strcmp( tokens[0], "HETATM" ) )
- (system->N)++;
+ {
+ ++system->N;
+ }
line[0] = 0;
}
- if ( ferror ( bgf ) )
+ if ( ferror( bgf ) )
{
- return FAILURE;
+ fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+ exit( INVALID_INPUT );
}
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "system->N: %d\n", system->N );
#endif
+
fclose( bgf );
/* memory allocations for atoms, atom maps, bond restrictions */
@@ -745,7 +728,7 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
{
/* check number of restrictions */
Check_Input_Range( token_cnt - 2, 0, MAX_RESTRICT,
- "CONECT line exceeds max restrictions allowed.\n" );
+ "CONECT line exceeds max restrictions allowed.\n" );
/* read bond restrictions */
if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[1]) ) )
@@ -777,16 +760,11 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
tokens[i][0] = 0;
}
}
- if ( ferror ( bgf ) )
+ if ( ferror( bgf ) )
{
- return FAILURE;
+ fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+ exit( INVALID_INPUT );
}
fclose( bgf );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "bgf file read\n" );
-#endif
-
- return SUCCESS;
}
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index 5e70721d8b52965170e7f098098c387be914be3b..61cec37a2563055994416eec53b49f045f4bb104 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -24,14 +24,12 @@
#include "mytypes.h"
+
// CUSTOM_BOXGEO: BOXGEO box_x box_y box_z angle1 angle2 angle3
#define CUSTOM_BOXGEO_FORMAT " %s %lf %lf %lf %lf %lf %lf"
// CUSTOM ATOM: serial element name x y z
#define CUSTOM_ATOM_FORMAT " %d %s %s %lf %lf %lf"
-char Read_Geo( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
-
/* PDB format :
http://www.rcsb.org/pdb/file_formats/pdb/pdbguide2.2/guide2.2_frame.html
@@ -117,13 +115,18 @@ COLUMNS DATA TYPE FIELD DEFINITION
#define BGF_CRYSTX_FORMAT "%8s%11s%11s%11s%11s%11s%11s"
-char Read_PDB( const char * const, reax_system*, control_params*,
+
+void Read_Geo( const char * const, reax_system*, control_params*,
simulation_data*, static_storage* );
-char Read_BGF( const char * const, reax_system*, control_params*,
+void Read_PDB( const char * const, reax_system*, control_params*,
simulation_data*, static_storage* );
-char Write_PDB( reax_system*, reax_list*, simulation_data*,
+void Read_BGF( const char * const, reax_system*, control_params*,
+ simulation_data*, static_storage* );
+
+void Write_PDB( reax_system*, reax_list*, simulation_data*,
control_params*, static_storage*, output_controls* );
+
#endif
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 9c04c3f287a5558fd82b230694aa5a6a9ef2e86a..b0749b06940cea6d5aa50a0493de2aa77d76fc05 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -26,7 +26,7 @@
#include "vector.h"
-int Estimate_GCell_Population( reax_system* system )
+static int Estimate_GCell_Population( reax_system* system )
{
int i, j, k, l;
int max_atoms;
@@ -66,7 +66,7 @@ int Estimate_GCell_Population( reax_system* system )
}
-void Allocate_Space_for_Grid( reax_system *system )
+static void Allocate_Space_for_Grid( reax_system *system )
{
int i, j, k, l;
grid *g;
@@ -172,7 +172,7 @@ void Allocate_Space_for_Grid( reax_system *system )
}
-void Deallocate_Grid_Space( grid *g )
+static void Deallocate_Grid_Space( grid *g )
{
int i, j, k;
@@ -216,7 +216,7 @@ void Deallocate_Grid_Space( grid *g )
}
-int Shift(int p, int dp, int dim, grid *g )
+int Shift( int p, int dp, int dim, grid *g )
{
int dim_len = 0;
int newp = p + dp;
@@ -248,8 +248,8 @@ int Shift(int p, int dp, int dim, grid *g )
/* finds the closest point between two grid cells denoted by c1 and c2.
periodic boundary conditions are taken into consideration as well. */
-void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
- int c2x, int c2y, int c2z, rvec closest_point )
+static void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
+ int c2x, int c2y, int c2z, rvec closest_point )
{
int i, d;
ivec c1 = { c1x, c1y, c1z };
@@ -294,7 +294,7 @@ void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
}
-void Find_Neighbor_GridCells( grid *g )
+static void Find_Neighbor_GridCells( grid *g )
{
int i, j, k;
int di, dj, dk;
@@ -543,7 +543,7 @@ static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
}
-void Copy_Storage( reax_system *system, static_storage *workspace,
+static void Copy_Storage( reax_system *system, static_storage *workspace,
control_params *control, int top, int old_id, int old_type, int *num_H,
real **v, real **s, real **t, int *orig_id, rvec *f_old )
{
@@ -578,7 +578,7 @@ void Copy_Storage( reax_system *system, static_storage *workspace,
}
-void Free_Storage( static_storage *workspace, control_params * control )
+static void Free_Storage( static_storage *workspace, control_params * control )
{
int i;
@@ -600,7 +600,7 @@ void Free_Storage( static_storage *workspace, control_params * control )
}
-void Assign_New_Storage( static_storage *workspace,
+static void Assign_New_Storage( static_storage *workspace,
real **v, real **s, real **t, int *orig_id, rvec *f_old )
{
workspace->v = v;
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index 6e83b740ed0a7f8035b45ad0423ef62e2ab0c850..af1cf0ace862dec750baddc7de7fb5cd89f1e2a3 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -25,19 +25,17 @@
#include "mytypes.h"
+int Shift( int, int, int, grid* );
+
void Setup_Grid( reax_system* );
void Update_Grid( reax_system* );
-void Finalize_Grid( reax_system* );
+void Bin_Atoms( reax_system*, static_storage* );
-int Shift( int, int, int, grid* );
+void Finalize_Grid( reax_system* );
void Cluster_Atoms( reax_system *, static_storage *, control_params * );
-void Bin_Atoms( reax_system*, static_storage* );
-
-void Reset_Marks( grid*, ivec*, int );
-
#endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index c809655fa4b9796f09122388657b05229f5f96ad..30e9ef4c26b8a30547272adf7e6b719a97d5d2a6 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -35,7 +35,7 @@
#include "vector.h"
-void Generate_Initial_Velocities( reax_system *system, real T )
+static void Generate_Initial_Velocities( reax_system *system, real T )
{
int i;
real scale, norm;
@@ -73,7 +73,7 @@ void Generate_Initial_Velocities( reax_system *system, real T )
}
-void Init_System( reax_system *system, control_params *control,
+static void Init_System( reax_system *system, control_params *control,
simulation_data *data )
{
int i;
@@ -128,7 +128,7 @@ void Init_System( reax_system *system, control_params *control,
}
-void Init_Simulation_Data( reax_system *system, control_params *control,
+static void Init_Simulation_Data( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control,
evolve_function *Evolve )
{
@@ -230,7 +230,7 @@ void Init_Simulation_Data( reax_system *system, control_params *control,
/* Initialize Taper params */
-void Init_Taper( control_params *control )
+static void Init_Taper( control_params *control )
{
real d1, d7;
real swa, swa2, swa3;
@@ -273,7 +273,7 @@ void Init_Taper( control_params *control )
}
-void Init_Workspace( reax_system *system, control_params *control,
+static void Init_Workspace( reax_system *system, control_params *control,
static_storage *workspace )
{
int i;
@@ -682,7 +682,7 @@ void Init_Workspace( reax_system *system, control_params *control,
}
-void Init_Lists( reax_system *system, control_params *control,
+static void Init_Lists( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
@@ -728,20 +728,12 @@ void Init_Lists( reax_system *system, control_params *control,
break;
}
- if ( Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz ) == FAILURE )
- {
- fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz );
/* TODO: better estimate for H_sp?
* If so, need to refactor Estimate_Storage_Sizes
* to use various cut-off distances as parameters
* (non-bonded, hydrogen, 3body, etc.) */
- if ( Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz ) == FAILURE )
- {
- fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
@@ -814,7 +806,7 @@ void Init_Lists( reax_system *system, control_params *control,
}
-void Init_Out_Controls( reax_system *system, control_params *control,
+static void Init_Out_Controls( reax_system *system, control_params *control,
static_storage *workspace, output_controls *out_control )
{
#define TEMP_SIZE (1000)
@@ -1105,7 +1097,7 @@ void Initialize( reax_system *system, control_params *control,
}
-void Finalize_System( reax_system *system, control_params *control,
+static void Finalize_System( reax_system *system, control_params *control,
simulation_data *data )
{
int i, j, k;
@@ -1147,13 +1139,13 @@ void Finalize_System( reax_system *system, control_params *control,
}
-void Finalize_Simulation_Data( reax_system *system, control_params *control,
+static void Finalize_Simulation_Data( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control )
{
}
-void Finalize_Workspace( reax_system *system, control_params *control,
+static void Finalize_Workspace( reax_system *system, control_params *control,
static_storage *workspace )
{
int i;
@@ -1369,7 +1361,7 @@ void Finalize_Workspace( reax_system *system, control_params *control,
}
-void Finalize_Lists( control_params *control, reax_list **lists )
+static void Finalize_Lists( control_params *control, reax_list **lists )
{
Delete_List( TYP_FAR_NEIGHBOR, &(*lists)[FAR_NBRS] );
if ( control->hb_cut > 0.0 )
@@ -1386,7 +1378,7 @@ void Finalize_Lists( control_params *control, reax_list **lists )
}
-void Finalize_Out_Controls( reax_system *system, control_params *control,
+static void Finalize_Out_Controls( reax_system *system, control_params *control,
static_storage *workspace, output_controls *out_control )
{
/* close trajectory file */
@@ -1523,7 +1515,6 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
}
#endif
-
#ifdef TEST_FORCES
/* close bond orders file */
if ( out_control->fbo )
@@ -1594,6 +1585,9 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
}
+/* Deallocate top-level data structures, close file handles, etc.
+ *
+ */
void Finalize( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control, const int output_enabled )
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index aff852c30006753fc38cd2890906a65892cf77ee..fbcc1667d79c9e504d53c37154106fd6587abf42 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -92,7 +92,6 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
}
-
void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params* control,
simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control, interaction_function *interaction_functions )
@@ -204,7 +203,6 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
}
-
/* uses Berendsen-type coupling for both T and P.
All box dimensions are scaled by the same amount,
there is no change in the angles between axes. */
@@ -220,6 +218,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
dt = control->dt;
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
+
#if defined(DEBUG_FOCUS)
//fprintf( out_control->prs,
// "tau_t: %g tau_p: %g dt/tau_t: %g dt/tau_p: %g\n",
@@ -291,9 +290,13 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
/* temperature scaler */
lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
if ( lambda < MIN_dT )
+ {
lambda = MIN_dT;
+ }
else if (lambda > MAX_dT )
+ {
lambda = MAX_dT;
+ }
lambda = SQRT( lambda );
/* Scale velocities and positions at t+dt */
@@ -446,13 +449,10 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
/************************************************/
#ifdef ANISOTROPIC
-
-void Velocity_Verlet_Nose_Hoover_NVT(reax_system* system,
- control_params* control,
- simulation_data *data,
- static_storage *workspace,
- reax_list **lists,
- output_controls *out_control )
+void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
+ control_params* control, simulation_data *data,
+ static_storage *workspace, reax_list **lists,
+ output_controls *out_control,
interaction_function *interaction_functions )
{
int i;
@@ -522,7 +522,6 @@ void Velocity_Verlet_Nose_Hoover_NVT(reax_system* system,
}
-
void Velocity_Verlet_Isotropic_NPT( reax_system* system,
control_params* control, simulation_data *data,
static_storage *workspace, reax_list **lists,
@@ -585,7 +584,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
//Update_Box_Isotropic( EXP( 3.0 * iso_bar->eps ), &(system->box) );
Update_Box_Isotropic( &(system->box), EXP( 3.0 * iso_bar->eps ) );
-
// Calculate new forces, f(t + dt)
Reset( system, control, data, workspace );
fprintf(out_control->log, "reset-");
@@ -604,7 +602,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
fprintf(out_control->log, "forces\n");
fflush( out_control->log );
-
// Compute iteration constants for each atom's velocity and for P_internal
// Compute kinetic energy for initial velocities of the iteration
P_int_const = E_kin = 0;
@@ -624,7 +621,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
rvec_Dot( system->atoms[i].v, system->atoms[i].v ) );
}
-
// Compute initial p_int
inv_3V = 1.0 / (3.0 * system->box.volume);
P_int = inv_3V * ( 2.0 * E_kin + P_int_const );
@@ -640,7 +636,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
v_xi_old = v_xi_new;
v_eps_old = v_eps_new;
-
for ( i = 0; i < system->N; ++i )
{
coef_v = 1.0 / (1.0 + 0.5 * dt * exp_deps *
@@ -655,13 +650,11 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
v_eps_new = coef_v_eps * ( iso_bar->v_eps +
0.5 * dt * ( iso_bar->a_eps + a_eps_new ) );
-
G_xi_new = control->Tau_T * ( 2.0 * E_kin +
SQR( v_eps_old ) / control->Tau_P -
(data->N_f + 1) * K_B * control->T );
v_xi_new = therm->v_xi + 0.5 * dt * ( therm->G_xi + G_xi_new );
-
E_kin = 0;
for ( i = 0; i < system->N; ++i )
E_kin += (0.5 * system->reaxprm.sbp[system->atoms[i].type].mass *
@@ -669,14 +662,12 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
P_int = inv_3V * ( 2.0 * E_kin + P_int_const );
-
fprintf( out_control->log,
"itr %d E_kin: %8.3f veps_n:%8.3f veps_o:%8.3f vxi_n:%8.3f vxi_o: %8.3f\n",
itr, E_kin, v_eps_new, v_eps_old, v_xi_new, v_xi_old );
}
while ( FABS(v_eps_new - v_eps_old) + FABS(v_xi_new - v_xi_old) > 2e-3 );
-
therm->v_xi_old = therm->v_xi;
therm->v_xi = v_xi_new;
therm->G_xi = G_xi_new;
@@ -696,6 +687,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
#endif
+
/* uses Berendsen-type coupling for both T and P.
All box dimensions are scaled by the same amount,
there is no change in the angles between axes. */
@@ -710,7 +702,6 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
rvec dx;
reax_atom *atom;
-
#if defined(DEBUG_FOCUS)
fprintf( stderr, "step%d\n", data->step );
#endif
@@ -764,9 +755,13 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
Compute_Kinetic_Energy( system, data );
lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
if ( lambda < MIN_dT )
+ {
lambda = MIN_dT;
+ }
else if (lambda > MAX_dT )
+ {
lambda = MAX_dT;
+ }
lambda = SQRT( lambda );
fprintf (stderr, "step:%d lambda -> %f \n", data->step, lambda);
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 11b7d9801d5867f42080116b7322e650ed0604aa..e0ea4d5f09c50d7a7f331332d4261e730460483c 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -56,11 +56,7 @@ static sparse_matrix * create_test_mat( void )
unsigned int i, n;
sparse_matrix *H_test;
- if ( Allocate_Matrix( &H_test, 3, 6 ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for test matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &H_test, 3, 6 );
//3x3, SPD, store lower half
i = 0;
@@ -173,27 +169,15 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
if ( *A_full == NULL )
{
- if ( Allocate_Matrix( A_full, A->n, 2 * A->m - A->n ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for full A. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
}
else if ( (*A_full)->m < 2 * A->m - A->n )
{
Deallocate_Matrix( *A_full );
- if ( Allocate_Matrix( A_full, A->n, 2 * A->m - A->n ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for full A. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
}
- if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for full A. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &A_t, A->n, A->m );
/* Set up the sparse matrix data structure for A. */
Transpose( A, A_t );
@@ -250,20 +234,12 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
if ( *A_spar_patt == NULL )
{
- if ( Allocate_Matrix( A_spar_patt, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "[SAI] Not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( A_spar_patt, A->n, A->m );
}
else if ( ((*A_spar_patt)->m) < (A->m) )
{
Deallocate_Matrix( *A_spar_patt );
- if ( Allocate_Matrix( A_spar_patt, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "[SAI] Not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( A_spar_patt, A->n, A->m );
}
@@ -271,7 +247,7 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
/* list: values from the matrix*/
/* left-right: search space of the quick-select */
- list = (real *) smalloc( sizeof(real) * (A->start[A->n]),"Sparse_Approx_Inverse::list" );
+ list = (real *) smalloc( sizeof(real) * (A->start[A->n]),"setup_sparse_approx_inverse::list" );
left = 0;
right = A->start[A->n] - 1;
@@ -366,12 +342,19 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
compute_full_sparse_matrix( A, A_full );
compute_full_sparse_matrix( *A_spar_patt, A_spar_patt_full );
- /* A_app_inv has the same sparsity pattern
- * * as A_spar_patt_full (omit non-zero values) */
- if ( Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m ) == FAILURE )
+ if ( *A_app_inv == NULL )
{
- fprintf( stderr, "not enough memory for approximate inverse matrix. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* A_app_inv has the same sparsity pattern
+ * * as A_spar_patt_full (omit non-zero values) */
+ Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
+ }
+ else if ( ((*A_app_inv)->m) < (A->m) )
+ {
+ Deallocate_Matrix( *A_app_inv );
+
+ /* A_app_inv has the same sparsity pattern
+ * * as A_spar_patt_full (omit non-zero values) */
+ Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
}
}
@@ -602,11 +585,7 @@ real SuperLU_Factorize( const sparse_matrix * const A,
"SuperLU_Factorize::Ltop" );
Utop = (unsigned int*) smalloc( (A->n + 1) * sizeof(unsigned int),
"SuperLU_Factorize::Utop" );
- if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &A_t, A->n, A->m );
/* Set up the sparse matrix data structure for A. */
Transpose( A, A_t );
@@ -1002,11 +981,7 @@ real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
D = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D" );
D_inv = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D_inv" );
Utop = (int*) smalloc( (A->n + 1) * sizeof(int), "ICHOL_PAR::Utop" );
- if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for ICHOL_PAR preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &DAD, A->n, A->m );
#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
@@ -1179,11 +1154,7 @@ real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
D = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D" );
D_inv = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D_inv" );
- if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "[ILU_PAR] Not enough memory for preconditioning matrices. Terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &DAD, A->n, A->m );
#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
@@ -1392,13 +1363,9 @@ real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
start = Get_Time( );
- if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE ||
- Allocate_Matrix( &L_temp, A->n, A->m ) == FAILURE ||
- Allocate_Matrix( &U_temp, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for ILUT_PAR preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &DAD, A->n, A->m );
+ Allocate_Matrix( &L_temp, A->n, A->m );
+ Allocate_Matrix( &U_temp, A->n, A->m );
D = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D" );
D_inv = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D_inv" );
@@ -1636,11 +1603,19 @@ real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
#if defined(HAVE_LAPACKE) || defined(HAVE_LAPACKE_MKL)
/* Compute M^{1} \approx A which minimizes
- * \sum_{j=1}^{N} ||e_j - Am_j||_2^2
- * where: e_j is the j-th column of the NxN identify matrix,
- * m_j is the j-th column of the NxN approximate sparse matrix M
+ * \sum_{j=1}^{N} ||m_j^T*A - e_j^T||_2^2
+ * where: e_j^T is the j-th row of the NxN identify matrix,
+ * m_j^T is the j-th row of the NxN approximate sparse matrix M
*
- * Internally, use LAPACKE to solve the least-squares problems
+ * Internally, use LAPACKE to solve the independent least-squares problems.
+ * Furthermore, internally solve the related problem
+ * \sum_{j=1}^{N} ||A*m_j - e_j||_2^2
+ * for j-th columns m_j and e_j, but store the transpose of M to solve
+ * the original problem. That is, for the problems
+ * min ||M_1*A - I||_F^2 and min ||A*M_2 - I||_F^2,
+ * it can be shown that if A = A^T, M_1 = M_2^T.
+ * Hence, we solve for M_2, and stores its transpose as the result
+ * (more efficient for for CSR matrices, row-major storage).
*
* A: symmetric, sparse matrix, stored in full CSR format
* A_spar_patt: sparse matrix used as template sparsity pattern
@@ -1676,10 +1651,10 @@ real sparse_approx_inverse( const sparse_matrix * const A,
shared(A_app_inv, stderr)
#endif
{
- X = (char *) smalloc( sizeof(char) * A->n, "Sparse_Approx_Inverse::X" );
- Y = (char *) smalloc( sizeof(char) * A->n, "Sparse_Approx_Inverse::Y" );
- pos_x = (int *) smalloc( sizeof(int) * A->n, "Sparse_Approx_Inverse::pos_x" );
- pos_y = (int *) smalloc( sizeof(int) * A->n, "Sparse_Approx_Inverse::pos_y" );
+ X = (char *) smalloc( sizeof(char) * A->n, "sparse_approx_inverse::X" );
+ Y = (char *) smalloc( sizeof(char) * A->n, "sparse_approx_inverse::Y" );
+ pos_x = (int *) smalloc( sizeof(int) * A->n, "sparse_approx_inverse::pos_x" );
+ pos_y = (int *) smalloc( sizeof(int) * A->n, "sparse_approx_inverse::pos_y" );
for ( i = 0; i < A->n; ++i )
{
@@ -1733,7 +1708,7 @@ real sparse_approx_inverse( const sparse_matrix * const A,
// allocate memory for NxM dense matrix
dense_matrix = (real *) smalloc( sizeof(real) * N * M,
- "Sparse_Approx_Inverse::dense_matrix" );
+ "sparse_approx_inverse::dense_matrix" );
// fill in the entries of dense matrix
for ( d_i = 0; d_i < M; ++d_i)
@@ -1758,7 +1733,7 @@ real sparse_approx_inverse( const sparse_matrix * const A,
/* create the right hand side of the linear equation
that is the full column of the identity matrix*/
e_j = (real *) smalloc( sizeof(real) * M,
- "Sparse_Approx_Inverse::e_j" );
+ "sparse_approx_inverse::e_j" );
for ( k = 0; k < M; ++k )
{
@@ -1797,8 +1772,8 @@ real sparse_approx_inverse( const sparse_matrix * const A,
}
//empty variables that will be used next iteration
- sfree( dense_matrix, "Sparse_Approx_Inverse::dense_matrix" );
- sfree( e_j, "Sparse_Approx_Inverse::e_j" );
+ sfree( dense_matrix, "sparse_approx_inverse::dense_matrix" );
+ sfree( e_j, "sparse_approx_inverse::e_j" );
for ( k = 0; k < A->n; ++k )
{
X[k] = 0;
@@ -1808,10 +1783,10 @@ real sparse_approx_inverse( const sparse_matrix * const A,
}
}
- sfree( pos_y, "Sparse_Approx_Inverse::pos_y" );
- sfree( pos_x, "Sparse_Approx_Inverse::pos_x" );
- sfree( Y, "Sparse_Approx_Inverse::Y" );
- sfree( X, "Sparse_Approx_Inverse::X" );
+ sfree( pos_y, "sparse_approx_inverse::pos_y" );
+ sfree( pos_x, "sparse_approx_inverse::pos_x" );
+ sfree( Y, "sparse_approx_inverse::Y" );
+ sfree( X, "sparse_approx_inverse::X" );
}
return Get_Timing_Info( start );
@@ -1961,11 +1936,7 @@ void Transpose_I( sparse_matrix * const A )
{
sparse_matrix * A_t;
- if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
- {
- fprintf( stderr, "not enough memory for transposing matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &A_t, A->n, A->m );
Transpose( A, A_t );
@@ -2490,11 +2461,7 @@ void permute_matrix( const static_storage * const workspace,
int i, pj, nr, nc;
sparse_matrix *LUtemp;
- if ( Allocate_Matrix( &LUtemp, LU->n, LU->m ) == FAILURE )
- {
- fprintf( stderr, "Not enough space for graph coloring (factor permutation). Terminating...\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( &LUtemp, LU->n, LU->m );
/* count nonzeros in each row of permuted factor (re-use color_top for counting) */
memset( workspace->color_top, 0, sizeof(unsigned int) * (LU->n + 1) );
@@ -2629,20 +2596,12 @@ void setup_graph_coloring( const control_params * const control,
{
if ( *H_p == NULL )
{
- if ( Allocate_Matrix( H_p, H->n, H->m ) == FAILURE )
- {
- fprintf( stderr, "[ERROR] not enough memory for graph coloring. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( H_p, H->n, H->m );
}
else if ( (*H_p)->m < H->m )
{
Deallocate_Matrix( *H_p );
- if ( Allocate_Matrix( H_p, H->n, H->m ) == FAILURE )
- {
- fprintf( stderr, "[ERROR] not enough memory for graph coloring. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ Allocate_Matrix( H_p, H->n, H->m );
}
compute_full_sparse_matrix( H, H_full );
@@ -3003,7 +2962,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
-/* generalized minimual residual method with restarting for sparse linear systems */
+/* Generalized minimual residual method with restarting and
+ * left preconditioning for sparse linear systems */
int GMRES( const static_storage * const workspace, const control_params * const control,
simulation_data * const data, const sparse_matrix * const H, const real * const b,
const real tol, real * const x, const int fresh_pre )
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index e6fbd2f4377e5edacc19f4c538101e1cc8e5efa4..fb7ddaed09ce0af278ee96c0adfeb61cf2bdd9bb 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -97,52 +97,28 @@ void Delete_List( int type, reax_list* l )
switch ( type )
{
case TYP_VOID:
- if ( l->select.v != NULL )
- {
- sfree( l->select.v, "Delete_List::l->select.v" );
- }
+ sfree( l->select.v, "Delete_List::l->select.v" );
break;
case TYP_THREE_BODY:
- if ( l->select.three_body_list != NULL )
- {
- sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
- }
+ sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
break;
case TYP_BOND:
- if ( l->select.bond_list != NULL )
- {
- sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
- }
+ sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
break;
case TYP_DBO:
- if ( l->select.dbo_list != NULL )
- {
- sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
- }
+ sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
break;
case TYP_DDELTA:
- if ( l->select.dDelta_list != NULL )
- {
- sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
- }
+ sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
break;
case TYP_FAR_NEIGHBOR:
- if ( l->select.far_nbr_list != NULL )
- {
- sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
- }
+ sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
break;
case TYP_NEAR_NEIGHBOR:
- if ( l->select.near_nbr_list != NULL )
- {
- sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
- }
+ sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
break;
case TYP_HBOND:
- if ( l->select.hbond_list != NULL )
- {
- sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
- }
+ sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
break;
default:
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index c10931436b73d78feb0be46d0b905028ff6e64e7..aba004d849c4457734af6690f2f5bc2d142cf8b3 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -52,7 +52,7 @@ void Make_Lookup_Table( real xmin, real xmax, int n,
/* Fills solution into x. Warning: will modify c and d! */
-void Tridiagonal_Solve( const real *a, const real *b,
+static void Tridiagonal_Solve( const real *a, const real *b,
real *c, real *d, real *x, unsigned int n)
{
int i;
@@ -77,7 +77,7 @@ void Tridiagonal_Solve( const real *a, const real *b,
}
-void Natural_Cubic_Spline( const real *h, const real *f,
+static void Natural_Cubic_Spline( const real *h, const real *f,
cubic_spline_coef *coef, unsigned int n )
{
int i;
@@ -151,8 +151,7 @@ void Natural_Cubic_Spline( const real *h, const real *f,
}
-
-void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
+static void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
cubic_spline_coef *coef, unsigned int n )
{
int i;
@@ -219,7 +218,7 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
}
-void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
+static void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
{
int i;
real base, dif;
@@ -502,7 +501,7 @@ void Finalize_LR_Lookup_Table( reax_system *system, control_params *control,
}
-int Lookup_Index_Of( real x, lookup_table* t )
+static int Lookup_Index_Of( real x, lookup_table* t )
{
return (int)( t->a * ( x - t->xmin ) );
}
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index 6fbd6ffafb921fc41ff17be571e1fb7892f6dac1..804226c89023ea0de605043f4835b7c3cfd12d50 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -25,14 +25,11 @@
#include "mytypes.h"
-/* Function pointer definitions */
typedef real (*lookup_function)(real);
void Make_Lookup_Table( real, real, int, lookup_function, lookup_table* );
-int Lookup_Index_Of( real, lookup_table* );
-
real Lookup( real, lookup_table* );
void Make_LR_Lookup_Table( reax_system*, control_params*,
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 61de424a6363a79a6ce51f993ebade2ca7d7abbd..a173ecca352c8460181947f859c9d17491f1c236 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -27,14 +27,14 @@
#include "config.h"
#endif
-#include "math.h"
-#include "random.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "sys/time.h"
-#include "time.h"
-#include "zlib.h"
+#include <math.h>
+#include <random.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/time.h>
+#include <time.h>
+#include <zlib.h>
#ifdef _OPENMP
#include <omp.h>
@@ -549,7 +549,6 @@ typedef struct
int periodic_boundaries;
int restrict_bonds;
int tabulate;
- ivec periodic_images;
real dt;
int reneighbor;
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index b63a7fca09a8249c2ab2cab71ffc7741f9df9200..1b9e38a61197f17a0fddc7873962e73b9c757eeb 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -259,7 +259,6 @@ int Estimate_NumNeighbors( reax_system *system, control_params *control,
#if defined DONE
-
void Choose_Neighbor_Finder( reax_system *system, control_params *control,
get_far_neighbors_function *Get_Far_Neighbors )
{
@@ -620,7 +619,6 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
}
-
void Generate_Neighbor_Lists( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
@@ -745,7 +743,4 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
system->N * far_nbrs->intrs_per_unit );
#endif
}
-
-
-
#endif
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index fde27ebe626616fcf8a2ab1842a4f41c81b695d7..4239433c4426eaf8d73c31995017ed8eaa74ee1c 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -24,9 +24,12 @@
#include "mytypes.h"
+
void Generate_Neighbor_Lists( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
int Estimate_NumNeighbors( reax_system*, control_params*,
- static_storage*, reax_list** );
+ static_storage*, reax_list** );
+
+
#endif
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 8c3ea3095354baddad4c763e69357de35132c5a6..447e7cc326a7b10e4f0c60aca94bfe16829c30be 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -506,7 +506,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
}
-int fn_qsort_intcmp( const void *a, const void *b )
+static int fn_qsort_intcmp( const void *a, const void *b )
{
return ( *(int *)a - * (int *)b);
}
@@ -577,36 +577,9 @@ void Output_Results( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
- int i, type_i;
- real e_pol, q, f_update;
+ real f_update;
real t_elapsed = 0;
- /* Compute Polarization Energy */
- e_pol = 0.0;
-
-#ifdef _OPENMP
- #pragma omp parallel for default(none) private(q, type_i) shared(system) \
- reduction(+: e_pol) schedule(static)
-#endif
- for ( i = 0; i < system->N; i++ )
- {
- q = system->atoms[i].q;
- type_i = system->atoms[i].type;
-
- e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
- (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
- KCALpMOL_to_EV;
- }
-
- data->E_Pol = e_pol;
-
- data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
- data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
- data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
-
-
- data->E_Tot = data->E_Pot + E_CONV * data->E_Kin;
-
/* output energies if it is the time */
if ( out_control->energy_update_freq > 0 &&
data->step % out_control->energy_update_freq == 0 )
@@ -724,28 +697,9 @@ void Output_Results( reax_system *system, control_params *control,
//t_elapsed = Get_Timing_Info( t_start );
//fprintf(stdout, "append_frame took %.6f seconds\n", t_elapsed );
}
-
- if ( IS_NAN_REAL(data->E_Pol) )
- {
- fprintf( stderr, "[ERROR] NaN detected for polarization energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
-
- if ( IS_NAN_REAL(data->E_Pot) )
- {
- fprintf( stderr, "[ERROR] NaN detected for potential energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
-
- if ( IS_NAN_REAL(data->E_Tot) )
- {
- fprintf( stderr, "[ERROR] NaN detected for total energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
}
-
void Print_Linear_System( reax_system *system, control_params *control,
static_storage *workspace, int step )
{
diff --git a/sPuReMD/src/print_utils.h b/sPuReMD/src/print_utils.h
index fe75d198d61efab104a25073f4b804ab9dba4300..de0fa4abf9fa009d4695fcf8c688acc0c00d4986 100644
--- a/sPuReMD/src/print_utils.h
+++ b/sPuReMD/src/print_utils.h
@@ -24,9 +24,6 @@
#include "mytypes.h"
-typedef void (*print_interaction)(reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls*);
-print_interaction Print_Interactions[NO_OF_INTERACTIONS];
char *Get_Element( reax_system*, int );
@@ -60,28 +57,39 @@ void Print_Sparse_Matrix2( sparse_matrix*, char*, char* );
void Print_Sparse_Matrix_Binary( sparse_matrix*, char* );
void Print_Bonds( reax_system*, reax_list*, char* );
+
void Print_Bond_List2( reax_system*, reax_list*, char* );
#ifdef TEST_FORCES
void Dummy_Printer( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Print_Bond_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Print_LonePair_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Print_OverUnderCoor_Forces(reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*);
+
void Print_Three_Body_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Print_Hydrogen_Bond_Forces(reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*);
+
void Print_Four_Body_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Print_vdW_Coulomb_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Compare_Total_Forces( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
+
void Init_Force_Test_Functions( );
#endif
+
#endif
diff --git a/sPuReMD/src/random.c b/sPuReMD/src/random.c
index 9b09e7526b7a8418470cbf8c1b45bd1940dcbfa9..4afe6aeac249822b2303a553bbf4938b08892f80 100644
--- a/sPuReMD/src/random.c
+++ b/sPuReMD/src/random.c
@@ -23,24 +23,26 @@
/* System random number generator used linear congruance method with
- large periodicity for generation of pseudo random number. function
- Random returns this random number appropriately scaled so that
- 0 <= Random(range) < range */
+ * large periodicity for generation of pseudo random number. function
+ * Random returns this random number appropriately scaled so that
+ * 0 <= Random(range) < range */
double Random(double range)
{
return (random() * range) / 2147483647L;
}
+
/* This function seeds the system pseudo random number generator with
- current time. Use this function once in the begining to initialize
- the system */
+ * current time. Use this function once in the begining to initialize
+ * the system */
void Randomize()
{
srandom(time(NULL));
}
+
/* GRandom return random number with gaussian distribution with mean
- and standard deviation "sigma" */
+ * and standard deviation "sigma" */
double GRandom(double mean, double sigma)
{
double v1 = Random(2.0) - 1.0;
diff --git a/sPuReMD/src/random.h b/sPuReMD/src/random.h
index a90b181dc78a9491944ffc4b666ba4201ae32e26..1fd1642f4e5a7044f898fcc003ded9d428e4d07b 100644
--- a/sPuReMD/src/random.h
+++ b/sPuReMD/src/random.h
@@ -22,21 +22,23 @@
#ifndef __RANDOM_H_
#define __RANDOM_H_
-#include <mytypes.h>
+#include "mytypes.h"
+
/* System random number generator used linear congruance method with
large periodicity for generation of pseudo random number. function
Random returns this random number appropriately scaled so that
0 <= Random(range) < range */
-double Random(double);
+double Random( double );
/* This function seeds the system pseudo random number generator with
current time. Use this function once in the begining to initialize
the system */
-void Randomize();
+void Randomize( );
/* GRandom return random number with gaussian distribution with mean
and standard deviation "sigma" */
-double GRandom(double, double);
+double GRandom( double, double );
+
#endif
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index 3edfe9072c0b635f7d3cfb97e555dd37d19d4a37..89f35972c4c8a56aec619ae2d38546a4b6d4500b 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -36,7 +36,7 @@ void Reset_Atoms( reax_system* system )
}
-void Reset_Pressures( simulation_data *data )
+static void Reset_Pressures( simulation_data *data )
{
rtensor_MakeZero( data->flex_bar.P );
data->iso_bar.P = 0.0;
@@ -115,7 +115,7 @@ void Reset_Workspace( reax_system *system, static_storage *workspace )
void Reset_Neighbor_Lists( reax_system *system, control_params *control,
- static_storage *workspace, reax_list **lists )
+ static_storage *workspace, reax_list **lists )
{
int i, tmp;
reax_list *bonds;
@@ -159,10 +159,6 @@ void Reset( reax_system *system, control_params *control,
Reset_Workspace( system, workspace );
Reset_Neighbor_Lists( system, control, workspace, lists );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "reset - ");
-#endif
}
diff --git a/sPuReMD/src/reset_utils.h b/sPuReMD/src/reset_utils.h
index 4355460b524403d28c4c7e9e6331377b0e6075a8..97598adbcfd27633b1c01d0129c0210b64bfef47 100644
--- a/sPuReMD/src/reset_utils.h
+++ b/sPuReMD/src/reset_utils.h
@@ -24,9 +24,8 @@
#include "mytypes.h"
-void Reset_Atoms( reax_system* );
-void Reset_Pressures( simulation_data* );
+void Reset_Atoms( reax_system* );
void Reset_Simulation_Data( simulation_data* );
@@ -42,10 +41,9 @@ void Reset_Neighbor_Lists( reax_system*, control_params*,
void Reset( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list** );
-//void Reset_Neighbor_Lists( reax_system*, static_storage*, reax_list** );
-
void Reset_Grid( grid* );
void Reset_Marks( grid*, ivec*, int );
+
#endif
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 4a9ff3102d5d86e33db27d71c19ed101a465c54f..cf22fb6a3b4e1bf4d1929babbd48e538213d5f75 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -92,7 +92,9 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
/* correction for C2 */
if ( system->reaxprm.gp.l[5] > 0.001 &&
!strcmp( system->reaxprm.sbp[type_i].name, "C" ) )
+ {
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
+ {
if ( i < bonds->select.bond_list[pj].nbr )
{
j = bonds->select.bond_list[pj].nbr;
@@ -128,9 +130,10 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
}
}
+ }
+ }
}
-
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
@@ -138,8 +141,13 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
/* over-coordination energy */
if ( sbp_i->mass > 21.0 )
+ {
dfvl = 0.0;
- else dfvl = 1.0; // only for 1st-row elements
+ }
+ else
+ {
+ dfvl = 1.0; // only for 1st-row elements
+ }
p_ovun2 = sbp_i->p_ovun2;
sum_ovun1 = 0;
@@ -178,7 +186,6 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
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;
p_ovun5 = sbp_i->p_ovun5;
@@ -190,7 +197,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
inv_exp_ovun8 = 1.0 / (1.0 + exp_ovun8);
data->E_Un += e_un =
- -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
+ -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
CEunder1 = inv_exp_ovun2n * ( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
p_ovun2 * e_un * exp_ovun2n);
@@ -211,7 +218,6 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
Add_dDelta( system, lists, i, CEunder3, workspace->f_un ); // UnCoor - 1st
#endif
-
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
pbond = &(bonds->select.bond_list[pj]);
@@ -220,23 +226,20 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
bo_ij = &(pbond->bo_data);
twbp = &(system->reaxprm.tbp[ type_i ][ type_j ]);
-
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->BO_pi + bo_ij->BO_pi2); // OvCoor - 3a
bo_ij->Cdbopi += CEover4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //OvCoor-3b
+ (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->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->BO_pi + bo_ij->BO_pi2); // UnCoor - 2a
bo_ij->Cdbopi += CEunder4 *
- (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
+ (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->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
#ifdef TEST_ENERGY
/* fprintf( out_control->eov, "%6d%23.15e%23.15e"
diff --git a/sPuReMD/src/single_body_interactions.h b/sPuReMD/src/single_body_interactions.h
index 8283c6243ad25902576e0eafce60d9266f502e21..f575453bde1e4dc3f5f9a502b2a710f30bf96de0 100644
--- a/sPuReMD/src/single_body_interactions.h
+++ b/sPuReMD/src/single_body_interactions.h
@@ -22,9 +22,12 @@
#ifndef __SINGLE_BODY_INTERACTIONS_H_
#define __SINGLE_BODY_INTERACTIONS_H_
-#include <mytypes.h>
+#include "mytypes.h"
+
void LonePair_OverUnder_Coordination_Energy( reax_system*, control_params*,
simulation_data*, static_storage*,
reax_list**, output_controls* );
+
+
#endif
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index f18e71d521d5e0195db3957bb29bfb760d4b94bb..cb3e7418a86f6fec84ae147d7c19b91e974da5d0 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -74,6 +74,8 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
rvec_ScaledAdd( system->atoms[i].v, -1., cross );
}
}
+
+ Compute_Total_Energy( system, data );
}
@@ -178,8 +180,6 @@ void* setup( const char * const geo_file, const char * const ffield_file,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
- //TODO: if errors detected, set handle to NULL to indicate failure
-
return (void*) spmd_handle;
}
@@ -235,11 +235,18 @@ int simulate( const void * const handle )
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
+ Compute_Total_Energy( spmd_handle->system, spmd_handle->data );
+
if ( spmd_handle->output_enabled == TRUE )
{
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
}
+ if ( spmd_handle->callback != NULL )
+ {
+ spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
+ spmd_handle->lists );
+ }
++spmd_handle->data->step;
//}
@@ -262,11 +269,13 @@ int simulate( const void * const handle )
{
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+
Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
}
steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
+
if ( steps && spmd_handle->out_control->restart_freq &&
steps % spmd_handle->out_control->restart_freq == 0 &&
spmd_handle->output_enabled == TRUE )
@@ -350,3 +359,21 @@ reax_atom* get_atoms( const void * const handle )
return atoms;
}
+
+
+int set_output_enabled( const void * const handle, const int enabled )
+{
+ int ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+ spmd_handle->output_enabled = enabled;
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 3a78242bd768e2674c727ac882f74b69b2a5e336..679aeaf3bdd3a9a14a44099994100d059fbee19e 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -22,11 +22,11 @@
#ifndef __SPUREMD_H_
#define __SPUREMD_H_
-#define SPUREMD_SUCCESS (0)
-#define SPUREMD_FAILURE (-1)
+#include "mytypes.h"
-#include "mytypes.h"
+#define SPUREMD_SUCCESS (0)
+#define SPUREMD_FAILURE (-1)
#ifdef __cplusplus
@@ -44,6 +44,8 @@ int cleanup( const void * const );
reax_atom* get_atoms( const void * const );
+int set_output_enabled( const void * const, const int );
+
#ifdef __cplusplus
}
#endif
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index fad302b4c08b46a8b9d92d18f9a2a24a73e66e2f..cff459fd897c8b46df06ff3d01ce93687365e37a 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "system_props.h"
+
#include "tool_box.h"
#include "vector.h"
@@ -209,6 +210,57 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
}
+void Compute_Total_Energy( reax_system* system, simulation_data* data )
+{
+ int i, type_i;
+ real e_pol, q;
+
+ /* Compute Polarization Energy */
+ e_pol = 0.0;
+
+#ifdef _OPENMP
+ #pragma omp parallel for default(none) private(q, type_i) shared(system) \
+ reduction(+: e_pol) schedule(static)
+#endif
+ for ( i = 0; i < system->N; i++ )
+ {
+ q = system->atoms[i].q;
+ type_i = system->atoms[i].type;
+
+ e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
+ (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
+ KCALpMOL_to_EV;
+ }
+
+ data->E_Pol = e_pol;
+
+ data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
+ data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
+ data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
+
+
+ data->E_Tot = data->E_Pot + E_CONV * data->E_Kin;
+
+ if ( IS_NAN_REAL(data->E_Pol) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for polarization energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+
+ if ( IS_NAN_REAL(data->E_Pot) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for potential energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+
+ if ( IS_NAN_REAL(data->E_Tot) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for total energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+}
+
+
/* IMPORTANT: This function assumes that current kinetic energy and
* the center of mass of the system is already computed before.
*
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index bf8f36e00f502fef36d695de684fecee52763b26..de65bf666b619732f56dc2b47049a57b268d7aa2 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -22,7 +22,7 @@
#ifndef __SYSTEM_PROP_H_
#define __SYSTEM_PROP_H_
-#include <mytypes.h>
+#include "mytypes.h"
void Temperature_Control( control_params*, simulation_data*, output_controls* );
@@ -33,9 +33,13 @@ void Compute_Center_of_Mass( reax_system*, simulation_data*, FILE* );
void Compute_Kinetic_Energy( reax_system*, simulation_data* );
-void Compute_Pressure( reax_system*, simulation_data*, static_storage* );
+void Compute_Total_Energy( reax_system*, simulation_data* );
void Compute_Pressure_Isotropic( reax_system*, control_params*, simulation_data*, output_controls* );
+void Compute_Pressure_Isotropic_Klein( reax_system*, simulation_data* );
+
+void Compute_Pressure( reax_system*, simulation_data*, static_storage* );
+
#endif
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 485be5fdf934f5e0d113a35569a44ed44a5fb14f..1ea18a872fc9d18d78c44507a6b2d5340c61d10d 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "three_body_interactions.h"
+
#include "bond_orders.h"
#include "list.h"
#include "lookup.h"
diff --git a/sPuReMD/src/three_body_interactions.h b/sPuReMD/src/three_body_interactions.h
index c654751e8de354d6cc158e00866626fc8acc519e..30fe5c1c5f4b7f0121ef7165a6fd6a4507a2d76d 100644
--- a/sPuReMD/src/three_body_interactions.h
+++ b/sPuReMD/src/three_body_interactions.h
@@ -24,14 +24,16 @@
#include "mytypes.h"
-void Three_Body_Interactions( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
-
-void Hydrogen_Bonds( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
void Calculate_Theta( rvec, real, rvec, real, real*, real* );
void Calculate_dCos_Theta( rvec, real, rvec, real, rvec*, rvec*, rvec* );
+void Three_Body_Interactions( reax_system*, control_params*, simulation_data*,
+ static_storage*, reax_list**, output_controls* );
+
+void Hydrogen_Bonds( reax_system*, control_params*, simulation_data*,
+ static_storage*, reax_list**, output_controls* );
+
+
#endif
diff --git a/sPuReMD/src/traj.h b/sPuReMD/src/traj.h
index e45921dffcf01f63a9453f9b0ddd3eb49fc0b1cc..738bc19333b244a586e8a4bb63d29c275f7915cf 100644
--- a/sPuReMD/src/traj.h
+++ b/sPuReMD/src/traj.h
@@ -23,7 +23,8 @@
#define __TRAJ_H__
#include "mytypes.h"
-#include "zlib.h"
+
+#include <zlib.h>
#define BLOCK_MARK "REAX_BLOCK_MARK "
@@ -74,24 +75,28 @@
#define SIZE_INFO_LINE3 "%-10d %-10d %-10d\n"
#define SIZE_INFO_LEN3 33
-enum ATOM_LINE_OPTS {OPT_NOATOM = 0, OPT_ATOM_BASIC = 4, OPT_ATOM_wF = 5,
- OPT_ATOM_wV = 6, OPT_ATOM_FULL = 7
- };
-enum BOND_LINE_OPTS {OPT_NOBOND, OPT_BOND_BASIC, OPT_BOND_FULL};
-enum ANGLE_LINE_OPTS {OPT_NOANGLE, OPT_ANGLE_BASIC};
-typedef struct
+enum ATOM_LINE_OPTS
{
- int no_of_sub_blocks;
- int size;
- char* buffer;
- struct block** sub_blocks;
-} block;
-
+ OPT_NOATOM = 0,
+ OPT_ATOM_BASIC = 4,
+ OPT_ATOM_wF = 5,
+ OPT_ATOM_wV = 6,
+ OPT_ATOM_FULL = 7,
+};
+
+enum BOND_LINE_OPTS
+{
+ OPT_NOBOND = 0,
+ OPT_BOND_BASIC = 1,
+ OPT_BOND_FULL = 2,
+};
-int Write_Block( gzFile, block* );
-int Read_Next_Block( gzFile, block*, int* );
-int Skip_Next_Block( gzFile, int*);
+enum ANGLE_LINE_OPTS
+{
+ OPT_NOANGLE = 0,
+ OPT_ANGLE_BASIC = 1,
+};
/*
@@ -141,12 +146,11 @@ int Skip_Next_Block( gzFile, int*);
No. of torsion entries (int)
Torsion info lines as per torsion format.
*/
-
-
int Write_Custom_Header( reax_system*, control_params*,
- static_storage*, output_controls* );
-int Write_xyz_Header ( reax_system*, control_params*,
- static_storage*, output_controls* );
+ static_storage*, output_controls* );
+
+int Write_xyz_Header( reax_system*, control_params*,
+ static_storage*, output_controls* );
/*
Write_Traj_Header( gzfile file,
@@ -157,34 +161,13 @@ int Write_xyz_Header ( reax_system*, control_params*,
*/
char Write_Traj_Header( FILE*, int, char**, char**, control_params* );
-
-/*
- Push_Traj_Frame(gzfile file,
- reax_system* system,
- control_params* control,
- simulation_data* data,
- static_storage* workspace,
- reax_list** lists,
- char** various flags);
-*/
-int Push_Traj_Frame( /*gzfile*/ FILE*, reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, char** );
-
-/*
- Append_Traj_Frame( gzfile file,
- reax_system* system,
- control_params* control,
- simulation_data* data,
- static_storage* workspace,
- reax_list** lists,
- char** various flags);
-*/
int Append_Custom_Frame( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
-int Append_xyz_Frame ( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
+int Append_xyz_Frame( reax_system*, control_params*, simulation_data*,
+ static_storage*, reax_list**, output_controls* );
void Read_Traj( output_controls*, char * );
+
#endif
diff --git a/sPuReMD/src/two_body_interactions.h b/sPuReMD/src/two_body_interactions.h
index b176b5ebca0f577dab2300d93dcad915d1de5f42..042ee246ce0171f5f2be170352804938acdcbc1a 100644
--- a/sPuReMD/src/two_body_interactions.h
+++ b/sPuReMD/src/two_body_interactions.h
@@ -22,13 +22,19 @@
#ifndef __TWO_BODY_INTERACTIONS_H_
#define __TWO_BODY_INTERACTIONS_H_
-#include <mytypes.h>
+#include "mytypes.h"
+
void Bond_Energy( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
+
void vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
+
void LR_vdW_Coulomb( reax_system*, control_params*, int, int, real, LR_data* );
+
void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+ static_storage*, reax_list**, output_controls* );
+
+
#endif
diff --git a/tools/driver.py b/tools/driver.py
index 12a4a0882981ef73d8b252dbb360505e830d35f4..59a8c4a8ba7b6eda19fe08b8dcb9e61fd18e38fc 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -227,6 +227,7 @@ class SimulationData(Structure):
class ReaxAtom(Structure):
_fields_ = [
+ ("type", c_int),
("name", c_char * 8),
("x", c_double * 3),
("v", c_double * 3),
@@ -264,12 +265,18 @@ if __name__ == '__main__':
setup_callback = lib.setup_callback
setup_callback.restype = c_int
+ set_output_enabled = lib.set_output_enabled
+ set_output_enabled.argtypes = [c_void_p, c_int]
+ set_output_enabled.restype = c_int
+
handle = setup(b"data/benchmarks/water/water_6540.pdb",
b"data/benchmarks/water/ffield.water",
b"environ/param.gpu.water")
ret = setup_callback(handle, CALLBACKFUNC(get_simulation_step_results))
+ ret = set_output_enabled(handle, c_int(0))
+
print("{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
ret = simulate(handle)