From edf4762f8efb49b2de0554b6ddffe62731b917dc Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 16 May 2018 13:37:57 -0400
Subject: [PATCH] PG-PuReMD: updates to bonds / hbonds list initialization for
no Coulombic interaction case.
---
PG-PuReMD/src/forces.c | 70 +++++++++++++++++++++++-------------------
1 file changed, 38 insertions(+), 32 deletions(-)
diff --git a/PG-PuReMD/src/forces.c b/PG-PuReMD/src/forces.c
index 491f888e..6ab553d6 100644
--- a/PG-PuReMD/src/forces.c
+++ b/PG-PuReMD/src/forces.c
@@ -79,11 +79,11 @@ void Dummy_Interaction( reax_system *system, control_params *control,
void Init_Force_Functions( control_params *control )
{
Interaction_Functions[0] = BO;
- Interaction_Functions[1] = Bonds; //Dummy_Interaction;
- Interaction_Functions[2] = Atom_Energy; //Dummy_Interaction;
- Interaction_Functions[2] = Atom_Energy; //Dummy_Interaction;
- Interaction_Functions[3] = Valence_Angles; //Dummy_Interaction;
- Interaction_Functions[4] = Torsion_Angles; //Dummy_Interaction;
+ Interaction_Functions[1] = Bonds;
+ Interaction_Functions[2] = Atom_Energy;
+ Interaction_Functions[2] = Atom_Energy;
+ Interaction_Functions[3] = Valence_Angles;
+ Interaction_Functions[4] = Torsion_Angles;
if ( control->hbond_cut > 0.0 )
{
Interaction_Functions[5] = Hydrogen_Bonds;
@@ -92,10 +92,10 @@ void Init_Force_Functions( control_params *control )
{
Interaction_Functions[5] = Dummy_Interaction;
}
- Interaction_Functions[6] = Dummy_Interaction; //empty
- Interaction_Functions[7] = Dummy_Interaction; //empty
- Interaction_Functions[8] = Dummy_Interaction; //empty
- Interaction_Functions[9] = Dummy_Interaction; //empty
+ Interaction_Functions[6] = Dummy_Interaction;
+ Interaction_Functions[7] = Dummy_Interaction;
+ Interaction_Functions[8] = Dummy_Interaction;
+ Interaction_Functions[9] = Dummy_Interaction;
}
@@ -105,13 +105,13 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
{
int i;
- /* Mark beginning of a new timestep in bonded energy files */
#if defined(TEST_ENERGY)
+ /* Mark beginning of a new timestep in bonded energy files */
Debug_Marker_Bonded( out_control, data->step );
#endif
/* Implement all force calls as function pointers */
- for( i = 0; i < NUM_INTRS; i++ )
+ for ( i = 0; i < NUM_INTRS; i++ )
{
(Interaction_Functions[i])( system, control, data, workspace, lists, out_control );
}
@@ -123,8 +123,8 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
reax_list **lists, output_controls *out_control,
mpi_datatypes *mpi_data )
{
- /* Mark beginning of a new timestep in nonbonded energy files */
#if defined(TEST_ENERGY)
+ /* Mark beginning of a new timestep in nonbonded energy files */
Debug_Marker_Nonbonded( out_control, data->step );
#endif
@@ -463,8 +463,8 @@ int 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, bond_list, control->bo_cut,
+ nbr_pj->d <= control->bond_cut
+ && BOp( workspace, bond_list, control->bo_cut,
i, btop_i, nbr_pj, sbp_i, sbp_j, twbp ) == TRUE )
{
++btop_i;
@@ -557,8 +557,12 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
type_i = atom_i->type;
start_i = Start_Index( i, far_nbr_list );
end_i = End_Index( i, far_nbr_list );
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
btop_i = End_Index( i, bond_list );
sbp_i = &system->reax_param.sbp[type_i];
+ ihb = NON_H_BONDING_ATOM;
+ ihb_top = -1;
if ( i < system->n )
{
@@ -571,13 +575,14 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
cutoff = control->bond_cut;
}
- ihb = NON_H_BONDING_ATOM;
- ihb_top = -1;
if ( local == TRUE && control->hbond_cut > 0.0 )
{
ihb = sbp_i->p_hbond;
+
if ( ihb == H_ATOM )
{
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
ihb_top = End_Index( atom_i->Hindex, hbond_list );
}
else
@@ -589,9 +594,9 @@ int Init_Forces_No_Charges( 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_nbr_list->far_nbr_list[pj] );
+ nbr_pj = &far_nbr_list->far_nbr_list[pj];
j = nbr_pj->nbr;
- atom_j = &(system->my_atoms[j]);
+ atom_j = &system->my_atoms[j];
if ( renbr == TRUE )
{
@@ -610,9 +615,10 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
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) )
+
+ if ( nbr_pj->d <= SQR( cutoff ) )
{
- nbr_pj->d = SQRT(nbr_pj->d);
+ nbr_pj->d = SQRT( nbr_pj->d );
flag = TRUE;
}
else
@@ -631,8 +637,9 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
if ( local == TRUE )
{
/* hydrogen bond lists */
- if ( control->hbond_cut > 0 && (ihb == H_ATOM || ihb == H_BONDING_ATOM) &&
- nbr_pj->d <= control->hbond_cut )
+ if ( control->hbond_cut > 0.0
+ && (ihb == H_ATOM || ihb == H_BONDING_ATOM)
+ && nbr_pj->d <= control->hbond_cut )
{
jhb = sbp_j->p_hbond;
@@ -643,6 +650,7 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
hbond_list->hbond_list[ihb_top].ptr = nbr_pj;
++ihb_top;
}
+ /* only add to list for local j */
else if ( j < system->n && ihb == H_BONDING_ATOM && jhb == H_ATOM )
{
jhb_top = End_Index( atom_j->Hindex, hbond_list );
@@ -656,9 +664,9 @@ int Init_Forces_No_Charges( 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, bond_list, control->bo_cut,
- i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) )
+ nbr_pj->d <= control->bond_cut
+ && BOp( workspace, bond_list, control->bo_cut,
+ i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) == TRUE )
{
++btop_i;
@@ -679,25 +687,23 @@ int Init_Forces_No_Charges( reax_system *system, control_params *control,
{
Set_End_Index( atom_i->Hindex, ihb_top, hbond_list );
}
+ }
- /* reallocation checks */
+ /* reallocation checks */
+ for ( i = 0; i < system->N; ++i )
+ {
if ( Num_Entries( i, bond_list ) > system->max_bonds[i] )
{
workspace->realloc.bonds = TRUE;
}
if ( ihb == H_ATOM
- && Num_Entries( atom_i->Hindex, hbond_list ) > system->max_bonds[atom_i->Hindex] )
+ && Num_Entries( atom_i->Hindex, hbond_list ) > system->max_hbonds[atom_i->Hindex] )
{
workspace->realloc.hbonds = TRUE;
}
}
-#if defined( DEBUG )
- Print_Bonds( system, bond_list, "debugbonds.out" );
- Print_Bond_List2( system, bond_list, "pbonds.out" );
-#endif
-
return (workspace->realloc.bonds == TRUE
|| workspace->realloc.hbonds == TRUE) ? FAILURE : SUCCESS;
}
--
GitLab