From 7851f9aae5335d510de2ddf3e84788f7cf36b070 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Mon, 15 Oct 2018 22:37:02 -0400
Subject: [PATCH] sPuReMD: refactoring to align code with PuReMD. Disable
SMALL_BOX_SUPPORT. Refactor to remove union-based polymorphic behavior for
reax_list. Other general code clean-ups.
---
sPuReMD/src/allocate.c | 22 +--
sPuReMD/src/analyze.c | 124 ++++++-------
sPuReMD/src/bond_orders.c | 100 +++++-----
sPuReMD/src/box.c | 21 ++-
sPuReMD/src/control.c | 36 ++--
sPuReMD/src/forces.c | 248 ++++++++++++++-----------
sPuReMD/src/four_body_interactions.c | 20 +-
sPuReMD/src/geo_tools.c | 64 ++++---
sPuReMD/src/init_md.c | 16 +-
sPuReMD/src/integrate.c | 11 +-
sPuReMD/src/list.c | 52 +++---
sPuReMD/src/lookup.c | 6 +-
sPuReMD/src/neighbors.c | 25 ++-
sPuReMD/src/print_utils.c | 66 +++----
sPuReMD/src/reax_types.h | 90 ++++++---
sPuReMD/src/reset_utils.c | 2 +-
sPuReMD/src/single_body_interactions.c | 20 +-
sPuReMD/src/spuremd.c | 18 +-
sPuReMD/src/system_props.c | 15 +-
sPuReMD/src/three_body_interactions.c | 78 ++++----
sPuReMD/src/tool_box.c | 75 +++-----
sPuReMD/src/tool_box.h | 4 +-
sPuReMD/src/traj.c | 48 ++---
sPuReMD/src/two_body_interactions.c | 22 +--
24 files changed, 613 insertions(+), 570 deletions(-)
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 4c7bca5c..50c5e410 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -31,24 +31,24 @@ void PreAllocate_Space( reax_system *system, control_params *control,
{
int i;
- system->atoms = (reax_atom*) scalloc( system->N,
- sizeof(reax_atom), "PreAllocate_Space::system->atoms" );
- workspace->orig_id = (int*) scalloc( system->N,
- sizeof(int), "PreAllocate_Space::workspace->orid_id" );
+ system->atoms = scalloc( system->N, sizeof(reax_atom),
+ "PreAllocate_Space::system->atoms" );
+ workspace->orig_id = scalloc( system->N, 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), "PreAllocate_Space::workspace->restricted_atoms" );
+ workspace->restricted = scalloc( system->N, sizeof(int),
+ "PreAllocate_Space::workspace->restricted_atoms" );
- workspace->restricted_list = (int**) scalloc( system->N,
- sizeof(int*), "PreAllocate_Space::workspace->restricted_list" );
+ workspace->restricted_list = scalloc( system->N, sizeof(int*),
+ "PreAllocate_Space::workspace->restricted_list" );
for ( i = 0; i < system->N; ++i )
{
- workspace->restricted_list[i] = (int*) scalloc( MAX_RESTRICT,
- sizeof(int), "PreAllocate_Space::workspace->restricted_list[i]" );
+ workspace->restricted_list[i] = scalloc( MAX_RESTRICT, sizeof(int),
+ "PreAllocate_Space::workspace->restricted_list[i]" );
}
}
}
@@ -258,7 +258,7 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
workspace->U = NULL;
}
- if ( control->hb_cut > 0.0 && realloc->hbonds > 0 )
+ if ( control->hbond_cut > 0.0 && realloc->hbonds > 0 )
{
Reallocate_HBonds_List( system->N, workspace->num_H, workspace->hbond_index,
lists[HBONDS] );
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index 68e94284..c0a3f816 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -68,13 +68,13 @@ static void Copy_Bond_List( reax_system *system, control_params *control,
// fprintf( stdout, "%d: ", i );
for ( j = Start_Index( i, new_bonds ); j < End_Index( i, new_bonds ); ++j )
- if ( new_bonds->select.bond_list[j].bo_data.BO >= control->bg_cut )
+ if ( new_bonds->bond_list[j].bo_data.BO >= control->bg_cut )
{
- // fprintf( stderr, "%d ", new_bonds->select.bond_list[j].nbr );
- old_bonds->select.bond_list[ top_old ].nbr =
- new_bonds->select.bond_list[j].nbr;
- old_bonds->select.bond_list[ top_old ].bo_data.BO =
- new_bonds->select.bond_list[j].bo_data.BO;
+ // fprintf( stderr, "%d ", new_bonds->bond_list[j].nbr );
+ old_bonds->bond_list[ top_old ].nbr =
+ new_bonds->bond_list[j].nbr;
+ old_bonds->bond_list[ top_old ].bo_data.BO =
+ new_bonds->bond_list[j].bo_data.BO;
top_old++;
}
@@ -95,14 +95,14 @@ static int Compare_Bond_Lists( int atom, control_params *control, reax_list **li
/*fprintf( stdout, "\n%d\nold_bonds:", atom );
for( oldp = Start_Index( atom, old_bonds );
oldp < End_Index( atom, old_bonds ); ++oldp )
- if( old_bonds->select.bond_list[oldp].bo_data.BO >= control->bg_cut )
- fprintf( stdout, "%5d", old_bonds->select.bond_list[oldp].nbr );
+ if( old_bonds->bond_list[oldp].bo_data.BO >= control->bg_cut )
+ fprintf( stdout, "%5d", old_bonds->bond_list[oldp].nbr );
fprintf( stdout, "\nnew_bonds:" );
for( newp = Start_Index( atom, new_bonds );
newp < End_Index( atom, new_bonds ); ++newp )
- if( new_bonds->select.bond_list[newp].bo_data.BO >= control->bg_cut )
- fprintf( stdout, "%5d", new_bonds->select.bond_list[newp].nbr );*/
+ if( new_bonds->bond_list[newp].bo_data.BO >= control->bg_cut )
+ fprintf( stdout, "%5d", new_bonds->bond_list[newp].nbr );*/
for ( oldp = Start_Index( atom, old_bonds ),
@@ -112,28 +112,28 @@ static int Compare_Bond_Lists( int atom, control_params *control, reax_list **li
newp = MIN( newp + 1, End_Index( atom, new_bonds ) ) )
{
while ( oldp < End_Index( atom, old_bonds )
- && old_bonds->select.bond_list[oldp].bo_data.BO < control->bg_cut )
+ && old_bonds->bond_list[oldp].bo_data.BO < control->bg_cut )
{
++oldp;
}
while ( newp < End_Index( atom, new_bonds )
- && new_bonds->select.bond_list[newp].bo_data.BO < control->bg_cut )
+ && new_bonds->bond_list[newp].bo_data.BO < control->bg_cut )
{
++newp;
}
/*fprintf( fout, "%d, oldp: %d - %d, newp: %d - %d",
- atom, oldp, old_bonds->select.bond_list[oldp].nbr,
- newp, new_bonds->select.bond_list[newp].nbr );*/
+ atom, oldp, old_bonds->bond_list[oldp].nbr,
+ newp, new_bonds->bond_list[newp].nbr );*/
if ( oldp < End_Index( atom, old_bonds ) )
{
/* there are some other bonds in the old list */
if ( newp < End_Index( atom, new_bonds ) )
{
- if ( old_bonds->select.bond_list[oldp].nbr !=
- new_bonds->select.bond_list[newp].nbr )
+ if ( old_bonds->bond_list[oldp].nbr !=
+ new_bonds->bond_list[newp].nbr )
{
//fprintf( fout, " --> case1, return 1\n" );
return 1;
@@ -187,10 +187,10 @@ static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
for ( i = start; i < end; ++i )
{
- if ( bonds->select.bond_list[i].bo_data.BO >= control->bg_cut &&
- !mark[bonds->select.bond_list[i].nbr] )
+ if ( bonds->bond_list[i].bo_data.BO >= control->bg_cut &&
+ !mark[bonds->bond_list[i].nbr] )
{
- Get_Molecule( bonds->select.bond_list[i].nbr, m, mark,
+ Get_Molecule( bonds->bond_list[i].nbr, m, mark,
system, control, bonds, print, fout );
}
}
@@ -399,11 +399,11 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
fprintf( fout, "%5d(%s) was connected to:", rev1, system->atoms[a1].name );
for ( i = Start_Index(a1, old_bonds); i < End_Index(a1, old_bonds); ++i )
{
- if ( old_bonds->select.bond_list[i].bo_data.BO >= control->bg_cut )
+ if ( old_bonds->bond_list[i].bo_data.BO >= control->bg_cut )
{
fprintf( fout, " %5d(%s)",
- workspace->orig_id[ old_bonds->select.bond_list[i].nbr ],
- system->atoms[ old_bonds->select.bond_list[i].nbr ].name );
+ workspace->orig_id[ old_bonds->bond_list[i].nbr ],
+ system->atoms[ old_bonds->bond_list[i].nbr ].name );
}
}
fprintf( fout, "\n" );
@@ -411,11 +411,11 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
fprintf( fout, "%5d(%s) was connected to:", rev2, system->atoms[a2].name );
for ( i = Start_Index(a2, old_bonds); i < End_Index(a2, old_bonds); ++i )
{
- if ( old_bonds->select.bond_list[i].bo_data.BO >= control->bg_cut )
+ if ( old_bonds->bond_list[i].bo_data.BO >= control->bg_cut )
{
fprintf( fout, " %5d(%s)",
- workspace->orig_id[ old_bonds->select.bond_list[i].nbr ],
- system->atoms[ old_bonds->select.bond_list[i].nbr ].name );
+ workspace->orig_id[ old_bonds->bond_list[i].nbr ],
+ system->atoms[ old_bonds->bond_list[i].nbr ].name );
}
}
fprintf( fout, "\n" );
@@ -433,69 +433,69 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
/* fprintf( fout, "\n%d\nold_bonds:", atom );
for( oldp = Start_Index( atom, old_bonds );
oldp < End_Index( atom, old_bonds ); ++oldp )
- if( old_bonds->select.bond_list[oldp].bo_data.BO >= control->bg_cut )
- fprintf( fout, "%5d", old_bonds->select.bond_list[oldp].nbr );
+ if( old_bonds->bond_list[oldp].bo_data.BO >= control->bg_cut )
+ fprintf( fout, "%5d", old_bonds->bond_list[oldp].nbr );
fprintf( fout, "\nnew_bonds:" );
for( newp = Start_Index( atom, new_bonds );
newp < End_Index( atom, new_bonds ); ++newp )
- if( new_bonds->select.bond_list[newp].bo_data.BO >= control->bg_cut )
- fprintf( fout, "%6d", new_bonds->select.bond_list[newp].nbr );
+ if( new_bonds->bond_list[newp].bo_data.BO >= control->bg_cut )
+ fprintf( fout, "%6d", new_bonds->bond_list[newp].nbr );
fprintf( fout, "\n" ); */
for ( oldp = Start_Index( atom, old_bonds );
oldp < End_Index( atom, old_bonds )
- && old_bonds->select.bond_list[oldp].nbr < atom; ++oldp )
+ && old_bonds->bond_list[oldp].nbr < atom; ++oldp )
;
for ( newp = Start_Index( atom, new_bonds );
newp < End_Index( atom, new_bonds )
- && new_bonds->select.bond_list[newp].nbr < atom; ++newp )
+ && new_bonds->bond_list[newp].nbr < atom; ++newp )
;
while ( oldp < End_Index( atom, old_bonds ) ||
newp < End_Index( atom, new_bonds ) )
{
while ( oldp < End_Index( atom, old_bonds ) &&
- old_bonds->select.bond_list[oldp].bo_data.BO < control->bg_cut )
+ old_bonds->bond_list[oldp].bo_data.BO < control->bg_cut )
{
++oldp;
}
while ( newp < End_Index( atom, new_bonds ) &&
- new_bonds->select.bond_list[newp].bo_data.BO < control->bg_cut )
+ new_bonds->bond_list[newp].bo_data.BO < control->bg_cut )
{
++newp;
}
/*fprintf( fout, "%d, oldp: %d - %d: %f newp: %d - %d: %f",
- atom, oldp, old_bonds->select.bond_list[oldp].nbr,
- old_bonds->select.bond_list[oldp].bo_data.BO,
- newp, new_bonds->select.bond_list[newp].nbr,
- new_bonds->select.bond_list[newp].bo_data.BO ); */
+ atom, oldp, old_bonds->bond_list[oldp].nbr,
+ old_bonds->bond_list[oldp].bo_data.BO,
+ newp, new_bonds->bond_list[newp].nbr,
+ new_bonds->bond_list[newp].bo_data.BO ); */
if ( oldp < End_Index( atom, old_bonds ) )
{
/* there are some more bonds in the old list */
if ( newp < End_Index( atom, new_bonds ) )
{
- if ( old_bonds->select.bond_list[oldp].nbr <
- new_bonds->select.bond_list[newp].nbr )
+ if ( old_bonds->bond_list[oldp].nbr <
+ new_bonds->bond_list[newp].nbr )
{
// fprintf( fout, "%5d-%5d bond broken\n",
- // atom, old_bonds->select.bond_list[oldp].nbr );
+ // atom, old_bonds->bond_list[oldp].nbr );
Report_Bond_Change( system, control, workspace, old_bonds, new_bonds,
- atom, old_bonds->select.bond_list[oldp].nbr, 0,
+ atom, old_bonds->bond_list[oldp].nbr, 0,
fout );
++oldp;
}
- else if ( old_bonds->select.bond_list[oldp].nbr >
- new_bonds->select.bond_list[newp].nbr )
+ else if ( old_bonds->bond_list[oldp].nbr >
+ new_bonds->bond_list[newp].nbr )
{
// fprintf( fout, "%5d-%5d bond formed\n",
- // atom, new_bonds->select.bond_list[newp].nbr );
+ // atom, new_bonds->bond_list[newp].nbr );
Report_Bond_Change( system, control, workspace, old_bonds, new_bonds,
- atom, new_bonds->select.bond_list[newp].nbr, 1,
+ atom, new_bonds->bond_list[newp].nbr, 1,
fout );
++newp;
}
@@ -509,13 +509,13 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
/* there is no other bond in the new list */
while ( oldp < End_Index( atom, old_bonds ) )
{
- if ( old_bonds->select.bond_list[oldp].bo_data.BO >= control->bg_cut )
+ if ( old_bonds->bond_list[oldp].bo_data.BO >= control->bg_cut )
{
// fprintf( fout, "%5d-%5d bond broken\n",
- // atom, old_bonds->select.bond_list[oldp].nbr );
+ // atom, old_bonds->bond_list[oldp].nbr );
Report_Bond_Change( system, control, workspace,
old_bonds, new_bonds, atom,
- old_bonds->select.bond_list[oldp].nbr, 0,
+ old_bonds->bond_list[oldp].nbr, 0,
fout );
}
++oldp;
@@ -529,13 +529,13 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
/* there is at least one other bond in the new list */
while ( newp < End_Index( atom, new_bonds ) )
{
- if ( new_bonds->select.bond_list[newp].bo_data.BO >= control->bg_cut )
+ if ( new_bonds->bond_list[newp].bo_data.BO >= control->bg_cut )
{
// fprintf( fout, "%5d-%5d bond formed\n",
- // atom, new_bonds->select.bond_list[newp].nbr );
+ // atom, new_bonds->bond_list[newp].nbr );
Report_Bond_Change( system, control, workspace,
old_bonds, new_bonds, atom,
- new_bonds->select.bond_list[newp].nbr, 1, fout );
+ new_bonds->bond_list[newp].nbr, 1, fout );
}
++newp;
}
@@ -568,8 +568,8 @@ static void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
end = End_Index( atom, bonds );
for ( i = start; i < end; ++i )
{
- nbr = bonds->select.bond_list[i].nbr;
- bo = bonds->select.bond_list[i].bo_data.BO;
+ nbr = bonds->bond_list[i].nbr;
+ bo = bonds->bond_list[i].bo_data.BO;
if ( bo >= control->bg_cut && !mark[nbr] )
Visit_Bonds( nbr, mark, type, system, control, bonds, ignore );
}
@@ -675,7 +675,7 @@ static void Analyze_Silica( reax_system *system, control_params *control,
for ( newp = Start_Index( atom, new_bonds );
newp < End_Index( atom, new_bonds ); ++newp )
{
- if ( new_bonds->select.bond_list[newp].bo_data.BO >= control->bg_cut )
+ if ( new_bonds->bond_list[newp].bo_data.BO >= control->bg_cut )
{
++coord;
}
@@ -726,40 +726,40 @@ static void Analyze_Silica( reax_system *system, control_params *control,
{
for ( pi = Start_Index(j, new_bonds); pi < End_Index(j, new_bonds); ++pi )
{
- if ( new_bonds->select.bond_list[pi].bo_data.BO >= control->bg_cut )
+ if ( new_bonds->bond_list[pi].bo_data.BO >= control->bg_cut )
{
- i = new_bonds->select.bond_list[pi].nbr;
+ i = new_bonds->bond_list[pi].nbr;
if (system->atoms[i].type == O_ATOM || system->atoms[i].type == SI_ATOM)
{
for ( pk = Start_Index( pi, thb_intrs );
pk < End_Index( pi, thb_intrs ); ++pk )
{
- k = thb_intrs->select.three_body_list[pk].thb;
- pk_j = thb_intrs->select.three_body_list[pk].pthb;
+ k = thb_intrs->three_body_list[pk].thb;
+ pk_j = thb_intrs->three_body_list[pk].pthb;
// get k's pointer on j's bond list
- if ( new_bonds->select.bond_list[pk_j].bo_data.BO >=
+ if ( new_bonds->bond_list[pk_j].bo_data.BO >=
control->bg_cut ) // physical j&k bond
{
/*fprintf( fout, "%5d(%d) %5d(%d) %5d(%d) %8.3f\n",
i, system->atoms[i].type, j, system->atoms[j].type,
k, system->atoms[k].type,
- thb_intrs->select.three_body_list[pk].theta );*/
+ thb_intrs->three_body_list[pk].theta );*/
if ( system->atoms[i].type == O_ATOM &&
system->atoms[j].type == SI_ATOM &&
system->atoms[k].type == O_ATOM )
{
O_SI_O_count++;
- O_SI_O += thb_intrs->select.three_body_list[pk].theta;
+ O_SI_O += thb_intrs->three_body_list[pk].theta;
}
else if ( system->atoms[i].type == SI_ATOM &&
system->atoms[j].type == O_ATOM &&
system->atoms[k].type == SI_ATOM )
{
SI_O_SI_count++;
- SI_O_SI += thb_intrs->select.three_body_list[pk].theta;
+ SI_O_SI += thb_intrs->three_body_list[pk].theta;
}
}
}
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index 0b02a6dd..dd85ba4e 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -44,14 +44,14 @@ void Get_dBO( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for ( k = start_pj; k < end_pj; ++k )
{
- rvec_Scale( v[dBOs->select.dbo_list[k].wrt],
- C, dBOs->select.dbo_list[k].dBO );
+ rvec_Scale( v[dBOs->dbo_list[k].wrt],
+ C, dBOs->dbo_list[k].dBO );
}
}
@@ -66,13 +66,13 @@ void Get_dBOpinpi2( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for ( k = start_pj; k < end_pj; ++k )
{
- dbo_k = &(dBOs->select.dbo_list[k]);
+ dbo_k = &dBOs->dbo_list[k];
rvec_Scale( vpi[dbo_k->wrt], Cpi, dbo_k->dBOpi );
rvec_Scale( vpi2[dbo_k->wrt], Cpi2, dbo_k->dBOpi2 );
}
@@ -88,7 +88,7 @@ void Add_dBO( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
@@ -96,8 +96,8 @@ void Add_dBO( reax_system *system, reax_list **lists,
for ( k = start_pj; k < end_pj; ++k )
{
- rvec_ScaledAdd( v[dBOs->select.dbo_list[k].wrt],
- C, dBOs->select.dbo_list[k].dBO );
+ rvec_ScaledAdd( v[dBOs->dbo_list[k].wrt],
+ C, dBOs->dbo_list[k].dBO );
}
}
@@ -112,13 +112,13 @@ void Add_dBOpinpi2( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for ( k = start_pj; k < end_pj; ++k )
{
- dbo_k = &(dBOs->select.dbo_list[k]);
+ dbo_k = &dBOs->dbo_list[k];
rvec_ScaledAdd( vpi[dbo_k->wrt], Cpi, dbo_k->dBOpi );
rvec_ScaledAdd( vpi2[dbo_k->wrt], Cpi2, dbo_k->dBOpi2 );
}
@@ -134,14 +134,14 @@ void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for ( k = start_pj; k < end_pj; ++k )
{
- rvec_ScaledAdd( system->atoms[dBOs->select.dbo_list[k].wrt].f,
- C, dBOs->select.dbo_list[k].dBO );
+ rvec_ScaledAdd( system->atoms[dBOs->dbo_list[k].wrt].f,
+ C, dBOs->dbo_list[k].dBO );
}
}
@@ -156,13 +156,13 @@ void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- pj = bonds->select.bond_list[pj].dbond_index;
+ pj = bonds->bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for ( k = start_pj; k < end_pj; ++k )
{
- dbo_k = &(dBOs->select.dbo_list[k]);
+ dbo_k = &dBOs->dbo_list[k];
rvec_ScaledAdd( system->atoms[dbo_k->wrt].f, Cpi, dbo_k->dBOpi );
rvec_ScaledAdd( system->atoms[dbo_k->wrt].f, Cpi2, dbo_k->dBOpi2 );
}
@@ -180,8 +180,8 @@ void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v
for ( k = start; k < end; ++k )
{
- rvec_ScaledAdd( v[dDeltas->select.dDelta_list[k].wrt],
- C, dDeltas->select.dDelta_list[k].dVal );
+ rvec_ScaledAdd( v[dDeltas->dDelta_list[k].wrt],
+ C, dDeltas->dDelta_list[k].dVal );
}
}
@@ -197,8 +197,8 @@ void Add_dDelta_to_Forces( reax_system *system, reax_list **lists, int i, real C
for ( k = start; k < end; ++k )
{
- rvec_ScaledAdd( system->atoms[dDeltas->select.dDelta_list[k].wrt].f,
- C, dDeltas->select.dDelta_list[k].dVal );
+ rvec_ScaledAdd( system->atoms[dDeltas->dDelta_list[k].wrt].f,
+ C, dDeltas->dDelta_list[k].dVal );
}
}
@@ -215,23 +215,23 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
bonds = lists[BONDS];
dBOs = lists[DBO];
- j = bonds->select.bond_list[pj].nbr;
- bo_ij = &bonds->select.bond_list[pj].bo_data;
+ j = bonds->bond_list[pj].nbr;
+ bo_ij = &bonds->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->dbo_list[ *top ];
for ( k = start_i; k < end_i; ++k )
{
- nbr_k = &bonds->select.bond_list[k];
+ nbr_k = &bonds->bond_list[k];
- for ( ; l < end_j && bonds->select.bond_list[l].nbr < nbr_k->nbr; ++l )
+ for ( ; l < end_j && bonds->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];
+ nbr_l = &bonds->bond_list[l];
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
@@ -271,10 +271,10 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
rvec_Scale( top_dbo->dBOpi, -bo_ij->C3dbopi, dBOp ); //dBOpi-3
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C3dbopi2, dBOp );//dBOpp-3
- if ( l < end_j && bonds->select.bond_list[l].nbr == nbr_k->nbr )
+ if ( l < end_j && bonds->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->bond_list[l];
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C3dbo, dBOp ); //dBO,3rd
@@ -309,7 +309,7 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
{
/* 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->bond_list[l];
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
@@ -360,10 +360,10 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
#endif
bonds = lists[BONDS];
- nbr_j = &bonds->select.bond_list[pj];
+ nbr_j = &bonds->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_ji = &bonds->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];
@@ -396,7 +396,7 @@ 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->bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
f_k = &workspace->f_local[tid * system->N + k];
@@ -457,7 +457,7 @@ 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->bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
f_k = &workspace->f_local[tid * system->N + k];
@@ -549,10 +549,10 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
/* Initializations */
bonds = lists[BONDS];
- nbr_j = &bonds->select.bond_list[pj];
+ nbr_j = &bonds->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_ji = &bonds->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];
@@ -581,7 +581,7 @@ 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->bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
f_k = &workspace->f_local[tid * system->N + k];
@@ -625,7 +625,7 @@ 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->bond_list[pk];
k = nbr_k->nbr;
#ifdef _OPENMP
f_k = &workspace->f_local[tid * system->N + k];
@@ -776,9 +776,9 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
type_j = system->atoms[j].type;
- bo_ij = &bonds->select.bond_list[pj].bo_data;
+ bo_ij = &bonds->bond_list[pj].bo_data;
if ( i < j )
{
@@ -812,7 +812,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->dbo_list[ top_dbo ];
/* compute dBO_ij/dr_i */
pdbo->wrt = i;
@@ -983,11 +983,11 @@ 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->dDelta_list[top_dDelta];
for ( pj = start_i; pj < end_i; ++pj )
{
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
if ( !rvec_isZero( workspace->dDelta[j] ) )
{
@@ -1002,7 +1002,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
end_j = End_Index(j, bonds);
for ( pk = start_j; pk < end_j; ++pk )
{
- k = bonds->select.bond_list[pk].nbr;
+ k = bonds->bond_list[pk].nbr;
if ( !rvec_isZero( workspace->dDelta[k] ) )
{
ptop_dDelta->wrt = k;
@@ -1018,10 +1018,10 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
/*for( pj=Start_Index(i,dDeltas); pj<End_Index(i,dDeltas); ++pj )
fprintf( stdout, "dDel: %d %d [%g %g %g]\n",
- i+1, dDeltas->select.dDelta_list[pj].wrt+1,
- dDeltas->select.dDelta_list[pj].dVal[0],
- dDeltas->select.dDelta_list[pj].dVal[1],
- dDeltas->select.dDelta_list[pj].dVal[2] );*/
+ i+1, dDeltas->dDelta_list[pj].wrt+1,
+ dDeltas->dDelta_list[pj].dVal[0],
+ dDeltas->dDelta_list[pj].dVal[1],
+ dDeltas->dDelta_list[pj].dVal[2] );*/
#endif
}
@@ -1047,7 +1047,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
type_j = system->atoms[j].type;
if ( type_j < 0 )
@@ -1064,10 +1064,10 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
{
/* We only need to update bond orders from bo_ji
everything else is set in uncorrected_bo calculations */
- sym_index = bonds->select.bond_list[pj].sym_index;
+ sym_index = bonds->bond_list[pj].sym_index;
- bo_ij = &(bonds->select.bond_list[ pj ].bo_data);
- bo_ji = &(bonds->select.bond_list[ sym_index ].bo_data);
+ bo_ij = &bonds->bond_list[ pj ].bo_data;
+ bo_ji = &bonds->bond_list[ sym_index ].bo_data;
bo_ij->BO = bo_ji->BO;
bo_ij->BO_s = bo_ji->BO_s;
bo_ij->BO_pi = bo_ji->BO_pi;
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index fef4145c..e142f24c 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -60,15 +60,16 @@ void Make_Consistent( simulation_box* box )
box->box_inv[2][2] = (box->box[0][0] * box->box[1][1] -
box->box[0][1] * box->box[1][0]) * one_vol;
- box->box_norms[0] = SQRT( SQR(box->box[0][0]) +
- SQR(box->box[0][1]) +
- SQR(box->box[0][2]) );
- box->box_norms[1] = SQRT( SQR(box->box[1][0]) +
- SQR(box->box[1][1]) +
- SQR(box->box[1][2]) );
- box->box_norms[2] = SQRT( SQR(box->box[2][0]) +
- SQR(box->box[2][1]) +
- SQR(box->box[2][2]) );
+ box->box_norms[0] = SQRT( SQR(box->box[0][0]) + SQR(box->box[0][1])
+ + SQR(box->box[0][2]) );
+ box->box_norms[1] = SQRT( SQR(box->box[1][0]) + SQR(box->box[1][1])
+ + SQR(box->box[1][2]) );
+ box->box_norms[2] = SQRT( SQR(box->box[2][0]) + SQR(box->box[2][1])
+ + SQR(box->box[2][2]) );
+
+ box->max[0] = box->min[0] + box->box_norms[0];
+ box->max[1] = box->min[1] + box->box_norms[1];
+ box->max[2] = box->min[2] + box->box_norms[2];
box->trans[0][0] = box->box[0][0] / box->box_norms[0];
box->trans[0][1] = box->box[1][0] / box->box_norms[0];
@@ -169,6 +170,8 @@ void Setup_Box( real a, real b, real c, real alpha, real beta, real gamma,
s_gamma = SIN( DEG2RAD(gamma) );
zi = (c_alpha - c_beta * c_gamma) / s_gamma;
+ rvec_MakeZero( box->min );
+
box->box[0][0] = a;
box->box[0][1] = 0.0;
box->box[0][2] = 0.0;
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index baacc07b..1c952338 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -57,13 +57,18 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->periodic_boundaries = 1;
control->reneighbor = 1;
- control->vlist_cut = 0;
- control->nbr_cut = 4.;
- control->r_cut = 10.;
- control->r_sp_cut = 10.;
- control->bo_cut = 0.01;
+
+ /* interaction cutoffs from force field global paramters */
+ control->bo_cut = 0.01 * system->reaxprm.gp.l[29];
+ control->nonb_low = system->reaxprm.gp.l[11];
+ control->nonb_cut = system->reaxprm.gp.l[12];
+
+ /* defaults values for other cutoffs */
+ control->vlist_cut = control->nonb_cut;
+ control->bond_cut = 5.0;
+ control->bg_cut = 0.3;
control->thb_cut = 0.001;
- control->hb_cut = 0.0;
+ control->hbond_cut = 0.0;
control->tabulate = 0;
@@ -125,7 +130,6 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->molec_anal = NO_ANALYSIS;
control->freq_molec_anal = 0;
- control->bg_cut = 0.3;
control->num_ignored = 0;
memset( control->ignore, 0, sizeof(int) * MAX_ATOM_TYPES );
@@ -223,12 +227,12 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
else if ( strncmp(tmp[0], "vlist_buffer", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
- control->vlist_cut = val;
+ control->vlist_cut = val + control->nonb_cut;
}
else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
- control->nbr_cut = val;
+ control->bond_cut = val;
}
else if ( strncmp(tmp[0], "thb_cutoff", MAX_LINE) == 0 )
{
@@ -238,7 +242,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
else if ( strncmp(tmp[0], "hbond_cutoff", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
- control->hb_cut = val;
+ control->hbond_cut = val;
}
else if ( strncmp(tmp[0], "charge_method", MAX_LINE) == 0 )
{
@@ -560,7 +564,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
if (ferror(fp))
{
- fprintf(stderr, "Error reading control file. Terminating.\n");
+ fprintf( stderr, "Error reading control file. Terminating.\n" );
exit( INVALID_INPUT );
}
@@ -574,14 +578,10 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->T = control->T_init;
}
- /* near neighbor and far neighbor cutoffs */
- control->bo_cut = 0.01 * system->reaxprm.gp.l[29];
- control->r_low = system->reaxprm.gp.l[11];
- control->r_cut = system->reaxprm.gp.l[12];
- control->r_sp_cut = control->r_cut * control->cm_domain_sparsity;
- control->vlist_cut += control->r_cut;
+ /* derived cutoffs */
+ control->nonb_sp_cut = control->nonb_cut * control->cm_domain_sparsity;
- system->g.cell_size = control->vlist_cut / 2.;
+ system->g.cell_size = control->vlist_cut / 2.0;
for ( i = 0; i < 3; ++i )
{
system->g.spread[i] = 2;
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index f28343f4..1f29e730 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -56,7 +56,7 @@ void Init_Bonded_Force_Functions( control_params *control )
//*/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 )
+ if ( control->hbond_cut > 0.0 )
{
control->intr_funcs[5] = &Hydrogen_Bonds; //*/Dummy_Interaction;
}
@@ -182,7 +182,7 @@ static void Compute_Total_Force( reax_system *system, control_params *control,
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- if ( i < bonds->select.bond_list[pj].nbr )
+ if ( i < bonds->bond_list[pj].nbr )
{
if ( control->ensemble == NVE || control->ensemble == nhNVT
|| control->ensemble == bNVT )
@@ -497,18 +497,18 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
{
- if ( far_nbrs->select.far_nbr_list[pj].d <= control->r_cut )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- j = far_nbrs->select.far_nbr_list[pj].nbr;
+ j = far_nbrs->far_nbr_list[pj].nbr;
xcut = ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
+ system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 )
/ 2.0;
- if ( far_nbrs->select.far_nbr_list[pj].d < xcut &&
- far_nbrs->select.far_nbr_list[pj].d > 0.001 )
+ if ( far_nbrs->far_nbr_list[pj].d < xcut &&
+ far_nbrs->far_nbr_list[pj].d > 0.001 )
{
- d = far_nbrs->select.far_nbr_list[pj].d / xcut;
+ d = far_nbrs->far_nbr_list[pj].d / xcut;
bond_softness = system->reaxprm.gp.l[34] * POW( d, 3.0 ) * POW( 1.0 - d, 6.0 );
H->j[*Htop] = system->N + j + 1;
@@ -631,49 +631,58 @@ static void Init_Forces( reax_system *system, control_params *control,
H_sp->start[i] = H_sp_top;
btop_i = End_Index( i, bonds );
sbp_i = &system->reaxprm.sbp[type_i];
- ihb = ihb_top = -1;
+ ihb = -1;
+ ihb_top = -1;
- if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
+ if ( control->hbond_cut > 0.0 )
{
- ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ ihb = sbp_i->p_hbond;
+
+ if ( ihb == 1 )
+ {
+ ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ }
}
+ /* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->select.far_nbr_list[pj];
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
atom_j = &system->atoms[j];
+ flag = FALSE;
+ flag_sp = FALSE;
- flag = 0;
- flag_sp = 0;
+ /* check if reneighboring step */
if ( (data->step - data->prev_steps) % control->reneighbor == 0 )
{
- if ( nbr_pj->d <= control->r_cut )
+ if ( nbr_pj->d <= control->nonb_cut )
{
- flag = 1;
- if ( nbr_pj->d <= control->r_sp_cut )
+ flag = TRUE;
+ if ( nbr_pj->d <= control->nonb_sp_cut )
{
- flag_sp = 1;
+ flag_sp = TRUE;
}
}
- else
- {
- flag = 0;
- flag_sp = 0;
- }
}
- else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &(system->box),
- nbr_pj->dvec)) <= SQR(control->r_cut))
+ else
{
- if ( nbr_pj->d <= SQR(control->r_sp_cut))
+ nbr_pj->d = Sq_Distance_on_T3( atom_i->x, atom_j->x,
+ &system->box, nbr_pj->dvec );
+
+ if ( nbr_pj->d <= SQR(control->nonb_cut) )
{
- flag_sp = 1;
+ if ( nbr_pj->d <= SQR(control->nonb_sp_cut) )
+ {
+ flag_sp = TRUE;
+ }
+
+ nbr_pj->d = SQRT( nbr_pj->d );
+ flag = TRUE;
}
- nbr_pj->d = SQRT( nbr_pj->d );
- flag = 1;
}
- if ( flag )
+ if ( flag == TRUE )
{
type_j = system->atoms[j].type;
r_ij = nbr_pj->d;
@@ -686,7 +695,7 @@ static void Init_Forces( reax_system *system, control_params *control,
++Htop;
/* H_sp matrix entry */
- if ( flag_sp )
+ if ( flag_sp == TRUE )
{
H_sp->j[H_sp_top] = j;
H_sp->val[H_sp_top] = H->val[Htop - 1];
@@ -694,32 +703,33 @@ static void Init_Forces( reax_system *system, control_params *control,
}
/* hydrogen bond lists */
- if ( control->hb_cut > 0 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hb_cut )
+ if ( control->hbond_cut > 0.0 && (ihb == 1 || ihb == 2)
+ && nbr_pj->d <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
+
if ( ihb == 1 && jhb == 2 )
{
- hbonds->select.hbond_list[ihb_top].nbr = j;
- hbonds->select.hbond_list[ihb_top].scl = 1;
- hbonds->select.hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].nbr = j;
+ hbonds->hbond_list[ihb_top].scl = 1;
+ hbonds->hbond_list[ihb_top].ptr = nbr_pj;
++ihb_top;
++num_hbonds;
}
else if ( ihb == 2 && jhb == 1 )
{
jhb_top = End_Index( workspace->hbond_index[j], hbonds );
- hbonds->select.hbond_list[jhb_top].nbr = i;
- hbonds->select.hbond_list[jhb_top].scl = -1;
- hbonds->select.hbond_list[jhb_top].ptr = nbr_pj;
+ hbonds->hbond_list[jhb_top].nbr = i;
+ hbonds->hbond_list[jhb_top].scl = -1;
+ hbonds->hbond_list[jhb_top].ptr = nbr_pj;
Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
++num_hbonds;
}
}
/* uncorrected bond orders */
- if ( far_nbrs->select.far_nbr_list[pj].d <= control->nbr_cut )
+ if ( nbr_pj->d <= control->bond_cut )
{
r2 = SQR( r_ij );
@@ -730,8 +740,8 @@ static void Init_Forces( reax_system *system, control_params *control,
}
else
{
- BO_s = 0.0;
C12 = 0.0;
+ BO_s = 0.0;
}
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
@@ -741,8 +751,8 @@ static void Init_Forces( reax_system *system, control_params *control,
}
else
{
- BO_pi = 0.0;
C34 = 0.0;
+ BO_pi = 0.0;
}
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
@@ -752,8 +762,8 @@ static void Init_Forces( reax_system *system, control_params *control,
}
else
{
- BO_pi2 = 0.0;
C56 = 0.0;
+ BO_pi2 = 0.0;
}
/* Initially BO values are the uncorrected ones, page 1 */
@@ -763,9 +773,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->bond_list[btop_i];
btop_j = End_Index( j, bonds );
- jbond = &bonds->select.bond_list[btop_j];
+ jbond = &bonds->bond_list[btop_j];
ibond->nbr = j;
jbond->nbr = i;
@@ -783,11 +793,15 @@ static void Init_Forces( reax_system *system, control_params *control,
Set_End_Index( j, btop_j + 1, bonds );
bo_ij = &ibond->bo_data;
+ bo_ij->BO = BO;
+ bo_ij->BO_s = BO_s;
+ bo_ij->BO_pi = BO_pi;
+ bo_ij->BO_pi2 = BO_pi2;
bo_ji = &jbond->bo_data;
- bo_ji->BO = bo_ij->BO = BO;
- bo_ji->BO_s = bo_ij->BO_s = BO_s;
- bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
- bo_ji->BO_pi2 = bo_ij->BO_pi2 = BO_pi2;
+ bo_ji->BO = BO;
+ bo_ji->BO_s = BO_s;
+ bo_ji->BO_pi = BO_pi;
+ bo_ji->BO_pi2 = BO_pi2;
/* Bond Order page2-3, derivative of total bond order prime */
Cln_BOp_s = twbp->p_bo2 * C12 / r2;
@@ -799,17 +813,16 @@ static void Init_Forces( reax_system *system, control_params *control,
rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi2,
- -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
+ -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
/* Only dBOp wrt. dr_i is stored here, note that
dBOp/dr_i = -dBOp/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dBOp,
- -(bo_ij->BO_s * Cln_BOp_s +
- bo_ij->BO_pi * Cln_BOp_pi +
- bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+ rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
+ + bo_ij->BO_pi * Cln_BOp_pi
+ + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
@@ -915,54 +928,64 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
H_sp->start[i] = H_sp_top;
btop_i = End_Index( i, bonds );
sbp_i = &system->reaxprm.sbp[type_i];
- ihb = ihb_top = -1;
+ ihb = -1;
+ ihb_top = -1;
- if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
+ if ( control->hbond_cut > 0.0 )
{
- ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ ihb = sbp_i->p_hbond;
+
+ if ( ihb == 1 )
+ {
+ ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ }
}
+ /* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->select.far_nbr_list[pj];
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
atom_j = &system->atoms[j];
+ flag = FALSE;
+ flag_sp = FALSE;
- flag = 0;
- flag_sp = 0;
+ /* check if reneighboring step */
if ( (data->step - data->prev_steps) % control->reneighbor == 0 )
{
- if (nbr_pj->d <= control->r_cut)
+ if (nbr_pj->d <= control->nonb_cut)
{
- flag = 1;
- if ( nbr_pj->d <= control->r_sp_cut )
+ flag = TRUE;
+ if ( nbr_pj->d <= control->nonb_sp_cut )
{
- flag_sp = 1;
+ flag_sp = TRUE;
}
}
- else
- {
- flag = 0;
- flag_sp = 0;
- }
}
else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &system->box,
- nbr_pj->dvec)) <= SQR(control->r_cut))
+ nbr_pj->dvec)) <= SQR(control->nonb_cut))
{
- if ( nbr_pj->d <= SQR(control->r_sp_cut))
+ nbr_pj->d = Sq_Distance_on_T3( atom_i->x, atom_j->x,
+ &system->box, nbr_pj->dvec );
+
+ if ( nbr_pj->d <= SQR(control->nonb_cut) )
{
- flag_sp = 1;
+ if ( nbr_pj->d <= SQR(control->nonb_sp_cut))
+ {
+ flag_sp = TRUE;
+ }
+
+ nbr_pj->d = SQRT( nbr_pj->d );
+ flag = TRUE;
}
- nbr_pj->d = SQRT( nbr_pj->d );
- flag = 1;
}
- if ( flag )
+ if ( flag == TRUE )
{
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,
@@ -970,7 +993,7 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
++Htop;
/* H_sp matrix entry */
- if ( flag_sp )
+ if ( flag_sp == TRUE )
{
H_sp->j[H_sp_top] = j;
H_sp->val[H_sp_top] = H->val[Htop - 1];
@@ -978,32 +1001,33 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
}
/* hydrogen bond lists */
- if ( control->hb_cut > 0 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hb_cut )
+ if ( control->hbond_cut > 0 && (ihb == 1 || ihb == 2)
+ && nbr_pj->d <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
+
if ( ihb == 1 && jhb == 2 )
{
- hbonds->select.hbond_list[ihb_top].nbr = j;
- hbonds->select.hbond_list[ihb_top].scl = 1;
- hbonds->select.hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].nbr = j;
+ hbonds->hbond_list[ihb_top].scl = 1;
+ hbonds->hbond_list[ihb_top].ptr = nbr_pj;
++ihb_top;
++num_hbonds;
}
else if ( ihb == 2 && jhb == 1 )
{
jhb_top = End_Index( workspace->hbond_index[j], hbonds );
- hbonds->select.hbond_list[jhb_top].nbr = i;
- hbonds->select.hbond_list[jhb_top].scl = -1;
- hbonds->select.hbond_list[jhb_top].ptr = nbr_pj;
+ hbonds->hbond_list[jhb_top].nbr = i;
+ hbonds->hbond_list[jhb_top].scl = -1;
+ hbonds->hbond_list[jhb_top].ptr = nbr_pj;
Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
++num_hbonds;
}
}
/* uncorrected bond orders */
- if ( far_nbrs->select.far_nbr_list[pj].d <= control->nbr_cut )
+ if ( nbr_pj->d <= control->bond_cut )
{
r2 = SQR( r_ij );
@@ -1014,8 +1038,8 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
}
else
{
- BO_s = 0.0;
C12 = 0.0;
+ BO_s = 0.0;
}
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
@@ -1025,8 +1049,8 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
}
else
{
- BO_pi = 0.0;
C34 = 0.0;
+ BO_pi = 0.0;
}
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
@@ -1036,8 +1060,8 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
}
else
{
- BO_pi2 = 0.0;
C56 = 0.0;
+ BO_pi2 = 0.0;
}
/* Initially BO values are the uncorrected ones, page 1 */
@@ -1047,9 +1071,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->bond_list[btop_i];
btop_j = End_Index( j, bonds );
- jbond = &bonds->select.bond_list[btop_j];
+ jbond = &bonds->bond_list[btop_j];
ibond->nbr = j;
jbond->nbr = i;
@@ -1066,12 +1090,16 @@ static void Init_Forces_Tab( 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_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;
- bo_ji->BO_pi2 = bo_ij->BO_pi2 = BO_pi2;
+ bo_ij = &ibond->bo_data;
+ bo_ij->BO = BO;
+ bo_ij->BO_s = BO_s;
+ bo_ij->BO_pi = BO_pi;
+ bo_ij->BO_pi2 = BO_pi2;
+ bo_ji = &jbond->bo_data;
+ bo_ji->BO = BO;
+ bo_ji->BO_s = BO_s;
+ bo_ji->BO_pi = BO_pi;
+ bo_ji->BO_pi2 = BO_pi2;
/* Bond Order page2-3, derivative of total bond order prime */
Cln_BOp_s = twbp->p_bo2 * C12 / r2;
@@ -1083,17 +1111,16 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
rvec_Scale( bo_ij->dln_BOp_pi2,
- -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
+ -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
/* Only dBOp wrt. dr_i is stored here, note that
dBOp/dr_i = -dBOp/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dBOp,
- -(bo_ij->BO_s * Cln_BOp_s +
- bo_ij->BO_pi * Cln_BOp_pi +
- bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+ rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
+ + bo_ij->BO_pi * Cln_BOp_pi
+ + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
@@ -1178,27 +1205,28 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
{
atom_i = &system->atoms[i];
type_i = atom_i->type;
- start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
+ start_i = Start_Index( i, far_nbrs );
+ end_i = End_Index( i, far_nbrs );
sbp_i = &system->reaxprm.sbp[type_i];
ihb = sbp_i->p_hbond;
+ /* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->select.far_nbr_list[pj];
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
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];
- if ( nbr_pj->d <= control->r_cut )
+ if ( nbr_pj->d <= control->nonb_cut )
{
++(*Htop);
/* hydrogen bond lists */
- if ( control->hb_cut > 0.1 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hb_cut )
+ if ( control->hbond_cut > 0.0 && (ihb == 1 || ihb == 2)
+ && nbr_pj->d <= control->hbond_cut )
{
jhb = sbp_j->p_hbond;
@@ -1213,7 +1241,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
}
/* uncorrected bond orders */
- if ( nbr_pj->d <= control->nbr_cut )
+ if ( nbr_pj->d <= control->bond_cut )
{
r_ij = nbr_pj->d;
@@ -1268,7 +1296,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
hb_top[i] = MAX( hb_top[i] * SAFE_HBONDS, MIN_HBONDS );
- *num_3body += SQR(bond_top[i]);
+ *num_3body += SQR( bond_top[i] );
bond_top[i] = MAX( bond_top[i] * 2, MIN_BONDS );
}
*num_3body *= SAFE_ZONE;
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index 23e24408..23896294 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -239,7 +239,7 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
for ( pk = start_j; pk < end_j; ++pk )
{
- pbond_jk = &( bonds->select.bond_list[pk] );
+ pbond_jk = &( bonds->bond_list[pk] );
k = pbond_jk->nbr;
bo_jk = &( pbond_jk->bo_data );
BOA_jk = bo_jk->BO - control->thb_cut;
@@ -280,10 +280,10 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
/* pick i up from j-k interaction where j is the centre atom */
for ( pi = start_pk; pi < end_pk; ++pi )
{
- p_ijk = &( thb_intrs->select.three_body_list[pi] );
+ p_ijk = &thb_intrs->three_body_list[pi];
pij = p_ijk->pthb; // pij is pointer to i on j's bond_list
- pbond_ij = &( bonds->select.bond_list[pij] );
- bo_ij = &( pbond_ij->bo_data );
+ pbond_ij = &bonds->bond_list[pij];
+ bo_ij = &pbond_ij->bo_data;
if ( bo_ij->BO > control->thb_cut/*0*/ )
{
@@ -321,15 +321,15 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
/* pick l up from j-k intr. where k is the centre */
for ( pl = start_pj; pl < end_pj; ++pl )
{
- p_jkl = &( thb_intrs->select.three_body_list[pl] );
+ p_jkl = &thb_intrs->three_body_list[pl];
l = p_jkl->thb;
plk = p_jkl->pthb; //pointer to l on k's bond_list!
- pbond_kl = &( bonds->select.bond_list[plk] );
- bo_kl = &( pbond_kl->bo_data );
+ pbond_kl = &bonds->bond_list[plk];
+ bo_kl = &pbond_kl->bo_data;
type_l = system->atoms[l].type;
- fbh = &(system->reaxprm.fbp[type_i][type_j][type_k][type_l]);
- fbp = &(system->reaxprm.fbp[type_i][type_j]
- [type_k][type_l].prm[0]);
+ fbh = &system->reaxprm.fbp[type_i][type_j][type_k][type_l];
+ fbp = &system->reaxprm.fbp[type_i][type_j]
+ [type_k][type_l].prm[0];
if ( i != l && fbh->cnt && bo_kl->BO > control->thb_cut &&
bo_ij->BO * bo_jk->BO * bo_kl->BO > control->thb_cut )
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 45e36abd..2e09f60b 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -56,7 +56,7 @@ static void Count_Geo_Atoms( FILE *geo, reax_system *system )
static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
{
- int ret;
+ int ret, cryst_len;
char *cryst;
char line[MAX_LINE + 1];
char descriptor[9];
@@ -73,10 +73,12 @@ static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
{
case PDB:
cryst = "CRYST1";
+ cryst_len = 6;
break;
default:
cryst = "BOX";
+ cryst_len = 3;
break;
}
@@ -84,7 +86,7 @@ static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
initialize the big box */
while ( fgets( line, MAX_LINE, geo ) )
{
- if ( strncmp( line, cryst, 6 ) == 0 )
+ if ( strncmp( line, cryst, cryst_len ) == 0 )
{
if ( geo_format == PDB )
{
@@ -100,6 +102,7 @@ static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
&system->box );
ret = SUCCESS;
+ break;
}
}
}
@@ -172,15 +175,14 @@ static void Count_PDB_Atoms( FILE *geo, reax_system *system )
char s_x[9], s_y[9], s_z[9];
rvec x;
- /* initialize variables */
- fseek( geo, 0, SEEK_SET ); /* set the pointer to the beginning of the file */
+ fseek( geo, 0, SEEK_SET );
system->N = 0;
/* increment number of atoms for each line denoting an atom desc */
while ( fgets( line, MAX_LINE, geo ) )
{
- if ( strncmp( line, "ATOM", 4 ) == 0 ||
- strncmp( line, "HETATM", 6 ) == 0 )
+ if ( strncmp( line, "ATOM", 4 ) == 0
+ || strncmp( line, "HETATM", 6 ) == 0 )
{
system->N++;
@@ -231,8 +233,8 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
/* read box information */
if ( Read_Box_Info( system, pdb, PDB ) == FAILURE )
{
- fprintf( stderr, "Read_Box_Info: no CRYST line in the pdb file!" );
- fprintf( stderr, "terminating...\n" );
+ fprintf( stderr, "[ERROR] Read_Box_Info: no CRYST line in the pdb file!" );
+ fprintf( stderr, " Terminating...\n" );
exit( INVALID_GEO );
}
@@ -289,7 +291,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
strncpy( &charge[0], s1 + 78, 2 );
charge[2] = 0;
}
- else if (strncmp(tmp[0], "HETATM", 6) == 0)
+ else if ( strncmp(tmp[0], "HETATM", 6) == 0 )
{
strncpy( &descriptor[0], s1, 6 );
descriptor[6] = 0;
@@ -324,36 +326,38 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
}
/* if the point is inside my_box, add it to my lists */
- Make_Point( strtod( &s_x[0], &endptr ),
- strtod( &s_y[0], &endptr ),
- strtod( &s_z[0], &endptr ), &x );
+ Make_Point( strtod( &s_x[0], &endptr ), strtod( &s_y[0], &endptr ),
+ strtod( &s_z[0], &endptr ), &x );
- Fit_to_Periodic_Box( &(system->box), &x );
+ Fit_to_Periodic_Box( &system->box, &x );
- /* store orig_id, type, name and coord info of the new atom */
- atom = &(system->atoms[top]);
- pdb_serial = (int) strtod( &serial[0], &endptr );
- workspace->orig_id[top] = pdb_serial;
+ if ( is_Inside_Box( &system->box, x ) )
+ {
+ /* store orig_id, type, name and coord info of the new atom */
+ atom = &system->atoms[top];
+ pdb_serial = (int) strtod( &serial[0], &endptr );
+ workspace->orig_id[top] = pdb_serial;
- Trim_Spaces( element, 9 );
- atom->type = Get_Atom_Type( &(system->reaxprm), element );
- strncpy( atom->name, atom_name, 8 );
+ Trim_Spaces( element, 9 );
+ atom->type = Get_Atom_Type( &system->reaxprm, element );
+ strncpy( atom->name, atom_name, 8 );
- rvec_Copy( atom->x, x );
- rvec_MakeZero( atom->v );
- rvec_MakeZero( atom->f );
- atom->q = 0;
+ rvec_Copy( atom->x, x );
+ rvec_MakeZero( atom->v );
+ rvec_MakeZero( atom->f );
+ atom->q = 0;
- top++;
+ top++;
+ }
c++;
}
/* IMPORTANT: We do not check for the soundness of restrictions here.
- When atom2 is on atom1's restricted list, and there is a restriction
- on atom2, then atom1 has to be on atom2's restricted list, too.
- However, we do not check if this is the case in the input file,
- this is upto the user. */
- else if (!strncmp( tmp[0], "CONECT", 6 ))
+ * When atom2 is on atom1's restricted list, and there is a restriction
+ * on atom2, then atom1 has to be on atom2's restricted list, too.
+ * However, we do not check if this is the case in the input file,
+ * this is upto the user. */
+ else if ( !strncmp( tmp[0], "CONECT", 6 ) )
{
if ( control->restrict_bonds )
{
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 47b4ecb4..af1d3517 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -152,6 +152,7 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
case nhNVT:
data->N_f = 3 * system->N + 1;
+ *Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
//control->Tau_T = 100 * data->N_f * K_B * control->T_final;
if ( !control->restart || (control->restart && control->random_vel) )
@@ -168,8 +169,6 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
data->N_f, data->therm.v_xi );
#endif
}
-
- *Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
break;
/* anisotropic NPT */
@@ -178,6 +177,7 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
exit( UNKNOWN_OPTION );
data->N_f = 3 * system->N + 9;
+ *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
if ( !control->restart )
{
@@ -188,8 +188,6 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
//data->inv_W = 1. / (data->N_f*K_B*control->T*SQR(control->Tau_P));
//Compute_Pressure( system, data, workspace );
}
-
- *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
break;
/* semi-isotropic NPT */
@@ -238,8 +236,8 @@ static void Init_Taper( control_params *control, static_storage *workspace )
real swa, swa2, swa3;
real swb, swb2, swb3;
- swa = control->r_low;
- swb = control->r_cut;
+ swa = control->nonb_low;
+ swb = control->nonb_cut;
if ( FABS( swa ) > 0.01 )
{
@@ -779,7 +777,7 @@ static void Init_Lists( reax_system *system, control_params *control,
#endif
workspace->num_H = 0;
- if ( control->hb_cut > 0.0 )
+ if ( control->hbond_cut > 0.0 )
{
/* init H indexes */
for ( i = 0; i < system->N; ++i )
@@ -797,7 +795,7 @@ static void Init_Lists( reax_system *system, control_params *control,
if ( workspace->num_H == 0 )
{
- control->hb_cut = 0.0;
+ control->hbond_cut = 0.0;
}
else
{
@@ -1439,7 +1437,7 @@ 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] );
- if ( control->hb_cut > 0.0 )
+ if ( control->hbond_cut > 0.0 )
{
Delete_List( TYP_HBOND, lists[HBONDS] );
}
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 5698379c..62d4fc64 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -57,11 +57,11 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
rvec_ScaledSum( dx, dt, system->atoms[i].v,
- -0.5 * dt_sqr * F_CONV * inv_m, system->atoms[i].f );
- Inc_on_T3( system->atoms[i].x, dx, &( system->box ) );
-
+ 0.5 * dt_sqr * -F_CONV * inv_m, system->atoms[i].f );
+// Inc_on_T3( system->atoms[i].x, dx, &system->box );
+ rvec_Add( system->atoms[i].x, dx );
rvec_ScaledAdd( system->atoms[i].v,
- -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ 0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
}
#if defined(DEBUG_FOCUS)
@@ -69,6 +69,7 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
#endif
Reallocate( system, control, workspace, lists, renbr );
+
Reset( system, control, data, workspace, lists );
if ( renbr )
@@ -83,7 +84,7 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
{
inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
rvec_ScaledAdd( system->atoms[i].v,
- -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ 0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
}
#if defined(DEBUG_FOCUS)
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index 85805eb9..02246843 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -35,48 +35,48 @@ void Make_List( int n, int total_intrs, int type, reax_list* l )
switch ( type )
{
case TYP_VOID:
- l->select.v = smalloc( l->total_intrs * sizeof(void),
- "Make_List::l->select.v" );
+ l->v = smalloc( l->total_intrs * sizeof(void),
+ "Make_List::l->v" );
break;
case TYP_THREE_BODY:
- l->select.three_body_list = smalloc( l->total_intrs * sizeof(three_body_interaction_data),
- "Make_List::l->select.three_body_list" );
+ l->three_body_list = smalloc( l->total_intrs * sizeof(three_body_interaction_data),
+ "Make_List::l->three_body_list" );
break;
case TYP_BOND:
- l->select.bond_list = smalloc( l->total_intrs * sizeof(bond_data),
- "Make_List::l->select.bond_list" );
+ l->bond_list = smalloc( l->total_intrs * sizeof(bond_data),
+ "Make_List::l->bond_list" );
break;
case TYP_DBO:
- l->select.dbo_list = smalloc( l->total_intrs * sizeof(dbond_data),
- "Make_List::l->select.dbo_list" );
+ l->dbo_list = smalloc( l->total_intrs * sizeof(dbond_data),
+ "Make_List::l->dbo_list" );
break;
case TYP_DDELTA:
- l->select.dDelta_list = smalloc( l->total_intrs * sizeof(dDelta_data),
- "Make_List::l->select.dDelta_list" );
+ l->dDelta_list = smalloc( l->total_intrs * sizeof(dDelta_data),
+ "Make_List::l->dDelta_list" );
break;
case TYP_FAR_NEIGHBOR:
- l->select.far_nbr_list = smalloc( l->total_intrs * sizeof(far_neighbor_data),
- "Make_List::l->select.far_nbr_list" );
+ l->far_nbr_list = smalloc( l->total_intrs * sizeof(far_neighbor_data),
+ "Make_List::l->far_nbr_list" );
break;
case TYP_NEAR_NEIGHBOR:
- l->select.near_nbr_list = smalloc( l->total_intrs * sizeof(near_neighbor_data),
- "Make_List::l->select.near_nbr_list" );
+ l->near_nbr_list = smalloc( l->total_intrs * sizeof(near_neighbor_data),
+ "Make_List::l->near_nbr_list" );
break;
case TYP_HBOND:
- l->select.hbond_list = smalloc( l->total_intrs * sizeof(hbond_data),
- "Make_List::l->select.hbond_list" );
+ l->hbond_list = smalloc( l->total_intrs * sizeof(hbond_data),
+ "Make_List::l->hbond_list" );
break;
default:
- l->select.v = smalloc( l->total_intrs * sizeof(void),
- "Make_List::l->select.v" );
+ l->v = smalloc( l->total_intrs * sizeof(void),
+ "Make_List::l->v" );
break;
}
}
@@ -90,35 +90,35 @@ void Delete_List( int type, reax_list* l )
switch ( type )
{
case TYP_VOID:
- sfree( l->select.v, "Delete_List::l->select.v" );
+ sfree( l->v, "Delete_List::l->v" );
break;
case TYP_THREE_BODY:
- sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
+ sfree( l->three_body_list, "Delete_List::l->three_body_list" );
break;
case TYP_BOND:
- sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
+ sfree( l->bond_list, "Delete_List::l->bond_list" );
break;
case TYP_DBO:
- sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
+ sfree( l->dbo_list, "Delete_List::l->dbo_list" );
break;
case TYP_DDELTA:
- sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
+ sfree( l->dDelta_list, "Delete_List::l->dDelta_list" );
break;
case TYP_FAR_NEIGHBOR:
- sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
+ sfree( l->far_nbr_list, "Delete_List::l->far_nbr_list" );
break;
case TYP_NEAR_NEIGHBOR:
- sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
+ sfree( l->near_nbr_list, "Delete_List::l->near_nbr_list" );
break;
case TYP_HBOND:
- sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
+ sfree( l->hbond_list, "Delete_List::l->hbond_list" );
break;
default:
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 6af97014..332b1e14 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -245,7 +245,7 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
v0_vdw = 0;
num_atom_types = system->reaxprm.num_atom_types;
- dr = control->r_cut / control->tabulate;
+ dr = control->nonb_cut / control->tabulate;
h = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::h" );
fh = scalloc( (control->tabulate + 1), sizeof(real),
@@ -292,10 +292,10 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
if ( existing_types[j] )
{
workspace->LR[i][j].xmin = 0;
- workspace->LR[i][j].xmax = control->r_cut;
+ workspace->LR[i][j].xmax = control->nonb_cut;
workspace->LR[i][j].n = control->tabulate + 1;
workspace->LR[i][j].dx = dr;
- workspace->LR[i][j].inv_dx = control->tabulate / control->r_cut;
+ workspace->LR[i][j].inv_dx = control->tabulate / control->nonb_cut;
workspace->LR[i][j].y =
smalloc( workspace->LR[i][j].n * sizeof(LR_data),
"Make_LR_Lookup_Table::LR[i][j].y" );
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 8a6d5a42..aab95c34 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -61,11 +61,10 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
far_neighbor_data nbr_data;
g = &system->g;
+ num_far = 0;
Bin_Atoms( system, workspace );
- num_far = 0;
-
/* first pick up a cell in the grid */
for ( i = 0; i < g->ncell[0]; i++ )
{
@@ -149,9 +148,9 @@ void Choose_Neighbor_Finder( reax_system *system, control_params *control,
if ( control->periodic_boundaries )
{
#ifdef SMALL_BOX_SUPPORT
- if ( system->box.box_norms[0] > 2.0 * control->vlist_cut
- && system->box.box_norms[1] > 2.0 * control->vlist_cut
- && system->box.box_norms[2] > 2.0 * control->vlist_cut )
+ if ( system->box.box_norms[0] >= 2.0 * control->vlist_cut
+ && system->box.box_norms[1] >= 2.0 * control->vlist_cut
+ && system->box.box_norms[2] >= 2.0 * control->vlist_cut )
{
*Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Big_Box;
}
@@ -245,7 +244,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
if ( atom1 > atom2 )
{
- nbr_data = &far_nbrs->select.far_nbr_list[num_far];
+ nbr_data = &far_nbrs->far_nbr_list[num_far];
count = Find_Far_Neighbors( system->atoms[atom1].x, system->atoms[atom2].x,
atom1, atom2, &system->box, control->vlist_cut, nbr_data );
@@ -281,7 +280,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
#if defined(DEBUG)
for ( i = 0; i < system->N; ++i )
{
- qsort( &far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ],
+ qsort( &far_nbrs->far_nbr_list[ Start_Index(i, far_nbrs) ],
Num_Entries(i, far_nbrs), sizeof(far_neighbor_data),
compare_far_nbrs );
}
@@ -304,7 +303,7 @@ static inline int is_Near_Neighbor( reax_list *near_nbrs, int atom1, int atom2 )
for ( i = Start_Index(atom1, near_nbrs); i < End_Index(atom1, near_nbrs); ++i )
{
- if ( near_nbrs->select.near_nbr_list[i].nbr == atom2 )
+ if ( near_nbrs->near_nbr_list[i].nbr == atom2 )
{
return FALSE;
}
@@ -408,7 +407,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
if ( atom1 >= atom2 )
{
- nbr_data = &far_nbrs->select.far_nbr_list[num_far];
+ nbr_data = &far_nbrs->far_nbr_list[num_far];
//top_near1 = End_Index( atom1, near_nbrs );
//Set_Start_Index( atom1, num_far, far_nbrs );
@@ -461,11 +460,11 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
system->atoms[ atom2 ].x, &system->box, control,
new_nbrs, &count );
- Set_Near_Neighbor( &near_nbrs->select.near_nbr_list[ top_near1 ], atom2,
+ Set_Near_Neighbor( &near_nbrs->near_nbr_list[ top_near1 ], atom2,
new_nbrs[c].d, 1.0, new_nbrs[c].dvec, new_nbrs[c].rel_box );
++top_near1;
- Set_Near_Neighbor( &near_nbrs->select.near_nbr_list[ top_near2 ], atom1,
+ Set_Near_Neighbor( &near_nbrs->near_nbr_list[ top_near2 ], atom1,
new_nbrs[c].d, -1.0, new_nbrs[c].dvec, new_nbrs[c].rel_box );
Set_End_Index( atom2, top_near2 + 1, near_nbrs );
}
@@ -504,7 +503,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
- qsort( &near_nbrs->select.near_nbr_list[ Start_Index(i, near_nbrs) ],
+ qsort( &near_nbrs->near_nbr_list[ Start_Index(i, near_nbrs) ],
Num_Entries(i, near_nbrs), sizeof(near_neighbor_data),
compare_near_nbrs );
}
@@ -512,7 +511,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
#ifdef TEST_ENERGY
// for ( i = 0; i < system->N; ++i )
// {
-// qsort( &far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ],
+// qsort( &far_nbrs->far_nbr_list[ Start_Index(i, far_nbrs) ],
// Num_Entries(i, far_nbrs), sizeof(far_neighbor_data), compare_far_nbrs );
// }
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 22aceadf..6f75f431 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -54,13 +54,13 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- bo_ij = &(bonds->select.bond_list[pj].bo_data);
+ bo_ij = &(bonds->bond_list[pj].bo_data);
fprintf( out_control->fbo, "%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
//workspace->orig_id[i],
- //workspace->orig_id[bonds->select.bond_list[pj].nbr],
+ //workspace->orig_id[bonds->bond_list[pj].nbr],
i + 1,
- bonds->select.bond_list[pj].nbr + 1,
- bonds->select.bond_list[pj].d,
+ bonds->bond_list[pj].nbr + 1,
+ bonds->bond_list[pj].d,
bo_ij->BO, bo_ij->BO_s, bo_ij->BO_pi, bo_ij->BO_pi2 );
}
}
@@ -75,29 +75,29 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
{
/*fprintf( out_control->fdbo, "%6d %6d\tstart: %6d\tend: %6d\n",
workspace->orig_id[i],
- workspace->orig_id[bonds->select.bond_list[pj].nbr],
+ workspace->orig_id[bonds->bond_list[pj].nbr],
Start_Index( pj, dBOs ), End_Index( pj, dBOs ) );*/
for ( pk = Start_Index(pj, dBOs); pk < End_Index(pj, dBOs); ++pk )
{
- dbo_k = &(dBOs->select.dbo_list[pk]);
+ dbo_k = &dBOs->dbo_list[pk];
//if( !rvec_isZero( dbo_k->dBO ) )
fprintf( out_control->fdbo, "%6d%6d%6d%23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
- workspace->orig_id[bonds->select.bond_list[pj].nbr],
+ workspace->orig_id[bonds->bond_list[pj].nbr],
workspace->orig_id[dbo_k->wrt],
dbo_k->dBO[0], dbo_k->dBO[1], dbo_k->dBO[2] );
fprintf( out_control->fdbo, "%6d%6d%6d%23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
- workspace->orig_id[bonds->select.bond_list[pj].nbr],
+ workspace->orig_id[bonds->bond_list[pj].nbr],
workspace->orig_id[dbo_k->wrt],
dbo_k->dBOpi[0], dbo_k->dBOpi[1], dbo_k->dBOpi[2] );
fprintf( out_control->fdbo, "%6d%6d%6d%23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
- workspace->orig_id[bonds->select.bond_list[pj].nbr],
+ workspace->orig_id[bonds->bond_list[pj].nbr],
workspace->orig_id[dbo_k->wrt],
dbo_k->dBOpi2[0], dbo_k->dBOpi2[1], dbo_k->dBOpi2[2] );
}
@@ -415,15 +415,15 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
for ( j = Start_Index(i, near_nbrs); j < End_Index(i, near_nbrs); ++j )
{
- id_j = workspace->orig_id[near_nbrs->select.near_nbr_list[j].nbr];
+ id_j = workspace->orig_id[near_nbrs->near_nbr_list[j].nbr];
// if( id_i < id_j )
fprintf( fout, "%6d%6d%23.15e%23.15e%23.15e%23.15e\n",
id_i, id_j,
- near_nbrs->select.near_nbr_list[j].d,
- near_nbrs->select.near_nbr_list[j].dvec[0],
- near_nbrs->select.near_nbr_list[j].dvec[1],
- near_nbrs->select.near_nbr_list[j].dvec[2] );
+ near_nbrs->near_nbr_list[j].d,
+ near_nbrs->near_nbr_list[j].dvec[0],
+ near_nbrs->near_nbr_list[j].dvec[1],
+ near_nbrs->near_nbr_list[j].dvec[2] );
}
}
@@ -449,15 +449,15 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
fprintf( fout, "%6d:", id_i);
for ( j = Start_Index(i, near_nbrs); j < End_Index(i, near_nbrs); ++j )
{
- id_j = workspace->orig_id[near_nbrs->select.near_nbr_list[j].nbr];
+ id_j = workspace->orig_id[near_nbrs->near_nbr_list[j].nbr];
fprintf( fout, "%6d", id_j);
/* fprintf( fout, "%6d%6d%23.15e%23.15e%23.15e%23.15e\n",
id_i, id_j,
- near_nbrs->select.near_nbr_list[j].d,
- near_nbrs->select.near_nbr_list[j].dvec[0],
- near_nbrs->select.near_nbr_list[j].dvec[1],
- near_nbrs->select.near_nbr_list[j].dvec[2] ); */
+ near_nbrs->near_nbr_list[j].d,
+ near_nbrs->near_nbr_list[j].dvec[0],
+ near_nbrs->near_nbr_list[j].dvec[1],
+ near_nbrs->near_nbr_list[j].dvec[2] ); */
}
fprintf( fout, "\n");
}
@@ -484,21 +484,21 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
for ( j = Start_Index(i, far_nbrs); j < End_Index(i, far_nbrs); ++j )
{
- id_j = workspace->orig_id[far_nbrs->select.far_nbr_list[j].nbr];
+ id_j = workspace->orig_id[far_nbrs->far_nbr_list[j].nbr];
fprintf( fout, "%6d%6d%23.15e%23.15e%23.15e%23.15e\n",
id_i, id_j,
- far_nbrs->select.far_nbr_list[j].d,
- far_nbrs->select.far_nbr_list[j].dvec[0],
- far_nbrs->select.far_nbr_list[j].dvec[1],
- far_nbrs->select.far_nbr_list[j].dvec[2] );
+ far_nbrs->far_nbr_list[j].d,
+ far_nbrs->far_nbr_list[j].dvec[0],
+ far_nbrs->far_nbr_list[j].dvec[1],
+ far_nbrs->far_nbr_list[j].dvec[2] );
fprintf( fout, "%6d%6d%23.15e%23.15e%23.15e%23.15e\n",
id_j, id_i,
- far_nbrs->select.far_nbr_list[j].d,
- -far_nbrs->select.far_nbr_list[j].dvec[0],
- -far_nbrs->select.far_nbr_list[j].dvec[1],
- -far_nbrs->select.far_nbr_list[j].dvec[2] );
+ far_nbrs->far_nbr_list[j].d,
+ -far_nbrs->far_nbr_list[j].dvec[0],
+ -far_nbrs->far_nbr_list[j].dvec[1],
+ -far_nbrs->far_nbr_list[j].dvec[2] );
}
}
@@ -533,7 +533,7 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
for ( j = Start_Index(i, far_nbrs); j < End_Index(i, far_nbrs); ++j )
{
- id_j = workspace->orig_id[far_nbrs->select.far_nbr_list[j].nbr];
+ id_j = workspace->orig_id[far_nbrs->far_nbr_list[j].nbr];
temp[num++] = id_j;
}
qsort(&temp, num, sizeof(int), fn_qsort_intcmp);
@@ -974,8 +974,8 @@ void Print_Bonds( reax_system *system, reax_list *bonds, char *fname )
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- pbond = &(bonds->select.bond_list[pj]);
- bo_ij = &(pbond->bo_data);
+ pbond = &bonds->bond_list[pj];
+ bo_ij = &pbond->bo_data;
//fprintf( f, "%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
// i+1, pbond->nbr+1, pbond->d,
// bo_ij->BO, bo_ij->BO_s, bo_ij->BO_pi, bo_ij->BO_pi2 );
@@ -1002,7 +1002,7 @@ void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
fprintf( f, "%6d:", id_i);
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- nbr = bonds->select.bond_list[pj].nbr;
+ nbr = bonds->bond_list[pj].nbr;
id_j = nbr + 1; //system->my_atoms[nbr].orig_id;
if ( id_i < id_j )
{
@@ -1010,7 +1010,7 @@ void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
}
}
- qsort(&temp, num, sizeof(int), fn_qsort_intcmp);
+ qsort( &temp, num, sizeof(int), fn_qsort_intcmp );
for ( j = 0; j < num; j++ )
{
fprintf( f, "%6d", temp[j] );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 87be5e87..3355a85c 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -46,7 +46,7 @@
#define REORDER_ATOMS
/* enables support for small simulation boxes (i.e. a simulation box with any
* dimension less than twice the Verlet list cutoff distance, vlist_cut) */
-#define SMALL_BOX_SUPPORT
+//#define SMALL_BOX_SUPPORT
#define SUCCESS (1)
#define FAILURE (0)
@@ -562,6 +562,8 @@ struct simulation_box
{
real volume;
+ rvec min;
+ rvec max;
rvec box_norms;
rvec side_prop;
rvec nbr_box_press[27];
@@ -632,22 +634,40 @@ struct control_params
* 4 : iNPT (Parrinello-Rehman-Nose-Hoover) isotropic
* 5 : aNPT (Parrinello-Rehman-Nose-Hoover) anisotropic */
int ensemble;
+ /* number of simulation time steps */
int nsteps;
+ /* */
int periodic_boundaries;
+ /* */
int restrict_bonds;
+ /* */
int tabulate;
+ /* simulation time step length (in fs) */
real dt;
-
+ /* number of simulation steps to elapse before
+ * recomputing the verlet lists */
int reneighbor;
+ /* cutoff (in Angstroms) used for for constructing the
+ * long range interaction Verlet list (a.k.a. far neighbor list) */
real vlist_cut;
- real nbr_cut;
- /* upper and lower taper */
- real r_cut;
- real r_sp_cut;
- real r_low;
+ /* bonded interaction cutoff (in Angstroms) */
+ real bond_cut;
+ /* nonbonded interaction cutoff (in Angstroms),
+ * also used as upper taper radius */
+ real nonb_cut;
+ /* shorter version of nonbonded interaction cutoff (in Angstroms),
+ * used for computing sparser charge matrix
+ * for electrostatic interactions */
+ real nonb_sp_cut;
+ /* shorter version of nonbonded interaction cutoff (in Angstroms),
+ * used for lower taper radius */
+ real nonb_low;
+ /* bond order cutoff (in Angstroms) */
real bo_cut;
+ /* three body (valence angle) cutoff (in Angstroms) */
real thb_cut;
- real hb_cut;
+ /* hydrogen bonding cutoff (in Angstroms) */
+ real hbond_cut;
real T_init;
real T_final;
@@ -1223,25 +1243,41 @@ struct reax_list
/* max. num. of interactions per atom */
int *max_intrs;
/* interaction list (polymorphic via union dereference) */
- union
- {
- /* typeless interaction list */
- void *v;
- /* three body interaction list */
- three_body_interaction_data *three_body_list;
- /* bond interaction list */
- bond_data *bond_list;
- /* bond interaction list */
- dbond_data *dbo_list;
- /* test forces interaction list */
- dDelta_data *dDelta_list;
- /* far neighbor interaction list */
- far_neighbor_data *far_nbr_list;
- /* near neighbor interaction list */
- near_neighbor_data *near_nbr_list;
- /* hydrogen bond interaction list */
- hbond_data *hbond_list;
- } select;
+// union
+// {
+// /* typeless interaction list */
+// void *v;
+// /* three body interaction list */
+// three_body_interaction_data *three_body_list;
+// /* bond interaction list */
+// bond_data *bond_list;
+// /* bond interaction list */
+// dbond_data *dbo_list;
+// /* test forces interaction list */
+// dDelta_data *dDelta_list;
+// /* far neighbor interaction list */
+// far_neighbor_data *far_nbr_list;
+// /* near neighbor interaction list */
+// near_neighbor_data *near_nbr_list;
+// /* hydrogen bond interaction list */
+// hbond_data *hbond_list;
+// } select;
+ /* typeless interaction list */
+ void *v;
+ /* three body interaction list */
+ three_body_interaction_data *three_body_list;
+ /* bond interaction list */
+ bond_data *bond_list;
+ /* bond interaction list */
+ dbond_data *dbo_list;
+ /* test forces interaction list */
+ dDelta_data *dDelta_list;
+ /* far neighbor interaction list */
+ far_neighbor_data *far_nbr_list;
+ /* near neighbor interaction list */
+ near_neighbor_data *near_nbr_list;
+ /* hydrogen bond interaction list */
+ hbond_data *hbond_list;
};
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index aa58850f..6a6a808b 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -130,7 +130,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
Set_End_Index( i, tmp, bonds );
}
- if ( control->hb_cut > 0 )
+ if ( control->hbond_cut > 0 )
{
for ( i = 0; i < system->N; ++i )
{
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 667c7ae3..192b45eb 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -95,15 +95,15 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- if ( i < bonds->select.bond_list[pj].nbr )
+ if ( i < bonds->bond_list[pj].nbr )
{
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
type_j = system->atoms[j].type;
if ( !strncmp( system->reaxprm.sbp[type_j].name, "C", 15 ) )
{
- twbp = &( system->reaxprm.tbp[type_i][type_j]);
- bo_ij = &( bonds->select.bond_list[pj].bo_data );
+ twbp = &system->reaxprm.tbp[type_i][type_j];
+ bo_ij = &bonds->bond_list[pj].bo_data;
Di = workspace->Delta[i];
vov3 = bo_ij->BO - Di - 0.040 * POW(Di, 4.);
@@ -155,10 +155,10 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
type_j = system->atoms[j].type;
- bo_ij = &(bonds->select.bond_list[pj].bo_data);
- twbp = &(system->reaxprm.tbp[ type_i ][ type_j ]);
+ bo_ij = &bonds->bond_list[pj].bo_data;
+ twbp = &system->reaxprm.tbp[ type_i ][ type_j ];
sum_ovun1 += twbp->p_ovun1 * twbp->De_s * bo_ij->BO;
sum_ovun2 += (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]) *
@@ -220,11 +220,11 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
- pbond = &(bonds->select.bond_list[pj]);
+ pbond = &bonds->bond_list[pj];
j = pbond->nbr;
type_j = system->atoms[j].type;
- bo_ij = &(pbond->bo_data);
- twbp = &(system->reaxprm.tbp[ type_i ][ type_j ]);
+ 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]) *
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 26724de8..ab7b7e4c 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -44,9 +44,8 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
rvec diff, cross;
/* remove rotational and translational velocity of the center of mass */
- if ( control->ensemble != NVE &&
- control->remove_CoM_vel &&
- data->step && data->step % control->remove_CoM_vel == 0 )
+ if ( control->ensemble != NVE && control->remove_CoM_vel
+ && data->step % control->remove_CoM_vel == 0 )
{
/* compute velocity of the center of mass */
Compute_Center_of_Mass( system, data, out_control->prs );
@@ -227,10 +226,16 @@ 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 )
{
+ if ( (spmd_handle->out_control->energy_update_freq > 0
+ && spmd_handle->data->step % spmd_handle->out_control->energy_update_freq == 0)
+ || (spmd_handle->out_control->write_steps > 0
+ && spmd_handle->data->step % spmd_handle->out_control->write_steps == 0) )
+ {
+ Compute_Total_Energy( spmd_handle->system, spmd_handle->data );
+ }
+
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
@@ -242,10 +247,9 @@ int simulate( const void * const handle )
spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
spmd_handle->lists );
}
- ++spmd_handle->data->step;
//}
- for ( ; spmd_handle->data->step <= spmd_handle->control->nsteps; spmd_handle->data->step++ )
+ for ( ++spmd_handle->data->step; spmd_handle->data->step <= spmd_handle->control->nsteps; spmd_handle->data->step++ )
{
if ( spmd_handle->control->T_mode )
{
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index a7bf0a2f..412843d4 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -30,7 +30,8 @@ void Temperature_Control( control_params *control, simulation_data *data,
{
real tmp;
- if ( control->T_mode == 1 ) // step-wise temperature control
+ /* step-wise temperature control */
+ if ( control->T_mode == 1 )
{
if ( (data->step - data->prev_steps) %
((int)(control->T_freq / control->dt)) == 0 )
@@ -45,7 +46,8 @@ void Temperature_Control( control_params *control, simulation_data *data,
}
}
}
- else if ( control->T_mode == 2 ) // constant slope control
+ /* constant slope control */
+ else if ( control->T_mode == 2 )
{
tmp = control->T_rate * control->dt / control->T_freq;
@@ -195,15 +197,12 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
rvec_Scale( p, m, system->atoms[i].v );
data->E_Kin += 0.5 * rvec_Dot( p, system->atoms[i].v );
-
- /* fprintf(stderr,"%d, %lf, %lf, %lf %lf\n",
- i,system->atoms[i].v[0], system->atoms[i].v[1], system->atoms[i].v[2],
- system->reaxprm.sbp[system->atoms[i].type].mass); */
}
- data->therm.T = (2. * data->E_Kin) / (data->N_f * K_B);
+ data->therm.T = (2.0 * data->E_Kin) / (data->N_f * K_B);
- if ( FABS(data->therm.T) < ALMOST_ZERO ) /* avoid T being an absolute zero! */
+ /* avoid T being an absolute zero! */
+ if ( FABS( data->therm.T ) < ALMOST_ZERO )
{
data->therm.T = ALMOST_ZERO;
}
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 41973c57..545d4b80 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -95,9 +95,9 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
total_bo = workspace->total_bond_order;
bonds = lists[BONDS];
- bond_list = bonds->select.bond_list;
+ bond_list = bonds->bond_list;
thb_intrs = lists[THREE_BODIES];
- thb_list = thb_intrs->select.three_body_list;
+ thb_list = thb_intrs->three_body_list;
/* global parameters used in these calculations */
p_pen2 = system->reaxprm.gp.l[19];
p_pen3 = system->reaxprm.gp.l[20];
@@ -158,9 +158,9 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
//#ifdef _OPENMP
-// f_j = &(workspace->f_local[tid * system->N + j]);
+// f_j = &workspace->f_local[tid * system->N + j];
//#else
- f_j = &(system->atoms[j].f);
+ f_j = &system->atoms[j].f;
//#endif
p_val3 = system->reaxprm.sbp[ type_j ].p_val3;
@@ -168,10 +168,11 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
SBOp = 0.0;
prod_SBO = 1.0;
+
for ( t = start_j; t < end_j; ++t )
{
- bo_jt = &(bond_list[t].bo_data);
- SBOp += (bo_jt->BO_pi + bo_jt->BO_pi2);
+ bo_jt = &bond_list[t].bo_data;
+ SBOp += bo_jt->BO_pi + bo_jt->BO_pi2;
temp = SQR( bo_jt->BO );
temp *= temp;
temp *= temp;
@@ -222,8 +223,8 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
for ( pi = start_j; pi < end_j; ++pi )
{
Set_Start_Index( pi, num_thb_intrs, thb_intrs );
- pbond_ij = &(bond_list[pi]);
- bo_ij = &(pbond_ij->bo_data);
+ pbond_ij = &bond_list[pi];
+ bo_ij = &pbond_ij->bo_data;
BOA_ij = bo_ij->BO - control->thb_cut;
if ( BOA_ij > 0.0 )
@@ -233,7 +234,7 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
//#ifdef _OPENMP
// f_i = &(workspace->f_local[tid * system->N + i]);
//#else
- f_i = &(system->atoms[i].f);
+ f_i = &system->atoms[i].f;
//#endif
/* first copy 3-body intrs from previously computed ones where i>k.
@@ -272,16 +273,16 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
/* and this is the second for loop mentioned above */
for ( pk = pi + 1; pk < end_j; ++pk )
{
- pbond_jk = &(bond_list[pk]);
- bo_jk = &(pbond_jk->bo_data);
+ pbond_jk = &bond_list[pk];
+ bo_jk = &pbond_jk->bo_data;
BOA_jk = bo_jk->BO - control->thb_cut;
k = pbond_jk->nbr;
type_k = system->atoms[k].type;
- p_ijk = &( thb_list[num_thb_intrs] );
+ p_ijk = &thb_list[num_thb_intrs];
//#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
//CHANGE ORIGINAL
@@ -312,10 +313,9 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
++num_thb_intrs;
- if ( BOA_jk > 0.0 &&
- (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut)/*0*/)
+ if ( BOA_jk > 0.0 && (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut) )
{
- thbh = &( system->reaxprm.thbp[type_i][type_j][type_k] );
+ thbh = &system->reaxprm.thbp[type_i][type_j][type_k];
/* if( workspace->orig_id[i] < workspace->orig_id[k] )
fprintf( stdout, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
@@ -330,7 +330,7 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
{
if ( FABS(thbh->prm[cnt].p_val1) > 0.001 )
{
- thbp = &( thbh->prm[cnt] );
+ thbp = &thbh->prm[cnt];
/* ANGLE ENERGY */
p_val1 = thbp->p_val1;
@@ -373,11 +373,11 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
CEval1 = Cf7ij * f7_jk * f8_Dj * expval12theta;
CEval2 = Cf7jk * f7_ij * f8_Dj * expval12theta;
CEval3 = Cf8j * f7_ij * f7_jk * expval12theta;
- CEval4 = -2.0 * p_val1 * p_val2 * f7_ij * f7_jk * f8_Dj *
- expval2theta * (theta_0 - theta);
+ CEval4 = -2.0 * p_val1 * p_val2 * f7_ij * f7_jk * f8_Dj
+ * expval2theta * (theta_0 - theta);
- Ctheta_0 = p_val10 * DEG2RAD(theta_00) *
- EXP( -p_val10 * (2.0 - SBO2) );
+ Ctheta_0 = p_val10 * DEG2RAD(theta_00)
+ * EXP( -p_val10 * (2.0 - SBO2) );
CEval5 = -CEval4 * Ctheta_0 * CSBO2;
CEval6 = CEval5 * dSBO1;
@@ -452,8 +452,8 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
for ( t = start_j; t < end_j; ++t )
{
- pbond_jt = &( bond_list[t] );
- bo_jt = &(pbond_jt->bo_data);
+ pbond_jt = &bond_list[t];
+ bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
@@ -465,7 +465,7 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
#ifdef _OPENMP
// #pragma omp atomic
#endif
- bo_jt->Cdbo += (CEval6 * pBOjt7);
+ bo_jt->Cdbo += CEval6 * pBOjt7;
#ifdef _OPENMP
// #pragma omp atomic
#endif
@@ -593,8 +593,8 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
for ( t = start_j; t < end_j; ++t )
{
- pbond_jt = &( bond_list[t] );
- bo_jt = &(pbond_jt->bo_data);
+ pbond_jt = &bond_list[t];
+ bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
@@ -708,9 +708,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
num_hb_intrs = 0;
#endif
bonds = lists[BONDS];
- bond_list = bonds->select.bond_list;
+ bond_list = bonds->bond_list;
hbonds = lists[HBONDS];
- hbond_list = hbonds->select.hbond_list;
+ hbond_list = hbonds->hbond_list;
/* loops below discover the Hydrogen bonds between i-j-k triplets.
here j is H atom and there has to be some bond between i and j.
@@ -732,17 +732,17 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
hb_start_j = Start_Index( workspace->hbond_index[j], hbonds );
hb_end_j = End_Index( workspace->hbond_index[j], hbonds );
#ifdef _OPENMP
- f_j = &(workspace->f_local[tid * system->N + j]);
+ f_j = &workspace->f_local[tid * system->N + j];
#else
- f_j = &(system->atoms[j].f);
+ f_j = &system->atoms[j].f;
#endif
top = 0;
for ( pi = start_j; pi < end_j; ++pi )
{
- pbond_ij = &( bond_list[pi] );
+ pbond_ij = &bond_list[pi];
i = pbond_ij->nbr;
- bo_ij = &(pbond_ij->bo_data);
+ bo_ij = &pbond_ij->bo_data;
type_i = system->atoms[i].type;
if ( system->reaxprm.sbp[type_i].p_hbond == 2 &&
@@ -761,9 +761,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
r_jk = nbr_jk->d;
rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
#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
for ( itr = 0; itr < top; ++itr )
@@ -774,14 +774,14 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
if ( i != k )
{
- bo_ij = &(pbond_ij->bo_data);
+ bo_ij = &pbond_ij->bo_data;
type_i = system->atoms[i].type;
r_ij = pbond_ij->d;
- hbp = &(system->reaxprm.hbp[ type_i ][ type_j ][ type_k ]);
+ hbp = &system->reaxprm.hbp[ type_i ][ type_j ][ type_k ];
#ifdef _OPENMP
- f_i = &(workspace->f_local[tid * system->N + i]);
+ f_i = &workspace->f_local[tid * system->N + i];
#else
- f_i = &(system->atoms[i].f);
+ f_i = &system->atoms[i].f;
#endif
#ifdef TEST_FORCES
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index c46f8615..bec55b60 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -74,66 +74,41 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec *p )
for ( i = 0; i < 3; ++i )
{
- //TODO: verify box boundary coordinates -- assuming orthogonal box pinned at origin
- if ( (*p)[i] < 0. )
+ if ( (*p)[i] < box->min[i] )
{
/* handle lower coords */
- while ( (*p)[i] < 0. )
+ while ( (*p)[i] < box->min[i] )
+ {
(*p)[i] += box->box_norms[i];
+ }
}
- else if ( (*p)[i] >= box->box_norms[i] )
+ else if ( (*p)[i] >= box->max[i] )
{
/* handle higher coords */
- while ( (*p)[i] >= box->box_norms[i] )
+ while ( (*p)[i] >= box->max[i] )
+ {
(*p)[i] -= box->box_norms[i];
+ }
}
-// if ( (*p)[i] < box->min[i] )
-// {
-// /* handle lower coords */
-// while ( (*p)[i] < box->min[i] )
-// (*p)[i] += box->box_norms[i];
-// }
-// else if ( (*p)[i] >= box->max[i] )
-// {
-// /* handle higher coords */
-// while ( (*p)[i] >= box->max[i] )
-// (*p)[i] -= box->box_norms[i];
-// }
}
}
-/* determine the touch point, tp, of a box to
- its neighbor denoted by the relative coordinate rl */
-/*
-static inline void Box_Touch_Point( simulation_box *box, ivec rl, rvec tp )
-{
- int d;
-
- for ( d = 0; d < 3; ++d )
- if ( rl[d] == -1 )
- tp[d] = box->min[d];
- else if ( rl[d] == 0 )
- tp[d] = NEG_INF - 1.;
- else
- tp[d] = box->max[d];
-}
-*/
-
-
/* determine whether point p is inside the box */
/* assumes orthogonal box */
-/*
-static inline int is_Inside_Box( simulation_box *box, rvec p )
+int is_Inside_Box( simulation_box *box, rvec p )
{
- if ( p[0] < box->min[0] || p[0] >= box->max[0] ||
- p[1] < box->min[1] || p[1] >= box->max[1] ||
- p[2] < box->min[2] || p[2] >= box->max[2] )
- return FALSE;
+ int ret = TRUE;
- return TRUE;
+ if ( p[0] < box->min[0] || p[0] >= box->max[0]
+ || p[1] < box->min[1] || p[1] >= box->max[1]
+ || p[2] < box->min[2] || p[2] >= box->max[2] )
+ {
+ ret = FALSE;
+ }
+
+ return ret;
}
-*/
/*
@@ -365,16 +340,13 @@ void Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
{
int i;
- *line = (char*) smalloc( sizeof(char) * MAX_LINE,
- "Allocate_Tokenizer_Space::*line" );
- *backup = (char*) smalloc( sizeof(char) * MAX_LINE,
- "Allocate_Tokenizer_Space::*backup" );
- *tokens = (char**) smalloc( sizeof(char*) * MAX_TOKENS,
- "Allocate_Tokenizer_Space::*tokens" );
+ *line = smalloc( sizeof(char) * MAX_LINE, "Allocate_Tokenizer_Space::*line" );
+ *backup = smalloc( sizeof(char) * MAX_LINE, "Allocate_Tokenizer_Space::*backup" );
+ *tokens = smalloc( sizeof(char*) * MAX_TOKENS, "Allocate_Tokenizer_Space::*tokens" );
for ( i = 0; i < MAX_TOKENS; i++ )
{
- (*tokens)[i] = (char*) smalloc(sizeof(char) * MAX_TOKEN_LEN,
+ (*tokens)[i] = smalloc( sizeof(char) * MAX_TOKEN_LEN,
"Allocate_Tokenizer_Space::(*tokens)[i]" );
}
}
@@ -404,7 +376,8 @@ int Tokenize( char* s, char*** tok )
strncpy( test, s, MAX_LINE );
- for ( word = strtok_r(test, sep, &saveptr); word; word = strtok_r(NULL, sep, &saveptr) )
+ for ( word = strtok_r(test, sep, &saveptr); word != NULL;
+ word = strtok_r(NULL, sep, &saveptr) )
{
strncpy( (*tok)[count], word, MAX_LINE );
count++;
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index e0931a18..247c516a 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -32,9 +32,7 @@ void Transform_to_UnitBox( rvec, simulation_box*, char, rvec );
void Fit_to_Periodic_Box( simulation_box*, rvec* );
-//void Box_Touch_Point( simulation_box*, ivec, rvec );
-
-//int is_Inside_Box( simulation_box*, rvec );
+int is_Inside_Box( simulation_box*, rvec );
//int iown_midpoint( simulation_box*, rvec, rvec );
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index 43cd6c70..d65e7827 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -53,12 +53,12 @@ int Write_Custom_Header( reax_system *system, control_params *control,
control->reposition_atoms,
control->restrict_bonds,
control->tabulate,
- control->nbr_cut,
- control->r_cut,
+ control->bond_cut,
+ control->nonb_cut,
control->bg_cut,
control->bo_cut,
control->thb_cut,
- control->hb_cut,
+ control->hbond_cut,
control->cm_solver_q_err,
control->T_init,
control->T_final,
@@ -230,8 +230,8 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
{
for ( j = Start_Index( i, bonds ); j < End_Index( i, bonds ); ++j )
{
- if ( i < bonds->select.bond_list[j].nbr &&
- bonds->select.bond_list[j].bo_data.BO >= control->bg_cut )
+ if ( i < bonds->bond_list[j].nbr &&
+ bonds->bond_list[j].bo_data.BO >= control->bg_cut )
{
++num_bonds;
}
@@ -259,19 +259,19 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
{
for ( pi = Start_Index(j, bonds); pi < End_Index(j, bonds); ++pi )
{
- if ( bonds->select.bond_list[pi].bo_data.BO >= control->bg_cut )
+ if ( bonds->bond_list[pi].bo_data.BO >= control->bg_cut )
{
// physical j&i bond
for ( pk = Start_Index( pi, thb_intrs );
pk < End_Index( pi, thb_intrs ); ++pk )
{
- if ( bonds->select.bond_list[pi].nbr <
- thb_intrs->select.three_body_list[pk].thb )
+ if ( bonds->bond_list[pi].nbr <
+ thb_intrs->three_body_list[pk].thb )
{
// get k's pointer on j's bond list
- pk_j = thb_intrs->select.three_body_list[pk].pthb;
+ pk_j = thb_intrs->three_body_list[pk].pthb;
- if ( bonds->select.bond_list[pk_j].bo_data.BO >= control->bg_cut )
+ if ( bonds->bond_list[pk_j].bo_data.BO >= control->bg_cut )
{
// physical j&k bond
++num_thb_intrs;
@@ -408,10 +408,10 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
{
for ( j = Start_Index( i, bonds ); j < End_Index( i, bonds ); ++j )
{
- if ( i < bonds->select.bond_list[j].nbr &&
- bonds->select.bond_list[j].bo_data.BO >= control->bg_cut )
+ if ( i < bonds->bond_list[j].nbr &&
+ bonds->bond_list[j].bo_data.BO >= control->bg_cut )
{
- bo_ij = &( bonds->select.bond_list[j] );
+ bo_ij = &bonds->bond_list[j];
out_control->write( out_control->trj, BOND_BASIC,
workspace->orig_id[i],
workspace->orig_id[bo_ij->nbr],
@@ -426,10 +426,10 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
{
for ( j = Start_Index( i, bonds ); j < End_Index( i, bonds ); ++j )
{
- if ( i < bonds->select.bond_list[j].nbr &&
- bonds->select.bond_list[j].bo_data.BO >= control->bg_cut )
+ if ( i < bonds->bond_list[j].nbr &&
+ bonds->bond_list[j].bo_data.BO >= control->bg_cut )
{
- bo_ij = &( bonds->select.bond_list[j] );
+ bo_ij = &bonds->bond_list[j];
out_control->write( out_control->trj, BOND_FULL,
workspace->orig_id[i],
workspace->orig_id[bo_ij->nbr],
@@ -454,26 +454,26 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
{
for ( pi = Start_Index(j, bonds); pi < End_Index(j, bonds); ++pi )
{
- if ( bonds->select.bond_list[pi].bo_data.BO >= control->bg_cut )
+ if ( bonds->bond_list[pi].bo_data.BO >= control->bg_cut )
{
// physical j&i bond
for ( pk = Start_Index( pi, thb_intrs );
pk < End_Index( pi, thb_intrs ); ++pk )
{
- if ( bonds->select.bond_list[pi].nbr <
- thb_intrs->select.three_body_list[pk].thb )
+ if ( bonds->bond_list[pi].nbr <
+ thb_intrs->three_body_list[pk].thb )
{
- pk_j = thb_intrs->select.three_body_list[pk].pthb;
+ pk_j = thb_intrs->three_body_list[pk].pthb;
// get k's pointer on j's bond list
- if ( bonds->select.bond_list[pk_j].bo_data.BO >= control->bg_cut )
+ if ( bonds->bond_list[pk_j].bo_data.BO >= control->bg_cut )
{
// physical j&k bond
out_control->write( out_control->trj, ANGLE_BASIC,
- workspace->orig_id[bonds->select.bond_list[pi].nbr],
+ workspace->orig_id[bonds->bond_list[pi].nbr],
workspace->orig_id[j],
- workspace->orig_id[thb_intrs->select.three_body_list[pk].thb],
- RAD2DEG(thb_intrs->select.three_body_list[pk].theta) );
+ workspace->orig_id[thb_intrs->three_body_list[pk].thb],
+ RAD2DEG(thb_intrs->three_body_list[pk].theta) );
}
}
}
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index 2a179cab..e6316f81 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -67,16 +67,16 @@ void Bond_Energy( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- if ( i < bonds->select.bond_list[pj].nbr )
+ if ( i < bonds->bond_list[pj].nbr )
{
/* set the pointers */
- j = bonds->select.bond_list[pj].nbr;
+ j = bonds->bond_list[pj].nbr;
type_i = system->atoms[i].type;
type_j = system->atoms[j].type;
- sbp_i = &( system->reaxprm.sbp[type_i] );
- sbp_j = &( system->reaxprm.sbp[type_j] );
- twbp = &( system->reaxprm.tbp[type_i][type_j] );
- bo_ij = &( bonds->select.bond_list[pj].bo_data );
+ sbp_i = &system->reaxprm.sbp[type_i];
+ sbp_j = &system->reaxprm.sbp[type_j];
+ twbp = &system->reaxprm.tbp[type_i][type_j];
+ bo_ij = &bonds->bond_list[pj].bo_data;
/* calculate the constants */
pow_BOs_be2 = POW( bo_ij->BO_s, twbp->p_be2 );
@@ -218,9 +218,9 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- if ( far_nbrs->select.far_nbr_list[pj].d <= control->r_cut )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- nbr_pj = &far_nbrs->select.far_nbr_list[pj];
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
r_ij = nbr_pj->d;
twbp = &system->reaxprm.tbp[ system->atoms[i].type ]
@@ -563,16 +563,16 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- if ( far_nbrs->select.far_nbr_list[pj].d <= control->r_cut )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
type_j = system->atoms[j].type;
r_ij = nbr_pj->d;
self_coef = (i == j) ? 0.5 : 1.0;
tmin = MIN( type_i, type_j );
tmax = MAX( type_i, type_j );
- t = &( workspace->LR[tmin][tmax] );
+ t = &workspace->LR[tmin][tmax];
/* Cubic Spline Interpolation */
r = (int)(r_ij * t->inv_dx);
--
GitLab