diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 907fb1a0b50859381961ec56eb8da6ab6af94f74..c959c7564a21d0f847d26d4b11b79941f518170c 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -158,7 +158,7 @@ static void Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *h
int i;
int *hb_top;
- hb_top = (int *) scalloc( n, sizeof(int), "Reallocate_HBonds_List::hb_top" );
+ hb_top = scalloc( n, sizeof(int), "Reallocate_HBonds_List::hb_top" );
for ( i = 0; i < n; ++i )
{
if ( h_index[i] >= 0 )
@@ -240,19 +240,19 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
reallocate_data *realloc;
grid *g;
- realloc = &(workspace->realloc);
- g = &(system->g);
+ realloc = &workspace->realloc;
+ g = &system->g;
if ( realloc->num_far > 0 && nbr_flag )
{
- Reallocate_Neighbor_List( &(*lists)[FAR_NBRS],
+ Reallocate_Neighbor_List( lists[FAR_NBRS],
system->N, realloc->num_far * SAFE_ZONE );
realloc->num_far = -1;
}
if ( realloc->Htop > 0 )
{
- Reallocate_Matrix( &(workspace->H), system->N_cm,
+ Reallocate_Matrix( &workspace->H, system->N_cm,
realloc->Htop * SAFE_ZONE );
realloc->Htop = -1;
@@ -265,30 +265,30 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
if ( control->hb_cut > 0.0 && realloc->hbonds > 0 )
{
Reallocate_HBonds_List( system->N, workspace->num_H, workspace->hbond_index,
- &(*lists)[HBONDS] );
+ lists[HBONDS] );
realloc->hbonds = -1;
}
num_bonds = est_3body = -1;
if ( realloc->bonds > 0 )
{
- Reallocate_Bonds_List( system->N, &(*lists)[BONDS], &num_bonds, &est_3body );
+ Reallocate_Bonds_List( system->N, lists[BONDS], &num_bonds, &est_3body );
realloc->bonds = -1;
realloc->num_3body = MAX( realloc->num_3body, est_3body );
}
if ( realloc->num_3body > 0 )
{
- Delete_List( TYP_THREE_BODY, &(*lists)[THREE_BODIES] );
+ Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
if ( num_bonds == -1 )
{
- num_bonds = (*lists)[BONDS].total_intrs;
+ num_bonds = lists[BONDS]->total_intrs;
}
realloc->num_3body *= SAFE_ZONE;
Make_List( num_bonds, realloc->num_3body,
- TYP_THREE_BODY, &(*lists)[THREE_BODIES] );
+ TYP_THREE_BODY, lists[THREE_BODIES] );
realloc->num_3body = -1;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "reallocating 3 bodies\n" );
@@ -314,8 +314,7 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
{
// reallocate g->atoms
sfree( g->atoms[i][j][k], "Reallocate::g->atoms[i][j][k]" );
- g->atoms[i][j][k] = (int*)
- scalloc( workspace->realloc.gcell_atoms, sizeof(int),
+ g->atoms[i][j][k] = scalloc( workspace->realloc.gcell_atoms, sizeof(int),
"Reallocate::g->atoms[i][j][k]" );
}
}
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index 10a401de7b0e6b959df0417e1b2ac72961c28fc2..ae5270d849b78b516a634d705b741f9b93c42d97 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -22,7 +22,7 @@
#ifndef __ALLOCATE_H_
#define __ALLOCATE_H_
-#include "mytypes.h"
+#include "reax_types.h"
void PreAllocate_Space( reax_system*, control_params*, static_storage* );
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index c2a6a6ba1e6bb24cf25f2882c3a0f1e5b9331f5b..eba15b977dfc1dba6ffc6133760ed093f95307b0 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -59,8 +59,8 @@ static void Copy_Bond_List( reax_system *system, control_params *control,
reax_list **lists )
{
int i, j, top_old;
- reax_list *new_bonds = &(*lists)[BONDS];
- reax_list *old_bonds = &(*lists)[OLD_BONDS];
+ reax_list *new_bonds = lists[BONDS];
+ reax_list *old_bonds = lists[OLD_BONDS];
for ( top_old = 0, i = 0; i < system->N; ++i )
{
@@ -89,8 +89,8 @@ static void Copy_Bond_List( reax_system *system, control_params *control,
static int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
{
int oldp, newp;
- reax_list *new_bonds = &(*lists)[BONDS];
- reax_list *old_bonds = &(*lists)[OLD_BONDS];
+ reax_list *new_bonds = lists[BONDS];
+ reax_list *old_bonds = lists[OLD_BONDS];
/*fprintf( stdout, "\n%d\nold_bonds:", atom );
for( oldp = Start_Index( atom, old_bonds );
@@ -224,8 +224,8 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
int *old_mark = workspace->old_mark;
int num_old, num_new;
char s[MAX_MOLECULE_SIZE * 10];
- reax_list *new_bonds = &(*lists)[BONDS];
- reax_list *old_bonds = &(*lists)[OLD_BONDS];
+ reax_list *new_bonds = lists[BONDS];
+ reax_list *old_bonds = lists[OLD_BONDS];
molecule old_molecules[20], new_molecules[20];
fprintf( fout, "molecular analysis @ %d\n", data->step );
@@ -578,8 +578,8 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
char fragments[MAX_FRAGMENT_TYPES][MAX_ATOM_TYPES];
int fragment_count[MAX_FRAGMENT_TYPES];
molecule m;
- reax_list *new_bonds = &(*lists)[BONDS];
-// reax_list *old_bonds = &(*lists)[OLD_BONDS];
+ reax_list *new_bonds = lists[BONDS];
+// reax_list *old_bonds = lists[OLD_BONDS];
/* fragment analysis */
fprintf( fout, "step%d fragments\n", data->step );
@@ -648,8 +648,8 @@ static void Analyze_Silica( reax_system *system, control_params *control,
int O_SI_O_count, SI_O_SI_count;
int si_coord[10], ox_coord[10];
real O_SI_O, SI_O_SI;
- reax_list *new_bonds = &(*lists)[BONDS];
- reax_list *thb_intrs = &(*lists)[THREE_BODIES];
+ reax_list *new_bonds = lists[BONDS];
+ reax_list *thb_intrs = lists[THREE_BODIES];
Analyze_Fragments( system, control, data, workspace, lists, fout, 0 );
diff --git a/sPuReMD/src/analyze.h b/sPuReMD/src/analyze.h
index c8328240315db071b996c4bf19a1bd10ab1fc097..dcc5a451a24a71f80ff275c2e3f26a19c5099574 100644
--- a/sPuReMD/src/analyze.h
+++ b/sPuReMD/src/analyze.h
@@ -22,7 +22,7 @@
#ifndef __ANALYZE_H_
#define __ANALYZE_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Analysis( reax_system*, control_params*, simulation_data*,
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index 977bcffc5749768ca39ef8aa28c623522a5926b4..98d26e046c246f0ea6259eb9ede779a3f3c84663 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -42,8 +42,8 @@ void Get_dBO( reax_system *system, reax_list **lists,
reax_list *dBOs;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -64,8 +64,8 @@ void Get_dBOpinpi2( reax_system *system, reax_list **lists,
dbond_data *dbo_k;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -86,8 +86,8 @@ void Add_dBO( reax_system *system, reax_list **lists,
reax_list *dBOs;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -110,8 +110,8 @@ void Add_dBOpinpi2( reax_system *system, reax_list **lists,
dbond_data *dbo_k;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -132,8 +132,8 @@ void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
reax_list *dBOs;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -154,8 +154,8 @@ void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
dbond_data *dbo_k;
int start_pj, end_pj, k;
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -174,9 +174,9 @@ void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v
reax_list *dDeltas;
int start, end, k;
- dDeltas = &((*lists)[DDELTA]);
- start = Start_Index(i, dDeltas);
- end = End_Index(i, dDeltas);
+ dDeltas = lists[DDELTA];
+ start = Start_Index( i, dDeltas );
+ end = End_Index( i, dDeltas );
for ( k = start; k < end; ++k )
{
@@ -191,9 +191,9 @@ void Add_dDelta_to_Forces( reax_system *system, reax_list **lists, int i, real C
reax_list *dDeltas;
int start, end, k;
- dDeltas = &((*lists)[DDELTA]);
- start = Start_Index(i, dDeltas);
- end = End_Index(i, dDeltas);
+ dDeltas = lists[DDELTA];
+ start = Start_Index( i, dDeltas );
+ end = End_Index( i, dDeltas );
for ( k = start; k < end; ++k )
{
@@ -213,26 +213,25 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
bond_order_data *bo_ij;
dbond_data *top_dbo;
- /* Initializations */
- bonds = &(*lists)[BONDS];
- dBOs = &(*lists)[DBO];
+ bonds = lists[BONDS];
+ dBOs = lists[DBO];
j = bonds->select.bond_list[pj].nbr;
- bo_ij = &(bonds->select.bond_list[pj].bo_data);
+ bo_ij = &bonds->select.bond_list[pj].bo_data;
start_i = Start_Index( i, bonds );
end_i = End_Index( i, bonds );
l = Start_Index( j, bonds );
end_j = End_Index( j, bonds );
- top_dbo = &(dBOs->select.dbo_list[ (*top) ]);
+ top_dbo = &dBOs->select.dbo_list[ *top ];
for ( k = start_i; k < end_i; ++k )
{
- nbr_k = &(bonds->select.bond_list[k]);
+ nbr_k = &bonds->select.bond_list[k];
for ( ; l < end_j && bonds->select.bond_list[l].nbr < nbr_k->nbr; ++l )
{
/* These are the neighbors of j which aren't in the neighbor_list of i
- Note that they might also include i! */
- nbr_l = &(bonds->select.bond_list[l]);
+ * Note that they might also include i! */
+ nbr_l = &bonds->select.bond_list[l];
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
@@ -275,7 +274,7 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
if ( l < end_j && bonds->select.bond_list[l].nbr == nbr_k->nbr )
{
/* This is a common neighbor of i and j. */
- nbr_l = &(bonds->select.bond_list[l]);
+ nbr_l = &bonds->select.bond_list[l];
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C3dbo, dBOp ); //dBO,3rd
@@ -302,14 +301,15 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C4dbopi2, dDeltap_self ); //4th
}
- ++(*top), ++top_dbo;
+ ++(*top);
+ ++top_dbo;
}
for ( ; l < end_j; ++l )
{
/* These are the remaining neighbors of j which are not in the
neighbor_list of i. Note that they might also include i!*/
- nbr_l = &(bonds->select.bond_list[l]);
+ nbr_l = &bonds->select.bond_list[l];
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
@@ -337,7 +337,8 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C3dbopi2, dDeltap_self );//3rd
}
- ++(*top), ++top_dbo;
+ ++(*top);
+ ++top_dbo;
}
}
#endif
@@ -358,18 +359,17 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
int tid = omp_get_thread_num( );
#endif
- /* Initializations */
- bonds = &(*lists)[BONDS];
- nbr_j = &(bonds->select.bond_list[pj]);
+ bonds = lists[BONDS];
+ nbr_j = &bonds->select.bond_list[pj];
j = nbr_j->nbr;
- bo_ij = &(nbr_j->bo_data);
- bo_ji = &(bonds->select.bond_list[ nbr_j->sym_index ].bo_data);
+ bo_ij = &nbr_j->bo_data;
+ bo_ji = &bonds->select.bond_list[ nbr_j->sym_index ].bo_data;
#ifdef _OPENMP
- f_i = &(workspace->f_local[tid * system->N + i]);
- f_j = &(workspace->f_local[tid * system->N + j]);
+ f_i = &workspace->f_local[tid * system->N + i];
+ f_j = &workspace->f_local[tid * system->N + j];
#else
- f_i = &(system->atoms[i].f);
- f_j = &(system->atoms[j].f);
+ f_i = &system->atoms[i].f;
+ f_j = &system->atoms[j].f;
#endif
coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo);
@@ -396,12 +396,12 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
************************************/
for ( pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk )
{
- nbr_k = &(bonds->select.bond_list[pk]);
+ nbr_k = &bonds->select.bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
- f_k = &(workspace->f_local[tid * system->N + k]);
+ f_k = &workspace->f_local[tid * system->N + k];
#else
- f_k = &(system->atoms[k].f);
+ f_k = &system->atoms[k].f;
#endif
rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp ); /*2nd,dBO*/
@@ -457,12 +457,12 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
***************************************************************************/
for ( pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk )
{
- nbr_k = &(bonds->select.bond_list[pk]);
+ nbr_k = &bonds->select.bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
- f_k = &(workspace->f_local[tid * system->N + k]);
+ f_k = &workspace->f_local[tid * system->N + k];
#else
- f_k = &(system->atoms[k].f);
+ f_k = &system->atoms[k].f;
#endif
rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp ); /*3rd,dBO*/
@@ -548,17 +548,17 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
#endif
/* Initializations */
- bonds = &(*lists)[BONDS];
- nbr_j = &(bonds->select.bond_list[pj]);
+ bonds = lists[BONDS];
+ nbr_j = &bonds->select.bond_list[pj];
j = nbr_j->nbr;
- bo_ij = &(nbr_j->bo_data);
- bo_ji = &(bonds->select.bond_list[ nbr_j->sym_index ].bo_data);
+ bo_ij = &nbr_j->bo_data;
+ bo_ji = &bonds->select.bond_list[ nbr_j->sym_index ].bo_data;
#ifdef _OPENMP
- f_i = &(workspace->f_local[tid * system->N + i]);
- f_j = &(workspace->f_local[tid * system->N + j]);
+ f_i = &workspace->f_local[tid * system->N + i];
+ f_j = &workspace->f_local[tid * system->N + j];
#else
- f_i = &(system->atoms[i].f);
- f_j = &(system->atoms[j].f);
+ f_i = &system->atoms[i].f;
+ f_j = &system->atoms[j].f;
#endif
coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo);
@@ -581,12 +581,12 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
for ( pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk )
{
- nbr_k = &(bonds->select.bond_list[pk]);
+ nbr_k = &bonds->select.bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
- f_k = &(workspace->f_local[tid * system->N + k]);
+ f_k = &workspace->f_local[tid * system->N + k];
#else
- f_k = &(system->atoms[k].f);
+ f_k = &system->atoms[k].f;
#endif
rvec_ScaledAdd( *f_k, -coef.C2dbo, nbr_k->bo_data.dBOp );
@@ -625,12 +625,12 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
for ( pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk )
{
- nbr_k = &(bonds->select.bond_list[pk]);
+ nbr_k = &bonds->select.bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
- f_k = &(workspace->f_local[tid * system->N + k]);
+ f_k = &workspace->f_local[tid * system->N + k];
#else
- f_k = &(system->atoms[k].f);
+ f_k = &system->atoms[k].f;
#endif
rvec_ScaledAdd( *f_k, -coef.C3dbo, nbr_k->bo_data.dBOp );
@@ -743,7 +743,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
p_lp1 = system->reaxprm.gp.l[15];
p_boc1 = system->reaxprm.gp.l[0];
p_boc2 = system->reaxprm.gp.l[1];
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
#ifdef _OPENMP
#pragma omp parallel default(shared)
@@ -773,8 +773,8 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
top_dbo = 0;
top_dDelta = 0;
- dDeltas = &(*lists)[DDELTA];
- dBOs = &(*lists)[DBO];
+ dDeltas = lists[DDELTA];
+ dBOs = lists[DBO];
#endif
/* Calculate Deltaprime, Deltaprime_boc values */
@@ -784,7 +784,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
- sbp_i = &(system->reaxprm.sbp[type_i]);
+ sbp_i = &system->reaxprm.sbp[type_i];
workspace->Deltap[i] = workspace->total_bond_order[i] - sbp_i->valency;
workspace->Deltap_boc[i] =
workspace->total_bond_order[i] - sbp_i->valency_val;
@@ -803,7 +803,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
- sbp_i = &(system->reaxprm.sbp[type_i]);
+ sbp_i = &system->reaxprm.sbp[type_i];
val_i = sbp_i->valency;
Deltap_i = workspace->Deltap[i];
Deltap_boc_i = workspace->Deltap_boc[i];
@@ -814,11 +814,11 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
{
j = bonds->select.bond_list[pj].nbr;
type_j = system->atoms[j].type;
- bo_ij = &( bonds->select.bond_list[pj].bo_data );
+ bo_ij = &bonds->select.bond_list[pj].bo_data;
if ( i < j )
{
- twbp = &( system->reaxprm.tbp[type_i][type_j] );
+ twbp = &system->reaxprm.tbp[type_i][type_j];
#ifdef TEST_FORCES
Set_Start_Index( pj, top_dbo, dBOs );
@@ -848,7 +848,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
bo_ij->C4dbopi2 = 0.0;
#ifdef TEST_FORCES
- pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
+ pdbo = &dBOs->select.dbo_list[ top_dbo ];
/* compute dBO_ij/dr_i */
pdbo->wrt = i;
@@ -910,7 +910,8 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
{
/* No overcoordination correction! */
f1 = 1.0;
- Cf1_ij = Cf1_ji = 0.0;
+ Cf1_ij = 0.0;
+ Cf1_ji = 0.0;
}
if ( twbp->v13cor >= 0.001 )
@@ -1018,7 +1019,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
#ifdef TEST_FORCES
Set_Start_Index( i, top_dDelta, dDeltas );
- ptop_dDelta = &( dDeltas->select.dDelta_list[top_dDelta] );
+ ptop_dDelta = &dDeltas->select.dDelta_list[top_dDelta];
for ( pj = start_i; pj < end_i; ++pj )
{
@@ -1029,7 +1030,8 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
ptop_dDelta->wrt = j;
rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[j] );
rvec_MakeZero( workspace->dDelta[j] );
- ++top_dDelta, ++ptop_dDelta;
+ ++top_dDelta;
+ ++ptop_dDelta;
}
start_j = Start_Index(j, bonds);
@@ -1042,7 +1044,8 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
ptop_dDelta->wrt = k;
rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[k] );
rvec_MakeZero( workspace->dDelta[k] );
- ++top_dDelta, ++ptop_dDelta;
+ ++top_dDelta;
+ ++ptop_dDelta;
}
}
}
@@ -1069,10 +1072,12 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
+
if ( type_i < 0 )
{
continue;
}
+
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
@@ -1080,6 +1085,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
{
j = bonds->select.bond_list[pj].nbr;
type_j = system->atoms[j].type;
+
if ( type_j < 0 )
{
continue;
@@ -1126,7 +1132,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( j = 0; j < system->N; ++j )
{
type_j = system->atoms[j].type;
- sbp_j = &(system->reaxprm.sbp[ type_j ]);
+ sbp_j = &system->reaxprm.sbp[ type_j ];
workspace->Delta[j] = workspace->total_bond_order[j] - sbp_j->valency;
workspace->Delta_e[j] = workspace->total_bond_order[j] - sbp_j->valency_e;
diff --git a/sPuReMD/src/bond_orders.h b/sPuReMD/src/bond_orders.h
index 37cd62428e92c822cb1453e122808db4d42b5ac8..d4c18b09849001bd84e70a1fb812f0fa2d85e8c2 100644
--- a/sPuReMD/src/bond_orders.h
+++ b/sPuReMD/src/bond_orders.h
@@ -22,7 +22,7 @@
#ifndef __BOND_ORDERS_H_
#define __BOND_ORDERS_H_
-#include "mytypes.h"
+#include "reax_types.h"
typedef struct
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index 8282c14cf972725dfb1c648322ddbd374955fb03..af3fbcce93d968b78d3ad3681e258ab829b51c3d 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -22,7 +22,7 @@
#ifndef __BOX_H__
#define __BOX_H__
-#include "mytypes.h"
+#include "reax_types.h"
/* Computes all the transformations,
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index 85ac44d32fca2a17955f1286a5546fd1a32a5074..72d3deb299a8964520454313ec1d12431011f3f7 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -22,7 +22,7 @@
#ifndef __CHARGES_H_
#define __CHARGES_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
diff --git a/sPuReMD/src/control.h b/sPuReMD/src/control.h
index da79fe527aa670fb62269a3b8b35ef5c7d09afe9..14fe54931dc25a61b30c0bc038ed7a496b673c19 100644
--- a/sPuReMD/src/control.h
+++ b/sPuReMD/src/control.h
@@ -22,7 +22,7 @@
#ifndef __CONTROL_H_
#define __CONTROL_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index ed0251a71735d9b7a23314c33c1bb34112efea76..f561da8ff57aa430f50e3e04b7b25d52ccd74ccd 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -63,7 +63,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
reax->gp.l = (real*) smalloc( sizeof(real) * n,
"Read_Force_Field::reax->gp-l" );
- /* see mytypes.h for mapping between l[i] and the lambdas used in ff */
+ /* see reax_types.h for mapping between l[i] and the lambdas used in ff */
for (i = 0; i < n; i++)
{
fgets(s, MAX_LINE, fp);
@@ -503,7 +503,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
}
/* 3-body parameters -
- * supports multi-well potentials (upto MAX_3BODY_PARAM in mytypes.h) */
+ * supports multi-well potentials (upto MAX_3BODY_PARAM in reax_types.h) */
/* clear entries first */
for ( i = 0; i < reax->num_atom_types; ++i )
{
@@ -572,7 +572,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
* correspond to any type of pair of atoms in 1 and 4
* position. However, explicit X-Y-Z-W takes precedence over the
* default description.
- * supports multi-well potentials (upto MAX_4BODY_PARAM in mytypes.h)
+ * supports multi-well potentials (upto MAX_4BODY_PARAM in reax_types.h)
* IMPORTANT: for now, directions on how to read multi-entries from ffield
* is not clear */
diff --git a/sPuReMD/src/ffield.h b/sPuReMD/src/ffield.h
index e2e88a10a6e7eee68402c85a5eba12eea2dc9023..611a705d7afc3bf823995b75efcf9f0e2778e23b 100644
--- a/sPuReMD/src/ffield.h
+++ b/sPuReMD/src/ffield.h
@@ -22,7 +22,7 @@
#ifndef __FFIELD_H_
#define __FFIELD_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Read_Force_Field( FILE*, reax_interaction* );
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 5b2f0a7d2a7ca29d85a74b86b1fda3db248a15d5..f28343f4e6199d376b8757b9c3e2df0b8b37033b 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -48,33 +48,32 @@ static void Dummy_Interaction( reax_system *system, control_params *control,
}
-void Init_Bonded_Force_Functions( control_params *control,
- interaction_function *interaction_functions )
+void Init_Bonded_Force_Functions( control_params *control )
{
- interaction_functions[0] = Calculate_Bond_Orders;
- interaction_functions[1] = Bond_Energy; //*/Dummy_Interaction;
- interaction_functions[2] = LonePair_OverUnder_Coordination_Energy;
+ control->intr_funcs[0] = &Calculate_Bond_Orders;
+ control->intr_funcs[1] = &Bond_Energy; //*/Dummy_Interaction;
+ control->intr_funcs[2] = &LonePair_OverUnder_Coordination_Energy;
//*/Dummy_Interaction;
- interaction_functions[3] = Three_Body_Interactions; //*/Dummy_Interaction;
- interaction_functions[4] = Four_Body_Interactions; //*/Dummy_Interaction;
+ control->intr_funcs[3] = &Three_Body_Interactions; //*/Dummy_Interaction;
+ control->intr_funcs[4] = &Four_Body_Interactions; //*/Dummy_Interaction;
if ( control->hb_cut > 0.0 )
{
- interaction_functions[5] = Hydrogen_Bonds; //*/Dummy_Interaction;
+ control->intr_funcs[5] = &Hydrogen_Bonds; //*/Dummy_Interaction;
}
else
{
- interaction_functions[5] = Dummy_Interaction;
+ control->intr_funcs[5] = &Dummy_Interaction;
}
- interaction_functions[6] = Dummy_Interaction; //empty
- interaction_functions[7] = Dummy_Interaction; //empty
- interaction_functions[8] = Dummy_Interaction; //empty
- interaction_functions[9] = Dummy_Interaction; //empty
+ control->intr_funcs[6] = &Dummy_Interaction; //empty
+ control->intr_funcs[7] = &Dummy_Interaction; //empty
+ control->intr_funcs[8] = &Dummy_Interaction; //empty
+ control->intr_funcs[9] = &Dummy_Interaction; //empty
}
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 )
+ output_controls *out_control )
{
int i;
@@ -111,7 +110,7 @@ static void Compute_Bonded_Forces( reax_system *system, control_params *control,
/* Implement all the function calls as function pointers */
for ( i = 0; i < NO_OF_INTERACTIONS; i++ )
{
- (interaction_functions[i])( system, control, data, workspace,
+ (control->intr_funcs[i])( system, control, data, workspace,
lists, out_control );
#ifdef TEST_FORCES
@@ -165,7 +164,7 @@ static void Compute_Total_Force( reax_system *system, control_params *control,
int i;
reax_list *bonds;
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
#ifdef _OPENMP
#pragma omp parallel default(shared)
@@ -218,8 +217,8 @@ static void Validate_Lists( static_storage *workspace, reax_list **lists, int st
int i, flag;
reax_list *bonds, *hbonds;
- bonds = &(*lists)[BONDS];
- hbonds = &(*lists)[HBONDS];
+ bonds = lists[BONDS];
+ hbonds = lists[HBONDS];
/* far neighbors */
if ( Htop > Hmax * DANGER_ZONE )
@@ -368,8 +367,8 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
static inline real Init_Charge_Matrix_Entry( reax_system *system,
- control_params *control, int i, int j,
- real r_ij, MATRIX_ENTRY_POSITION pos )
+ control_params *control, static_storage *workspace,
+ int i, int j, real r_ij, MATRIX_ENTRY_POSITION pos )
{
real Tap, dr3gamij_1, dr3gamij_3, ret;
@@ -383,13 +382,13 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
switch ( pos )
{
case OFF_DIAGONAL:
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
+ Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
+ Tap = Tap * r_ij + workspace->Tap[5];
+ Tap = Tap * r_ij + workspace->Tap[4];
+ Tap = Tap * r_ij + workspace->Tap[3];
+ Tap = Tap * r_ij + workspace->Tap[2];
+ Tap = Tap * r_ij + workspace->Tap[1];
+ Tap = Tap * r_ij + workspace->Tap[0];
/* shielding */
dr3gamij_1 = ( r_ij * r_ij * r_ij +
@@ -610,9 +609,9 @@ static void Init_Forces( reax_system *system, control_params *control,
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
- far_nbrs = &(*lists)[FAR_NBRS];
- bonds = &(*lists)[BONDS];
- hbonds = &(*lists)[HBONDS];
+ far_nbrs = lists[FAR_NBRS];
+ bonds = lists[BONDS];
+ hbonds = lists[HBONDS];
H = workspace->H;
H_sp = workspace->H_sp;
Htop = 0;
@@ -624,14 +623,14 @@ static void Init_Forces( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
- atom_i = &(system->atoms[i]);
- type_i = atom_i->type;
+ atom_i = &system->atoms[i];
+ type_i = atom_i->type;
start_i = Start_Index( i, far_nbrs );
end_i = End_Index( i, far_nbrs );
H->start[i] = Htop;
H_sp->start[i] = H_sp_top;
btop_i = End_Index( i, bonds );
- sbp_i = &(system->reaxprm.sbp[type_i]);
+ sbp_i = &system->reaxprm.sbp[type_i];
ihb = ihb_top = -1;
if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
@@ -641,9 +640,9 @@ static void Init_Forces( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
+ nbr_pj = &far_nbrs->select.far_nbr_list[pj];
j = nbr_pj->nbr;
- atom_j = &(system->atoms[j]);
+ atom_j = &system->atoms[j];
flag = 0;
flag_sp = 0;
@@ -678,12 +677,12 @@ static void Init_Forces( reax_system *system, control_params *control,
{
type_j = system->atoms[j].type;
r_ij = nbr_pj->d;
- sbp_j = &(system->reaxprm.sbp[type_j]);
- twbp = &(system->reaxprm.tbp[type_i][type_j]);
+ sbp_j = &system->reaxprm.sbp[type_j];
+ twbp = &system->reaxprm.tbp[type_i][type_j];
H->j[Htop] = j;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j,
- r_ij, OFF_DIAGONAL );
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, j, r_ij, OFF_DIAGONAL );
++Htop;
/* H_sp matrix entry */
@@ -764,9 +763,9 @@ static void Init_Forces( reax_system *system, control_params *control,
{
num_bonds += 2;
/****** bonds i-j and j-i ******/
- ibond = &( bonds->select.bond_list[btop_i] );
+ ibond = &bonds->select.bond_list[btop_i];
btop_j = End_Index( j, bonds );
- jbond = &(bonds->select.bond_list[btop_j]);
+ jbond = &bonds->select.bond_list[btop_j];
ibond->nbr = j;
jbond->nbr = i;
@@ -783,8 +782,8 @@ static void Init_Forces( reax_system *system, control_params *control,
++btop_i;
Set_End_Index( j, btop_j + 1, bonds );
- bo_ij = &( ibond->bo_data );
- bo_ji = &( jbond->bo_data );
+ bo_ij = &ibond->bo_data;
+ bo_ji = &jbond->bo_data;
bo_ji->BO = bo_ij->BO = BO;
bo_ji->BO_s = bo_ij->BO_s = BO_s;
bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
@@ -837,8 +836,8 @@ static void Init_Forces( reax_system *system, control_params *control,
/* diagonal entry */
H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, i,
- r_ij, DIAGONAL );
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control, workspace,
+ i, i, r_ij, DIAGONAL );
++Htop;
H_sp->j[H_sp_top] = i;
@@ -894,9 +893,9 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
- far_nbrs = &(*lists)[FAR_NBRS];
- bonds = &(*lists)[BONDS];
- hbonds = &(*lists)[HBONDS];
+ far_nbrs = lists[FAR_NBRS];
+ bonds = lists[BONDS];
+ hbonds = lists[HBONDS];
H = workspace->H;
H_sp = workspace->H_sp;
Htop = 0;
@@ -908,14 +907,14 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
- atom_i = &(system->atoms[i]);
- type_i = atom_i->type;
+ atom_i = &system->atoms[i];
+ type_i = atom_i->type;
start_i = Start_Index( i, far_nbrs );
end_i = End_Index( i, far_nbrs );
H->start[i] = Htop;
H_sp->start[i] = H_sp_top;
btop_i = End_Index( i, bonds );
- sbp_i = &(system->reaxprm.sbp[type_i]);
+ sbp_i = &system->reaxprm.sbp[type_i];
ihb = ihb_top = -1;
if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
@@ -925,9 +924,9 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
+ nbr_pj = &far_nbrs->select.far_nbr_list[pj];
j = nbr_pj->nbr;
- atom_j = &(system->atoms[j]);
+ atom_j = &system->atoms[j];
flag = 0;
flag_sp = 0;
@@ -947,7 +946,7 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
flag_sp = 0;
}
}
- else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &(system->box),
+ else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &system->box,
nbr_pj->dvec)) <= SQR(control->r_cut))
{
if ( nbr_pj->d <= SQR(control->r_sp_cut))
@@ -966,8 +965,8 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
twbp = &(system->reaxprm.tbp[type_i][type_j]);
H->j[Htop] = j;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j,
- r_ij, OFF_DIAGONAL );
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, j, r_ij, OFF_DIAGONAL );
++Htop;
/* H_sp matrix entry */
@@ -1048,9 +1047,9 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
{
num_bonds += 2;
/****** bonds i-j and j-i ******/
- ibond = &( bonds->select.bond_list[btop_i] );
+ ibond = &bonds->select.bond_list[btop_i];
btop_j = End_Index( j, bonds );
- jbond = &(bonds->select.bond_list[btop_j]);
+ jbond = &bonds->select.bond_list[btop_j];
ibond->nbr = j;
jbond->nbr = i;
@@ -1121,8 +1120,8 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
/* diagonal entry */
H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, i,
- r_ij, DIAGONAL );
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, i, r_ij, DIAGONAL );
++Htop;
H_sp->j[H_sp_top] = i;
@@ -1173,25 +1172,25 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
- far_nbrs = *lists + FAR_NBRS;
+ far_nbrs = lists[FAR_NBRS];
for ( i = 0; i < system->N; ++i )
{
- atom_i = &(system->atoms[i]);
- type_i = atom_i->type;
+ atom_i = &system->atoms[i];
+ type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
- sbp_i = &(system->reaxprm.sbp[type_i]);
+ end_i = End_Index(i, far_nbrs);
+ sbp_i = &system->reaxprm.sbp[type_i];
ihb = sbp_i->p_hbond;
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
+ nbr_pj = &far_nbrs->select.far_nbr_list[pj];
j = nbr_pj->nbr;
- atom_j = &(system->atoms[j]);
+ atom_j = &system->atoms[j];
type_j = atom_j->type;
- sbp_j = &(system->reaxprm.sbp[type_j]);
- twbp = &(system->reaxprm.tbp[type_i][type_j]);
+ sbp_j = &system->reaxprm.sbp[type_j];
+ twbp = &system->reaxprm.tbp[type_i][type_j];
if ( nbr_pj->d <= control->r_cut )
{
@@ -1202,10 +1201,15 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
nbr_pj->d <= control->hb_cut )
{
jhb = sbp_j->p_hbond;
+
if ( ihb == 1 && jhb == 2 )
+ {
++hb_top[i];
+ }
else if ( ihb == 2 && jhb == 1 )
+ {
++hb_top[j];
+ }
}
/* uncorrected bond orders */
@@ -1218,21 +1222,33 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
C12 = twbp->p_bo1 * POW( r_ij / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + control->bo_cut) * EXP( C12 );
}
- else BO_s = C12 = 0.0;
+ else
+ {
+ C12 = 0.0;
+ BO_s = 0.0;
+ }
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0)
{
C34 = twbp->p_bo3 * POW( r_ij / twbp->r_p, twbp->p_bo4 );
BO_pi = EXP( C34 );
}
- else BO_pi = C34 = 0.0;
+ else
+ {
+ C34 = 0.0;
+ BO_pi = 0.0;
+ }
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0)
{
C56 = twbp->p_bo5 * POW( r_ij / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = EXP( C56 );
}
- else BO_pi2 = C56 = 0.0;
+ else
+ {
+ C56 = 0.0;
+ BO_pi2 = 0.0;
+ }
/* Initially BO values are the uncorrected ones, page 1 */
BO = BO_s + BO_pi + BO_pi2;
@@ -1261,8 +1277,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
void Compute_Forces( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
- reax_list** lists, output_controls *out_control,
- interaction_function *interaction_functions )
+ reax_list** lists, output_controls *out_control )
{
real t_start, t_elapsed;
@@ -1279,8 +1294,7 @@ void Compute_Forces( reax_system *system, control_params *control,
data->timing.init_forces += t_elapsed;
t_start = Get_Time( );
- Compute_Bonded_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Bonded_Forces( system, control, data, workspace, lists, out_control );
t_elapsed = Get_Timing_Info( t_start );
data->timing.bonded += t_elapsed;
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index 01c789b90ffe299c2d6f8ecde4641d3de8ad66f1..78587a509d3594523231173d12e42dffe5c68c26 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -22,15 +22,13 @@
#ifndef __FORCES_H_
#define __FORCES_H_
-#include "mytypes.h"
+#include "reax_types.h"
-void Init_Bonded_Force_Functions( control_params*,
- interaction_function* );
+void Init_Bonded_Force_Functions( control_params* );
void Compute_Forces( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ static_storage*, reax_list**, output_controls* );
void Estimate_Storage_Sizes( reax_system*, control_params*, reax_list**,
int*, int*, int*, int* );
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index 2444c43f2ee0bbafba924c34a75d9af6a1e00f3e..23e244084a7a6a90ea131ee9e6c01060f3b43f23 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -176,8 +176,8 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
p_tor3 = system->reaxprm.gp.l[24];
p_tor4 = system->reaxprm.gp.l[25];
p_cot2 = system->reaxprm.gp.l[27];
- bonds = &(*lists)[BONDS];
- thb_intrs = &(*lists)[THREE_BODIES];
+ bonds = lists[BONDS];
+ thb_intrs = lists[THREE_BODIES];
e_tor_total = 0.0;
e_con_total = 0.0;
diff --git a/sPuReMD/src/four_body_interactions.h b/sPuReMD/src/four_body_interactions.h
index e7a8e64fc55616b3853f2ee3ab9a5d1fb064ac31..11b886ec27ca079b3f44f0d045980ac6fdb78964 100644
--- a/sPuReMD/src/four_body_interactions.h
+++ b/sPuReMD/src/four_body_interactions.h
@@ -22,7 +22,7 @@
#ifndef __FOUR_BODY_INTERACTIONS_H_
#define __FOUR_BODY_INTERACTIONS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Four_Body_Interactions( reax_system*, control_params*, simulation_data*,
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index 61cec37a2563055994416eec53b49f045f4bb104..880edafaa667f9ae3018fc547018d342af9b5bbd 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -22,7 +22,7 @@
#ifndef __GEO_TOOLS_H_
#define __GEO_TOOLS_H_
-#include "mytypes.h"
+#include "reax_types.h"
// CUSTOM_BOXGEO: BOXGEO box_x box_y box_z angle1 angle2 angle3
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index af1cf0ace862dec750baddc7de7fb5cd89f1e2a3..b834249c8383812d7dbacdc609a6ca80cc386033 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -22,7 +22,7 @@
#ifndef __GRID_H_
#define __GRID_H_
-#include "mytypes.h"
+#include "reax_types.h"
int Shift( int, int, int, grid* );
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 6670dfc06f23b20b3f2bbf718d26d288dcef611c..cbe31041b30db93d6497574e2d5096b43a209a79 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -229,7 +229,7 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
/* Initialize Taper params */
-static void Init_Taper( control_params *control )
+static void Init_Taper( control_params *control, static_storage *workspace )
{
real d1, d7;
real swa, swa2, swa3;
@@ -260,14 +260,14 @@ static void Init_Taper( control_params *control )
swb2 = SQR( swb );
swb3 = CUBE( swb );
- control->Tap7 = 20.0 / d7;
- control->Tap6 = -70.0 * (swa + swb) / d7;
- control->Tap5 = 84.0 * (swa2 + 3.0 * swa * swb + swb2) / d7;
- control->Tap4 = -35.0 * (swa3 + 9.0 * swa2 * swb + 9.0 * swa * swb2 + swb3 ) / d7;
- control->Tap3 = 140.0 * (swa3 * swb + 3.0 * swa2 * swb2 + swa * swb3 ) / d7;
- control->Tap2 = -210.0 * (swa3 * swb2 + swa2 * swb3) / d7;
- control->Tap1 = 140.0 * swa3 * swb3 / d7;
- control->Tap0 = (-35.0 * swa3 * swb2 * swb2 + 21.0 * swa2 * swb3 * swb2 +
+ workspace->Tap[7] = 20.0 / d7;
+ workspace->Tap[6] = -70.0 * (swa + swb) / d7;
+ workspace->Tap[5] = 84.0 * (swa2 + 3.0 * swa * swb + swb2) / d7;
+ workspace->Tap[4] = -35.0 * (swa3 + 9.0 * swa2 * swb + 9.0 * swa * swb2 + swb3 ) / d7;
+ workspace->Tap[3] = 140.0 * (swa3 * swb + 3.0 * swa2 * swb2 + swa * swb3 ) / d7;
+ workspace->Tap[2] = -210.0 * (swa3 * swb2 + swa2 * swb3) / d7;
+ workspace->Tap[1] = 140.0 * swa3 * swb3 / d7;
+ workspace->Tap[0] = (-35.0 * swa3 * swb2 * swb2 + 21.0 * swa2 * swb3 * swb2 +
7.0 * swa * swb3 * swb3 + swb3 * swb3 * swb ) / d7;
}
@@ -692,7 +692,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
Reset_Workspace( system, workspace );
/* Initialize Taper function */
- Init_Taper( control );
+ Init_Taper( control, workspace );
}
@@ -708,7 +708,7 @@ static void Init_Lists( reax_system *system, control_params *control,
num_nbrs = Estimate_NumNeighbors( system, control, workspace, lists );
- Make_List( system->N, num_nbrs, TYP_FAR_NEIGHBOR, &(*lists)[FAR_NBRS] );
+ Make_List( system->N, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "memory allocated: far_nbrs = %ldMB\n",
@@ -742,12 +742,12 @@ static void Init_Lists( reax_system *system, control_params *control,
break;
}
- Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz );
+ 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.) */
- Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz );
+ Allocate_Matrix( &workspace->H_sp, system->N_cm, max_nnz );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
@@ -779,7 +779,7 @@ static void Init_Lists( reax_system *system, control_params *control,
else
{
Allocate_HBond_List( system->N, workspace->num_H, workspace->hbond_index,
- hb_top, &(*lists)[HBONDS] );
+ hb_top, lists[HBONDS] );
}
#if defined(DEBUG_FOCUS)
@@ -791,7 +791,7 @@ static void Init_Lists( reax_system *system, control_params *control,
}
/* bonds list */
- Allocate_Bond_List( system->N, bond_top, &(*lists)[BONDS] );
+ Allocate_Bond_List( system->N, bond_top, lists[BONDS] );
num_bonds = bond_top[system->N - 1];
#if defined(DEBUG_FOCUS)
@@ -801,7 +801,7 @@ static void Init_Lists( reax_system *system, control_params *control,
#endif
/* 3bodies list */
- Make_List( num_bonds, num_3body, TYP_THREE_BODY, &(*lists)[THREE_BODIES] );
+ Make_List( num_bonds, num_3body, TYP_THREE_BODY, lists[THREE_BODIES] );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "estimated storage - num_3body: %d\n", num_3body );
@@ -810,9 +810,9 @@ static void Init_Lists( reax_system *system, control_params *control,
#endif
#ifdef TEST_FORCES
- Make_List( system->N, num_bonds * 8, TYP_DDELTA, &(*lists)[DDELTA] );
+ Make_List( system->N, num_bonds * 8, TYP_DDELTA, lists[DDELTA] );
- Make_List( num_bonds, num_bonds * MAX_BONDS * 3, TYP_DBO, &(*lists)[DBO] );
+ Make_List( num_bonds, num_bonds * MAX_BONDS * 3, TYP_DBO, lists[DBO] );
#endif
sfree( hb_top, "Init_Lists::hb_top" );
@@ -1079,7 +1079,7 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
void Initialize( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control, evolve_function *Evolve,
- interaction_function *interaction_functions, const int output_enabled )
+ const int output_enabled )
{
#if defined(DEBUG)
real start, end;
@@ -1111,7 +1111,7 @@ void Initialize( reax_system *system, control_params *control,
}
/* These are done in forces.c, only forces.c can see all those functions */
- Init_Bonded_Force_Functions( control, interaction_functions );
+ Init_Bonded_Force_Functions( control );
#ifdef TEST_FORCES
Init_Force_Test_Functions( );
@@ -1412,17 +1412,17 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
static void Finalize_Lists( control_params *control, reax_list **lists )
{
- Delete_List( TYP_FAR_NEIGHBOR, &(*lists)[FAR_NBRS] );
+ Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
if ( control->hb_cut > 0.0 )
{
- Delete_List( TYP_HBOND, &(*lists)[HBONDS] );
+ Delete_List( TYP_HBOND, lists[HBONDS] );
}
- Delete_List( TYP_BOND, &(*lists)[BONDS] );
- Delete_List( TYP_THREE_BODY, &(*lists)[THREE_BODIES] );
+ Delete_List( TYP_BOND, lists[BONDS] );
+ Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
#ifdef TEST_FORCES
- Delete_List( TYP_DDELTA, &(*lists)[DDELTA] );
- Delete_List( TYP_DBO, &(*lists)[DBO] );
+ Delete_List( TYP_DDELTA, lists[DDELTA] );
+ Delete_List( TYP_DBO, lists[DBO] );
#endif
}
diff --git a/sPuReMD/src/init_md.h b/sPuReMD/src/init_md.h
index aff4d01d1cbed94e6e271b50ac962ac377891a4e..98d5cc90902422075a705c71de6508df5c8a36e2 100644
--- a/sPuReMD/src/init_md.h
+++ b/sPuReMD/src/init_md.h
@@ -22,12 +22,12 @@
#ifndef __INIT_MD_H_
#define __INIT_MD_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Initialize( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*, evolve_function*,
- interaction_function *, const int );
+ const int );
void Finalize( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*, const int );
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index af2d87bd2676d8dd705c5db5be254e846e5b43e3..cc89869966bc47ef8a139db7ab440f2f8e3c65eb 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -37,8 +37,7 @@
void Velocity_Verlet_NVE(reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control,
- interaction_function *interaction_functions )
+ reax_list **lists, output_controls *out_control )
{
int i, steps, renbr;
real inv_m, dt, dt_sqr;
@@ -76,8 +75,7 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
Generate_Neighbor_Lists( system, control, data, workspace,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Forces( system, control, data, workspace, lists, out_control );
for ( i = 0; i < system->N; i++ )
{
@@ -94,7 +92,7 @@ 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 )
+ output_controls *out_control )
{
int i, itr, steps, renbr;
real inv_m, coef_v, dt, dt_sqr;
@@ -139,8 +137,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
lists, out_control );
}
/* Calculate Forces at time (t + dt) */
- Compute_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Forces( system, control, data, workspace, lists, out_control );
/* Compute iteration constants for each atom's velocity */
for ( i = 0; i < system->N; ++i )
@@ -208,8 +205,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
there is no change in the angles between axes. */
void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
control_params* control, simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control,
- interaction_function *interaction_functions )
+ reax_list **lists, output_controls *out_control )
{
int i, steps, renbr;
real inv_m, dt, lambda, mu;
@@ -243,8 +239,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
Generate_Neighbor_Lists( system, control, data, workspace,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Forces( system, control, data, workspace, lists, out_control );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
@@ -308,8 +303,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
there is no change in the angles between axes. */
void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
control_params* control, simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control,
- interaction_function *interaction_functions )
+ reax_list **lists, output_controls *out_control )
{
int i, d, steps, renbr;
real dt, inv_m, lambda;
@@ -349,8 +343,7 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
Generate_Neighbor_Lists( system, control, data, workspace,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Forces( system, control, data, workspace, lists, out_control );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
@@ -419,8 +412,7 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
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 )
+ output_controls *out_control )
{
int i;
real inv_m;
@@ -461,8 +453,7 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
/* Compute_Charges( system, control, workspace, out_control );
fprintf(out_control->log,"qeq-"); fflush( out_control->log ); */
- Compute_Forces( system, control, data, workspace, lists, out_control,
- interaction_functions );
+ Compute_Forces( system, control, data, workspace, lists, out_control );
fprintf(out_control->log, "forces\n");
fflush( out_control->log );
@@ -565,7 +556,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
fprintf(out_control->log,"qeq-"); fflush( out_control->log ); */
Compute_Forces( system, control, data, workspace, lists,
- out_control, interaction_functions );
+ out_control );
fprintf(out_control->log, "forces\n");
fflush( out_control->log );
@@ -661,8 +652,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
void Velocity_Verlet_Berendsen_NVT( reax_system* system,
control_params* control, simulation_data *data,
static_storage *workspace, reax_list **lists,
- output_controls *out_control,
- interaction_function *interaction_functions )
+ output_controls *out_control )
{
int i, steps, renbr;
real inv_m, dt, lambda;
@@ -707,7 +697,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
}
Compute_Forces( system, control, data, workspace,
- lists, out_control, interaction_functions );
+ lists, out_control );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
diff --git a/sPuReMD/src/integrate.h b/sPuReMD/src/integrate.h
index 8b847601c012dda8d6192b1d116a5531870e1f3e..d5f0f3ae65853e8030ae441f4f88625d07382793 100644
--- a/sPuReMD/src/integrate.h
+++ b/sPuReMD/src/integrate.h
@@ -23,40 +23,32 @@
#define __INTEGRATE_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Velocity_Verlet_NVE( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Nose_Hoover_NVT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Flexible_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Isotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*,
- interaction_function* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
void Velocity_Verlet_Berendsen_NVT( reax_system* , control_params* ,
- simulation_data *, static_storage *, reax_list **, output_controls *,
- interaction_function* );
+ simulation_data *, static_storage *, reax_list **, output_controls * );
#endif
diff --git a/sPuReMD/src/lin_alg.h b/sPuReMD/src/lin_alg.h
index 696fdc9bd5401b9cf1442b9499b8c99e4940048c..429db79cc649198c7da8b157315ffa287e19ddb9 100644
--- a/sPuReMD/src/lin_alg.h
+++ b/sPuReMD/src/lin_alg.h
@@ -22,7 +22,7 @@
#ifndef __LIN_ALG_H_
#define __LIN_ALG_H_
-#include "mytypes.h"
+#include "reax_types.h"
typedef enum
{
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index fb7ddaed09ce0af278ee96c0adfeb61cf2bdd9bb..85805eb9998ee8038f02877766d6e5b7cef955d4 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -29,60 +29,53 @@ void Make_List( int n, int total_intrs, int type, reax_list* l )
l->n = n;
l->total_intrs = total_intrs;
- l->index = (int*) smalloc( n * sizeof(int), "Make_List::l->index" );
- l->end_index = (int*) smalloc( n * sizeof(int), "Make_List::l->end_index" );
+ l->index = smalloc( n * sizeof(int), "Make_List::l->index" );
+ l->end_index = smalloc( n * sizeof(int), "Make_List::l->end_index" );
switch ( type )
{
case TYP_VOID:
- l->select.v = (void *) smalloc( l->total_intrs * sizeof(void),
+ l->select.v = smalloc( l->total_intrs * sizeof(void),
"Make_List::l->select.v" );
break;
case TYP_THREE_BODY:
- l->select.three_body_list = (three_body_interaction_data*)
- smalloc( l->total_intrs * sizeof(three_body_interaction_data),
- "Make_List::l->select.three_body_list" );
+ l->select.three_body_list = smalloc( l->total_intrs * sizeof(three_body_interaction_data),
+ "Make_List::l->select.three_body_list" );
break;
case TYP_BOND:
- l->select.bond_list = (bond_data*)
- smalloc( l->total_intrs * sizeof(bond_data),
- "Make_List::l->select.bond_list" );
+ l->select.bond_list = smalloc( l->total_intrs * sizeof(bond_data),
+ "Make_List::l->select.bond_list" );
break;
case TYP_DBO:
- l->select.dbo_list = (dbond_data*)
- smalloc( l->total_intrs * sizeof(dbond_data),
- "Make_List::l->select.dbo_list" );
+ l->select.dbo_list = smalloc( l->total_intrs * sizeof(dbond_data),
+ "Make_List::l->select.dbo_list" );
break;
case TYP_DDELTA:
- l->select.dDelta_list = (dDelta_data*)
- smalloc( l->total_intrs * sizeof(dDelta_data),
- "Make_List::l->select.dDelta_list" );
+ l->select.dDelta_list = smalloc( l->total_intrs * sizeof(dDelta_data),
+ "Make_List::l->select.dDelta_list" );
break;
case TYP_FAR_NEIGHBOR:
- l->select.far_nbr_list = (far_neighbor_data*)
- smalloc( l->total_intrs * sizeof(far_neighbor_data),
- "Make_List::l->select.far_nbr_list" );
+ l->select.far_nbr_list = smalloc( l->total_intrs * sizeof(far_neighbor_data),
+ "Make_List::l->select.far_nbr_list" );
break;
case TYP_NEAR_NEIGHBOR:
- l->select.near_nbr_list = (near_neighbor_data*)
- smalloc( l->total_intrs * sizeof(near_neighbor_data),
- "Make_List::l->select.near_nbr_list" );
+ l->select.near_nbr_list = smalloc( l->total_intrs * sizeof(near_neighbor_data),
+ "Make_List::l->select.near_nbr_list" );
break;
case TYP_HBOND:
- l->select.hbond_list = (hbond_data*)
- smalloc( l->total_intrs * sizeof(hbond_data),
- "Make_List::l->select.hbond_list" );
+ l->select.hbond_list = smalloc( l->total_intrs * sizeof(hbond_data),
+ "Make_List::l->select.hbond_list" );
break;
default:
- l->select.v = (void *) smalloc( l->total_intrs * sizeof(void),
+ l->select.v = smalloc( l->total_intrs * sizeof(void),
"Make_List::l->select.v" );
break;
}
@@ -99,24 +92,31 @@ void Delete_List( int type, reax_list* l )
case TYP_VOID:
sfree( l->select.v, "Delete_List::l->select.v" );
break;
+
case TYP_THREE_BODY:
sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
break;
+
case TYP_BOND:
sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
break;
+
case TYP_DBO:
sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
break;
+
case TYP_DDELTA:
sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
break;
+
case TYP_FAR_NEIGHBOR:
sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
break;
+
case TYP_NEAR_NEIGHBOR:
sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
break;
+
case TYP_HBOND:
sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
break;
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index 9cd53cbf65291a35e285b249d0eabe98afc2850f..66464e97ed6d5e76875ee7af71878302135edf95 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -22,7 +22,7 @@
#ifndef __LIST_H_
#define __LIST_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Make_List( int, int, int, reax_list* );
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index aba004d849c4457734af6690f2f5bc2d142cf8b3..267c3542e2b0acc899282998c1454d202a7762c8 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -342,8 +342,8 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
for ( r = 1; r <= control->tabulate; ++r )
{
- LR_vdW_Coulomb( system, control, i, j, r * dr,
- &(workspace->LR[i][j].y[r]) );
+ LR_vdW_Coulomb( system, control, workspace,
+ i, j, r * dr, &workspace->LR[i][j].y[r] );
h[r] = workspace->LR[i][j].dx;
fh[r] = workspace->LR[i][j].y[r].H;
fvdw[r] = workspace->LR[i][j].y[r].e_vdW;
@@ -410,7 +410,7 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
if( existing_types[j] ) {
for( r = 1; r <= 100; ++r ) {
rand_dist = (real)rand()/RAND_MAX * control->r_cut;
- LR_vdW_Coulomb( system, control, i, j, rand_dist, &y );
+ LR_vdW_Coulomb( system, control, workspace, i, j, rand_dist, &y );
LR_Lookup( &(workspace->LR[i][j]), rand_dist, &y_spline );
evdw_abserr = FABS(y.e_vdW - y_spline.e_vdW);
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index 804226c89023ea0de605043f4835b7c3cfd12d50..8a63a1a9e3f72f61549b1420d9818ad00e6bec09 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -22,7 +22,7 @@
#ifndef __LOOKUP_H_
#define __LOOKUP_H_
-#include "mytypes.h"
+#include "reax_types.h"
typedef real (*lookup_function)(real);
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 1b9e38a61197f17a0fddc7873962e73b9c757eeb..441cb23838c736978031f093140c2102a7426c43 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -71,7 +71,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
t_start = Get_Time( );
// fprintf( stderr, "\n\tentered nbrs - " );
g = &( system->g );
- far_nbrs = &(*lists)[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
Bin_Atoms( system, workspace );
// fprintf( stderr, "atoms sorted - " );
@@ -387,7 +387,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
far_neighbor_data new_nbrs[125];
g = &( system->g );
- far_nbrs = &(*lists)[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
// fprintf( stderr, "\n\tentered nbrs - " );
if ( control->ensemble == iNPT || control->ensemble == sNPT ||
@@ -636,7 +636,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
far_neighbor_data new_nbrs[125];
g = &( system->g );
- far_nbrs = &(*lists)[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
// fprintf( stderr, "\n\tentered nbrs - " );
if ( control->ensemble == iNPT ||
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index 4239433c4426eaf8d73c31995017ed8eaa74ee1c..b2d2c0527bb2a5d1c2392487a3735d731a42168d 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -22,7 +22,7 @@
#ifndef __NEIGHBORS_H_
#define __NEIGHBORS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Generate_Neighbor_Lists( reax_system*, control_params*, simulation_data*,
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index d099c0be0d3bff33667042caaea0f27a97e1ac78..65e6246756e67c3594511d13a6c42c566be368b4 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -42,8 +42,8 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
{
int i, pj, pk;
bond_order_data *bo_ij;
- reax_list *bonds = &(*lists)[BONDS];
- reax_list *dBOs = &(*lists)[DBO];
+ reax_list *bonds = lists[BONDS];
+ reax_list *dBOs = lists[DBO];
dbond_data *dbo_k;
/* bond orders */
@@ -404,7 +404,7 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
int i, j, id_i, id_j;
char fname[MAX_STR];
FILE *fout;
- reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
+ reax_list *near_nbrs = lists[NEAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.near_nbrs", MAX_STR - 11, control->sim_name );
fout = sfopen( fname, "w" );
@@ -438,7 +438,7 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
int i, j, id_i, id_j;
char fname[MAX_STR];
FILE *fout;
- reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
+ reax_list *near_nbrs = lists[NEAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.near_nbrs_lgj", MAX_STR - 15, control->sim_name );
fout = sfopen( fname, "w" );
@@ -473,7 +473,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
int i, j, id_i, id_j;
char fname[MAX_STR];
FILE *fout;
- reax_list *far_nbrs = &(*lists)[FAR_NBRS];
+ reax_list *far_nbrs = lists[FAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.far_nbrs", MAX_STR - 10, control->sim_name );
fout = sfopen( fname, "w" );
@@ -518,7 +518,7 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
int i, j, id_i, id_j;
char fname[MAX_STR];
FILE *fout;
- reax_list *far_nbrs = &(*lists)[FAR_NBRS];
+ reax_list *far_nbrs = lists[FAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.far_nbrs_lgj", MAX_STR - 14, control->sim_name );
fout = sfopen( fname, "w" );
@@ -687,7 +687,7 @@ void Output_Results( reax_system *system, control_params *control,
out_control->append_traj_frame( system, control, data,
workspace, lists, out_control );
- //Write_PDB( system, &(*lists)[BONDS], data, control, workspace, out_control );
+ //Write_PDB( system, lists[BONDS], data, control, workspace, out_control );
//t_elapsed = Get_Timing_Info( t_start );
//fprintf(stdout, "append_frame took %.6f seconds\n", t_elapsed );
}
diff --git a/sPuReMD/src/print_utils.h b/sPuReMD/src/print_utils.h
index de0fa4abf9fa009d4695fcf8c688acc0c00d4986..4d67eef01181701be83e9ab7103d5d9e8ed90ede 100644
--- a/sPuReMD/src/print_utils.h
+++ b/sPuReMD/src/print_utils.h
@@ -22,7 +22,7 @@
#ifndef __PRINT_UTILS_H_
#define __PRINT_UTILS_H_
-#include "mytypes.h"
+#include "reax_types.h"
char *Get_Element( reax_system*, int );
diff --git a/sPuReMD/src/random.h b/sPuReMD/src/random.h
index cf8517ab151c3b8628a0c9715476f06db3cfc428..b9efb54438f26fad645a2c4e21e74cbf1348044e 100644
--- a/sPuReMD/src/random.h
+++ b/sPuReMD/src/random.h
@@ -23,7 +23,7 @@
#define __RANDOM_H_
/* Includes <stdlib.h> for random( ) */
-#include "mytypes.h"
+#include "reax_types.h"
/* This function seeds the system pseudo random number generator with
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/reax_types.h
similarity index 88%
rename from sPuReMD/src/mytypes.h
rename to sPuReMD/src/reax_types.h
index 754e0c93b030ef976f41b350f08dc3e9e34c9de6..6793dfb9ebc8d657014dca1384a6dae049487ffd 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/reax_types.h
@@ -19,8 +19,8 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
-#ifndef __MYTYPES_H_
-#define __MYTYPES_H_
+#ifndef __REAX_TYPES_H_
+#define __REAX_TYPES_H_
#if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_))
#define __CONFIG_H_
@@ -263,6 +263,59 @@ typedef int ivec[3];
typedef real rtensor[3][3];
+/* struct declarations, see definitions below for comments */
+typedef struct global_parameters global_parameters;
+typedef struct single_body_parameters single_body_parameters;
+typedef struct two_body_parameters two_body_parameters;
+typedef struct three_body_parameters three_body_parameters;
+typedef struct three_body_header three_body_header;
+typedef struct hbond_parameters hbond_parameters;
+typedef struct four_body_parameters four_body_parameters;
+typedef struct four_body_header four_body_header;
+typedef struct reax_interaction reax_interaction;
+typedef struct reax_atom reax_atom;
+typedef struct simulation_box simulation_box;
+typedef struct grid grid;
+typedef struct reax_system reax_system;
+typedef struct control_params control_params;
+typedef struct thermostat thermostat;
+typedef struct isotropic_barostat isotropic_barostat;
+typedef struct flexible_barostat flexible_barostat;
+typedef struct reax_timing reax_timing;
+typedef struct simulation_data simulation_data;
+typedef struct three_body_interaction_data three_body_interaction_data;
+typedef struct near_neighbor_data near_neighbor_data;
+typedef struct far_neighbor_data far_neighbor_data;
+typedef struct hbond_data hbond_data;
+typedef struct dDelta_data dDelta_data;
+typedef struct dbond_data dbond_data;
+typedef struct bond_order_data bond_order_data;
+typedef struct bond_data bond_data;
+typedef struct sparse_matrix sparse_matrix;
+typedef struct reallocate_data reallocate_data;
+typedef struct LR_data LR_data;
+typedef struct cubic_spline_coef cubic_spline_coef;
+typedef struct lookup_table lookup_table;
+typedef struct LR_lookup_table LR_lookup_table;
+typedef struct static_storage static_storage;
+typedef struct reax_list reax_list;
+typedef struct output_controls output_controls;
+typedef struct spuremd_handle spuremd_handle;
+
+
+/* function pointer definitions */
+/**/
+typedef void (*interaction_function)( reax_system*, control_params*,
+ simulation_data*, static_storage*, reax_list**, output_controls* );
+/**/
+typedef void (*evolve_function)( reax_system*, control_params*,
+ simulation_data*, static_storage*,
+ reax_list**, output_controls* );
+/**/
+typedef void (*callback_function)( reax_atom*, simulation_data*, reax_list** );
+
+
+/* struct definitions */
/* Force field global parameters mapping
* (contained in section 1 of file):
*
@@ -305,15 +358,15 @@ typedef real rtensor[3][3];
* l[36] = N/A
* l[37] = version number
* l[38] = p_coa3 */
-typedef struct
+struct global_parameters
{
int n_global;
real* l;
int vdw_type;
-} global_parameters;
+};
-typedef struct
+struct single_body_parameters
{
/* Line one in field file */
/* Two character atom name */
@@ -359,11 +412,11 @@ typedef struct
real rcore2;
real ecore2;
real acore2;
-} single_body_parameters;
+};
/* Two Body Parameters */
-typedef struct
+struct two_body_parameters
{
/* Bond Order parameters */
real p_bo1;
@@ -405,11 +458,11 @@ typedef struct
real v13cor;
real ovc;
-} two_body_parameters;
+};
/* 3-body parameters */
-typedef struct
+struct three_body_parameters
{
/* valence angle */
real theta_00;
@@ -423,28 +476,28 @@ typedef struct
/* 3-body conjugation */
real p_coa1;
-} three_body_parameters;
+};
-typedef struct
+struct three_body_header
{
int cnt;
three_body_parameters prm[MAX_3BODY_PARAM];
-} three_body_header;
+};
/* hydrogen-bond parameters */
-typedef struct
+struct hbond_parameters
{
real r0_hb;
real p_hb1;
real p_hb2;
real p_hb3;
-} hbond_parameters;
+};
/* 4-body parameters */
-typedef struct
+struct four_body_parameters
{
real V1;
real V2;
@@ -455,17 +508,17 @@ typedef struct
/* 4-body conjugation */
real p_cot1;
-} four_body_parameters;
+};
-typedef struct
+struct four_body_header
{
int cnt;
four_body_parameters prm[MAX_4BODY_PARAM];
-} four_body_header;
+};
-typedef struct
+struct reax_interaction
{
int num_atom_types;
global_parameters gp;
@@ -474,10 +527,10 @@ typedef struct
three_body_header ***thbp;
hbond_parameters ***hbp;
four_body_header ****fbp;
-} reax_interaction;
+};
-typedef struct
+struct reax_atom
{
/* Type of this atom */
int type;
@@ -491,10 +544,10 @@ typedef struct
rvec f;
/* Charge on the atom */
real q;
-} reax_atom;
+};
-typedef struct
+struct simulation_box
{
real volume;
@@ -508,10 +561,10 @@ typedef struct
rtensor trans;
rtensor trans_inv;
rtensor g;
-} simulation_box;
+};
-typedef struct
+struct grid
{
int max_atoms;
int max_nbrs;
@@ -530,10 +583,10 @@ typedef struct
int *** end;
ivec **** nbrs;
rvec **** nbrs_cp;
-} grid;
+};
-typedef struct
+struct reax_system
{
/* number of atoms */
int N;
@@ -547,11 +600,11 @@ typedef struct
simulation_box box;
/* grid structure used for binning atoms and tracking neighboring bins */
grid g;
-} reax_system;
+};
/* Simulation control parameters not related to the system */
-typedef struct
+struct control_params
{
char sim_name[MAX_STR];
char restart_from[MAX_STR];
@@ -583,14 +636,6 @@ typedef struct
real bo_cut;
real thb_cut;
real hb_cut;
- real Tap7;
- real Tap6;
- real Tap5;
- real Tap4;
- real Tap3;
- real Tap2;
- real Tap1;
- real Tap0;
real T_init;
real T_final;
@@ -664,38 +709,44 @@ typedef struct
/* num. of iterations used to apply preconditioner via
* Jacobi relaxation scheme (truncated Neumann series) */
unsigned int cm_solver_pre_app_jacobi_iters;
-
+ /**/
int molec_anal;
+ /**/
int freq_molec_anal;
+ /**/
real bg_cut;
+ /**/
int num_ignored;
+ /**/
int ignore[MAX_ATOM_TYPES];
-
+ /**/
int num_threads;
-} control_params;
+ /* function pointers for bonded interactions */
+ interaction_function intr_funcs[NO_OF_INTERACTIONS];
+};
-typedef struct
+struct thermostat
{
real T;
real xi;
real v_xi;
real v_xi_old;
real G_xi;
-} thermostat;
+};
-typedef struct
+struct isotropic_barostat
{
real P;
real eps;
real v_eps;
real v_eps_old;
real a_eps;
-} isotropic_barostat;
+};
-typedef struct
+struct flexible_barostat
{
rtensor P;
real P_scalar;
@@ -710,10 +761,10 @@ typedef struct
rtensor v_g0_old;
rtensor a_g0;
-} flexible_barostat;
+};
-typedef struct
+struct reax_timing
{
/* start time of event */
real start;
@@ -751,10 +802,10 @@ typedef struct
real cm_solver_tri_solve;
/* num. of retries in main sim. loop */
int num_retries;
-} reax_timing;
+};
-typedef struct
+struct simulation_data
{
int step;
int prev_steps;
@@ -808,10 +859,10 @@ typedef struct
rvec tot_press;
reax_timing timing;
-} simulation_data;
+};
-typedef struct
+struct three_body_interaction_data
{
int thb;
/* pointer to the third body on the central atom's nbrlist */
@@ -821,54 +872,54 @@ typedef struct
rvec dcos_di;
rvec dcos_dj;
rvec dcos_dk;
-} three_body_interaction_data;
+};
-typedef struct
+struct near_neighbor_data
{
int nbr;
ivec rel_box;
// rvec ext_factor;
real d;
rvec dvec;
-} near_neighbor_data;
+};
-typedef struct
+struct far_neighbor_data
{
int nbr;
ivec rel_box;
// rvec ext_factor;
real d;
rvec dvec;
-} far_neighbor_data;
+};
-typedef struct
+struct hbond_data
{
int nbr;
int scl;
far_neighbor_data *ptr;
-} hbond_data;
+};
-typedef struct
+struct dDelta_data
{
int wrt;
rvec dVal;
-} dDelta_data;
+};
-typedef struct
+struct dbond_data
{
int wrt;
rvec dBO;
rvec dBOpi;
rvec dBOpi2;
-} dbond_data;
+};
-typedef struct
+struct bond_order_data
{
real BO;
real BO_s;
@@ -892,10 +943,10 @@ typedef struct
rvec dln_BOp_s;
rvec dln_BOp_pi;
rvec dln_BOp_pi2;
-} bond_order_data;
+};
-typedef struct
+struct bond_data
{
int nbr;
int sym_index;
@@ -905,7 +956,7 @@ typedef struct
real d;
rvec dvec;
bond_order_data bo_data;
-} bond_data;
+};
/* compressed row storage (crs) format
@@ -917,17 +968,17 @@ typedef struct
* start: row pointer (last element contains ACTUAL NNZ)
* j: column index for corresponding matrix entry
* val: matrix entry */
-typedef struct
+struct sparse_matrix
{
unsigned int n;
unsigned int m;
unsigned int *start;
unsigned int *j;
real *val;
-} sparse_matrix;
+};
-typedef struct
+struct reallocate_data
{
int num_far;
int Htop;
@@ -937,29 +988,29 @@ typedef struct
int num_bonds;
int num_3body;
int gcell_atoms;
-} reallocate_data;
+};
-typedef struct
+struct LR_data
{
real H;
real e_vdW;
real CEvd;
real e_ele;
real CEclmb;
-} LR_data;
+};
-typedef struct
+struct cubic_spline_coef
{
real a;
real b;
real c;
real d;
-} cubic_spline_coef;
+};
-typedef struct
+struct lookup_table
{
real xmin;
real xmax;
@@ -972,10 +1023,10 @@ typedef struct
real c;
real *y;
-} lookup_table;
+};
-typedef struct
+struct LR_lookup_table
{
real xmin;
real xmax;
@@ -992,10 +1043,10 @@ typedef struct
cubic_spline_coef *CEvd;
cubic_spline_coef *ele;
cubic_spline_coef *CEclmb;
-} LR_lookup_table;
+};
-typedef struct
+struct static_storage
{
/* bond order related storage */
real *total_bond_order;
@@ -1112,6 +1163,9 @@ typedef struct
real *CdDelta; // coefficient of dDelta for force calculations
+ /* Taper */
+ real Tap[8];
+
int *mark;
int *old_mark; // storage for analysis
rvec *x_old;
@@ -1150,11 +1204,11 @@ typedef struct
/* Calculated on the fly in bond_orders.c */
rvec *dDelta;
#endif
-} static_storage;
+};
/* interaction lists */
-typedef struct
+struct reax_list
{
/* num. entries in list */
int n;
@@ -1186,10 +1240,10 @@ typedef struct
/* hydrogen bond interaction list */
hbond_data *hbond_list;
} select;
-} reax_list;
+};
-typedef struct
+struct output_controls
{
FILE *trj;
FILE *out;
@@ -1249,22 +1303,11 @@ typedef struct
FILE *ftot;
FILE *ftot2;
#endif
-} output_controls;
-
-
-/* Function pointer definitions */
-typedef void (*interaction_function)(reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls*);
-
-typedef void (*evolve_function)(reax_system*, control_params*,
- simulation_data*, static_storage*,
- reax_list**, output_controls*, interaction_function*);
-
-typedef void (*callback_function)(reax_atom*, simulation_data*, reax_list*);
+};
/* Handle for working with an instance of the sPuReMD library */
-typedef struct
+struct spuremd_handle
{
/* System info. struct pointer */
reax_system *system;
@@ -1275,16 +1318,14 @@ typedef struct
/* Internal workspace struct pointer */
static_storage *workspace;
/* Reax interaction list struct pointer */
- reax_list *lists;
+ reax_list **lists;
/* Output controls struct pointer */
output_controls *out_control;
/* TRUE if file I/O for simulation output enabled, FALSE otherwise */
int output_enabled;
- /* */
- interaction_function interaction_functions[NO_OF_INTERACTIONS];
/* Callback for getting simulation state at the end of each time step */
callback_function callback;
-} spuremd_handle;
+};
#endif
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index 89f35972c4c8a56aec619ae2d38546a4b6d4500b..a97ae246ef63e68ec426196192127d2436a783e6 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -121,8 +121,8 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
reax_list *bonds;
reax_list *hbonds;
- bonds = &(*lists)[BONDS];
- hbonds = &(*lists)[HBONDS];
+ bonds = lists[BONDS];
+ hbonds = lists[HBONDS];
for ( i = 0; i < system->N; ++i )
{
diff --git a/sPuReMD/src/reset_utils.h b/sPuReMD/src/reset_utils.h
index 97598adbcfd27633b1c01d0129c0210b64bfef47..fd2d6439ffd055ce940bd5e161341305b77b8b9c 100644
--- a/sPuReMD/src/reset_utils.h
+++ b/sPuReMD/src/reset_utils.h
@@ -22,7 +22,7 @@
#ifndef __RESET_UTILS_H_
#define __RESET_UTILS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Reset_Atoms( reax_system* );
diff --git a/sPuReMD/src/restart.h b/sPuReMD/src/restart.h
index 363bb58126162af424567e5baea66142cb391287..bf54a680b6818e9c23e58cdd4a805783eebb6319 100644
--- a/sPuReMD/src/restart.h
+++ b/sPuReMD/src/restart.h
@@ -22,7 +22,7 @@
#ifndef __RESTART_H_
#define __RESTART_H_
-#include "mytypes.h"
+#include "reax_types.h"
typedef struct
{
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 42b97041c0a5f8274b479dc0fb78e51d8cb66cb6..667c7ae3ce78c07e389c3396ba17332c34281f57 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -48,7 +48,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
bond_order_data *bo_ij;
reax_list *bonds;
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
/* Initialize parameters */
p_lp1 = system->reaxprm.gp.l[15];
p_lp3 = system->reaxprm.gp.l[5];
diff --git a/sPuReMD/src/single_body_interactions.h b/sPuReMD/src/single_body_interactions.h
index f575453bde1e4dc3f5f9a502b2a710f30bf96de0..587ed91200b672126518739ab0514e0c9ee7b237 100644
--- a/sPuReMD/src/single_body_interactions.h
+++ b/sPuReMD/src/single_body_interactions.h
@@ -22,7 +22,7 @@
#ifndef __SINGLE_BODY_INTERACTIONS_H_
#define __SINGLE_BODY_INTERACTIONS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void LonePair_OverUnder_Coordination_Energy( reax_system*, control_params*,
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index ffc631cece01911cd5408ca979412ecc031fca0e..26724de8e378d008d3c803824b0470d5abb48be2 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -137,6 +137,7 @@ static void Read_System( const char * const geo_file,
void* setup( const char * const geo_file, const char * const ffield_file,
const char * const control_file )
{
+ int i;
spuremd_handle *spmd_handle;
/* top-level allocation */
@@ -144,17 +145,23 @@ void* setup( const char * const geo_file, const char * const ffield_file,
"setup::spmd_handle" );
/* second-level allocations */
- spmd_handle->system = (reax_system*) smalloc( sizeof(reax_system),
+ spmd_handle->system = smalloc( sizeof(reax_system),
"Setup::spmd_handle->system" );
- spmd_handle->control = (control_params*) smalloc( sizeof(control_params),
+ spmd_handle->control = smalloc( sizeof(control_params),
"Setup::spmd_handle->control" );
- spmd_handle->data = (simulation_data*) smalloc( sizeof(simulation_data),
+ spmd_handle->data = smalloc( sizeof(simulation_data),
"Setup::spmd_handle->data" );
- spmd_handle->workspace = (static_storage*) smalloc( sizeof(static_storage),
+ spmd_handle->workspace = smalloc( sizeof(static_storage),
"Setup::spmd_handle->workspace" );
- spmd_handle->lists = (reax_list*) smalloc( sizeof(reax_list) * LIST_N,
+ spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
"Setup::spmd_handle->lists" );
- spmd_handle->out_control = (output_controls*) smalloc( sizeof(output_controls),
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ spmd_handle->lists[i] = smalloc( sizeof(reax_list),
+ "Setup::spmd_handle->lists[i]" );
+// spmd_handle->lists[i]->allocated = FALSE;
+ }
+ spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
spmd_handle->output_enabled = TRUE;
@@ -202,22 +209,21 @@ int simulate( const void * const handle )
spmd_handle = (spuremd_handle*) handle;
Initialize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists,
- spmd_handle->out_control, &Evolve, spmd_handle->interaction_functions,
+ spmd_handle->workspace, spmd_handle->lists,
+ spmd_handle->out_control, &Evolve,
spmd_handle->output_enabled );
/* compute f_0 */
//if( control.restart == 0 ) {
Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists );
+ spmd_handle->workspace, spmd_handle->lists );
Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
//fprintf( stderr, "total: %.2f secs\n", data.timing.nbrs);
Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control,
- spmd_handle->interaction_functions );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
@@ -226,7 +232,7 @@ int simulate( const void * const handle )
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 );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
Check_Energy( spmd_handle->data );
@@ -248,16 +254,15 @@ int simulate( const void * const handle )
}
Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control,
- spmd_handle->interaction_functions );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
Post_Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
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 );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
Check_Energy( spmd_handle->data );
@@ -265,7 +270,7 @@ int simulate( const void * const handle )
if ( spmd_handle->output_enabled == TRUE )
{
Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
@@ -287,7 +292,7 @@ int simulate( const void * const handle )
if ( spmd_handle->out_control->write_steps > 0 && spmd_handle->output_enabled == TRUE )
{
- Write_PDB( spmd_handle->system, &(spmd_handle->lists[BONDS]), spmd_handle->data,
+ Write_PDB( spmd_handle->system, spmd_handle->lists[BONDS], spmd_handle->data,
spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
}
@@ -308,7 +313,7 @@ int simulate( const void * const handle )
int cleanup( const void * const handle )
{
- int ret;
+ int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
@@ -318,10 +323,14 @@ int cleanup( const void * const handle )
spmd_handle = (spuremd_handle*) handle;
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->output_enabled );
sfree( spmd_handle->out_control, "cleanup::spmd_handle->out_control" );
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ sfree( spmd_handle->lists[i], "cleanup::spmd_handle->lists[i]" );
+ }
sfree( spmd_handle->lists, "cleanup::spmd_handle->lists" );
sfree( spmd_handle->workspace, "cleanup::spmd_handle->workspace" );
sfree( spmd_handle->data, "cleanup::spmd_handle->data" );
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 679aeaf3bdd3a9a14a44099994100d059fbee19e..43d9523827c0225c1aaeb0c2206dc2ae46ff7650 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -22,7 +22,7 @@
#ifndef __SPUREMD_H_
#define __SPUREMD_H_
-#include "mytypes.h"
+#include "reax_types.h"
#define SPUREMD_SUCCESS (0)
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 263c72cecd8c886add75ff18803e3e223219f8ec..50a81b1b04171fae64cf3ac87061179312270dba 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 "reax_types.h"
void Temperature_Control( control_params*, simulation_data*, output_controls* );
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 8e8fc825f09202e6a179275bf358f34efe618f96..41973c5703c0ab6e97ab6c14245364ad2d342957 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -94,9 +94,9 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
real e_ang_total, e_pen_total, e_coa_total;
total_bo = workspace->total_bond_order;
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
bond_list = bonds->select.bond_list;
- thb_intrs = &(*lists)[THREE_BODIES];
+ thb_intrs = lists[THREE_BODIES];
thb_list = thb_intrs->select.three_body_list;
/* global parameters used in these calculations */
p_pen2 = system->reaxprm.gp.l[19];
@@ -707,9 +707,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
#ifdef TEST_FORCES
num_hb_intrs = 0;
#endif
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
bond_list = bonds->select.bond_list;
- hbonds = &(*lists)[HBONDS];
+ hbonds = lists[HBONDS];
hbond_list = hbonds->select.hbond_list;
/* loops below discover the Hydrogen bonds between i-j-k triplets.
diff --git a/sPuReMD/src/three_body_interactions.h b/sPuReMD/src/three_body_interactions.h
index 30fe5c1c5f4b7f0121ef7165a6fd6a4507a2d76d..1aa39b517346d9b6fe5fc28017dd2ff86cacd795 100644
--- a/sPuReMD/src/three_body_interactions.h
+++ b/sPuReMD/src/three_body_interactions.h
@@ -22,7 +22,7 @@
#ifndef __THREE_BODY_INTERACTIONS_H_
#define __THREE_BODY_INTERACTIONS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Calculate_Theta( rvec, real, rvec, real, real*, real* );
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index a5d181e853157387bf74ad7b1564e60d2a6c0d0f..e0931a186d64a2857b80480a86eff208d9e46642 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -22,7 +22,7 @@
#ifndef __TOOL_BOX_H_
#define __TOOL_BOX_H_
-#include "mytypes.h"
+#include "reax_types.h"
/* from box.h */
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index 35e893584f524d5bab66ffc88a120aa61c87fd3b..2acb1bc5f718de880c32e7d4ac809c2c7136e523 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -180,8 +180,8 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
int frame_globals_len, num_bonds, num_thb_intrs;
real P;
char buffer[SIZE];
- reax_list *bonds = &(*lists)[BONDS];
- reax_list *thb_intrs = &(*lists)[THREE_BODIES];
+ reax_list *bonds = lists[BONDS];
+ reax_list *thb_intrs = lists[THREE_BODIES];
bond_data *bo_ij;
/* IMPORTANT: This whole part will go to init_trj after finalized! */
diff --git a/sPuReMD/src/traj.h b/sPuReMD/src/traj.h
index f4d074c8fdd8242737bb38693adad6c03c7638f5..da63f4549181a064a0dc8da373cf0f2fa9f1cf83 100644
--- a/sPuReMD/src/traj.h
+++ b/sPuReMD/src/traj.h
@@ -22,7 +22,7 @@
#ifndef __TRAJ_H__
#define __TRAJ_H__
-#include "mytypes.h"
+#include "reax_types.h"
#define BLOCK_MARK "REAX_BLOCK_MARK "
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index 8fb702c3eef0e19de3402e921883c6232f92c0db..eca3d5a050b6007a77b389b9924e2103a6b0af0b 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -35,7 +35,7 @@ void Bond_Energy( reax_system *system, control_params *control,
real gp3, gp4, gp7, gp10, gp37, ebond_total;
reax_list *bonds;
- bonds = &(*lists)[BONDS];
+ bonds = lists[BONDS];
gp3 = system->reaxprm.gp.l[3];
gp4 = system->reaxprm.gp.l[4];
gp7 = system->reaxprm.gp.l[7];
@@ -177,7 +177,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
p_vdW1 = system->reaxprm.gp.l[28];
p_vdW1i = 1.0 / p_vdW1;
- far_nbrs = &(*lists)[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
e_vdW_total = 0.0;
e_ele_total = 0.0;
@@ -229,20 +229,20 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
/* Calculate Taper and its derivative */
// Tap = nbr_pj->Tap; -- precomputed during compte_H
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
- dTap = 7 * control->Tap7 * r_ij + 6 * control->Tap6;
- dTap = dTap * r_ij + 5 * control->Tap5;
- dTap = dTap * r_ij + 4 * control->Tap4;
- dTap = dTap * r_ij + 3 * control->Tap3;
- dTap = dTap * r_ij + 2 * control->Tap2;
- dTap += control->Tap1 / r_ij;
+ Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
+ Tap = Tap * r_ij + workspace->Tap[5];
+ Tap = Tap * r_ij + workspace->Tap[4];
+ Tap = Tap * r_ij + workspace->Tap[3];
+ Tap = Tap * r_ij + workspace->Tap[2];
+ Tap = Tap * r_ij + workspace->Tap[1];
+ Tap = Tap * r_ij + workspace->Tap[0];
+
+ dTap = 7 * workspace->Tap[7] * r_ij + 6 * workspace->Tap[6];
+ dTap = dTap * r_ij + 5 * workspace->Tap[5];
+ dTap = dTap * r_ij + 4 * workspace->Tap[4];
+ dTap = dTap * r_ij + 3 * workspace->Tap[3];
+ dTap = dTap * r_ij + 2 * workspace->Tap[2];
+ dTap += workspace->Tap[1] / r_ij;
/* vdWaals Calculations */
if ( system->reaxprm.gp.vdw_type == 1 || system->reaxprm.gp.vdw_type == 3 )
@@ -404,7 +404,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
void LR_vdW_Coulomb( reax_system *system, control_params *control,
- int i, int j, real r_ij, LR_data *lr )
+ static_storage *workspace, int i, int j, real r_ij, LR_data *lr )
{
real p_vdW1 = system->reaxprm.gp.l[28];
real p_vdW1i = 1.0 / p_vdW1;
@@ -415,25 +415,25 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
real e_core, de_core;
two_body_parameters *twbp;
- twbp = &(system->reaxprm.tbp[i][j]);
+ twbp = &system->reaxprm.tbp[i][j];
e_core = 0;
de_core = 0;
/* calculate taper and its derivative */
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
- dTap = 7 * control->Tap7 * r_ij + 6 * control->Tap6;
- dTap = dTap * r_ij + 5 * control->Tap5;
- dTap = dTap * r_ij + 4 * control->Tap4;
- dTap = dTap * r_ij + 3 * control->Tap3;
- dTap = dTap * r_ij + 2 * control->Tap2;
- dTap += control->Tap1 / r_ij;
+ Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
+ Tap = Tap * r_ij + workspace->Tap[5];
+ Tap = Tap * r_ij + workspace->Tap[4];
+ Tap = Tap * r_ij + workspace->Tap[3];
+ Tap = Tap * r_ij + workspace->Tap[2];
+ Tap = Tap * r_ij + workspace->Tap[1];
+ Tap = Tap * r_ij + workspace->Tap[0];
+
+ dTap = 7 * workspace->Tap[7] * r_ij + 6 * workspace->Tap[6];
+ dTap = dTap * r_ij + 5 * workspace->Tap[5];
+ dTap = dTap * r_ij + 4 * workspace->Tap[4];
+ dTap = dTap * r_ij + 3 * workspace->Tap[3];
+ dTap = dTap * r_ij + 2 * workspace->Tap[2];
+ dTap += workspace->Tap[1] / r_ij;
/* vdWaals calculations */
@@ -528,7 +528,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
reax_list *far_nbrs;
real e_vdW_total, e_ele_total;
- far_nbrs = &(*lists)[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
steps = data->step - data->prev_steps;
update_freq = out_control->energy_update_freq;
update_energies = update_freq > 0 && steps % update_freq == 0;
diff --git a/sPuReMD/src/two_body_interactions.h b/sPuReMD/src/two_body_interactions.h
index 042ee246ce0171f5f2be170352804938acdcbc1a..eaf73aaad93dbc423952a6050161aac7057839f9 100644
--- a/sPuReMD/src/two_body_interactions.h
+++ b/sPuReMD/src/two_body_interactions.h
@@ -22,7 +22,7 @@
#ifndef __TWO_BODY_INTERACTIONS_H_
#define __TWO_BODY_INTERACTIONS_H_
-#include "mytypes.h"
+#include "reax_types.h"
void Bond_Energy( reax_system*, control_params*, simulation_data*,
@@ -31,7 +31,8 @@ void Bond_Energy( reax_system*, control_params*, simulation_data*,
void vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
-void LR_vdW_Coulomb( reax_system*, control_params*, int, int, real, LR_data* );
+void LR_vdW_Coulomb( reax_system*, control_params*, static_storage *,
+ int, int, real, LR_data* );
void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
diff --git a/sPuReMD/src/vector.h b/sPuReMD/src/vector.h
index 5550a2d7354362ffcfe24bbdbd94bc77af0afc13..1c804295c2d8f3e4180ee44824abe2b36ee42925 100644
--- a/sPuReMD/src/vector.h
+++ b/sPuReMD/src/vector.h
@@ -22,7 +22,7 @@
#ifndef __VECTOR_H_
#define __VECTOR_H_
-#include "mytypes.h"
+#include "reax_types.h"
#include "random.h"
diff --git a/tools/driver.py b/tools/driver.py
index 93163deaec3d723284d3d61d73d57d3d23df3fb3..1367be5ad25e65f51ac607c032837a612335e231 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -393,7 +393,7 @@ if __name__ == '__main__':
get_atoms.restype = POINTER(ReaxAtom)
CALLBACKFUNC = CFUNCTYPE(None, POINTER(ReaxAtom),
- POINTER(SimulationData), POINTER(ReaxList))
+ POINTER(SimulationData), POINTER(POINTER(ReaxList)))
setup_callback = lib.setup_callback
setup_callback.restype = c_int