From ef2bcd8a82cd900d19655d948c3a749736e8049c Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Sun, 23 Dec 2018 00:35:42 -0500
Subject: [PATCH] PuReMD-old: change from array of structs to struct of arrays
for far nbr list (far_neighbor_data).
---
PuReMD/src/bond_orders.c | 28 +++--
PuReMD/src/bond_orders.h | 2 +-
PuReMD/src/forces.c | 213 ++++++++++++++++++------------------
PuReMD/src/hydrogen_bonds.c | 24 ++--
PuReMD/src/io_tools.c | 25 +++--
PuReMD/src/list.c | 61 ++++++-----
PuReMD/src/neighbors.c | 12 +-
PuReMD/src/nonbonded.c | 105 ++++++++++--------
PuReMD/src/parallelreax.c | 1 -
PuReMD/src/reax_types.h | 21 +++-
10 files changed, 261 insertions(+), 231 deletions(-)
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index da07f839..f01993b8 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -663,11 +663,10 @@ void Add_dBond_to_Forces( int i, int pj,
int BOp( storage *workspace, reax_list *bonds, real bo_cut,
- int i, int btop_i, far_neighbor_data *nbr_pj,
+ int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec,
int far_nbr_list_format, single_body_parameters *sbp_i,
single_body_parameters *sbp_j, two_body_parameters *twbp )
{
- int j;
real r2, C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
@@ -677,26 +676,25 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
bond_data *jbond;
bond_order_data *bo_ji;
- j = nbr_pj->nbr;
- r2 = SQR(nbr_pj->d);
+ r2 = SQR(d);
if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 )
{
- C12 = twbp->p_bo1 * pow( nbr_pj->d / twbp->r_s, twbp->p_bo2 );
+ C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + bo_cut) * exp( C12 );
}
else BO_s = C12 = 0.0;
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
{
- C34 = twbp->p_bo3 * pow( nbr_pj->d / twbp->r_p, twbp->p_bo4 );
+ C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 );
BO_pi = exp( C34 );
}
else BO_pi = C34 = 0.0;
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
{
- C56 = twbp->p_bo5 * pow( nbr_pj->d / twbp->r_pp, twbp->p_bo6 );
+ C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = exp( C56 );
}
else BO_pi2 = C56 = 0.0;
@@ -707,25 +705,25 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
if ( BO >= bo_cut )
{
/****** bonds i-j and j-i ******/
- ibond = &( bonds->bond_list[btop_i] );
+ ibond = &bonds->bond_list[btop_i];
if ( far_nbr_list_format == HALF_LIST )
{
btop_j = End_Index( j, bonds );
- jbond = &(bonds->bond_list[btop_j]);
+ jbond = &bonds->bond_list[btop_j];
}
ibond->nbr = j;
- ibond->d = nbr_pj->d;
- rvec_Copy( ibond->dvec, nbr_pj->dvec );
- ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
+ ibond->d = d;
+ rvec_Copy( ibond->dvec, *dvec );
+ ivec_Copy( ibond->rel_box, *rel_box );
ibond->dbond_index = btop_i;
if ( far_nbr_list_format == HALF_LIST )
{
ibond->sym_index = btop_j;
jbond->nbr = i;
- jbond->d = nbr_pj->d;
- rvec_Scale( jbond->dvec, -1.0, nbr_pj->dvec );
- ivec_Scale( jbond->rel_box, -1.0, nbr_pj->rel_box );
+ jbond->d = d;
+ rvec_Scale( jbond->dvec, -1.0, *dvec );
+ ivec_Scale( jbond->rel_box, -1.0, *rel_box );
jbond->dbond_index = btop_i;
jbond->sym_index = btop_i;
diff --git a/PuReMD/src/bond_orders.h b/PuReMD/src/bond_orders.h
index a692fccb..fe142fb0 100644
--- a/PuReMD/src/bond_orders.h
+++ b/PuReMD/src/bond_orders.h
@@ -52,7 +52,7 @@ void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
void Add_dBond_to_Forces( int, int, storage*, reax_list** );
void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
storage*, reax_list** );
-int BOp( storage*, reax_list*, real, int, int, far_neighbor_data*,
+int BOp( storage*, reax_list*, real, int, int, int, ivec*, real, rvec*,
int, single_body_parameters*, single_body_parameters*,
two_body_parameters* );
void BO( reax_system*, control_params*, simulation_data*,
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 0f836ff4..d9e4d993 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -340,33 +340,30 @@ static void Init_Distance( reax_system *system, control_params *control,
int start_i, end_i;
int renbr;
reax_list *far_nbrs;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
far_nbrs = lists[FAR_NBRS];
-
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
for ( i = 0; i < system->N; ++i )
{
atom_i = &system->my_atoms[i];
- 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 );
/* update i-j distance for non-reneighboring steps */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
if ( !renbr )
{
- nbr_pj->dvec[0] = atom_j->x[0] - atom_i->x[0];
- nbr_pj->dvec[1] = atom_j->x[1] - atom_i->x[1];
- nbr_pj->dvec[2] = atom_j->x[2] - atom_i->x[2];
- nbr_pj->d = rvec_Norm_Sqr( nbr_pj->dvec );
- nbr_pj->d = sqrt( nbr_pj->d );
+ far_nbrs->far_nbr_list.dvec[pj][0] = atom_j->x[0] - atom_i->x[0];
+ far_nbrs->far_nbr_list.dvec[pj][1] = atom_j->x[1] - atom_i->x[1];
+ far_nbrs->far_nbr_list.dvec[pj][2] = atom_j->x[2] - atom_i->x[2];
+ far_nbrs->far_nbr_list.d[pj] = rvec_Norm_Sqr( far_nbrs->far_nbr_list.dvec[pj] );
+ far_nbrs->far_nbr_list.d[pj] = sqrt( far_nbrs->far_nbr_list.d[pj] );
}
}
@@ -389,7 +386,6 @@ static void Init_CM_Half_NT( reax_system *system, control_params *control,
reax_list *far_nbrs;
single_body_parameters *sbp_i;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int mark[6];
int total_cnt[6];
@@ -482,11 +478,10 @@ static void Init_CM_Half_NT( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -498,7 +493,7 @@ static void Init_CM_Half_NT( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
+ r_ij = far_nbrs->far_nbr_list.d[pj];
twbp = &system->reax_param.tbp[type_i][type_j];
if ( local == 1 )
@@ -615,7 +610,6 @@ static void Init_CM_Full_NT( reax_system *system, control_params *control,
reax_list *far_nbrs;
single_body_parameters *sbp_i;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int mark[6];
int total_cnt[6];
@@ -708,11 +702,10 @@ static void Init_CM_Full_NT( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -724,7 +717,7 @@ static void Init_CM_Full_NT( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
+ r_ij = far_nbrs->far_nbr_list.d[pj];
twbp = &system->reax_param.tbp[type_i][type_j];
if ( local == 1 )
@@ -842,7 +835,6 @@ static void Init_CM_Half_FS( reax_system *system, control_params *control,
reax_list *far_nbrs;
single_body_parameters *sbp_i;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
far_nbrs = lists[FAR_NBRS];
@@ -881,11 +873,10 @@ static void Init_CM_Half_FS( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -897,7 +888,7 @@ static void Init_CM_Half_FS( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
+ r_ij = far_nbrs->far_nbr_list.d[pj];
twbp = &system->reax_param.tbp[type_i][type_j];
if ( local )
@@ -958,7 +949,6 @@ static void Init_CM_Full_FS( reax_system *system, control_params *control,
reax_list *far_nbrs;
single_body_parameters *sbp_i;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
far_nbrs = lists[FAR_NBRS];
@@ -997,11 +987,10 @@ static void Init_CM_Full_FS( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -1013,7 +1002,7 @@ static void Init_CM_Full_FS( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
+ r_ij = far_nbrs->far_nbr_list.d[pj];
twbp = &system->reax_param.tbp[type_i][type_j];
if ( local )
@@ -1072,7 +1061,6 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int jhb_top;
@@ -1140,11 +1128,10 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
/* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -1164,7 +1151,7 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
/* hydrogen bond lists */
if ( control->hbond_cut > 0
&& (ihb == 1 || ihb == 2)
- && nbr_pj->d <= control->hbond_cut )
+ && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
@@ -1173,7 +1160,7 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
{
hbonds->hbond_list[ihb_top].nbr = j;
hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].ptr = pj;
++ihb_top;
++num_hbonds;
}
@@ -1183,7 +1170,7 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
jhb_top = End_Index( atom_j->Hindex, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = nbr_pj;
+ hbonds->hbond_list[jhb_top].ptr = pj;
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
@@ -1192,9 +1179,11 @@ static void Init_Bond_Half( reax_system *system, control_params *control,
/* uncorrected bond orders */
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- nbr_pj->d <= control->bond_cut
+ far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
&& BOp( workspace, bonds, control->bo_cut,
- i, btop_i, nbr_pj, far_nbrs->format,
+ i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+ &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+ &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
@@ -1257,7 +1246,6 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int start_j, end_j;
int btop_j;
@@ -1324,11 +1312,10 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
/* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
atom_j = &system->my_atoms[j];
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
flag = 1;
}
@@ -1348,7 +1335,7 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
/* hydrogen bond lists */
if ( control->hbond_cut > 0
&& (ihb == 1 || ihb == 2)
- && nbr_pj->d <= control->hbond_cut )
+ && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
@@ -1357,7 +1344,7 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
{
hbonds->hbond_list[ihb_top].nbr = j;
hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].ptr = pj;
++ihb_top;
++num_hbonds;
}
@@ -1366,9 +1353,11 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
/* uncorrected bond orders */
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- nbr_pj->d <= control->bond_cut
+ far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
&& BOp( workspace, bonds, control->bo_cut,
- i, btop_i, nbr_pj, far_nbrs->format,
+ i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+ &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+ &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
@@ -1523,7 +1512,6 @@ void Init_Forces( reax_system *system, control_params *control,
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int jhb_top;
int start_j, end_j;
@@ -1674,24 +1662,26 @@ void Init_Forces( reax_system *system, control_params *control,
// update i-j distance - check if j is within cutoff
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->far_nbr_list[pj] );
- j = nbr_pj->nbr;
- atom_j = &(system->my_atoms[j]);
+ j = far_nbrs->far_nbr_list.nbr[pj];
+ atom_j = &system->my_atoms[j];
+
if ( renbr )
{
- if (nbr_pj->d <= cutoff)
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
flag = 1;
- else flag = 0;
+ else
+ flag = 0;
}
else
{
- nbr_pj->dvec[0] = atom_j->x[0] - atom_i->x[0];
- nbr_pj->dvec[1] = atom_j->x[1] - atom_i->x[1];
- nbr_pj->dvec[2] = atom_j->x[2] - atom_i->x[2];
- nbr_pj->d = rvec_Norm_Sqr( nbr_pj->dvec );
- if ( nbr_pj->d <= SQR(cutoff) )
+ far_nbrs->far_nbr_list.dvec[pj][0] = atom_j->x[0] - atom_i->x[0];
+ far_nbrs->far_nbr_list.dvec[pj][1] = atom_j->x[1] - atom_i->x[1];
+ far_nbrs->far_nbr_list.dvec[pj][2] = atom_j->x[2] - atom_i->x[2];
+ far_nbrs->far_nbr_list.d[pj] = rvec_Norm_Sqr( far_nbrs->far_nbr_list.dvec[pj] );
+
+ if ( far_nbrs->far_nbr_list.d[pj] <= SQR(cutoff) )
{
- nbr_pj->d = sqrt(nbr_pj->d);
+ far_nbrs->far_nbr_list.d[pj] = sqrt( far_nbrs->far_nbr_list.d[pj] );
flag = 1;
}
else
@@ -1703,9 +1693,9 @@ void Init_Forces( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
- sbp_j = &(system->reax_param.sbp[type_j]);
- twbp = &(system->reax_param.tbp[type_i][type_j]);
+ r_ij = far_nbrs->far_nbr_list.d[pj];
+ sbp_j = &system->reax_param.sbp[type_j];
+ twbp = &system->reax_param.tbp[type_i][type_j];
if ( local == 1 )
{
@@ -1756,8 +1746,9 @@ void Init_Forces( reax_system *system, control_params *control,
#endif
// hydrogen bond lists
- if ( control->hbond_cut > 0 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hbond_cut )
+ if ( control->hbond_cut > 0.0
+ && (ihb == 1 || ihb == 2)
+ && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
@@ -1765,7 +1756,7 @@ void Init_Forces( reax_system *system, control_params *control,
{
hbonds->hbond_list[ihb_top].nbr = j;
hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].ptr = pj;
++ihb_top;
++num_hbonds;
}
@@ -1776,7 +1767,7 @@ void Init_Forces( reax_system *system, control_params *control,
jhb_top = End_Index( atom_j->Hindex, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = nbr_pj;
+ hbonds->hbond_list[jhb_top].ptr = pj;
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
@@ -1822,9 +1813,11 @@ void Init_Forces( reax_system *system, control_params *control,
// uncorrected bond orders
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- nbr_pj->d <= control->bond_cut &&
- BOp( workspace, bonds, control->bo_cut,
- i, btop_i, nbr_pj, far_nbrs->format,
+ far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
+ && BOp( workspace, bonds, control->bo_cut,
+ i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+ &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+ &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
@@ -1942,7 +1935,6 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
int jhb_top;
int start_j, end_j;
@@ -2021,25 +2013,30 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
/* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->far_nbr_list[pj] );
- j = nbr_pj->nbr;
- atom_j = &(system->my_atoms[j]);
+ j = far_nbrs->far_nbr_list.nbr[pj];
+ atom_j = &system->my_atoms[j];
if ( renbr )
{
- if ( nbr_pj->d <= cutoff )
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
+ {
flag = 1;
- else flag = 0;
+ }
+ else
+ {
+ flag = 0;
+ }
}
else
{
- nbr_pj->dvec[0] = atom_j->x[0] - atom_i->x[0];
- nbr_pj->dvec[1] = atom_j->x[1] - atom_i->x[1];
- nbr_pj->dvec[2] = atom_j->x[2] - atom_i->x[2];
- nbr_pj->d = rvec_Norm_Sqr( nbr_pj->dvec );
- if ( nbr_pj->d <= SQR(cutoff) )
+ far_nbrs->far_nbr_list.dvec[pj][0] = atom_j->x[0] - atom_i->x[0];
+ far_nbrs->far_nbr_list.dvec[pj][1] = atom_j->x[1] - atom_i->x[1];
+ far_nbrs->far_nbr_list.dvec[pj][2] = atom_j->x[2] - atom_i->x[2];
+ far_nbrs->far_nbr_list.d[pj] = rvec_Norm_Sqr( far_nbrs->far_nbr_list.dvec[pj] );
+
+ if ( far_nbrs->far_nbr_list.d[pj] <= SQR(cutoff) )
{
- nbr_pj->d = sqrt(nbr_pj->d);
+ far_nbrs->far_nbr_list.d[pj] = sqrt( far_nbrs->far_nbr_list.d[pj] );
flag = 1;
}
else
@@ -2051,15 +2048,16 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
if ( flag )
{
type_j = atom_j->type;
- r_ij = nbr_pj->d;
- sbp_j = &(system->reax_param.sbp[type_j]);
- twbp = &(system->reax_param.tbp[type_i][type_j]);
+ r_ij = far_nbrs->far_nbr_list.d[pj];
+ sbp_j = &system->reax_param.sbp[type_j];
+ twbp = &system->reax_param.tbp[type_i][type_j];
if ( local )
{
/* hydrogen bond lists */
- if ( control->hbond_cut > 0 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hbond_cut )
+ if ( control->hbond_cut > 0.0
+ && (ihb == 1 || ihb == 2)
+ && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
{
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
@@ -2067,7 +2065,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
{
hbonds->hbond_list[ihb_top].nbr = j;
hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
+ hbonds->hbond_list[ihb_top].ptr = pj;
++ihb_top;
++num_hbonds;
}
@@ -2078,7 +2076,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
jhb_top = End_Index( atom_j->Hindex, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = nbr_pj;
+ hbonds->hbond_list[jhb_top].ptr = pj;
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
@@ -2088,9 +2086,11 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
/* uncorrected bond orders */
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- nbr_pj->d <= control->bond_cut &&
- BOp( workspace, bonds, control->bo_cut,
- i, btop_i, nbr_pj, far_nbrs->format,
+ far_nbrs->far_nbr_list.d[pj] <= control->bond_cut
+ && BOp( workspace, bonds, control->bo_cut,
+ i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+ &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+ &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
@@ -2155,14 +2155,14 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
#endif
Validate_Lists( system, workspace, lists, data->step,
- system->n, system->N, system->numH, comm );
+ system->n, system->N, system->numH, comm );
}
void Estimate_Storages( reax_system *system, control_params *control,
- reax_list **lists, int *Htop, int *hb_top,
- int *bond_top, int *num_3body, MPI_Comm comm,
- int *matrix_dim, int cm_format )
+ reax_list **lists, int *Htop, int *hb_top,
+ int *bond_top, int *num_3body, MPI_Comm comm,
+ int *matrix_dim, int cm_format )
{
int i, j, pj;
int start_i, end_i;
@@ -2176,7 +2176,6 @@ void Estimate_Storages( reax_system *system, control_params *control,
reax_list *far_nbrs;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
far_nbrs = lists[FAR_NBRS];
@@ -2224,8 +2223,8 @@ void Estimate_Storages( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &( far_nbrs->far_nbr_list[pj] );
- j = nbr_pj->nbr;
+ j = far_nbrs->far_nbr_list.nbr[pj];
+
#if !defined(NEUTRAL_TERRITORY)
if ( far_nbrs->format == HALF_LIST )
#endif
@@ -2233,12 +2232,12 @@ void Estimate_Storages( reax_system *system, control_params *control,
atom_j = &system->my_atoms[j];
}
- if (nbr_pj->d <= cutoff)
+ if ( far_nbrs->far_nbr_list.d[pj] <= cutoff )
{
type_j = system->my_atoms[j].type;
- r_ij = nbr_pj->d;
- sbp_j = &(system->reax_param.sbp[type_j]);
- twbp = &(system->reax_param.tbp[type_i][type_j]);
+ r_ij = far_nbrs->far_nbr_list.d[pj];
+ sbp_j = &system->reax_param.sbp[type_j];
+ twbp = &system->reax_param.tbp[type_i][type_j];
if ( local == 1 )
{
@@ -2257,10 +2256,12 @@ void Estimate_Storages( reax_system *system, control_params *control,
#endif
/* hydrogen bond lists */
- if ( control->hbond_cut > 0.1 && (ihb == 1 || ihb == 2) &&
- nbr_pj->d <= control->hbond_cut )
+ if ( control->hbond_cut > 0.1
+ && (ihb == 1 || ihb == 2)
+ && far_nbrs->far_nbr_list.d[pj] <= control->hbond_cut )
{
jhb = sbp_j->p_hbond;
+
if ( ihb == 1 && jhb == 2 )
{
++hb_top[i];
@@ -2286,7 +2287,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
#endif
/* uncorrected bond orders */
- if ( nbr_pj->d <= control->bond_cut )
+ if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
{
r2 = SQR(r_ij);
diff --git a/PuReMD/src/hydrogen_bonds.c b/PuReMD/src/hydrogen_bonds.c
index cba267ff..699e7bc1 100644
--- a/PuReMD/src/hydrogen_bonds.c
+++ b/PuReMD/src/hydrogen_bonds.c
@@ -39,12 +39,13 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control )
{
- int i, j, k, pi, pk;
- int type_i, type_j, type_k;
- int start_j, end_j, hb_start_j, hb_end_j;
- int hblist[MAX_BONDS];
- int itr, top;
- int num_hb_intrs = 0;
+ int i, j, k, pi, pk;
+ int type_i, type_j, type_k;
+ int start_j, end_j, hb_start_j, hb_end_j;
+ int hblist[MAX_BONDS];
+ int itr, top;
+ int num_hb_intrs = 0;
+ int nbr_jk;
ivec rel_jk;
real r_ij, r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
@@ -54,11 +55,11 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
hbond_parameters *hbp;
bond_order_data *bo_ij;
bond_data *pbond_ij;
- far_neighbor_data *nbr_jk;
- reax_list *bonds, *hbonds;
+ reax_list *far_nbrs, *bonds, *hbonds;
bond_data *bond_list;
hbond_data *hbond_list;
+ far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
bond_list = bonds->bond_list;
hbonds = lists[HBONDS];
@@ -102,8 +103,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
k = hbond_list[pk].nbr;
type_k = system->my_atoms[k].type;
nbr_jk = hbond_list[pk].ptr;
- r_jk = nbr_jk->d;
- rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
+ r_jk = far_nbrs->far_nbr_list.d[nbr_jk];
+ rvec_Scale( dvec_jk, hbond_list[pk].scl, far_nbrs->far_nbr_list.dvec[nbr_jk] );
for ( itr = 0; itr < top; ++itr )
{
@@ -174,7 +175,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
rvec_ScaledAdd( workspace->f[j], +CEhb2, dcos_theta_dj );
- ivec_Scale( rel_jk, hbond_list[pk].scl, nbr_jk->rel_box );
+ ivec_Scale( rel_jk, hbond_list[pk].scl,
+ far_nbrs->far_nbr_list.rel_box[nbr_jk] );
rvec_Scale( force, +CEhb2, dcos_theta_dk );
rvec_Add( workspace->f[k], force );
rvec_iMultiply( ext_press, rel_jk, force );
diff --git a/PuReMD/src/io_tools.c b/PuReMD/src/io_tools.c
index 308086ea..2b58ae63 100644
--- a/PuReMD/src/io_tools.c
+++ b/PuReMD/src/io_tools.c
@@ -682,20 +682,20 @@ void Print_Far_Neighbors( reax_system *system, reax_list **lists,
for ( j = Start_Index(i, far_nbrs); j < End_Index(i, far_nbrs); ++j )
{
- nbr = far_nbrs->far_nbr_list[j].nbr;
+ nbr = far_nbrs->far_nbr_list.nbr[j];
id_j = system->my_atoms[nbr].orig_id;
fprintf( fout, "%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
- id_i, id_j, 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] );
+ id_i, id_j, far_nbrs->far_nbr_list.d[j],
+ far_nbrs->far_nbr_list.dvec[j][0],
+ far_nbrs->far_nbr_list.dvec[j][1],
+ far_nbrs->far_nbr_list.dvec[j][2] );
fprintf( fout, "%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
- id_j, id_i, 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] );
+ id_j, id_i, far_nbrs->far_nbr_list.d[j],
+ -far_nbrs->far_nbr_list.dvec[j][0],
+ -far_nbrs->far_nbr_list.dvec[j][1],
+ -far_nbrs->far_nbr_list.dvec[j][2] );
}
}
@@ -865,6 +865,7 @@ void Print_HBonds( reax_system *system, reax_list **lists,
char fname[MAX_STR];
hbond_data *phbond;
FILE *fout;
+ reax_list *far_nbrs = lists[FAR_NBRS];
reax_list *hbonds = lists[HBONDS];
sprintf( fname, "%s.hbonds.%d.%d", control->sim_name, step, system->my_rank );
@@ -877,7 +878,9 @@ void Print_HBonds( reax_system *system, reax_list **lists,
phbond = &hbonds->hbond_list[pj];
fprintf( fout, "%8d%8d %24.15e %24.15e %24.15e\n", i, phbond->nbr,
- phbond->ptr->dvec[0], phbond->ptr->dvec[1], phbond->ptr->dvec[2] );
+ far_nbrs->far_nbr_list.dvec[phbond->ptr][0],
+ far_nbrs->far_nbr_list.dvec[phbond->ptr][1],
+ far_nbrs->far_nbr_list.dvec[phbond->ptr][2] );
// fprintf( fout, "%8d%8d %8d %8d\n", i, phbond->nbr,
// phbond->scl, phbond->sym_index );
}
@@ -1024,7 +1027,7 @@ void Print_Far_Neighbors_List_Adj_Format( reax_system *system,
for ( pj = Start_Index(i, list); pj < End_Index(i, list); ++pj )
{
- nbr = list->far_nbr_list[pj].nbr;
+ nbr = list->far_nbr_list.nbr[pj];
id_j = system->my_atoms[nbr].orig_id;
intrs[cnt++] = id_j;
}
diff --git a/PuReMD/src/list.c b/PuReMD/src/list.c
index 5bfba249..922f42cc 100644
--- a/PuReMD/src/list.c
+++ b/PuReMD/src/list.c
@@ -38,8 +38,8 @@ int Make_List( int n, int num_intrs, int type, int format,
l->n = n;
l->num_intrs = num_intrs;
- l->index = (int*) smalloc( n * sizeof(int), "list:index", comm );
- l->end_index = (int*) smalloc( n * sizeof(int), "list:end_index", comm );
+ l->index = smalloc( n * sizeof(int), "Make_List:index", comm );
+ l->end_index = smalloc( n * sizeof(int), "Make_List:end_index", comm );
l->type = type;
l->format = format;
@@ -51,42 +51,48 @@ int Make_List( int n, int num_intrs, int type, int format,
switch ( l->type )
{
case TYP_VOID:
- l->v = (void*) smalloc(l->num_intrs * sizeof(void*), "list:v", comm);
+ l->v = smalloc( l->num_intrs * sizeof(void*),
+ "Make_List:v", comm );
break;
case TYP_THREE_BODY:
- l->three_body_list = (three_body_interaction_data*)
- smalloc( l->num_intrs * sizeof(three_body_interaction_data),
- "list:three_bodies", comm );
+ l->three_body_list = smalloc( l->num_intrs * sizeof(three_body_interaction_data),
+ "Make_List:three_bodies", comm );
break;
case TYP_BOND:
- l->bond_list = (bond_data*)
- smalloc( l->num_intrs * sizeof(bond_data), "list:bonds", comm );
+ l->bond_list = smalloc( l->num_intrs * sizeof(bond_data),
+ "Make_List:bonds", comm );
break;
case TYP_DBO:
- l->dbo_list = (dbond_data*)
- smalloc( l->num_intrs * sizeof(dbond_data), "list:dbonds", comm );
+ l->dbo_list = smalloc( l->num_intrs * sizeof(dbond_data),
+ "Make_List:dbonds", comm );
break;
case TYP_DDELTA:
- l->dDelta_list = (dDelta_data*)
- smalloc( l->num_intrs * sizeof(dDelta_data), "list:dDeltas", comm );
+ l->dDelta_list = smalloc( l->num_intrs * sizeof(dDelta_data),
+ "Make_List:dDeltas", comm );
break;
case TYP_FAR_NEIGHBOR:
- l->far_nbr_list = (far_neighbor_data*)
- smalloc(l->num_intrs * sizeof(far_neighbor_data), "list:far_nbrs", comm);
+ l->far_nbr_list.nbr = smalloc( l->num_intrs * sizeof(int),
+ "Make_List:far_nbr_list.nbr", comm );
+ l->far_nbr_list.rel_box = smalloc( l->num_intrs * sizeof(ivec),
+ "Make_List:far_nbr_list.rel_box", comm );
+ l->far_nbr_list.d = smalloc( l->num_intrs * sizeof(real),
+ "Make_List:far_nbr_list.d", comm );
+ l->far_nbr_list.dvec = smalloc( l->num_intrs * sizeof(rvec),
+ "Make_List:far_nbr_list.dvec", comm );
break;
case TYP_HBOND:
- l->hbond_list = (hbond_data*)
- smalloc( l->num_intrs * sizeof(hbond_data), "list:hbonds", comm );
+ l->hbond_list = smalloc( l->num_intrs * sizeof(hbond_data),
+ "Make_List:hbonds", comm );
break;
default:
- fprintf( stderr, "ERROR: no %d list type defined!\n", l->type );
+ fprintf( stderr, "[ERROR]: no %d list type defined!\n", l->type );
MPI_Abort( comm, INVALID_INPUT );
}
@@ -100,31 +106,34 @@ void Delete_List( reax_list *l, MPI_Comm comm )
return;
l->allocated = 0;
- sfree( l->index, "list:index" );
- sfree( l->end_index, "list:end_index" );
+ sfree( l->index, "Delete_List:index" );
+ sfree( l->end_index, "Delete_List:end_index" );
switch (l->type)
{
case TYP_VOID:
- sfree( l->v, "list:v" );
+ sfree( l->v, "Delete_List:v" );
break;
case TYP_HBOND:
- sfree( l->hbond_list, "list:hbonds" );
+ sfree( l->hbond_list, "Delete_List:hbonds" );
break;
case TYP_FAR_NEIGHBOR:
- sfree( l->far_nbr_list, "list:far_nbrs" );
+ sfree( l->far_nbr_list.nbr, "Delete_List:far_nbr_list.nbr" );
+ sfree( l->far_nbr_list.rel_box, "Delete_List:far_nbr_list.rel_box" );
+ sfree( l->far_nbr_list.d, "Delete_List:far_nbr_list.d" );
+ sfree( l->far_nbr_list.dvec, "Delete_List:far_nbr_list.dvec" );
break;
case TYP_BOND:
- sfree( l->bond_list, "list:bonds" );
+ sfree( l->bond_list, "Delete_List:bonds" );
break;
case TYP_DBO:
- sfree( l->dbo_list, "list:dbos" );
+ sfree( l->dbo_list, "Delete_List:dbos" );
break;
case TYP_DDELTA:
- sfree( l->dDelta_list, "list:dDeltas" );
+ sfree( l->dDelta_list, "Delete_List:dDeltas" );
break;
case TYP_THREE_BODY:
- sfree( l->three_body_list, "list:three_bodies" );
+ sfree( l->three_body_list, "Delete_List:three_bodies" );
break;
default:
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index e59269ab..a0701698 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -74,7 +74,6 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
grid *g;
grid_cell *gci, *gcj;
reax_list *far_nbrs;
- far_neighbor_data *nbr_data;
reax_atom *atom1, *atom2;
#if defined(LOG_PERFORMANCE)
@@ -140,12 +139,11 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
d = rvec_Norm_Sqr( dvec );
if ( d <= cutoff )
{
- nbr_data = &(far_nbrs->far_nbr_list[num_far]);
- nbr_data->nbr = m;
- nbr_data->d = sqrt(d);
- rvec_Copy( nbr_data->dvec, dvec );
- ivec_ScaledSum( nbr_data->rel_box,
- 1, gcj->rel_box, -1, gci->rel_box );
+ far_nbrs->far_nbr_list.nbr[num_far] = m;
+ far_nbrs->far_nbr_list.d[num_far] = sqrt(d);
+ rvec_Copy( far_nbrs->far_nbr_list.dvec[num_far], dvec );
+ ivec_ScaledSum( far_nbrs->far_nbr_list.rel_box[num_far],
+ 1, gcj->rel_box, -1, gci->rel_box );
++num_far;
}
}
diff --git a/PuReMD/src/nonbonded.c b/PuReMD/src/nonbonded.c
index 4f07103f..174bfada 100644
--- a/PuReMD/src/nonbonded.c
+++ b/PuReMD/src/nonbonded.c
@@ -47,7 +47,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
real e_ele, e_vdW, e_core;
rvec temp, ext_press;
two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
reax_list *far_nbrs;
// rtensor temp_rtensor, total_rtensor;
@@ -60,27 +59,25 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
for ( i = 0; i < natoms; ++i )
{
- 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 );
orig_i = system->my_atoms[i].orig_id;
//fprintf( stderr, "i:%d, start_i: %d, end_i: %d\n", i, start_i, end_i );
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &(far_nbrs->far_nbr_list[pj]);
- j = nbr_pj->nbr;
- orig_j = system->my_atoms[j].orig_id;
+ j = far_nbrs->far_nbr_list.nbr[pj];
+ orig_j = system->my_atoms[j].orig_id;
- if ( nbr_pj->d <= control->nonb_cut
+ if ( far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut
&& ((far_nbrs->format == HALF_LIST && (j < natoms || orig_i < orig_j))
|| (far_nbrs->format == FULL_LIST && orig_i < orig_j)) )
{
- r_ij = nbr_pj->d;
- twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
- [ system->my_atoms[j].type ]);
+ r_ij = far_nbrs->far_nbr_list.d[pj];
+ twbp = &system->reax_param.tbp[
+ system->my_atoms[i].type ][ system->my_atoms[j].type ];
/* Calculate Taper and its derivative */
- // Tap = nbr_pj->Tap; -- precomputed during compte_H
Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
Tap = Tap * r_ij + workspace->Tap[5];
Tap = Tap * r_ij + workspace->Tap[4];
@@ -96,12 +93,13 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
dTap = dTap * r_ij + 2 * workspace->Tap[2];
dTap += workspace->Tap[1] / r_ij;
- /*vdWaals Calculations*/
- if (system->reax_param.gp.vdw_type == 1 || system->reax_param.gp.vdw_type == 3)
+ /* vdWaals Calculations */
+ if ( system->reax_param.gp.vdw_type == 1
+ || system->reax_param.gp.vdw_type == 3 )
{
// shielding
- powr_vdW1 = pow(r_ij, p_vdW1);
- powgi_vdW1 = pow( 1.0 / twbp->gamma_w, p_vdW1);
+ powr_vdW1 = pow( r_ij, p_vdW1 );
+ powgi_vdW1 = pow( 1.0 / twbp->gamma_w, p_vdW1 );
fn13 = pow( powr_vdW1 + powgi_vdW1, p_vdW1i );
exp1 = exp( twbp->alpha * (1.0 - fn13 / twbp->r_vdW) );
@@ -110,11 +108,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_vdW = twbp->D * (exp1 - 2.0 * exp2);
data->my_en.e_vdW += Tap * e_vdW;
- dfn13 = pow( powr_vdW1 + powgi_vdW1, p_vdW1i - 1.0) *
- pow(r_ij, p_vdW1 - 2.0);
+ dfn13 = pow( powr_vdW1 + powgi_vdW1, p_vdW1i - 1.0 )
+ * pow( r_ij, p_vdW1 - 2.0 );
- CEvd = dTap * e_vdW -
- Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) * dfn13;
+ CEvd = dTap * e_vdW - Tap * twbp->D
+ * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) * dfn13;
}
else // no shielding
{
@@ -156,24 +154,30 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
if ( control->virial == 0 )
{
- rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
+ rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb),
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb),
+ far_nbrs->far_nbr_list.dvec[pj] );
}
else /* NPT, iNPT or sNPT */
{
/* for pressure coupling, terms not related to bond order
derivatives are added directly into pressure vector/tensor */
- rvec_Scale( temp, CEvd + CEclmb, nbr_pj->dvec );
+ rvec_Scale( temp, CEvd + CEclmb,
+ far_nbrs->far_nbr_list.dvec[pj] );
rvec_ScaledAdd( workspace->f[i], -1., temp );
rvec_Add( workspace->f[j], temp );
- rvec_iMultiply( ext_press, nbr_pj->rel_box, temp );
+ rvec_iMultiply( ext_press,
+ far_nbrs->far_nbr_list.rel_box[pj], temp );
rvec_Add( data->my_ext_press, ext_press );
// fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)
// force(%f %f %f) ext_press (%12.6f %12.6f %12.6f)\n",
- // i, j, nbr_pj->rel_box[0], nbr_pj->rel_box[1], nbr_pj->rel_box[2],
+ // i, j, far_nbrs->far_nbr_list.rel_box[pj][0],
+ // far_nbrs->far_nbr_list.rel_box[pj][1],
+ // far_nbrs->far_nbr_list.rel_box[pj][2],
// temp[0], temp[1], temp[2],
// data->ext_press[0], data->ext_press[1], data->ext_press[2] );
}
@@ -195,10 +199,14 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_ele, data->my_en.e_ele );
#endif
#ifdef TEST_FORCES
- rvec_ScaledAdd( workspace->f_vdw[i], -CEvd, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_vdw[j], +CEvd, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_ele[i], -CEclmb, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_ele[j], +CEclmb, nbr_pj->dvec );
+ rvec_ScaledAdd( workspace->f_vdw[i], -CEvd,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_vdw[j], +CEvd,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_ele[i], -CEclmb,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_ele[j], +CEclmb,
+ far_nbrs->far_nbr_list.dvec[pj] );
#endif
}
}
@@ -227,7 +235,6 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
real e_vdW, e_ele;
real CEvd, CEclmb;
rvec temp, ext_press;
- far_neighbor_data *nbr_pj;
reax_list *far_nbrs;
LR_lookup_table *t;
@@ -247,21 +254,19 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
for ( pj = start_i; pj < end_i; ++pj )
{
- nbr_pj = &(far_nbrs->far_nbr_list[pj]);
- j = nbr_pj->nbr;
- orig_j = system->my_atoms[j].orig_id;
+ j = far_nbrs->far_nbr_list.nbr[pj];
+ orig_j = system->my_atoms[j].orig_id;
- if ( nbr_pj->d <= control->nonb_cut
+ if ( far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut
&& ((far_nbrs->format == HALF_LIST && (j < natoms || orig_i < orig_j))
|| (far_nbrs->format == FULL_LIST && orig_i < orig_j)) )
{
- j = nbr_pj->nbr;
type_j = system->my_atoms[j].type;
- r_ij = nbr_pj->d;
- tmin = MIN( type_i, type_j );
- tmax = MAX( type_i, type_j );
- t = &( LR[tmin][tmax] );
- // table = &( LR[type_i][type_j] );
+ r_ij = far_nbrs->far_nbr_list.d[pj];
+ tmin = MIN( type_i, type_j );
+ tmax = MAX( type_i, type_j );
+ t = &LR[tmin][tmax];
+ // table = &LR[type_i][type_j];
/* Cubic Spline Interpolation */
r = (int)(r_ij * t->inv_dx);
@@ -292,19 +297,21 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
if ( control->virial == 0 )
{
- rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
+ rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb),
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb),
+ far_nbrs->far_nbr_list.dvec[pj] );
}
else // NPT, iNPT or sNPT
{
/* for pressure coupling, terms not related to bond order derivatives
are added directly into pressure vector/tensor */
- rvec_Scale( temp, CEvd + CEclmb, nbr_pj->dvec );
+ rvec_Scale( temp, CEvd + CEclmb, far_nbrs->far_nbr_list.dvec[pj] );
rvec_ScaledAdd( workspace->f[i], -1., temp );
rvec_Add( workspace->f[j], temp );
- rvec_iMultiply( ext_press, nbr_pj->rel_box, temp );
+ rvec_iMultiply( ext_press, far_nbrs->far_nbr_list.rel_box[pj], temp );
rvec_Add( data->my_ext_press, ext_press );
}
@@ -320,10 +327,14 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_ele, data->my_en.e_ele );
#endif
#ifdef TEST_FORCES
- rvec_ScaledAdd( workspace->f_vdw[i], -CEvd, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_vdw[j], +CEvd, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_ele[i], -CEclmb, nbr_pj->dvec );
- rvec_ScaledAdd( workspace->f_ele[j], +CEclmb, nbr_pj->dvec );
+ rvec_ScaledAdd( workspace->f_vdw[i], -CEvd,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_vdw[j], +CEvd,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_ele[i], -CEclmb,
+ far_nbrs->far_nbr_list.dvec[pj] );
+ rvec_ScaledAdd( workspace->f_ele[j], +CEclmb,
+ far_nbrs->far_nbr_list.dvec[pj] );
#endif
}
}
diff --git a/PuReMD/src/parallelreax.c b/PuReMD/src/parallelreax.c
index 0ff49df1..e633a3d3 100644
--- a/PuReMD/src/parallelreax.c
+++ b/PuReMD/src/parallelreax.c
@@ -160,7 +160,6 @@ int main( int argc, char* argv[] )
lists[i]->bond_list = NULL;
lists[i]->dbo_list = NULL;
lists[i]->dDelta_list = NULL;
- lists[i]->far_nbr_list = NULL;
lists[i]->hbond_list = NULL;
}
out_control = (output_controls *)
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index b99f95d3..ab8a7b50 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -897,10 +897,16 @@ typedef struct
typedef struct
{
- int nbr;
- ivec rel_box;
- real d;
- rvec dvec;
+ /* neighbor atom IDs */
+ int *nbr;
+ /* set of three integers which deterimine if the neighbor
+ * atom is a non-periodic neighbor (all zeros) or a periodic
+ * neighbor and which perioidic image this neighbor comes from */
+ ivec *rel_box;
+ /* distance to the neighboring atom */
+ real *d;
+ /* difference between positions of this atom and its neighboring atom */
+ rvec *dvec;
} far_neighbor_data;
@@ -917,9 +923,12 @@ typedef struct
typedef struct
{
+ /* neighbor atom ID */
int nbr;
+ /* ??? */
int scl;
- far_neighbor_data *ptr;
+ /* position of neighbor in far neighbor list */
+ int ptr;
} hbond_data;
@@ -1127,7 +1136,7 @@ typedef struct
bond_data *bond_list;
dbond_data *dbo_list;
dDelta_data *dDelta_list;
- far_neighbor_data *far_nbr_list;
+ far_neighbor_data far_nbr_list;
#if defined(NEUTRAL_TERRITORY)
nt_neighbor_data *nt_nbr_list;
#endif
--
GitLab