From 3f6e1de4336579c9947e9ea538861d73c0da4918 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 6 Sep 2017 23:32:46 -0400
Subject: [PATCH] sPuReMD: further OpenMP work on force computations and
interaction lists. Minor data structure scope changes. Fix numereous
non-critical memory leaks. Add tear-down code. Other general refactoring.
---
sPuReMD/src/allocate.c | 8 +-
sPuReMD/src/analyze.c | 36 +-
sPuReMD/src/bond_orders.c | 160 ++-
sPuReMD/src/box.c | 69 +-
sPuReMD/src/charges.c | 108 +-
sPuReMD/src/control.c | 25 +-
sPuReMD/src/ffield.c | 45 +-
sPuReMD/src/forces.c | 124 ++-
sPuReMD/src/four_body_interactions.c | 82 +-
sPuReMD/src/geo_tools.c | 37 +-
sPuReMD/src/grid.c | 15 +-
sPuReMD/src/grid.h | 4 +
sPuReMD/src/init_md.c | 403 ++++++-
sPuReMD/src/init_md.h | 7 +-
sPuReMD/src/integrate.c | 4 +-
sPuReMD/src/list.c | 73 +-
sPuReMD/src/list.h | 18 +-
sPuReMD/src/lookup.c | 26 +-
sPuReMD/src/lookup.h | 10 +-
sPuReMD/src/mytypes.h | 504 +++++----
sPuReMD/src/neighbors.c | 19 +-
sPuReMD/src/print_utils.c | 277 +++--
sPuReMD/src/reset_utils.c | 39 +-
sPuReMD/src/single_body_interactions.c | 17 +-
sPuReMD/src/system_props.c | 54 +-
sPuReMD/src/testmd.c | 16 +-
sPuReMD/src/three_body_interactions.c | 1356 +++++++++++++-----------
sPuReMD/src/tool_box.c | 15 +
sPuReMD/src/tool_box.h | 1 +
sPuReMD/src/two_body_interactions.c | 59 +-
sPuReMD/src/vector.c | 57 +-
sPuReMD/src/vector.h | 2 +
32 files changed, 2338 insertions(+), 1332 deletions(-)
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index ab075f94..6da598bd 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -56,7 +56,7 @@ int PreAllocate_Space( reax_system *system, control_params *control,
void Reallocate_Neighbor_List( list *far_nbrs, int n, int num_intrs )
{
- Delete_List( far_nbrs );
+ Delete_List( TYP_FAR_NEIGHBOR, far_nbrs );
if (!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs ))
{
@@ -188,7 +188,7 @@ int Reallocate_HBonds_List( int n, int num_h, int *h_index, list *hbonds )
}
}
- Delete_List( hbonds );
+ Delete_List( TYP_HBOND, hbonds );
Allocate_HBond_List( n, num_h, h_index, hb_top, hbonds );
@@ -249,7 +249,7 @@ int Reallocate_Bonds_List( int n, list *bonds, int *num_bonds, int *est_3body )
bond_top[i] = MAX( Num_Entries( i, bonds ) * 2, MIN_BONDS );
}
- Delete_List( bonds );
+ Delete_List( TYP_BOND, bonds );
Allocate_Bond_List( n, bond_top, bonds );
*num_bonds = bond_top[n - 1];
@@ -306,7 +306,7 @@ void Reallocate( reax_system *system, static_storage *workspace, list **lists,
if ( realloc->num_3body > 0 )
{
- Delete_List( (*lists) + THREE_BODIES );
+ Delete_List( TYP_THREE_BODY, (*lists) + THREE_BODIES );
if ( num_bonds == -1 )
num_bonds = ((*lists) + BONDS)->num_intrs;
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index b9f5c33a..e78f533d 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -20,12 +20,42 @@
----------------------------------------------------------------------*/
#include "analyze.h"
+
#include "box.h"
#include "list.h"
#include "vector.h"
+
#define MAX_FRAGMENT_TYPES 100
+
+enum atoms
+{
+ C_ATOM = 0,
+ H_ATOM = 1,
+ O_ATOM = 2,
+ N_ATOM = 3,
+ S_ATOM = 4,
+ SI_ATOM = 5,
+ GE_ATOM = 6,
+ X_ATOM = 7,
+};
+
+enum molecule_type
+{
+ UNKNOWN = 0,
+ WATER = 1,
+};
+
+
+typedef struct
+{
+ int atom_count;
+ int atom_list[MAX_MOLECULE_SIZE];
+ int mtypes[MAX_ATOM_TYPES];
+} molecule;
+
+
// copy bond list into old bond list
void Copy_Bond_List( reax_system *system, control_params *control,
list **lists )
@@ -772,9 +802,9 @@ void Calculate_Drift( reax_system *system, control_params *control,
Distance_on_T3_Gen( workspace->x_old[i], system->atoms[i].x,
&(system->box), driftvec );
- if ( fabs( driftvec[0] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
- fabs( driftvec[1] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
- fabs( driftvec[2] ) >= system->box.box_norms[0] / 2.0 - 2.0 )
+ if ( FABS( driftvec[0] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
+ FABS( driftvec[1] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
+ FABS( driftvec[2] ) >= system->box.box_norms[0] / 2.0 - 2.0 )
{
/* the atom has moved almost half the box size.
exclude it from further drift computations as it might have an
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index dd902994..686cf3b0 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -29,7 +29,7 @@
inline real Cf45( real p1, real p2 )
{
return -EXP(-p2 / 2) /
- ( SQR( EXP(-p1 / 2) + EXP(p1 / 2) ) * (EXP(-p2 / 2) + EXP(p2 / 2)) );
+ ( SQR( EXP(-p1 / 2) + EXP(p1 / 2) ) * (EXP(-p2 / 2) + EXP(p2 / 2)) );
}
@@ -188,7 +188,6 @@ void Add_dDelta_to_Forces( reax_system *system, list **lists, int i, real C )
void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
int *top )
{
- /* Initializations */
int j, k, l, start_i, end_i, end_j;
rvec dDeltap_self, dBOp;
list *bonds, *dBOs;
@@ -196,34 +195,20 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
bond_order_data *bo_ij;
dbond_data *top_dbo;
+ /* Initializations */
bonds = (*lists) + BONDS;
dBOs = (*lists) + DBO;
-
j = bonds->select.bond_list[pj].nbr;
bo_ij = &(bonds->select.bond_list[pj].bo_data);
-
- /*rvec due_j[1000], due_i[1000];
- rvec due_j_pi[1000], due_i_pi[1000];
-
- memset(due_j, 0, sizeof(rvec)*1000 );
- memset(due_i, 0, sizeof(rvec)*1000 );
- memset(due_j_pi, 0, sizeof(rvec)*1000 );
- memset(due_i_pi, 0, sizeof(rvec)*1000 );*/
-
- //fprintf( stderr,"dbo %d-%d\n",workspace->orig_id[i],workspace->orig_id[j] );
-
- start_i = Start_Index(i, bonds);
- end_i = End_Index(i, bonds);
-
- l = Start_Index(j, bonds);
- end_j = End_Index(j, bonds);
-
+ 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) ]);
for ( k = start_i; k < end_i; ++k )
{
nbr_k = &(bonds->select.bond_list[k]);
- //fprintf( stderr, "\tnbr_k = %d\n", workspace->orig_id[nbr_k->nbr] );
for ( ; l < end_j && bonds->select.bond_list[l].nbr < nbr_k->nbr; ++l )
{
@@ -232,12 +217,10 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
nbr_l = &(bonds->select.bond_list[l]);
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
- //fprintf( stderr,"\t\tnbr_l = %d\n",workspace->orig_id[nbr_l->nbr] );
rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, dBOp ); // dBO, 3rd
rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, dBOp ); // dBOpi, 4th
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, dBOp );// dBOpipi, 4th
- //rvec_ScaledAdd(due_j[top_dbo->wrt],-bo_ij->BO*bo_ij->A2_ji, dBOp);
if ( nbr_l->nbr == i )
{
@@ -248,24 +231,19 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
rvec_ScaledAdd( top_dbo->dBO, bo_ij->C2dbo, dDeltap_self ); //2nd
/* dBOpi */
- rvec_ScaledAdd(top_dbo->dBOpi, bo_ij->C1dbopi, bo_ij->dln_BOp_pi); //1
- rvec_ScaledAdd(top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp); //2nd
- rvec_ScaledAdd(top_dbo->dBOpi, bo_ij->C3dbopi, dDeltap_self); //3rd
+ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C1dbopi, bo_ij->dln_BOp_pi ); //1
+ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp ); //2nd
+ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C3dbopi, dDeltap_self ); //3rd
/* dBOpp, 1st */
- rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2);
- rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp); //2nd
- rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C3dbopi2, dDeltap_self); //3rd
-
- /* do the adjustments on i */
- //rvec_ScaledAdd( due_i[i],
- //bo_ij->A0_ij + bo_ij->BO * bo_ij->A1_ij, bo_ij->dBOp );//1st,dBO
- //rvec_ScaledAdd( due_i[i], bo_ij->BO * bo_ij->A2_ij,
- //dDeltap_self ); //2nd, dBO
+ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 );
+ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp ); //2nd
+ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C3dbopi2, dDeltap_self ); //3rd
+
}
- //rvec_Add( workspace->dDelta[nbr_l->nbr], top_dbo->dBO );
- ++(*top), ++top_dbo;
+ ++(*top);
+ ++top_dbo;
}
/* Now we are processing neighbor k of i. */
@@ -275,11 +253,6 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
rvec_Scale( top_dbo->dBO, -bo_ij->C2dbo, dBOp ); //dBO-2
rvec_Scale( top_dbo->dBOpi, -bo_ij->C3dbopi, dBOp ); //dBOpi-3
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C3dbopi2, dBOp );//dBOpp-3
- //rvec_ScaledAdd(due_i[top_dbo->wrt],-bo_ij->BO*bo_ij->A2_ij,dBOp);//dBO-2
-
- // fprintf( stderr, "\tnbr_k = %d, nbr_l = %d, l = %d, end_j = %d\n",
- // workspace->orig_id[nbr_k->nbr],
- // workspace->orig_id[bonds->select.bond_list[l].nbr], l, end_j );
if ( l < end_j && bonds->select.bond_list[l].nbr == nbr_k->nbr )
{
@@ -291,9 +264,6 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C4dbopi, dBOp ); //dBOpi,4th
rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C4dbopi2, dBOp );//dBOpp.4th
++l;
-
- //rvec_ScaledAdd( due_j[top_dbo->wrt], -bo_ij->BO * bo_ij->A2_ji,
- //nbr_l->bo_data.dBOp ); //3rd, dBO
}
else if ( k == pj )
{
@@ -304,22 +274,16 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
rvec_ScaledAdd( top_dbo->dBO, bo_ij->C3dbo, dDeltap_self );// 3rd, dBO
/* dBOpi, 1st */
- rvec_ScaledAdd(top_dbo->dBOpi, -bo_ij->C1dbopi, bo_ij->dln_BOp_pi);
- rvec_ScaledAdd(top_dbo->dBOpi, -bo_ij->C2dbopi, bo_ij->dBOp); //2nd
+ rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C1dbopi, bo_ij->dln_BOp_pi );
+ rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C2dbopi, bo_ij->dBOp ); //2nd
rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C4dbopi, dDeltap_self ); //4th
/* dBOpi2, 1st */
- rvec_ScaledAdd(top_dbo->dBOpi2, -bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 );
- rvec_ScaledAdd(top_dbo->dBOpi2, -bo_ij->C2dbopi2, bo_ij->dBOp ); //2nd
- rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C4dbopi2, dDeltap_self ); //4th
-
- //rvec_ScaledAdd( due_j[j], -(bo_ij->A0_ij + bo_ij->BO*bo_ij->A1_ij),
- //bo_ij->dBOp ); //1st, dBO
- //rvec_ScaledAdd( due_j[j], bo_ij->BO * bo_ij->A2_ji,
- //workspace->dDeltap_self[j] ); //3rd, dBO
+ rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 );
+ rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C2dbopi2, bo_ij->dBOp ); //2nd
+ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C4dbopi2, dDeltap_self ); //4th
}
- // rvec_Add( workspace->dDelta[nbr_k->nbr], top_dbo->dBO );
++(*top), ++top_dbo;
}
@@ -330,15 +294,11 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
nbr_l = &(bonds->select.bond_list[l]);
top_dbo->wrt = nbr_l->nbr;
rvec_Copy( dBOp, nbr_l->bo_data.dBOp );
- //fprintf( stderr,"\tl=%d, nbr_l=%d\n",l,workspace->orig_id[nbr_l->nbr] );
rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, dBOp ); //3rd, dBO
rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, dBOp ); //4th, dBOpi
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, dBOp );//4th, dBOpp
- // rvec_ScaledAdd( due_j[top_dbo->wrt], -bo_ij->BO * bo_ij->A2_ji,
- // nbr_l->bo_data.dBOp );
-
if ( nbr_l->nbr == i )
{
/* do the adjustments on i */
@@ -357,23 +317,10 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2);
rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp ); //2nd
rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C3dbopi2, dDeltap_self );//3rd
-
- //rvec_ScaledAdd( due_i[i], bo_ij->A0_ij + bo_ij->BO * bo_ij->A1_ij,
- //bo_ij->dBOp ); /*1st, dBO*/
- //rvec_ScaledAdd( due_i[i], bo_ij->BO * bo_ij->A2_ij,
- //dDeltap_self ); /*2nd, dBO*/
}
- // rvec_Add( workspace->dDelta[nbr_l->nbr], top_dbo->dBO );
++(*top), ++top_dbo;
}
-
- /*for( k = 0; k < 21; ++k ){
- fprintf( stderr, "%d %d %d, due_i:[%g %g %g]\n",
- i+1, j+1, k+1, due_i[k][0], due_i[k][1], due_i[k][2] );
- fprintf( stderr, "%d %d %d, due_j:[%g %g %g]\n",
- i+1, j+1, k+1, due_j[k][0], due_j[k][1], due_j[k][2] );
- }*/
}
#endif
@@ -447,7 +394,9 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
rvec_Add( *f_k, temp );
/* pressure */
rvec_iMultiply( ext_press, nbr_k->rel_box, temp );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Add_dBond_to_Forces_NPT_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -483,7 +432,6 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
rvec_Add( *f_i, temp );
/* ext pressure due to i dropped, counting force on j only will be enough */
-
/****************************************************************************
* forces and pressure related to atom j *
* first neighbors of atom j *
@@ -510,7 +458,9 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
{
ivec_Sum( rel_box, nbr_k->rel_box, nbr_j->rel_box );//k's rel_box wrt i
rvec_iMultiply( ext_press, rel_box, temp );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Add_dBond_to_Forces_NPT_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -547,7 +497,9 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
rvec_Add( *f_j, temp );
/* pressure */
rvec_iMultiply( ext_press, nbr_j->rel_box, temp );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Add_dBond_to_Forces_NPT_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -773,7 +725,9 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
p_boc2 = system->reaxprm.gp.l[1];
bonds = (*lists) + BONDS;
+#ifdef _OPENMP
#pragma omp parallel default(shared)
+#endif
{
int i, j, pj, type_i, type_j;
int start_i, end_i;
@@ -791,15 +745,22 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
single_body_parameters *sbp_i, *sbp_j;
#if defined(TEST_FORCES)
int k, pk, start_j, end_j;
- int top_dbo = 0, top_dDelta = 0;
+ int top_dbo, top_dDelta;
dbond_data *pdbo;
dDelta_data *ptop_dDelta;
- list *dDeltas = (*lists) + DDELTA;
- list *dBOs = (*lists) + DBO;
+ list *dDeltas;
+ list *dBOs;
+
+ top_dbo = 0;
+ top_dDelta = 0;
+ dDeltas = (*lists) + DDELTA;
+ dBOs = (*lists) + DBO;
#endif
/* Calculate Deltaprime, Deltaprime_boc values */
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
@@ -807,14 +768,18 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
workspace->Deltap[i] = workspace->total_bond_order[i] - sbp_i->valency;
workspace->Deltap_boc[i] =
workspace->total_bond_order[i] - sbp_i->valency_val;
- workspace->total_bond_order[i] = 0;
+ workspace->total_bond_order[i] = 0.0;
}
/* wait until initialization complete */
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* Corrected Bond Order calculations */
+#ifdef _OPENMP
#pragma omp for schedule(guided)
+#endif
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
@@ -841,6 +806,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
workspace->reverse_map[i], workspace->reverse_map[j],
twbp->ovc, twbp->v13cor, bo_ij->BO ); */
#endif
+
if ( twbp->ovc < 0.001 && twbp->v13cor < 0.001 )
{
/* There is no correction to bond orders nor to derivatives of
@@ -912,15 +878,13 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
//Cf1_ij = -Cf1A_ij * p_boc1 * exp_p1i +
// Cf1B_ij * exp_p2i / ( exp_p2i + exp_p2j );
Cf1_ij = 0.50 * ( -p_boc1 * exp_p1i / u1_ij -
- ((val_i + f2) / SQR(u1_ij)) *
- ( -p_boc1 * exp_p1i +
- exp_p2i / ( exp_p2i + exp_p2j ) ) +
- -p_boc1 * exp_p1i / u1_ji -
- ((val_j + f2) / SQR(u1_ji)) * ( -p_boc1 * exp_p1i +
- exp_p2i / ( exp_p2i + exp_p2j ) ));
+ ((val_i + f2) / SQR(u1_ij)) * ( -p_boc1 * exp_p1i +
+ exp_p2i / ( exp_p2i + exp_p2j ) ) + -p_boc1 * exp_p1i / u1_ji -
+ ((val_j + f2) / SQR(u1_ji)) * ( -p_boc1 * exp_p1i +
+ exp_p2i / ( exp_p2i + exp_p2j ) ));
Cf1_ji = -Cf1A_ij * p_boc1 * exp_p1j +
- Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j );
+ Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j );
}
else
{
@@ -960,17 +924,17 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
/* Bond Order page 10, derivative of total bond order */
A0_ij = f1 * f4f5;
A1_ij = -2 * twbp->p_boc3 * twbp->p_boc4 * bo_ij->BO *
- (Cf45_ij + Cf45_ji);
+ (Cf45_ij + Cf45_ji);
A2_ij = Cf1_ij / f1 + twbp->p_boc3 * Cf45_ij;
A2_ji = Cf1_ji / f1 + twbp->p_boc3 * Cf45_ji;
A3_ij = A2_ij + Cf1_ij / f1;
A3_ji = A2_ji + Cf1_ji / f1;
/* find corrected bond order values and their deriv coefs */
- bo_ij->BO = bo_ij->BO * A0_ij;
+ bo_ij->BO = bo_ij->BO * A0_ij;
bo_ij->BO_pi = bo_ij->BO_pi * A0_ij * f1;
bo_ij->BO_pi2 = bo_ij->BO_pi2 * A0_ij * f1;
- bo_ij->BO_s = bo_ij->BO - ( bo_ij->BO_pi + bo_ij->BO_pi2 );
+ bo_ij->BO_s = bo_ij->BO - ( bo_ij->BO_pi + bo_ij->BO_pi2 );
bo_ij->C1dbo = A0_ij + bo_ij->BO * A1_ij;
bo_ij->C2dbo = bo_ij->BO * A2_ij;
@@ -1074,9 +1038,13 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
}
/* wait for bo_ij to be updated */
+#ifdef _OPENMP
#pragma omp barrier
+#endif
+#ifdef _OPENMP
#pragma omp for schedule(guided)
+#endif
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
@@ -1124,11 +1092,15 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
}
/* need to wait for total_bond_order to be accumulated */
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* Calculate some helper variables that are used at many places
throughout force calculations */
+#ifdef _OPENMP
#pragma omp for schedule(guided)
+#endif
for ( j = 0; j < system->N; ++j )
{
type_j = system->atoms[j].type;
@@ -1139,7 +1111,7 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
workspace->Delta_boc[j] = workspace->total_bond_order[j] -
sbp_j->valency_boc;
- workspace->vlpex[j] = workspace->Delta_e[j] -
+ workspace->vlpex[j] = workspace->Delta_e[j] -
2.0 * (int)(workspace->Delta_e[j] / 2.0);
explp1 = EXP(-p_lp1 * SQR(2.0 + workspace->vlpex[j]));
workspace->nlp[j] = explp1 - (int)(workspace->Delta_e[j] / 2.0);
@@ -1148,14 +1120,14 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
/* Adri uses different dDelta_lp values than the ones in notes... */
workspace->dDelta_lp[j] = workspace->Clp[j];
//workspace->dDelta_lp[j] = workspace->Clp[j] + (0.5-workspace->Clp[j]) *
- //((fabs(workspace->Delta_e[j]/2.0 -
+ //((FABS(workspace->Delta_e[j]/2.0 -
// (int)(workspace->Delta_e[j]/2.0)) < 0.1) ? 1 : 0 );
if ( sbp_j->mass > 21.0 )
{
workspace->nlp_temp[j] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
- workspace->dDelta_lp_temp[j] = 0.;
+ workspace->dDelta_lp_temp[j] = 0.0;
}
else
{
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index 6c2fda0b..3cddaf2d 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -190,9 +190,13 @@ void Update_Box( rtensor box_tensor, simulation_box* box )
{
int i, j;
- for (i = 0; i < 3; i++)
- for (j = 0; j < 3; j++)
+ for ( i = 0; i < 3; i++ )
+ {
+ for ( j = 0; j < 3; j++ )
+ {
box->box[i][j] = box_tensor[i][j];
+ }
+ }
Make_Consistent( box );
}
@@ -240,15 +244,20 @@ void Inc_on_T3( rvec x, rvec dx, simulation_box *box )
{
tmp = x[i] + dx[i];
if ( tmp <= -box->box_norms[i] || tmp >= box->box_norms[i] )
+ {
tmp = FMOD( tmp, box->box_norms[i] );
+ }
- if ( tmp < 0 ) tmp += box->box_norms[i];
+ if ( tmp < 0 )
+ {
+ tmp += box->box_norms[i];
+ }
x[i] = tmp;
}
}
-real Sq_Distance_on_T3(rvec x1, rvec x2, simulation_box* box, rvec r)
+real Sq_Distance_on_T3( rvec x1, rvec x2, simulation_box* box, rvec r )
{
real norm = 0.0;
real d, tmp;
@@ -262,9 +271,13 @@ real Sq_Distance_on_T3(rvec x1, rvec x2, simulation_box* box, rvec r)
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
if (x2[i] > x1[i])
+ {
d -= box->box_norms[i];
+ }
else
+ {
d += box->box_norms[i];
+ }
r[i] = d;
norm += SQR(d);
@@ -323,7 +336,9 @@ real Metric_Product( rvec x1, rvec x2, simulation_box* box )
{
tmp = 0.0;
for ( j = 0; j < 3; j++ )
+ {
tmp += box->g[i][j] * x2[j];
+ }
dist += x1[i] * tmp;
}
@@ -332,7 +347,7 @@ real Metric_Product( rvec x1, rvec x2, simulation_box* box )
int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
- real cutoff, far_neighbor_data *data )
+ real cutoff, far_neighbor_data *data )
{
real norm_sqr, d, tmp;
int i;
@@ -368,9 +383,9 @@ int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
}
}
- if ( norm_sqr <= SQR(cutoff) )
+ if ( norm_sqr <= SQR( cutoff ) )
{
- data->d = sqrt(norm_sqr);
+ data->d = SQRT( norm_sqr );
return 1;
}
@@ -382,13 +397,11 @@ int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
If so, this neighborhood is added to the list of far neighbors.
Periodic boundary conditions do not apply. */
void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
- control_params *control,
- far_neighbor_data *new_nbrs, int *count )
+ control_params *control, far_neighbor_data *new_nbrs, int *count )
{
real norm_sqr;
rvec_ScaledSum( new_nbrs[0].dvec, 1.0, x2, -1.0, x1 );
-
norm_sqr = rvec_Norm_Sqr( new_nbrs[0].dvec );
if ( norm_sqr <= SQR( control->vlist_cut ) )
@@ -399,7 +412,10 @@ void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
ivec_MakeZero( new_nbrs[0].rel_box );
// rvec_MakeZero( new_nbrs[0].ext_factor );
}
- else *count = 0;
+ else
+ {
+ *count = 0;
+ }
}
@@ -408,9 +424,7 @@ void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
If the periodic distance between x1 and x2 is than vlist_cut, this
neighborhood is added to the list of far neighbors. */
void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
- control_params *control,
- far_neighbor_data *periodic_nbrs,
- int *count )
+ control_params *control, far_neighbor_data *periodic_nbrs, int *count )
{
real norm_sqr, d, tmp;
int i;
@@ -456,7 +470,10 @@ void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
*count = 1;
periodic_nbrs[0].d = SQRT( norm_sqr );
}
- else *count = 0;
+ else
+ {
+ *count = 0;
+ }
}
@@ -468,9 +485,7 @@ void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
periodic images of x2 that are two boxs away!!!
*/
void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box,
- control_params *control,
- far_neighbor_data *periodic_nbrs,
- int *count )
+ control_params *control, far_neighbor_data *periodic_nbrs, int *count )
{
int i, j, k;
int imax, jmax, kmax;
@@ -489,13 +504,16 @@ void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box
for ( i = -imax; i <= imax; ++i )
- if (fabs(d_i = ((x2[0] + i * box->box_norms[0]) - x1[0])) <= control->vlist_cut)
+ {
+ if (FABS(d_i = ((x2[0] + i * box->box_norms[0]) - x1[0])) <= control->vlist_cut)
{
for ( j = -jmax; j <= jmax; ++j )
- if (fabs(d_j = ((x2[1] + j * box->box_norms[1]) - x1[1])) <= control->vlist_cut)
+ {
+ if (FABS(d_j = ((x2[1] + j * box->box_norms[1]) - x1[1])) <= control->vlist_cut)
{
for ( k = -kmax; k <= kmax; ++k )
- if (fabs(d_k = ((x2[2] + k * box->box_norms[2]) - x1[2])) <= control->vlist_cut)
+ {
+ if (FABS(d_k = ((x2[2] + k * box->box_norms[2]) - x1[2])) <= control->vlist_cut)
{
sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
if ( sqr_norm <= SQR(control->vlist_cut) )
@@ -533,8 +551,11 @@ void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box
++(*count);
}
}
+ }
}
+ }
}
+ }
}
@@ -586,7 +607,9 @@ void Print_Box( simulation_box* box, FILE *out )
{
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
+ {
fprintf( out, "%8.3f ", box->box[i][j] );
+ }
fprintf( out, "}" );
}
fprintf( out, "}\n" );
@@ -600,7 +623,9 @@ void Print_Box( simulation_box* box, FILE *out )
{
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
+ {
fprintf( out, "%8.3f ", box->trans[i][j] );
+ }
fprintf( out, "}" );
}
fprintf( out, "}\n" );
@@ -610,7 +635,9 @@ void Print_Box( simulation_box* box, FILE *out )
{
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
+ {
fprintf( out, "%8.3f ", box->trans_inv[i][j] );
+ }
fprintf( out, "}" );
}
fprintf( out, "}\n" );
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index f81fba69..07641a69 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -107,7 +107,9 @@ static void Sort_Matrix_Rows( sparse_matrix * const A )
unsigned int i, j, si, ei;
sparse_matrix_entry *temp;
+#ifdef _OPENMP
// #pragma omp parallel default(none) private(i, j, si, ei, temp) shared(stderr)
+#endif
{
if ( ( temp = (sparse_matrix_entry *) malloc( A->n * sizeof(sparse_matrix_entry)) ) == NULL )
{
@@ -116,7 +118,9 @@ static void Sort_Matrix_Rows( sparse_matrix * const A )
}
/* sort each row of A using column indices */
+#ifdef _OPENMP
// #pragma omp for schedule(guided)
+#endif
for ( i = 0; i < A->n; ++i )
{
si = A->start[i];
@@ -153,7 +157,9 @@ static void Calculate_Droptol( const sparse_matrix * const A,
unsigned int tid;
#endif
+#ifdef _OPENMP
#pragma omp parallel default(none) private(i, j, k, val, tid), shared(droptol_local, stderr)
+#endif
{
#ifdef _OPENMP
tid = omp_get_thread_num();
@@ -185,10 +191,14 @@ static void Calculate_Droptol( const sparse_matrix * const A,
#endif
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* calculate sqaure of the norm of each row */
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < A->n; ++i )
{
for ( k = A->start[i]; k < A->start[i + 1] - 1; ++k )
@@ -214,9 +224,9 @@ static void Calculate_Droptol( const sparse_matrix * const A,
#endif
}
+#ifdef _OPENMP
#pragma omp barrier
-#ifdef _OPENMP
#pragma omp for schedule(static)
for ( i = 0; i < A->n; ++i )
{
@@ -226,13 +236,15 @@ static void Calculate_Droptol( const sparse_matrix * const A,
droptol[i] += droptol_local[k * A->n + i];
}
}
-#endif
#pragma omp barrier
+#endif
/* calculate local droptol for each row */
//fprintf( stderr, "droptol: " );
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < A->n; ++i )
{
//fprintf( stderr, "%f-->", droptol[i] );
@@ -246,19 +258,20 @@ static void Calculate_Droptol( const sparse_matrix * const A,
static int Estimate_LU_Fill( const sparse_matrix * const A, const real * const droptol )
{
- int i, j, pj;
+ int i, pj;
int fillin;
real val;
fillin = 0;
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) private(i, j, pj, val) reduction(+: fillin)
+ default(none) private(i, pj, val) reduction(+: fillin)
+#endif
for ( i = 0; i < A->n; ++i )
{
for ( pj = A->start[i]; pj < A->start[i + 1] - 1; ++pj )
{
- j = A->j[pj];
val = A->val[pj];
if ( FABS(val) > droptol[i] )
@@ -304,7 +317,7 @@ static real SuperLU_Factorize( const sparse_matrix * const A,
#ifdef _OPENMP
//TODO: set as global parameter and use
#pragma omp parallel \
- default(none) shared(nprocs)
+ default(none) shared(nprocs)
{
#pragma omp master
{
@@ -578,8 +591,10 @@ static real diag_pre_comp( const sparse_matrix * const H, real * const Hdia_inv
start = Get_Time( );
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
default(none) private(i)
+#endif
for ( i = 0; i < H->n; ++i )
{
if ( H->val[H->start[i + 1] - 1] != 0.0 )
@@ -749,6 +764,7 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
* Edmond Chow and Aftab Patel
* Fine-Grained Parallel Incomplete LU Factorization
* SIAM J. Sci. Comp. */
+#if defined(TESTING)
static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
sparse_matrix * const U_t, sparse_matrix * const U )
{
@@ -768,8 +784,10 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
exit( INSUFFICIENT_MEMORY );
}
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(D_inv, D) private(i)
+ default(none) shared(D_inv, D) private(i)
+#endif
for ( i = 0; i < A->n; ++i )
{
D_inv[i] = SQRT( A->val[A->start[i + 1] - 1] );
@@ -780,10 +798,12 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
memset( Utop, 0, sizeof(unsigned int) * (A->n + 1) );
/* to get convergence, A must have unit diagonal, so apply
- * transformation DAD, where D = D(1./sqrt(D(A))) */
+ * transformation DAD, where D = D(1./SQRT(D(A))) */
memcpy( DAD->start, A->start, sizeof(int) * (A->n + 1) );
+#ifdef _OPENMP
#pragma omp parallel for schedule(guided) \
- default(none) shared(DAD, D_inv, D) private(i, pj)
+ default(none) shared(DAD, D_inv, D) private(i, pj)
+#endif
for ( i = 0; i < A->n; ++i )
{
/* non-diagonals */
@@ -806,8 +826,10 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
for ( i = 0; i < sweeps; ++i )
{
/* for each nonzero */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, stderr) private(sum, ei_x, ei_y, k) firstprivate(x, y)
+ default(none) shared(DAD, stderr) private(sum, ei_x, ei_y, k) firstprivate(x, y)
+#endif
for ( j = 0; j < A->start[A->n]; ++j )
{
sum = ZERO;
@@ -879,8 +901,10 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
/* apply inverse transformation D^{-1}U^{T},
* since DAD \approx U^{T}U, so
* D^{-1}DADD^{-1} = A \approx D^{-1}U^{T}UD^{-1} */
+#ifdef _OPENMP
#pragma omp parallel for schedule(guided) \
- default(none) shared(D_inv) private(i, pj)
+ default(none) shared(D_inv) private(i, pj)
+#endif
for ( i = 0; i < A->n; ++i )
{
for ( pj = A->start[i]; pj < A->start[i + 1]; ++pj )
@@ -907,6 +931,7 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
return Get_Timing_Info( start );
}
+#endif
/* Fine-grained (parallel) incomplete LU factorization
@@ -936,8 +961,10 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
exit( INSUFFICIENT_MEMORY );
}
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(D, D_inv) private(i)
+ default(none) shared(D, D_inv) private(i)
+#endif
for ( i = 0; i < A->n; ++i )
{
D_inv[i] = SQRT( FABS( A->val[A->start[i + 1] - 1] ) );
@@ -946,10 +973,12 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
}
/* to get convergence, A must have unit diagonal, so apply
- * transformation DAD, where D = D(1./sqrt(abs(D(A)))) */
+ * transformation DAD, where D = D(1./SQRT(abs(D(A)))) */
memcpy( DAD->start, A->start, sizeof(int) * (A->n + 1) );
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, D) private(i, pj)
+ default(none) shared(DAD, D) private(i, pj)
+#endif
for ( i = 0; i < A->n; ++i )
{
/* non-diagonals */
@@ -974,7 +1003,9 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
memcpy( U->val, DAD->val, sizeof(real) * (DAD->start[DAD->n]) );
/* L has unit diagonal, by convention */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) default(none) private(i)
+#endif
for ( i = 0; i < A->n; ++i )
{
L->val[L->start[i + 1] - 1] = 1.0;
@@ -983,8 +1014,10 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
for ( i = 0; i < sweeps; ++i )
{
/* for each nonzero in L */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD) private(j, k, x, y, ei_x, ei_y, sum)
+ default(none) shared(DAD) private(j, k, x, y, ei_x, ei_y, sum)
+#endif
for ( j = 0; j < DAD->start[DAD->n]; ++j )
{
sum = ZERO;
@@ -1033,8 +1066,10 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
}
}
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD) private(j, k, x, y, ei_x, ei_y, sum)
+ default(none) shared(DAD) private(j, k, x, y, ei_x, ei_y, sum)
+#endif
for ( j = 0; j < DAD->start[DAD->n]; ++j )
{
sum = ZERO;
@@ -1084,8 +1119,10 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
/* apply inverse transformation:
* since DAD \approx LU, then
* D^{-1}DADD^{-1} = A \approx D^{-1}LUD^{-1} */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, D_inv) private(i, pj)
+ default(none) shared(DAD, D_inv) private(i, pj)
+#endif
for ( i = 0; i < DAD->n; ++i )
{
for ( pj = DAD->start[i]; pj < DAD->start[i + 1]; ++pj )
@@ -1146,8 +1183,10 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
exit( INSUFFICIENT_MEMORY );
}
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(D, D_inv) private(i)
+ default(none) shared(D, D_inv) private(i)
+#endif
for ( i = 0; i < A->n; ++i )
{
D_inv[i] = SQRT( FABS( A->val[A->start[i + 1] - 1] ) );
@@ -1155,10 +1194,12 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
}
/* to get convergence, A must have unit diagonal, so apply
- * transformation DAD, where D = D(1./sqrt(D(A))) */
+ * transformation DAD, where D = D(1./SQRT(D(A))) */
memcpy( DAD->start, A->start, sizeof(int) * (A->n + 1) );
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, D) private(i, pj)
+ default(none) shared(DAD, D) private(i, pj)
+#endif
for ( i = 0; i < A->n; ++i )
{
/* non-diagonals */
@@ -1183,8 +1224,10 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
memcpy( U_temp->val, DAD->val, sizeof(real) * (DAD->start[DAD->n]) );
/* L has unit diagonal, by convention */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) private(i) shared(L_temp)
+ default(none) private(i) shared(L_temp)
+#endif
for ( i = 0; i < A->n; ++i )
{
L_temp->val[L_temp->start[i + 1] - 1] = 1.0;
@@ -1193,8 +1236,10 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
for ( i = 0; i < sweeps; ++i )
{
/* for each nonzero in L */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, L_temp, U_temp) private(j, k, x, y, ei_x, ei_y, sum)
+ default(none) shared(DAD, L_temp, U_temp) private(j, k, x, y, ei_x, ei_y, sum)
+#endif
for ( j = 0; j < DAD->start[DAD->n]; ++j )
{
sum = ZERO;
@@ -1243,8 +1288,10 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
}
}
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, L_temp, U_temp) private(j, k, x, y, ei_x, ei_y, sum)
+ default(none) shared(DAD, L_temp, U_temp) private(j, k, x, y, ei_x, ei_y, sum)
+#endif
for ( j = 0; j < DAD->start[DAD->n]; ++j )
{
sum = ZERO;
@@ -1294,8 +1341,10 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
/* apply inverse transformation:
* since DAD \approx LU, then
* D^{-1}DADD^{-1} = A \approx D^{-1}LUD^{-1} */
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) shared(DAD, L_temp, U_temp, D_inv) private(i, pj)
+ default(none) shared(DAD, L_temp, U_temp, D_inv) private(i, pj)
+#endif
for ( i = 0; i < DAD->n; ++i )
{
for ( pj = DAD->start[i]; pj < DAD->start[i + 1]; ++pj )
@@ -1374,8 +1423,10 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
/* extrapolation for s & t */
//TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
- default(none) private(i, s_tmp, t_tmp)
+ default(none) private(i, s_tmp, t_tmp)
+#endif
for ( i = 0; i < system->N_cm; ++i )
{
// no extrapolation
@@ -1426,8 +1477,10 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
/* extrapolation for s */
//TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
+#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
default(none) private(i, s_tmp)
+#endif
for ( i = 0; i < system->N_cm; ++i )
{
// no extrapolation
@@ -2724,10 +2777,9 @@ static void ACKS2( reax_system * const system, control_params * const control,
void Compute_Charges( reax_system * const system, control_params * const control,
- simulation_data * const data, static_storage * const workspace,
- const list * const far_nbrs, const output_controls * const out_control )
+ simulation_data * const data, static_storage * const workspace,
+ const list * const far_nbrs, const output_controls * const out_control )
{
- int i;
#if defined(DEBUG_FOCUS)
char fname[200];
FILE * fp;
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 16083739..ad2d3d21 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -95,9 +95,9 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->P[0] = 0.000101325;
control->P[1] = 0.000101325;
control->P[2] = 0.000101325;
- control->Tau_P[0] = 500.0;
- control->Tau_P[1] = 500.0;
- control->Tau_P[2] = 500.0;
+ control->Tau_P[0] = 500.0;
+ control->Tau_P[1] = 500.0;
+ control->Tau_P[2] = 500.0;
control->Tau_PT = 500.0;
control->compressibility = 1.0;
control->press_mode = 0;
@@ -139,8 +139,10 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
/* memory allocations */
s = (char*) malloc(sizeof(char) * MAX_LINE);
tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
- for (i = 0; i < MAX_TOKENS; i++)
+ for ( i = 0; i < MAX_TOKENS; i++ )
+ {
tmp[i] = (char*) malloc(sizeof(char) * MAX_LINE);
+ }
/* read control parameters file */
while (fgets(s, MAX_LINE, fp))
@@ -319,7 +321,9 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->T_init = val;
if ( control->T_init < 0.001 )
+ {
control->T_init = 0.001;
+ }
}
else if ( strcmp(tmp[0], "temp_final") == 0 )
{
@@ -327,12 +331,15 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->T_final = val;
if ( control->T_final < 0.1 )
+ {
control->T_final = 0.1;
+ }
}
else if ( strcmp(tmp[0], "t_mass") == 0 )
{
val = atof(tmp[1]);
- control->Tau_T = val * 1.e-3; // convert t_mass from fs to ps
+ /* convert t_mass from fs to ps */
+ control->Tau_T = val * 1.e-3;
}
else if ( strcmp(tmp[0], "t_mode") == 0 )
{
@@ -549,9 +556,13 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
/* determine target T */
if ( control->T_mode == 0 )
+ {
control->T = control->T_final;
- else control->T = control->T_init;
-
+ }
+ else
+ {
+ control->T = control->T_init;
+ }
/* near neighbor and far neighbor cutoffs */
control->bo_cut = 0.01 * system->reaxprm.gp.l[29];
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 0fa59d4e..a39aacf1 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -33,11 +33,12 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
int c, i, j, k, l, m, n, o, p, cnt;
real val;
- s = (char*) malloc(sizeof(char) * MAX_LINE);
- tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
+ s = (char*) malloc( sizeof(char) * MAX_LINE );
+ tmp = (char**) malloc( sizeof(char*) * MAX_TOKENS );
for (i = 0; i < MAX_TOKENS; i++)
- tmp[i] = (char*) malloc(sizeof(char) * MAX_TOKEN_LEN);
-
+ {
+ tmp[i] = (char*) malloc( sizeof(char) * MAX_TOKEN_LEN );
+ }
/* reading first header comment */
fgets( s, MAX_LINE, fp );
@@ -56,7 +57,7 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
}
reax->gp.n_global = n;
- reax->gp.l = (real*) malloc(sizeof(real) * n);
+ reax->gp.l = (real*) malloc( sizeof(real) * n );
/* see mytypes.h for mapping between l[i] and the lambdas used in ff */
for (i = 0; i < n; i++)
@@ -397,55 +398,55 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
(reax->sbp[j].r_pi_pi + reax->sbp[i].r_pi_pi);
reax->tbp[i][j].p_boc3 =
- sqrt(reax->sbp[i].b_o_132 *
+ SQRT(reax->sbp[i].b_o_132 *
reax->sbp[j].b_o_132);
reax->tbp[j][i].p_boc3 =
- sqrt(reax->sbp[j].b_o_132 *
+ SQRT(reax->sbp[j].b_o_132 *
reax->sbp[i].b_o_132);
reax->tbp[i][j].p_boc4 =
- sqrt(reax->sbp[i].b_o_131 *
+ SQRT(reax->sbp[i].b_o_131 *
reax->sbp[j].b_o_131);
reax->tbp[j][i].p_boc4 =
- sqrt(reax->sbp[j].b_o_131 *
+ SQRT(reax->sbp[j].b_o_131 *
reax->sbp[i].b_o_131);
reax->tbp[i][j].p_boc5 =
- sqrt(reax->sbp[i].b_o_133 *
+ SQRT(reax->sbp[i].b_o_133 *
reax->sbp[j].b_o_133);
reax->tbp[j][i].p_boc5 =
- sqrt(reax->sbp[j].b_o_133 *
+ SQRT(reax->sbp[j].b_o_133 *
reax->sbp[i].b_o_133);
reax->tbp[i][j].D =
- sqrt(reax->sbp[i].epsilon *
+ SQRT(reax->sbp[i].epsilon *
reax->sbp[j].epsilon);
reax->tbp[j][i].D =
- sqrt(reax->sbp[j].epsilon *
+ SQRT(reax->sbp[j].epsilon *
reax->sbp[i].epsilon);
reax->tbp[i][j].alpha =
- sqrt(reax->sbp[i].alpha *
+ SQRT(reax->sbp[i].alpha *
reax->sbp[j].alpha);
reax->tbp[j][i].alpha =
- sqrt(reax->sbp[j].alpha *
+ SQRT(reax->sbp[j].alpha *
reax->sbp[i].alpha);
reax->tbp[i][j].r_vdW =
- 2.0 * sqrt(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);
+ 2.0 * SQRT(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);
reax->tbp[j][i].r_vdW =
- 2.0 * sqrt(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);
+ 2.0 * SQRT(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);
reax->tbp[i][j].gamma_w =
- sqrt(reax->sbp[i].gamma_w *
+ SQRT(reax->sbp[i].gamma_w *
reax->sbp[j].gamma_w);
reax->tbp[j][i].gamma_w =
- sqrt(reax->sbp[j].gamma_w *
+ SQRT(reax->sbp[j].gamma_w *
reax->sbp[i].gamma_w);
reax->tbp[i][j].gamma =
@@ -717,7 +718,9 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
/* deallocate helper storage */
for ( i = 0; i < MAX_TOKENS; i++ )
+ {
free( tmp[i] );
+ }
free( tmp );
free( s );
@@ -727,7 +730,9 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
for ( j = 0; j < reax->num_atom_types; j++ )
{
for ( k = 0; k < reax->num_atom_types; k++ )
+ {
free( tor_flag[i][j][k] );
+ }
free( tor_flag[i][j] );
}
@@ -735,6 +740,8 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
free( tor_flag[i] );
}
+ free( tor_flag );
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "force field read\n" );
#endif
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 0d046f37..b5c1a830 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -34,6 +34,10 @@
#include "vector.h"
+/* File scope variables */
+static interaction_function Interaction_Functions[NO_OF_INTERACTIONS];
+
+
typedef enum
{
DIAGONAL = 0,
@@ -179,14 +183,18 @@ void Compute_Total_Force( reax_system *system, control_params *control,
bonds = (*lists) + BONDS;
+#ifdef _OPENMP
#pragma omp parallel default(shared)
+#endif
{
int pj;
#ifdef _OPENMP
int j;
#endif
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < system->N; ++i )
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
@@ -402,7 +410,9 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
control_params *control, int i, int j,
real r_ij, MATRIX_ENTRY_POSITION pos )
{
- real Tap, gamij, dr3gamij_1, dr3gamij_3, ret = 0.0;
+ real Tap, gamij, dr3gamij_1, dr3gamij_3, ret;
+
+ ret = 0.0;
switch ( control->charge_method )
{
@@ -701,7 +711,6 @@ void Init_Forces( reax_system *system, control_params *control,
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
- real p_boc1, p_boc2;
sparse_matrix *H, *H_sp;
list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
@@ -714,7 +723,6 @@ void Init_Forces( reax_system *system, control_params *control,
far_nbrs = *lists + FAR_NBRS;
bonds = *lists + BONDS;
hbonds = *lists + HBONDS;
-
H = workspace->H;
H_sp = workspace->H_sp;
Htop = 0;
@@ -722,8 +730,6 @@ void Init_Forces( reax_system *system, control_params *control,
num_bonds = 0;
num_hbonds = 0;
btop_i = btop_j = 0;
- p_boc1 = system->reaxprm.gp.l[0];
- p_boc2 = system->reaxprm.gp.l[1];
for ( i = 0; i < system->N; ++i )
{
@@ -825,28 +831,40 @@ void Init_Forces( reax_system *system, control_params *control,
/* uncorrected bond orders */
if ( far_nbrs->select.far_nbr_list[pj].d <= control->nbr_cut )
{
- r2 = SQR(r_ij);
+ r2 = SQR( r_ij );
if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0)
{
C12 = twbp->p_bo1 * POW( r_ij / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + control->bo_cut) * EXP( C12 );
}
- else BO_s = C12 = 0.0;
+ else
+ {
+ BO_s = 0.0;
+ C12 = 0.0;
+ }
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0)
{
C34 = twbp->p_bo3 * POW( r_ij / twbp->r_p, twbp->p_bo4 );
BO_pi = EXP( C34 );
}
- else BO_pi = C34 = 0.0;
+ else
+ {
+ BO_pi = 0.0;
+ C34 = 0.0;
+ }
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0)
{
C56 = twbp->p_bo5 * POW( r_ij / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = EXP( C56 );
}
- else BO_pi2 = C56 = 0.0;
+ else
+ {
+ BO_pi2 = 0.0;
+ C56 = 0.0;
+ }
/* Initially BO values are the uncorrected ones, page 1 */
BO = BO_s + BO_pi + BO_pi2;
@@ -888,13 +906,13 @@ void Init_Forces( reax_system *system, control_params *control,
/* Only dln_BOp_xx wrt. dr_i is stored here, note that
dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
- 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);
- 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 );
+ 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 );
+ 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 */
@@ -913,8 +931,12 @@ void Init_Forces( reax_system *system, control_params *control,
bo_ji->BO -= control->bo_cut;
workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
- bo_ij->Cdbo = bo_ij->Cdbopi = bo_ij->Cdbopi2 = 0.0;
- bo_ji->Cdbo = bo_ji->Cdbopi = bo_ji->Cdbopi2 = 0.0;
+ bo_ij->Cdbo = 0.0;
+ bo_ij->Cdbopi = 0.0;
+ bo_ij->Cdbopi2 = 0.0;
+ bo_ji->Cdbo = 0.0;
+ bo_ji->Cdbopi = 0.0;
+ bo_ji->Cdbopi2 = 0.0;
/*fprintf( stderr, "%d %d %g %g %g\n",
i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/
@@ -953,7 +975,7 @@ void Init_Forces( reax_system *system, control_params *control,
/* diagonal entry */
H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j,
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, i,
r_ij, DIAGONAL );
++Htop;
@@ -966,8 +988,6 @@ void Init_Forces( reax_system *system, control_params *control,
{
Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
}
- //fprintf( stderr, "%d bonds start: %d, end: %d\n",
- // i, Start_Index( i, bonds ), End_Index( i, bonds ) );
}
Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
@@ -977,9 +997,6 @@ void Init_Forces( reax_system *system, control_params *control,
H->start[system->N_cm] = Htop;
H_sp->start[system->N_cm] = H_sp_top;
-// printf("Htop = %d\n", Htop);
-// printf("H_sp_top = %d\n", H_sp_top);
-
/* validate lists - decide if reallocation is required! */
Validate_Lists( workspace, lists,
data->step, system->N, H->m, Htop, num_bonds, num_hbonds );
@@ -1006,7 +1023,6 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
- real p_boc1, p_boc2;
sparse_matrix *H, *H_sp;
list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
@@ -1027,8 +1043,6 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
num_bonds = 0;
num_hbonds = 0;
btop_i = btop_j = 0;
- p_boc1 = system->reaxprm.gp.l[0];
- p_boc2 = system->reaxprm.gp.l[1];
for ( i = 0; i < system->N; ++i )
{
@@ -1075,7 +1089,7 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
{
flag_sp = 1;
}
- nbr_pj->d = sqrt(nbr_pj->d);
+ nbr_pj->d = SQRT( nbr_pj->d );
flag = 1;
}
@@ -1127,28 +1141,40 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
/* uncorrected bond orders */
if ( far_nbrs->select.far_nbr_list[pj].d <= control->nbr_cut )
{
- r2 = SQR(r_ij);
+ r2 = SQR( r_ij );
if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0)
{
C12 = twbp->p_bo1 * POW( r_ij / twbp->r_s, twbp->p_bo2 );
BO_s = (1.0 + control->bo_cut) * EXP( C12 );
}
- else BO_s = C12 = 0.0;
+ else
+ {
+ BO_s = 0.0;
+ C12 = 0.0;
+ }
if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0)
{
C34 = twbp->p_bo3 * POW( r_ij / twbp->r_p, twbp->p_bo4 );
BO_pi = EXP( C34 );
}
- else BO_pi = C34 = 0.0;
+ else
+ {
+ BO_pi = 0.0;
+ C34 = 0.0;
+ }
if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0)
{
C56 = twbp->p_bo5 * POW( r_ij / twbp->r_pp, twbp->p_bo6 );
BO_pi2 = EXP( C56 );
}
- else BO_pi2 = C56 = 0.0;
+ else
+ {
+ BO_pi2 = 0.0;
+ C56 = 0.0;
+ }
/* Initially BO values are the uncorrected ones, page 1 */
BO = BO_s + BO_pi + BO_pi2;
@@ -1191,13 +1217,13 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
/* Only dln_BOp_xx wrt. dr_i is stored here, note that
dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
- 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);
- 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 );
+ 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 );
+ 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 */
@@ -1216,8 +1242,12 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
bo_ji->BO -= control->bo_cut;
workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
- bo_ij->Cdbo = bo_ij->Cdbopi = bo_ij->Cdbopi2 = 0.0;
- bo_ji->Cdbo = bo_ji->Cdbopi = bo_ji->Cdbopi2 = 0.0;
+ bo_ij->Cdbo = 0.0;
+ bo_ij->Cdbopi = 0.0;
+ bo_ij->Cdbopi2 = 0.0;
+ bo_ji->Cdbo = 0.0;
+ bo_ji->Cdbopi = 0.0;
+ bo_ji->Cdbopi2 = 0.0;
Set_End_Index( j, btop_j + 1, bonds );
}
@@ -1237,15 +1267,17 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
Set_End_Index( i, btop_i, bonds );
if ( ihb == 1 )
+ {
Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
+ }
}
// mark the end of j list
H->start[i] = Htop;
H_sp->start[i] = H_sp_top;
/* validate lists - decide if reallocation is required! */
- Validate_Lists( workspace, lists,
- data->step, system->N, H->m, Htop, num_bonds, num_hbonds );
+ Validate_Lists( workspace, lists, data->step, system->N,
+ H->m, Htop, num_bonds, num_hbonds );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "step%d: Htop = %d, num_bonds = %d, num_hbonds = %d\n",
@@ -1264,10 +1296,9 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
int start_i, end_i;
int type_i, type_j;
int ihb, jhb;
- real r_ij, r2;
+ real r_ij;
real C12, C34, C56;
real BO, BO_s, BO_pi, BO_pi2;
- real p_boc1, p_boc2;
list *far_nbrs;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
@@ -1275,8 +1306,6 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
reax_atom *atom_i, *atom_j;
far_nbrs = *lists + FAR_NBRS;
- p_boc1 = system->reaxprm.gp.l[0];
- p_boc2 = system->reaxprm.gp.l[1];
for ( i = 0; i < system->N; ++i )
{
@@ -1315,7 +1344,6 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
if ( nbr_pj->d <= control->nbr_cut )
{
r_ij = nbr_pj->d;
- r2 = SQR(r_ij);
if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0)
{
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index aaaa701f..e0010e8c 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -169,11 +169,13 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
e_tor_total = 0.0;
e_con_total = 0.0;
+#ifdef _OPENMP
#pragma omp parallel default(shared) reduction(+: e_tor_total, e_con_total)
+#endif
{
int i, j, k, l, pi, pj, pk, pl, pij, plk;
int type_i, type_j, type_k, type_l;
- int start_j, end_j, start_k, end_k;
+ int start_j, end_j;
int start_pj, end_pj, start_pk, end_pk;
#ifdef TEST_FORCES
int num_frb_intrs = 0;
@@ -208,9 +210,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec *f_i, *f_j, *f_k, *f_l;
#ifdef _OPENMP
int tid = omp_get_thread_num( );
-#endif
#pragma omp for schedule(static)
+#endif
for ( j = 0; j < system->N; ++j )
{
type_j = system->atoms[j].type;
@@ -238,11 +240,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
/* see if there are any 3-body interactions involving j&k
where j is the central atom. Otherwise there is no point in
trying to form a 4-body interaction out of this neighborhood */
- if ( j < k && bo_jk->BO > control->thb_cut/*0*/ &&
+ if ( j < k && bo_jk->BO > control->thb_cut &&
Num_Entries(pk, thb_intrs) )
{
- start_k = Start_Index(k, bonds);
- end_k = End_Index(k, bonds);
pj = pbond_jk->sym_index; // pj points to j on k's list
/* do the same check as above: are there any 3-body interactions
@@ -265,7 +265,6 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
exp_tor34_inv = 1.0 / (1.0 + exp_tor3_DjDk + exp_tor4_DjDk);
f11_DjDk = (2.0 + exp_tor3_DjDk) * exp_tor34_inv;
-
/* pick i up from j-k interaction where j is the centre atom */
for ( pi = start_pk; pi < end_pk; ++pi )
{
@@ -320,8 +319,8 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
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/*0*/ &&
- bo_ij->BO * bo_jk->BO * bo_kl->BO > control->thb_cut/*0*/ )
+ if ( i != l && fbh->cnt && bo_kl->BO > control->thb_cut &&
+ bo_ij->BO * bo_jk->BO * bo_kl->BO > control->thb_cut )
{
#ifdef _OPENMP
f_l = &(workspace->f_local[tid * system->N + l]);
@@ -369,7 +368,7 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
/* end omega calculations */
/* torsion energy */
- exp_tor1 = EXP(fbp->p_tor1 * SQR(2. - bo_jk->BO_pi - f11_DjDk));
+ exp_tor1 = EXP( fbp->p_tor1 * SQR(2. - bo_jk->BO_pi - f11_DjDk) );
exp_tor2_kl = EXP( -p_tor2 * BOA_kl );
exp_cot2_kl = EXP( -p_cot2 * SQR(BOA_kl - 1.5) );
fn10 = (1.0 - exp_tor2_ij) * (1.0 - exp_tor2_jk) *
@@ -386,32 +385,32 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
e_tor_total += e_tor;
dfn11 = (-p_tor3 * exp_tor3_DjDk +
- (p_tor3 * exp_tor3_DjDk - p_tor4 * exp_tor4_DjDk) *
- (2. + exp_tor3_DjDk) * exp_tor34_inv) * exp_tor34_inv;
+ (p_tor3 * exp_tor3_DjDk - p_tor4 * exp_tor4_DjDk) *
+ (2. + exp_tor3_DjDk) * exp_tor34_inv) * exp_tor34_inv;
CEtors1 = sin_ijk * sin_jkl * CV;
CEtors2 = -fn10 * 2.0 * fbp->p_tor1 * fbp->V2 * exp_tor1 *
- (2.0 - bo_jk->BO_pi - f11_DjDk) * (1.0 - SQR(cos_omega)) *
- sin_ijk * sin_jkl;
+ (2.0 - bo_jk->BO_pi - f11_DjDk) * (1.0 - SQR(cos_omega)) *
+ sin_ijk * sin_jkl;
CEtors3 = CEtors2 * dfn11;
CEtors4 = CEtors1 * p_tor2 * exp_tor2_ij *
- (1.0 - exp_tor2_jk) * (1.0 - exp_tor2_kl);
+ (1.0 - exp_tor2_jk) * (1.0 - exp_tor2_kl);
CEtors5 = CEtors1 * p_tor2 * exp_tor2_jk *
- (1.0 - exp_tor2_ij) * (1.0 - exp_tor2_kl);
+ (1.0 - exp_tor2_ij) * (1.0 - exp_tor2_kl);
CEtors6 = CEtors1 * p_tor2 * exp_tor2_kl *
- (1.0 - exp_tor2_ij) * (1.0 - exp_tor2_jk);
+ (1.0 - exp_tor2_ij) * (1.0 - exp_tor2_jk);
cmn = -fn10 * CV;
CEtors7 = cmn * sin_jkl * tan_ijk_i;
CEtors8 = cmn * sin_ijk * tan_jkl_i;
CEtors9 = fn10 * sin_ijk * sin_jkl *
- (0.5 * fbp->V1 - 2.0 * fbp->V2 * exp_tor1 * cos_omega +
- 1.5 * fbp->V3 * (cos2omega + 2. * SQR(cos_omega)));
+ (0.5 * fbp->V1 - 2.0 * fbp->V2 * exp_tor1 * cos_omega +
+ 1.5 * fbp->V3 * (cos2omega + 2. * SQR(cos_omega)));
//cmn = -fn10 * CV;
//CEtors7 = cmn * sin_jkl * cos_ijk;
//CEtors8 = cmn * sin_ijk * cos_jkl;
@@ -420,7 +419,6 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
// fbp->V3 * (6*SQR(cos_omega) - 1.50));
/* end of torsion energy */
-
/* 4-body conjugation energy */
fn12 = exp_cot2_ij * exp_cot2_jk * exp_cot2_kl;
e_con = fbp->p_cot1 * fn12 *
@@ -459,11 +457,29 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
// sin_jkl, cos_jkl, tan_jkl_i );
/* forces */
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
bo_jk->Cdbopi += CEtors2;
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
workspace->CdDelta[j] += CEtors3;
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
workspace->CdDelta[k] += CEtors3;
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
bo_ij->Cdbo += (CEtors4 + CEconj1);
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
bo_jk->Cdbo += (CEtors5 + CEconj2);
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
bo_kl->Cdbo += (CEtors6 + CEconj3);
if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT )
@@ -502,7 +518,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors7 + CEconj4, p_ijk->dcos_dk );
rvec_Add( *f_i, force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -512,7 +530,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors7 + CEconj4, p_ijk->dcos_di );
rvec_Add( *f_k, force );
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -523,7 +543,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors8 + CEconj5, p_jkl->dcos_dj );
rvec_Add( *f_k, force );
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -531,7 +553,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors8 + CEconj5, p_jkl->dcos_dk );
rvec_Add( *f_l, force );
rvec_iMultiply( ext_press, rel_box_jl, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -540,7 +564,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors9 + CEconj6, dcos_omega_di );
rvec_Add( *f_i, force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -550,7 +576,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors9 + CEconj6, dcos_omega_dk );
rvec_Add( *f_k, force );
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -558,7 +586,9 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
rvec_Scale( force, CEtors9 + CEconj6, dcos_omega_dl );
rvec_Add( *f_l, force );
rvec_iMultiply( ext_press, rel_box_jl, force );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Four_Body_Interactions_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index f996af5e..442f1e0c 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -214,21 +214,21 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace )
{
- FILE *pdb;
+ FILE *pdb;
char **tmp;
- char *s, *s1;
- char descriptor[9], serial[9];
- char atom_name[9], res_name[9], res_seq[9];
- char s_x[9], s_y[9], s_z[9];
- char occupancy[9], temp_factor[9];
- char seg_id[9], element[9], charge[9];
- char alt_loc, chain_id, icode;
- char *endptr = NULL;
- int i, c, c1, pdb_serial, top;
- rvec x;
+ char *s, *s1;
+ char descriptor[9], serial[9];
+ char atom_name[9], res_name[9], res_seq[9];
+ char s_x[9], s_y[9], s_z[9];
+ char occupancy[9], temp_factor[9];
+ char seg_id[9], element[9], charge[9];
+ char alt_loc, chain_id, icode;
+ char *endptr = NULL;
+ int i, c, c1, pdb_serial, top;
+ rvec x;
reax_atom *atom;
/* open pdb file */
@@ -266,8 +266,9 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
#if defined(DEBUG_FOCUS)
fprintf( stderr, "starting to read the pdb file\n" );
#endif
+
fseek( pdb, 0, SEEK_SET );
- c = 0;
+ c = 0;
c1 = 0;
top = 0;
s[0] = 0;
@@ -424,7 +425,9 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
/* clear previous input line */
s[0] = 0;
for ( i = 0; i < c1; ++i )
+ {
tmp[i][0] = 0;
+ }
}
if ( ferror( pdb ) )
{
@@ -433,6 +436,8 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
fclose( pdb );
+ Deallocate_Tokenizer_Space( &s, &s1, &tmp );
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "finished reading the pdb file\n" );
#endif
@@ -483,7 +488,7 @@ char Write_PDB( reax_system* system, list* bonds, simulation_data *data,
(system->box.box_norms[2] * system->box.box_norms[1]) );
/*open pdb and write header*/
- sprintf(fname, "%s-%d.pdb", control->sim_name, data->step);
+ sprintf( fname, "%s-%d.pdb", control->sim_name, data->step );
pdb = fopen(fname, "w");
fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
@@ -537,8 +542,8 @@ char Write_PDB( reax_system* system, list* bonds, simulation_data *data,
}
*/
- free(buffer);
- free(line);
+ free( buffer );
+ free( line );
return SUCCESS;
}
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index c422a48a..d48644f6 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -69,7 +69,8 @@ void Allocate_Space_for_Grid( reax_system *system )
grid *g;
g = &(system->g);
- g->max_nbrs = (2 * g->spread[0] + 1) * (2 * g->spread[1] + 1) * (2 * g->spread[2] + 1) + 3;
+ g->max_nbrs = (2 * g->spread[0] + 1)
+ * (2 * g->spread[1] + 1) * (2 * g->spread[2] + 1) + 3;
/* allocate space for the new grid */
g->atoms = (int****) calloc( g->ncell[0], sizeof( int*** ));
@@ -156,6 +157,8 @@ void Deallocate_Grid_Space( grid *g )
free( g->atoms[i][j] );
free( g->top[i][j] );
+ free( g->start[i][j] );
+ free( g->end[i][j] );
free( g->mark[i][j] );
free( g->nbrs[i][j] );
free( g->nbrs_cp[i][j] );
@@ -163,6 +166,8 @@ void Deallocate_Grid_Space( grid *g )
free( g->atoms[i] );
free( g->top[i] );
+ free( g->start[i] );
+ free( g->end[i] );
free( g->mark[i] );
free( g->nbrs[i] );
free( g->nbrs_cp[i] );
@@ -170,6 +175,8 @@ void Deallocate_Grid_Space( grid *g )
free( g->atoms );
free( g->top );
+ free( g->start );
+ free( g->end );
free( g->mark );
free( g->nbrs );
free( g->nbrs_cp );
@@ -484,6 +491,12 @@ void Bin_Atoms( reax_system* system, static_storage *workspace )
}
+void Finalize_Grid( reax_system* system )
+{
+ Deallocate_Grid_Space( &( system->g ) );
+}
+
+
static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
{
dest->type = src->type;
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index cab82501..6e83b740 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -24,10 +24,13 @@
#include "mytypes.h"
+
void Setup_Grid( reax_system* );
void Update_Grid( reax_system* );
+void Finalize_Grid( reax_system* );
+
int Shift( int, int, int, grid* );
void Cluster_Atoms( reax_system *, static_storage *, control_params * );
@@ -36,4 +39,5 @@ void Bin_Atoms( reax_system*, static_storage* );
void Reset_Marks( grid*, ivec*, int );
+
#endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 6ebd65ae..26121468 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -41,11 +41,13 @@ void Generate_Initial_Velocities( reax_system *system, real T )
int i;
real scale, norm;
-
if ( T <= 0.1 )
{
- for (i = 0; i < system->N; i++)
+ for ( i = 0; i < system->N; i++ )
+ {
rvec_MakeZero( system->atoms[i].v );
+ }
+
#if defined(DEBUG)
fprintf( stderr, "no random velocities...\n" );
#endif
@@ -79,7 +81,9 @@ void Init_System( reax_system *system, control_params *control,
rvec dx;
if ( !control->restart )
+ {
Reset_Atoms( system );
+ }
Compute_Total_Mass( system, data );
Compute_Center_of_Mass( system, data, stderr );
@@ -117,7 +121,9 @@ void Init_System( reax_system *system, control_params *control,
/* Initialize velocities so that desired init T can be attained */
if ( !control->restart || (control->restart && control->random_vel) )
+ {
Generate_Initial_Velocities( system, control->T_init );
+ }
Setup_Grid( system );
}
@@ -131,7 +137,10 @@ void Init_Simulation_Data( reax_system *system, control_params *control,
Reset_Simulation_Data( data );
if ( !control->restart )
- data->step = data->prev_steps = 0;
+ {
+ data->step = 0;
+ data->prev_steps = 0;
+ }
switch ( control->ensemble )
{
@@ -328,8 +337,6 @@ void Init_Workspace( reax_system *system, control_params *control,
workspace->b_t = (real *) calloc( system->N_cm, sizeof( real ) );
workspace->b_prc = (real *) calloc( system->N_cm * 2, sizeof( real ) );
workspace->b_prm = (real *) calloc( system->N_cm * 2, sizeof( real ) );
- //TODO: check if unused
- //workspace->s_t = (real *) calloc( cm_lin_sys_size * 2, sizeof( real ) );
workspace->s = (real**) calloc( 5, sizeof( real* ) );
workspace->t = (real**) calloc( 5, sizeof( real* ) );
for ( i = 0; i < 5; ++i )
@@ -637,8 +644,8 @@ void Init_Lists( reax_system *system, control_params *control,
}
-void Init_Out_Controls(reax_system *system, control_params *control,
- static_storage *workspace, output_controls *out_control)
+void Init_Out_Controls( reax_system *system, control_params *control,
+ static_storage *workspace, output_controls *out_control )
{
char temp[1000];
@@ -867,7 +874,9 @@ void Initialize( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, list **lists,
output_controls *out_control, evolve_function *Evolve )
{
+#if defined(DEBUG)
real start, end;
+#endif
#ifdef _OPENMP
#pragma omp parallel default(shared)
@@ -898,14 +907,392 @@ void Initialize( reax_system *system, control_params *control,
if ( control->tabulate )
{
+#if defined(DEBUG)
start = Get_Time( );
+#endif
+
Make_LR_Lookup_Table( system, control );
+
+#if defined(DEBUG)
end = Get_Timing_Info( start );
- //fprintf( stderr, "Time for LR Lookup Table calculation is %f \n", end );
+ fprintf( stderr, "Time for LR Lookup Table calculation is %f \n", end );
+#endif
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "data structures have been initialized...\n" );
#endif
}
+
+
+void Finalize_System( reax_system *system, control_params *control,
+ simulation_data *data )
+{
+ int i, j, k;
+ reax_interaction *reax;
+
+ reax = &( system->reaxprm );
+
+ Finalize_Grid( system );
+
+ free( reax->gp.l );
+
+ for ( i = 0; i < reax->num_atom_types; i++ )
+ {
+ for ( j = 0; j < reax->num_atom_types; j++ )
+ {
+ for ( k = 0; k < reax->num_atom_types; k++ )
+ {
+ free( reax->fbp[i][j][k] );
+ }
+
+ free( reax->thbp[i][j] );
+ free( reax->hbp[i][j] );
+ free( reax->fbp[i][j] );
+ }
+
+ free( reax->tbp[i] );
+ free( reax->thbp[i] );
+ free( reax->hbp[i] );
+ free( reax->fbp[i] );
+ }
+
+ free( reax->sbp );
+ free( reax->tbp );
+ free( reax->thbp );
+ free( reax->hbp );
+ free( reax->fbp );
+
+ free( system->atoms );
+}
+
+
+void Finalize_Simulation_Data( reax_system *system, control_params *control,
+ simulation_data *data, output_controls *out_control )
+{
+}
+
+
+void Finalize_Workspace( reax_system *system, control_params *control,
+ static_storage *workspace )
+{
+ int i;
+
+ free( workspace->hbond_index );
+ free( workspace->total_bond_order );
+ free( workspace->Deltap );
+ free( workspace->Deltap_boc );
+ free( workspace->dDeltap_self );
+ free( workspace->Delta );
+ free( workspace->Delta_lp );
+ free( workspace->Delta_lp_temp );
+ free( workspace->dDelta_lp );
+ free( workspace->dDelta_lp_temp );
+ free( workspace->Delta_e );
+ free( workspace->Delta_boc );
+ free( workspace->nlp );
+ free( workspace->nlp_temp );
+ free( workspace->Clp );
+ free( workspace->CdDelta );
+ free( workspace->vlpex );
+
+ Deallocate_Matrix( workspace->H );
+ Deallocate_Matrix( workspace->H_sp );
+ if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+ control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+ control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+ {
+ Deallocate_Matrix( workspace->L );
+ Deallocate_Matrix( workspace->U );
+ }
+
+ for ( i = 0; i < 5; ++i )
+ {
+ free( workspace->s[i] );
+ free( workspace->t[i] );
+ }
+
+ free( workspace->Hdia_inv );
+ if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+ control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+ {
+ free( workspace->droptol );
+ }
+ //TODO: check if unused
+ //free( workspace->w );
+ //TODO: check if unused
+ free( workspace->b );
+ free( workspace->b_s );
+ free( workspace->b_t );
+ free( workspace->b_prc );
+ free( workspace->b_prm );
+ free( workspace->s );
+ free( workspace->t );
+
+ switch ( control->cm_solver_type )
+ {
+ /* GMRES storage */
+ case GMRES_S:
+ case GMRES_H_S:
+ for ( i = 0; i < control->cm_solver_restart + 1; ++i )
+ {
+ free( workspace->h[i] );
+ free( workspace->rn[i] );
+ free( workspace->v[i] );
+ }
+
+ free( workspace->y );
+ free( workspace->z );
+ free( workspace->g );
+ free( workspace->h );
+ free( workspace->hs );
+ free( workspace->hc );
+ free( workspace->rn );
+ free( workspace->v );
+
+ free( workspace->r );
+ free( workspace->d );
+ free( workspace->q );
+ free( workspace->p );
+ break;
+
+ /* CG storage */
+ case CG_S:
+ free( workspace->r );
+ free( workspace->d );
+ free( workspace->q );
+ free( workspace->p );
+ break;
+
+ case SDM_S:
+ free( workspace->r );
+ free( workspace->d );
+ free( workspace->q );
+ break;
+
+ default:
+ fprintf( stderr, "Unknown charge method linear solver type. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+
+ /* integrator storage */
+ free( workspace->a );
+ free( workspace->f_old );
+ free( workspace->v_const );
+
+#ifdef _OPENMP
+ free( workspace->f_local );
+#endif
+
+ /* storage for analysis */
+ if ( control->molec_anal || control->diffusion_coef )
+ {
+ free( workspace->mark );
+ free( workspace->old_mark );
+ }
+ else
+ {
+ free( workspace->mark );
+ }
+
+ if ( control->diffusion_coef )
+ {
+ free( workspace->x_old );
+ }
+ else
+ {
+ free( workspace->x_old );
+ }
+
+ free( workspace->orig_id );
+
+ /* space for keeping restriction info, if any */
+ if ( control->restrict_bonds )
+ {
+ for ( i = 0; i < system->N; ++i )
+ {
+ free( workspace->restricted_list[i] );
+ }
+
+ free( workspace->restricted );
+ free( workspace->restricted_list );
+ }
+
+#ifdef TEST_FORCES
+ free( workspace->dDelta );
+ free( workspace->f_ele );
+ free( workspace->f_vdw );
+ free( workspace->f_bo );
+ free( workspace->f_be );
+ free( workspace->f_lp );
+ free( workspace->f_ov );
+ free( workspace->f_un );
+ free( workspace->f_ang );
+ free( workspace->f_coa );
+ free( workspace->f_pen );
+ free( workspace->f_hb );
+ free( workspace->f_tor );
+ free( workspace->f_con );
+#endif
+}
+
+
+void Finalize_Lists( list **lists )
+{
+ Delete_List( TYP_FAR_NEIGHBOR, (*lists) + FAR_NBRS );
+ Delete_List( TYP_HBOND, (*lists) + HBONDS );
+ Delete_List( TYP_BOND, (*lists) + BONDS );
+ Delete_List( TYP_THREE_BODY, (*lists) + THREE_BODIES );
+
+#ifdef TEST_FORCES
+ Delete_List( TYP_DDELTA, (*lists) + DDELTA );
+ Delete_List( TYP_DBO, (*lists) + DBO );
+#endif
+}
+
+
+void Finalize_Out_Controls( reax_system *system, control_params *control,
+ static_storage *workspace, output_controls *out_control )
+{
+ /* close trajectory file */
+ if ( out_control->write_steps > 0 )
+ {
+ fclose( out_control->trj );
+ }
+
+ if ( out_control->energy_update_freq > 0 )
+ {
+ /* close out file */
+ fclose( out_control->out );
+
+ /* close potentials file */
+ fclose( out_control->pot );
+
+ /* close log file */
+ fclose( out_control->log );
+ }
+
+ /* close pressure file */
+ if ( control->ensemble == NPT ||
+ control->ensemble == iNPT ||
+ control->ensemble == sNPT )
+ {
+ fclose( out_control->prs );
+ }
+
+ /* close molecular analysis file */
+ if ( control->molec_anal )
+ {
+ fclose( out_control->mol );
+ }
+
+ /* close electric dipole moment analysis file */
+ if ( control->dipole_anal )
+ {
+ fclose( out_control->dpl );
+ }
+
+ /* close diffusion coef analysis file */
+ if ( control->diffusion_coef )
+ {
+ fclose( out_control->drft );
+ }
+
+
+#ifdef TEST_ENERGY
+ /* close bond energy file */
+ fclose( out_control->ebond );
+
+ /* close lone-pair energy file */
+ fclose( out_control->elp );
+
+ /* close overcoordination energy file */
+ fclose( out_control->eov );
+
+ /* close undercoordination energy file */
+ fclose( out_control->eun );
+
+ /* close angle energy file */
+ fclose( out_control->eval );
+
+ /* close penalty energy file */
+ fclose( out_control->epen );
+
+ /* close coalition energy file */
+ fclose( out_control->ecoa );
+
+ /* close hydrogen bond energy file */
+ fclose( out_control->ehb );
+
+ /* close torsion energy file */
+ fclose( out_control->etor );
+
+ /* close conjugation energy file */
+ fclose( out_control->econ );
+
+ /* close vdWaals energy file */
+ fclose( out_control->evdw );
+
+ /* close coulomb energy file */
+ fclose( out_control->ecou );
+#endif
+
+
+#ifdef TEST_FORCES
+ /* close bond orders file */
+ fclose( out_control->fbo );
+
+ /* close bond orders derivatives file */
+ fclose( out_control->fdbo );
+
+ /* close bond forces file */
+ fclose( out_control->fbond );
+
+ /* close lone-pair forces file */
+ fclose( out_control->flp );
+
+ /* close overcoordination forces file */
+ fclose( out_control->fatom );
+
+ /* close angle forces file */
+ fclose( out_control->f3body );
+
+ /* close hydrogen bond forces file */
+ fclose( out_control->fhb );
+
+ /* close torsion forces file */
+ fclose( out_control->f4body );
+
+ /* close nonbonded forces file */
+ fclose( out_control->fnonb );
+
+ /* close total force file */
+ fclose( out_control->ftot );
+
+ /* close coulomb forces file */
+ fclose( out_control->ftot2 );
+#endif
+}
+
+
+void Finalize( reax_system *system, control_params *control,
+ simulation_data *data, static_storage *workspace, list **lists,
+ output_controls *out_control )
+{
+ if ( control->tabulate )
+ {
+// Finalize_LR_Lookup_Table( system, control );
+ }
+
+ Finalize_Out_Controls( system, control, workspace, out_control );
+
+ Finalize_Lists( lists );
+
+ Finalize_Workspace( system, control, workspace );
+
+ Finalize_Simulation_Data( system, control, data, out_control );
+
+ Finalize_System( system, control, data );
+}
diff --git a/sPuReMD/src/init_md.h b/sPuReMD/src/init_md.h
index 3fc59053..3d311893 100644
--- a/sPuReMD/src/init_md.h
+++ b/sPuReMD/src/init_md.h
@@ -24,7 +24,12 @@
#include "mytypes.h"
+
void Initialize( reax_system*, control_params*, simulation_data*,
- static_storage*, list**, output_controls*, evolve_function* );
+ static_storage*, list**, output_controls*, evolve_function* );
+
+void Finalize( reax_system*, control_params*, simulation_data*,
+ static_storage*, list**, output_controls* );
+
#endif
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index ea400d4a..7f42a368 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -190,7 +190,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
itr, G_xi_new, v_xi_new, v_xi_old );
#endif
}
- while ( fabs(v_xi_new - v_xi_old ) > 1e-5 );
+ while ( FABS(v_xi_new - v_xi_old ) > 1e-5 );
therm->v_xi_old = therm->v_xi;
therm->v_xi = v_xi_new;
@@ -672,7 +672,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
"itr %d E_kin: %8.3f veps_n:%8.3f veps_o:%8.3f vxi_n:%8.3f vxi_o: %8.3f\n",
itr, E_kin, v_eps_new, v_eps_old, v_xi_new, v_xi_old );
}
- while ( fabs(v_eps_new - v_eps_old) + fabs(v_xi_new - v_xi_old) > 2e-3 );
+ while ( FABS(v_eps_new - v_eps_old) + FABS(v_xi_new - v_xi_old) > 2e-3 );
therm->v_xi_old = therm->v_xi;
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index f9044a6c..e02b22fb 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -21,7 +21,8 @@
#include "list.h"
-char Make_List(int n, int num_intrs, int type, list* l)
+
+char Make_List( int n, int num_intrs, int type, list* l )
{
char success = 1;
@@ -31,10 +32,16 @@ char Make_List(int n, int num_intrs, int type, list* l)
l->index = (int*) malloc( n * sizeof(int) );
l->end_index = (int*) malloc( n * sizeof(int) );
- if (l->index == NULL) success = 0;
- if (l->end_index == NULL) success = 0;
+ if ( l->index == NULL )
+ {
+ success = 0;
+ }
+ if ( l->end_index == NULL )
+ {
+ success = 0;
+ }
- switch (type)
+ switch ( type )
{
case TYP_VOID:
l->select.v = (void *) malloc(l->num_intrs * sizeof(void));
@@ -86,7 +93,6 @@ char Make_List(int n, int num_intrs, int type, list* l)
default:
l->select.v = (void *) malloc(l->num_intrs * sizeof(void));
if (l->select.v == NULL) success = 0;
- l->type = TYP_VOID;
break;
}
@@ -94,46 +100,66 @@ char Make_List(int n, int num_intrs, int type, list* l)
}
-void Delete_List(list* l)
+void Delete_List( int type, list* l )
{
if ( l->index != NULL )
- free(l->index);
+ {
+ free( l->index );
+ }
if ( l->end_index != NULL )
- free(l->end_index);
+ {
+ free( l->end_index );
+ }
- switch (l->type)
+ switch ( type )
{
case TYP_VOID:
if ( l->select.v != NULL )
- free(l->select.v);
+ {
+ free( l->select.v );
+ }
break;
case TYP_THREE_BODY:
if ( l->select.three_body_list != NULL )
- free(l->select.three_body_list);
+ {
+ free( l->select.three_body_list );
+ }
break;
case TYP_BOND:
if ( l->select.bond_list != NULL )
- free(l->select.bond_list);
+ {
+ free( l->select.bond_list );
+ }
break;
case TYP_DBO:
if ( l->select.dbo_list != NULL )
- free(l->select.dbo_list);
+ {
+ free( l->select.dbo_list );
+ }
break;
case TYP_DDELTA:
if ( l->select.dDelta_list != NULL )
- free(l->select.dDelta_list);
+ {
+ free( l->select.dDelta_list );
+ }
break;
case TYP_FAR_NEIGHBOR:
if ( l->select.far_nbr_list != NULL )
- free(l->select.far_nbr_list);
+ {
+ free( l->select.far_nbr_list );
+ }
break;
case TYP_NEAR_NEIGHBOR:
if ( l->select.near_nbr_list != NULL )
- free(l->select.near_nbr_list);
+ {
+ free( l->select.near_nbr_list );
+ }
break;
case TYP_HBOND:
if ( l->select.hbond_list != NULL )
- free(l->select.hbond_list);
+ {
+ free( l->select.hbond_list );
+ }
break;
default:
@@ -143,27 +169,32 @@ void Delete_List(list* l)
}
-inline int Num_Entries(int i, list* l)
+
+inline int Num_Entries( int i, list* l )
{
return l->end_index[i] - l->index[i];
}
-inline int Start_Index(int i, list *l )
+
+inline int Start_Index( int i, list *l )
{
return l->index[i];
}
+
inline int End_Index( int i, list *l )
{
return l->end_index[i];
}
-inline void Set_Start_Index(int i, int val, list *l)
+
+inline void Set_Start_Index( int i, int val, list *l )
{
l->index[i] = val;
}
-inline void Set_End_Index(int i, int val, list *l)
+
+inline void Set_End_Index( int i, int val, list *l )
{
l->end_index[i] = val;
}
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index e3ecc584..35a4b1de 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -24,14 +24,20 @@
#include "mytypes.h"
+
char Make_List( int, int, int, list* );
-void Delete_List( list* );
-int Num_Entries(int, list*);
-int Start_Index( int, list* );
-int End_Index( int, list* );
+void Delete_List( int, list* );
+
+int Num_Entries( int, list* );
+
+int Start_Index( int, list* );
+
+int End_Index( int, list* );
+
+void Set_Start_Index( int, int, list* );
+
+void Set_End_Index( int, int, list* );
-void Set_Start_Index(int, int, list*);
-void Set_End_Index(int, int, list*);
#endif
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 2ea39e3c..973ba5fa 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -20,10 +20,12 @@
----------------------------------------------------------------------*/
#include "lookup.h"
+
#include "two_body_interactions.h"
-void Make_Lookup_Table(real xmin, real xmax, int n,
- lookup_function f, lookup_table* t)
+
+void Make_Lookup_Table( real xmin, real xmax, int n,
+ lookup_function f, lookup_table* t )
{
int i;
@@ -36,7 +38,9 @@ void Make_Lookup_Table(real xmin, real xmax, int n,
t->y = (real*) malloc(n * sizeof(real));
for (i = 0; i < n; i++)
+ {
t->y[i] = f(i * t->dx + t->xmin);
+ }
// fprintf(stdout,"dx = %lf\n",t->dx);
// for(i=0; i < n; i++)
@@ -363,20 +367,20 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
LR_vdW_Coulomb( system, control, i, j, rand_dist, &y );
LR_Lookup( &(LR[i][j]), rand_dist, &y_spline );
- evdw_abserr = fabs(y.e_vdW - y_spline.e_vdW);
- evdw_relerr = fabs(evdw_abserr / y.e_vdW);
- fvdw_abserr = fabs(y.CEvd - y_spline.CEvd);
- fvdw_relerr = fabs(fvdw_abserr / y.CEvd);
- eele_abserr = fabs(y.e_ele - y_spline.e_ele);
- eele_relerr = fabs(eele_abserr / y.e_ele);
- fele_abserr = fabs(y.CEclmb - y_spline.CEclmb);
- fele_relerr = fabs(fele_abserr / y.CEclmb);
+ evdw_abserr = FABS(y.e_vdW - y_spline.e_vdW);
+ evdw_relerr = FABS(evdw_abserr / y.e_vdW);
+ fvdw_abserr = FABS(y.CEvd - y_spline.CEvd);
+ fvdw_relerr = FABS(fvdw_abserr / y.CEvd);
+ eele_abserr = FABS(y.e_ele - y_spline.e_ele);
+ eele_relerr = FABS(eele_abserr / y.e_ele);
+ fele_abserr = FABS(y.CEclmb - y_spline.CEclmb);
+ fele_relerr = FABS(fele_abserr / y.CEclmb);
if( evdw_relerr > 1e-10 || eele_relerr > 1e-10 ){
fprintf( stderr, "rand_dist = %24.15e\n", rand_dist );
fprintf( stderr, "%24.15e %24.15e %24.15e %24.15e\n",
y.H, y_spline.H,
- fabs(y.H-y_spline.H), fabs((y.H-y_spline.H)/y.H) );
+ FABS(y.H-y_spline.H), FABS((y.H-y_spline.H)/y.H) );
fprintf( stderr, "%24.15e %24.15e %24.15e %24.15e\n",
y.e_vdW, y_spline.e_vdW, evdw_abserr, evdw_relerr );
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index a0b9e516..bb1ba468 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -24,10 +24,18 @@
#include "mytypes.h"
+
+/* Function pointer definitions */
+typedef real (*lookup_function)(real);
+
+
void Make_Lookup_Table( real, real, int, lookup_function, lookup_table* );
-int Lookup_Index_Of( real, lookup_table* );
+
+int Lookup_Index_Of( real, lookup_table* );
+
real Lookup( real, lookup_table* );
void Make_LR_Lookup_Table( reax_system*, control_params* );
+
#endif
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 90fc139d..91f411e2 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -23,8 +23,8 @@
#define __MYTYPES_H_
#if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_))
-#define __CONFIG_H_
-#include "config.h"
+ #define __CONFIG_H_
+ #include "config.h"
#endif
#include "math.h"
@@ -131,27 +131,42 @@
#define LOOSE_ZONE 0.75
-typedef double real;
-typedef real rvec[3];
-typedef int ivec[3];
-typedef real rtensor[3][3];
-
/* config params */
enum ensemble
{
- NVE = 0, NVT = 1, NPT = 2, sNPT = 3, iNPT = 4, ensNR = 5, bNVT = 6,
+ NVE = 0,
+ NVT = 1,
+ NPT = 2,
+ sNPT = 3,
+ iNPT = 4,
+ ensNR = 5,
+ bNVT = 6,
};
enum interaction_list_offets
{
- FAR_NBRS = 0, NEAR_NBRS = 1, THREE_BODIES = 2, BONDS = 3, OLD_BONDS = 4,
- HBONDS = 5, DBO = 6, DDELTA = 7, LIST_N = 8,
+ FAR_NBRS = 0,
+ NEAR_NBRS = 1,
+ THREE_BODIES = 2,
+ BONDS = 3,
+ OLD_BONDS = 4,
+ HBONDS = 5,
+ DBO = 6,
+ DDELTA = 7,
+ LIST_N = 8,
};
enum interaction_type
{
- TYP_VOID = 0, TYP_THREE_BODY = 1, TYP_BOND = 2, TYP_HBOND = 3, TYP_DBO = 4,
- TYP_DDELTA = 5, TYP_FAR_NEIGHBOR = 6, TYP_NEAR_NEIGHBOR = 7, TYP_N = 8,
+ TYP_VOID = 0,
+ TYP_THREE_BODY = 1,
+ TYP_BOND = 2,
+ TYP_HBOND = 3,
+ TYP_DBO = 4,
+ TYP_DDELTA = 5,
+ TYP_FAR_NEIGHBOR = 6,
+ TYP_NEAR_NEIGHBOR = 7,
+ TYP_N = 8,
};
enum errors
@@ -168,96 +183,112 @@ enum errors
RUNTIME_ERROR = -19,
};
-enum atoms
-{
- C_ATOM = 0, H_ATOM = 1, O_ATOM = 2, N_ATOM = 3,
- S_ATOM = 4, SI_ATOM = 5, GE_ATOM = 6, X_ATOM = 7,
-};
-
-enum molecule_type
-{
- UNKNOWN = 0, WATER = 1,
-};
-
enum molecular_analysis_type
{
- NO_ANALYSIS = 0, FRAGMENTS = 1, REACTIONS = 2, NUM_ANALYSIS = 3,
+ NO_ANALYSIS = 0,
+ FRAGMENTS = 1,
+ REACTIONS = 2,
+ NUM_ANALYSIS = 3,
};
enum restart_format
{
- WRITE_ASCII = 0, WRITE_BINARY = 1, RF_N = 2,
+ WRITE_ASCII = 0,
+ WRITE_BINARY = 1,
+ RF_N = 2,
};
enum geo_formats
{
- CUSTOM = 0, PDB = 1, BGF = 2, ASCII_RESTART = 3, BINARY_RESTART = 4, GF_N = 5,
+ CUSTOM = 0,
+ PDB = 1,
+ BGF = 2,
+ ASCII_RESTART = 3,
+ BINARY_RESTART = 4,
+ GF_N = 5,
};
enum charge_method
{
- QEQ_CM = 0, EE_CM = 1, ACKS2_CM = 2,
+ QEQ_CM = 0,
+ EE_CM = 1,
+ ACKS2_CM = 2,
};
enum solver
{
- GMRES_S = 0, GMRES_H_S = 1, CG_S = 2, SDM_S = 3,
+ GMRES_S = 0,
+ GMRES_H_S = 1,
+ CG_S = 2,
+ SDM_S = 3,
};
enum pre_comp
{
- NONE_PC = 0, DIAG_PC = 1, ICHOLT_PC = 2, ILU_PAR_PC = 3, ILUT_PAR_PC = 4, ILU_SUPERLU_MT_PC = 5,
+ NONE_PC = 0,
+ DIAG_PC = 1,
+ ICHOLT_PC = 2,
+ ILU_PAR_PC = 3,
+ ILUT_PAR_PC = 4,
+ ILU_SUPERLU_MT_PC = 5,
};
enum pre_app
{
- TRI_SOLVE_PA = 0, TRI_SOLVE_LEVEL_SCHED_PA = 1, TRI_SOLVE_GC_PA = 2, JACOBI_ITER_PA = 3,
+ TRI_SOLVE_PA = 0,
+ TRI_SOLVE_LEVEL_SCHED_PA = 1,
+ TRI_SOLVE_GC_PA = 2,
+ JACOBI_ITER_PA = 3,
};
-/* Force field global params mapping */
-/*
-l[0] = p_boc1
-l[1] = p_boc2
-l[2] = p_coa2
-l[3] = N/A
-l[4] = N/A
-l[5] = N/A
-l[6] = p_ovun6
-l[7] = N/A
-l[8] = p_ovun7
-l[9] = p_ovun8
-l[10] = N/A
-l[11] = N/A
-l[12] = N/A
-l[13] = N/A
-l[14] = p_val6
-l[15] = p_lp1
-l[16] = p_val9
-l[17] = p_val10
-l[18] = N/A
-l[19] = p_pen2
-l[20] = p_pen3
-l[21] = p_pen4
-l[22] = N/A
-l[23] = p_tor2
-l[24] = p_tor3
-l[25] = p_tor4
-l[26] = N/A
-l[27] = p_cot2
-l[28] = p_vdW1
-l[29] = v_par30
-l[30] = p_coa4
-l[31] = p_ovun4
-l[32] = p_ovun3
-l[33] = p_val8
-l[34] = ACKS2 bond softness
-l[35] = N/A
-l[36] = N/A
-l[37] = version number
-l[38] = p_coa3
-*/
+typedef double real;
+typedef real rvec[3];
+typedef int ivec[3];
+typedef real rtensor[3][3];
+
+/* Force field global params mapping:
+ *
+ * l[0] = p_boc1
+ * l[1] = p_boc2
+ * l[2] = p_coa2
+ * l[3] = N/A
+ * l[4] = N/A
+ * l[5] = N/A
+ * l[6] = p_ovun6
+ * l[7] = N/A
+ * l[8] = p_ovun7
+ * l[9] = p_ovun8
+ * l[10] = N/A
+ * l[11] = N/A
+ * l[12] = N/A
+ * l[13] = N/A
+ * l[14] = p_val6
+ * l[15] = p_lp1
+ * l[16] = p_val9
+ * l[17] = p_val10
+ * l[18] = N/A
+ * l[19] = p_pen2
+ * l[20] = p_pen3
+ * l[21] = p_pen4
+ * l[22] = N/A
+ * l[23] = p_tor2
+ * l[24] = p_tor3
+ * l[25] = p_tor4
+ * l[26] = N/A
+ * l[27] = p_cot2
+ * l[28] = p_vdW1
+ * l[29] = v_par30
+ * l[30] = p_coa4
+ * l[31] = p_ovun4
+ * l[32] = p_ovun3
+ * l[33] = p_val8
+ * l[34] = ACKS2 bond softness
+ * l[35] = N/A
+ * l[36] = N/A
+ * l[37] = version number
+ * l[38] = p_coa3 */
typedef struct
{
int n_global;
@@ -266,7 +297,6 @@ typedef struct
} global_parameters;
-
typedef struct
{
/* Line one in field file */
@@ -289,8 +319,9 @@ typedef struct
real p_ovun5;
real chi;
real eta;
- int p_hbond; /* Determines whether this type of atom participates in H_bonds.
- It is 1 for donor H, 2 for acceptors (O,S,N), 0 for others*/
+ /* Determines whether this type of atom participates in H_bonds.
+ * It is 1 for donor H, 2 for acceptors (O,S,N), 0 for others*/
+ int p_hbond;
/* Line three in field file */
real r_pi_pi;
@@ -298,7 +329,8 @@ typedef struct
real b_o_131;
real b_o_132;
real b_o_133;
- real b_s_acks2; /* bond softness for ACKS2 */
+ /* bond softness for ACKS2 */
+ real b_s_acks2;
/* Line four in the field file */
real p_ovun2;
@@ -311,18 +343,29 @@ typedef struct
} single_body_parameters;
-
/* Two Body Parameters */
typedef struct
{
/* Bond Order parameters */
- real p_bo1, p_bo2, p_bo3, p_bo4, p_bo5, p_bo6;
- real r_s, r_p, r_pp; /* r_o distances in BO formula */
- real p_boc3, p_boc4, p_boc5;
+ real p_bo1;
+ real p_bo2;
+ real p_bo3;
+ real p_bo4;
+ real p_bo5;
+ real p_bo6;
+ real r_s;
+ real r_p;
+ real r_pp; /* r_o distances in BO formula */
+ real p_boc3;
+ real p_boc4;
+ real p_boc5;
/* Bond Energy parameters */
- real p_be1, p_be2;
- real De_s, De_p, De_pp;
+ real p_be1;
+ real p_be2;
+ real De_s;
+ real De_p;
+ real De_pp;
/* Over/Under coordination parameters */
real p_ovun1;
@@ -332,7 +375,9 @@ typedef struct
real alpha;
real r_vdW;
real gamma_w;
- real rcore, ecore, acore;
+ real rcore;
+ real ecore;
+ real acore;
/* electrostatic parameters */
real gamma; // note: this parameter is gamma^-3 and not gamma.
@@ -341,13 +386,15 @@ typedef struct
} two_body_parameters;
-
/* 3-body parameters */
typedef struct
{
/* valence angle */
real theta_00;
- real p_val1, p_val2, p_val4, p_val7;
+ real p_val1;
+ real p_val2;
+ real p_val4;
+ real p_val7;
/* penalty */
real p_pen1;
@@ -364,19 +411,22 @@ typedef struct
} three_body_header;
-
/* hydrogen-bond parameters */
typedef struct
{
- real r0_hb, p_hb1, p_hb2, p_hb3;
+ real r0_hb;
+ real p_hb1;
+ real p_hb2;
+ real p_hb3;
} hbond_parameters;
-
/* 4-body parameters */
typedef struct
{
- real V1, V2, V3;
+ real V1;
+ real V2;
+ real V3;
/* torsion angle */
real p_tor1;
@@ -407,14 +457,18 @@ typedef struct
typedef struct
{
- int type; /* Type of this atom */
+ /* Type of this atom */
+ int type;
+ /**/
char name[8];
-
- rvec x; // position
- rvec v; // velocity
- rvec f; // force
-
- real q; /* Charge on the atom */
+ /* position */
+ rvec x;
+ /* velocity */
+ rvec v;
+ /* force */
+ rvec f;
+ /* Charge on the atom */
+ real q;
} reax_atom;
@@ -425,7 +479,6 @@ typedef struct
rvec box_norms;
rvec side_prop;
rvec nbr_box_press[27];
- // rvec lower_end;
rtensor box, box_inv, old_box;
rtensor trans, trans_inv;
@@ -435,9 +488,9 @@ typedef struct
typedef struct
{
- int max_atoms;
- int max_nbrs;
- int total;
+ int max_atoms;
+ int max_nbrs;
+ int total;
real cell_size;
ivec spread;
@@ -446,10 +499,10 @@ typedef struct
rvec inv_len;
int**** atoms;
- int*** top;
- int*** mark;
- int*** start;
- int*** end;
+ int*** top;
+ int*** mark;
+ int*** start;
+ int*** end;
ivec**** nbrs;
rvec**** nbrs_cp;
} grid;
@@ -499,21 +552,34 @@ typedef struct
int reneighbor;
real vlist_cut;
real nbr_cut;
- real r_cut, r_sp_cut, r_low; // upper and lower taper
+ real r_cut;
+ real r_sp_cut;
+ real r_low; // upper and lower taper
real bo_cut;
real thb_cut;
real hb_cut;
- real Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
- int max_far_nbrs;
-
- real T_init, T_final, T;
+ real Tap7;
+ real Tap6;
+ real Tap5;
+ real Tap4;
+ real Tap3;
+ real Tap2;
+ real Tap1;
+ real Tap0;
+ int max_far_nbrs;
+
+ real T_init;
+ real T_final;
+ real T;
real Tau_T;
- int T_mode;
- real T_rate, T_freq;
+ int T_mode;
+ real T_rate;
+ real T_freq;
real Tau_PT;
- rvec P, Tau_P;
- int press_mode;
+ rvec P;
+ rvec Tau_P;
+ int press_mode;
real compressibility;
int remove_CoM_vel;
@@ -674,9 +740,13 @@ typedef struct
typedef struct
{
int thb;
- int pthb; /* pointer to the third body on the central atom's nbrlist */
- real theta, cos_theta;
- rvec dcos_di, dcos_dj, dcos_dk;
+ /* pointer to the third body on the central atom's nbrlist */
+ int pthb;
+ real theta;
+ real cos_theta;
+ rvec dcos_di;
+ rvec dcos_dj;
+ rvec dcos_dk;
} three_body_interaction_data;
@@ -684,7 +754,7 @@ typedef struct
{
int nbr;
ivec rel_box;
- // rvec ext_factor;
+// rvec ext_factor;
real d;
rvec dvec;
} near_neighbor_data;
@@ -694,10 +764,9 @@ typedef struct
{
int nbr;
ivec rel_box;
- // rvec ext_factor;
+// rvec ext_factor;
real d;
rvec dvec;
- // real H; //, Tap, inv_dr3gamij_1, inv_dr3gamij_3;
} far_neighbor_data;
@@ -719,26 +788,46 @@ typedef struct
typedef struct
{
int wrt;
- rvec dBO, dBOpi, dBOpi2;
+ rvec dBO;
+ rvec dBOpi;
+ rvec dBOpi2;
} dbond_data;
+
typedef struct
{
- real BO, BO_s, BO_pi, BO_pi2;
- real Cdbo, Cdbopi, Cdbopi2;
- real C1dbo, C2dbo, C3dbo;
- real C1dbopi, C2dbopi, C3dbopi, C4dbopi;
- real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
- rvec dBOp, dln_BOp_s, dln_BOp_pi, dln_BOp_pi2;
+ real BO;
+ real BO_s;
+ real BO_pi;
+ real BO_pi2;
+ real Cdbo;
+ real Cdbopi;
+ real Cdbopi2;
+ real C1dbo;
+ real C2dbo;
+ real C3dbo;
+ real C1dbopi;
+ real C2dbopi;
+ real C3dbopi;
+ real C4dbopi;
+ real C1dbopi2;
+ real C2dbopi2;
+ real C3dbopi2;
+ real C4dbopi2;
+ rvec dBOp;
+ rvec dln_BOp_s;
+ rvec dln_BOp_pi;
+ rvec dln_BOp_pi2;
} bond_order_data;
+
typedef struct
{
int nbr;
int sym_index;
int dbond_index;
ivec rel_box;
- // rvec ext_factor;
+// rvec ext_factor;
real d;
rvec dvec;
bond_order_data bo_data;
@@ -753,11 +842,11 @@ typedef struct
* n: number of rows
* start: row pointer (last element contains ACTUAL NNZ)
* j: column index for corresponding matrix entry
- * val: matrix entry
- * */
+ * val: matrix entry */
typedef struct
{
- unsigned int n, m;
+ unsigned int n;
+ unsigned int m;
unsigned int *start;
unsigned int *j;
real *val;
@@ -781,37 +870,63 @@ typedef struct
{
/* bond order related storage */
real *total_bond_order;
- real *Deltap, *Deltap_boc;
- real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc;
- real *dDelta_lp, *dDelta_lp_temp;
- real *nlp, *nlp_temp, *Clp, *vlpex;
+ real *Deltap;
+ real *Deltap_boc;
+ real *Delta;
+ real *Delta_lp;
+ real *Delta_lp_temp;
+ real *Delta_e;
+ real *Delta_boc;
+ real *dDelta_lp;
+ real *dDelta_lp_temp;
+ real *nlp;
+ real *nlp_temp;
+ real *Clp;
+ real *vlpex;
rvec *dDeltap_self;
/* charge method storage */
- sparse_matrix *H, *H_sp, *L, *U;
+ sparse_matrix *H;
+ sparse_matrix *H_sp;
+ sparse_matrix *L;
+ sparse_matrix *U;
real *droptol;
real *w;
real *Hdia_inv;
- real *b, *b_s, *b_t, *b_prc, *b_prm;
- real **s, **t;
- real *s_t; //, *s_old, *t_old, *s_oldest, *t_oldest;
+ real *b;
+ real *b_s;
+ real *b_t;
+ real *b_prc;
+ real *b_prm;
+ real **s;
+ real **t;
/* GMRES related storage */
- real *y, *z, *g;
- real *hc, *hs;
- real **h, **rn, **v;
+ real *y;
+ real *z;
+ real *g;
+ real *hc;
+ real *hs;
+ real **h;
+ real **rn;
+ real **v;
/* CG related storage */
- real *r, *d, *q, *p;
- int s_dims, t_dims;
+ real *r;
+ real *d;
+ real *q;
+ real *p;
int num_H;
int *hbond_index; // for hydrogen bonds
- rvec *v_const, *f_old, *a; // used in integrators
+ rvec *v_const;
+ rvec *f_old;
+ rvec *a; // used in integrators
real *CdDelta; // coefficient of dDelta for force calculations
- int *mark, *old_mark; // storage for analysis
+ int *mark;
+ int *old_mark; // storage for analysis
rvec *x_old;
/* storage space for bond restrictions */
@@ -841,7 +956,8 @@ typedef struct
rvec *f_hb;
rvec *f_tor;
rvec *f_con;
- rvec *dDelta; /* Calculated on the fly in bond_orders.c */
+ /* Calculated on the fly in bond_orders.c */
+ rvec *dDelta;
#endif
} static_storage;
@@ -853,7 +969,6 @@ typedef struct
int num_intrs;
int *index;
int *end_index;
- int type;
union
{
void *v;
@@ -874,45 +989,54 @@ typedef struct
FILE *out;
FILE *pot;
FILE *log;
- FILE *mol, *ign;
+ FILE *mol;
+ FILE *ign;
FILE *dpl;
FILE *drft;
FILE *pdb;
FILE *prs;
- int write_steps;
- int traj_compress;
- int traj_format;
+ int write_steps;
+ int traj_compress;
+ int traj_format;
char traj_title[81];
- int atom_format;
- int bond_info;
- int angle_info;
+ int atom_format;
+ int bond_info;
+ int angle_info;
- int restart_format;
- int restart_freq;
- int debug_level;
- int energy_update_freq;
+ int restart_format;
+ int restart_freq;
+ int debug_level;
+ int energy_update_freq;
- // trajectory output functions
+ /* trajectory output function pointer definitions */
int (* write_header)( reax_system*, control_params*, static_storage*, void* );
int (* append_traj_frame)(reax_system*, control_params*,
- simulation_data*, static_storage*, list **, void* );
+ simulation_data*, static_storage*, list **, void* );
int (* write)( FILE *, const char *, ... );
#ifdef TEST_ENERGY
FILE *ebond;
- FILE *elp, *eov, *eun;
- FILE *eval, *epen, *ecoa;
+ FILE *elp;
+ FILE *eov;
+ FILE *eun;
+ FILE *eval;
+ FILE *epen;
+ FILE *ecoa;
FILE *ehb;
- FILE *etor, *econ;
- FILE *evdw, *ecou;
+ FILE *etor;
+ FILE *econ;
+ FILE *evdw;
+ FILE *ecou;
#endif
FILE *ftot;
#ifdef TEST_FORCES
- FILE *fbo, *fdbo;
+ FILE *fbo;
+ FILE *fdbo;
FILE *fbond;
- FILE *flp, *fatom;
+ FILE *flp;
+ FILE *fatom;
FILE *f3body;
FILE *fhb;
FILE *f4body;
@@ -922,33 +1046,32 @@ typedef struct
} output_controls;
-typedef struct
-{
- int atom_count;
- int atom_list[MAX_MOLECULE_SIZE];
- int mtypes[MAX_ATOM_TYPES];
-} molecule;
-
-
typedef struct
{
real H;
- real e_vdW, CEvd;
- real e_ele, CEclmb;
+ real e_vdW;
+ real CEvd;
+ real e_ele;
+ real CEclmb;
} LR_data;
-
typedef struct
{
- real a, b, c, d;
+ real a;
+ real b;
+ real c;
+ real d;
} cubic_spline_coef;
+
typedef struct
{
- real xmin, xmax;
+ real xmin;
+ real xmax;
int n;
- real dx, inv_dx;
+ real dx;
+ real inv_dx;
real a;
real m;
@@ -960,36 +1083,35 @@ typedef struct
typedef struct
{
- real xmin, xmax;
+ real xmin;
+ real xmax;
int n;
- real dx, inv_dx;
+ real dx;
+ real inv_dx;
real a;
real m;
real c;
LR_data *y;
cubic_spline_coef *H;
- cubic_spline_coef *vdW, *CEvd;
- cubic_spline_coef *ele, *CEclmb;
+ cubic_spline_coef *vdW;
+ cubic_spline_coef *CEvd;
+ cubic_spline_coef *ele;
+ cubic_spline_coef *CEclmb;
} LR_lookup_table;
+/* Function pointer definitions */
typedef void (*interaction_function)(reax_system*, control_params*,
simulation_data*, static_storage*, list**, output_controls*);
-interaction_function Interaction_Functions[NO_OF_INTERACTIONS];
-
typedef void (*evolve_function)(reax_system*, control_params*,
simulation_data*, static_storage*,
list**, output_controls*);
-typedef real (*lookup_function)(real);
-lookup_table Exp, Sqrt, Cube_Root, Four_Third_Root, Cos, Sin, ACos;
+/* Global variables */
LR_lookup_table **LR;
-typedef void (*get_far_neighbors_function)(rvec, rvec, simulation_box*,
- control_params*, far_neighbor_data*, int*);
-
#endif
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 1208c826..1963bd14 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "neighbors.h"
+
#include "box.h"
#include "grid.h"
#include "list.h"
@@ -28,6 +29,10 @@
#include "vector.h"
+/* Function pointer definitions */
+typedef void (*get_far_neighbors_function)(rvec, rvec, simulation_box*,
+ control_params*, far_neighbor_data*, int*);
+
static inline real DistSqr_to_CP( rvec cp, rvec x )
{
@@ -47,14 +52,14 @@ static inline real DistSqr_to_CP( rvec cp, rvec x )
void Generate_Neighbor_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
- int i, j, k, l, m, itr;
- int x, y, z;
- int atom1, atom2, max;
- int num_far;
- int *nbr_atoms;
+ int i, j, k, l, m, itr;
+ int x, y, z;
+ int atom1, atom2, max;
+ int num_far;
+ int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
grid *g;
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 06b7d63d..9ff0f571 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "print_utils.h"
+
#include "list.h"
#include "geo_tools.h"
#include "system_props.h"
@@ -28,15 +29,15 @@
#ifdef TEST_FORCES
void Dummy_Printer( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
}
void Print_Bond_Orders( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i, pj, pk;
bond_order_data *bo_ij;
@@ -44,12 +45,12 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
list *dBOs = (*lists) + DBO;
dbond_data *dbo_k;
-
/* bond orders */
fprintf( out_control->fbo, "%6s%6s%12s%12s%12s%12s%12s\n",
"atom1", "atom2", "r_ij", "total_bo", "bo_s", "bo_p", "bo_pp" );
for ( i = 0; i < system->N; ++i )
+ {
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
bo_ij = &(bonds->select.bond_list[pj].bo_data);
@@ -61,12 +62,14 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
bonds->select.bond_list[pj].d,
bo_ij->BO, bo_ij->BO_s, bo_ij->BO_pi, bo_ij->BO_pi2 );
}
+ }
/* derivatives of bond orders */
/* fprintf( out_control->fbo, "%6s%6s%10s%10s%10s%10s\n",
"atom1", "atom2", "total_bo", "bo_s", "bo_p", "bo_pp"\n ); */
for ( i = 0; i < system->N; ++i )
+ {
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
/*fprintf( out_control->fdbo, "%6d %6d\tstart: %6d\tend: %6d\n",
@@ -98,14 +101,15 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
dbo_k->dBOpi2[0], dbo_k->dBOpi2[1], dbo_k->dBOpi2[2] );
}
}
+ }
- fflush(out_control->fdbo);
+ fflush( out_control->fdbo );
}
void Print_Bond_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i;
@@ -113,15 +117,17 @@ void Print_Bond_Forces( reax_system *system, control_params *control,
fprintf( out_control->fbond, "%6s\t%s\n", "atom", "fbond" );
for ( i = 0; i < system->N; ++i )
+ {
fprintf(out_control->fbond, "%6d %23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
workspace->f_be[i][0], workspace->f_be[i][1], workspace->f_be[i][2]);
+ }
}
void Print_LonePair_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i;
@@ -129,18 +135,19 @@ void Print_LonePair_Forces( reax_system *system, control_params *control,
fprintf( out_control->flp, "%6s\t%s\n", "atom", "f_lonepair" );
for ( i = 0; i < system->N; ++i )
+ {
fprintf(out_control->flp, "%6d %23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
workspace->f_lp[i][0], workspace->f_lp[i][1], workspace->f_lp[i][2]);
+ }
- fflush(out_control->flp);
+ fflush( out_control->flp );
}
void Print_OverUnderCoor_Forces( reax_system *system, control_params *control,
- simulation_data *data,
- static_storage *workspace, list **lists,
- output_controls *out_control )
+ simulation_data *data, static_storage *workspace, list **lists,
+ output_controls *out_control )
{
int i;
@@ -151,11 +158,14 @@ void Print_OverUnderCoor_Forces( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
if ( rvec_isZero( workspace->f_un[i] ) )
+ {
fprintf( out_control->fatom,
"%6d %23.15e%23.15e%23.15e 0 0 0\n",
workspace->orig_id[i], workspace->f_ov[i][0],
workspace->f_ov[i][1], workspace->f_ov[i][2] );
+ }
else
+ {
fprintf( out_control->fatom,
"%6d %23.15e%23.15e%23.15e %23.15e%23.15e%23.15e"\
"%23.15e%23.15e%23.15e\n",
@@ -167,15 +177,16 @@ void Print_OverUnderCoor_Forces( reax_system *system, control_params *control,
workspace->f_ov[i][2],
workspace->f_un[i][0], workspace->f_un[i][1],
workspace->f_un[i][2] );
+ }
}
- fflush(out_control->fatom);
+ fflush( out_control->fatom );
}
void Print_Three_Body_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int j;
@@ -186,10 +197,13 @@ void Print_Three_Body_Forces( reax_system *system, control_params *control,
for ( j = 0; j < system->N; ++j )
{
if ( rvec_isZero(workspace->f_pen[j]) && rvec_isZero(workspace->f_coa[j]) )
+ {
fprintf( out_control->f3body, "%6d %23.15e%23.15e%23.15e 0 0 0 0 0 0\n",
workspace->orig_id[j], workspace->f_ang[j][0],
workspace->f_ang[j][1], workspace->f_ang[j][2] );
+ }
else if ( rvec_isZero(workspace->f_coa[j]) )
+ {
fprintf( out_control->f3body,
"%6d %23.15e%23.15e%23.15e %23.15e%23.15e%23.15e "\
"%23.15e%23.15e%23.15e\n",
@@ -201,6 +215,7 @@ void Print_Three_Body_Forces( reax_system *system, control_params *control,
workspace->f_ang[j][2],
workspace->f_pen[j][0], workspace->f_pen[j][1],
workspace->f_pen[j][2] );
+ }
else
{
fprintf( out_control->f3body, "%6d %23.15e%23.15e%23.15e ",
@@ -224,14 +239,13 @@ void Print_Three_Body_Forces( reax_system *system, control_params *control,
}
}
- fflush(out_control->f3body);
+ fflush( out_control->f3body );
}
void Print_Hydrogen_Bond_Forces( reax_system *system, control_params *control,
- simulation_data *data,
- static_storage *workspace, list **lists,
- output_controls *out_control )
+ simulation_data *data, static_storage *workspace, list **lists,
+ output_controls *out_control )
{
int j;
@@ -239,17 +253,19 @@ void Print_Hydrogen_Bond_Forces( reax_system *system, control_params *control,
fprintf( out_control->fhb, "%6s\t%-38s\n", "atom", "f_hb" );
for ( j = 0; j < system->N; ++j )
+ {
fprintf(out_control->fhb, "%6d\t[%23.15e%23.15e%23.15e]\n",
workspace->orig_id[j],
workspace->f_hb[j][0], workspace->f_hb[j][1], workspace->f_hb[j][2]);
+ }
fflush(out_control->fhb);
}
void Print_Four_Body_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int j;
@@ -259,6 +275,7 @@ void Print_Four_Body_Forces( reax_system *system, control_params *control,
for ( j = 0; j < system->N; ++j )
{
if ( !rvec_isZero( workspace->f_con[j] ) )
+ {
fprintf( out_control->f4body,
"%6d %23.15e%23.15e%23.15e %23.15e%23.15e%23.15e "\
"%23.15e%23.15e%23.15e\n",
@@ -270,20 +287,23 @@ void Print_Four_Body_Forces( reax_system *system, control_params *control,
workspace->f_tor[j][2],
workspace->f_con[j][0], workspace->f_con[j][1],
workspace->f_con[j][2] );
+ }
else
+ {
fprintf( out_control->f4body,
"%6d %23.15e%23.15e%23.15e 0 0 0\n",
workspace->orig_id[j], workspace->f_tor[j][0],
workspace->f_tor[j][1], workspace->f_tor[j][2] );
+ }
}
- fflush(out_control->f4body);
+ fflush( out_control->f4body );
}
void Print_vdW_Coulomb_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i;
@@ -292,7 +312,9 @@ void Print_vdW_Coulomb_Forces( reax_system *system, control_params *control,
"atom", "nonbonded total", "f_vdw", "f_ele" );
for ( i = 0; i < system->N; ++i )
+ {
if ( !rvec_isZero(workspace->f_ele[i]) )
+ {
fprintf(out_control->fnonb,
"%6d %23.15e%23.15e%23.15e %23.15e%23.15e%23.15e "\
"%23.15e%23.15e%23.15e\n",
@@ -304,19 +326,23 @@ void Print_vdW_Coulomb_Forces( reax_system *system, control_params *control,
workspace->f_vdw[i][2],
workspace->f_ele[i][0], workspace->f_ele[i][1],
workspace->f_ele[i][2] );
+ }
else
+ {
fprintf(out_control->fnonb,
"%6d %23.15e%23.15e%23.15e 0 0 0\n",
workspace->orig_id[i], workspace->f_vdw[i][0],
workspace->f_vdw[i][1], workspace->f_vdw[i][2] );
+ }
+ }
- fflush(out_control->fnonb);
+ fflush( out_control->fnonb );
}
void Compare_Total_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i;
@@ -325,6 +351,7 @@ void Compare_Total_Forces( reax_system *system, control_params *control,
"atom", "f_total", "test_force total" );
for ( i = 0; i < system->N; ++i )
+ {
fprintf( out_control->ftot2,
"%6d %23.15e%23.15e%23.15e vs %23.15e%23.15e%23.15e\n",
workspace->orig_id[i],
@@ -347,8 +374,9 @@ void Compare_Total_Forces( reax_system *system, control_params *control,
workspace->f_coa[i][2] + workspace->f_hb[i][2] +
workspace->f_tor[i][2] + workspace->f_con[i][2] +
workspace->f_vdw[i][2] + workspace->f_ele[i][2] );
+ }
- fflush(out_control->ftot2);
+ fflush( out_control->ftot2 );
}
@@ -370,18 +398,16 @@ void Init_Force_Test_Functions( )
/* near nbrs contain both i-j, j-i nbrhood info */
void Print_Near_Neighbors( reax_system *system, control_params *control,
- static_storage *workspace, list **lists )
+ static_storage *workspace, list **lists )
{
- int i, j, id_i, id_j;
- char fname[MAX_STR];
+ int i, j, id_i, id_j;
+ char fname[MAX_STR];
FILE *fout;
list *near_nbrs = &((*lists)[NEAR_NBRS]);
sprintf( fname, "%s.near_nbrs", control->sim_name );
fout = fopen( fname, "w" );
- fprintf( fout, "hello:!\n" );
-
for ( i = 0; i < system->N; ++i )
{
id_i = workspace->orig_id[i];
@@ -406,10 +432,10 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
/* near nbrs contain both i-j, j-i nbrhood info */
void Print_Near_Neighbors2( reax_system *system, control_params *control,
- static_storage *workspace, list **lists )
+ static_storage *workspace, list **lists )
{
- int i, j, id_i, id_j;
- char fname[MAX_STR];
+ int i, j, id_i, id_j;
+ char fname[MAX_STR];
FILE *fout;
list *near_nbrs = &((*lists)[NEAR_NBRS]);
@@ -441,10 +467,10 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
/* far nbrs contain only i-j nbrhood info, no j-i. */
void Print_Far_Neighbors( reax_system *system, control_params *control,
- static_storage *workspace, list **lists )
+ static_storage *workspace, list **lists )
{
- int i, j, id_i, id_j;
- char fname[MAX_STR];
+ int i, j, id_i, id_j;
+ char fname[MAX_STR];
FILE *fout;
list *far_nbrs = &((*lists)[FAR_NBRS]);
@@ -486,10 +512,10 @@ int fn_qsort_intcmp( const void *a, const void *b )
void Print_Far_Neighbors2( reax_system *system, control_params *control,
- static_storage *workspace, list **lists )
+ static_storage *workspace, list **lists )
{
- int i, j, id_i, id_j;
- char fname[MAX_STR];
+ int i, j, id_i, id_j;
+ char fname[MAX_STR];
FILE *fout;
list *far_nbrs = &((*lists)[FAR_NBRS]);
@@ -514,13 +540,14 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
fprintf(fout, "%6d", temp[j]);
fprintf( fout, "\n");
}
+
fclose( fout );
}
void Print_Total_Force( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
int i;
#if !defined(TEST_FORCES)
@@ -530,51 +557,54 @@ void Print_Total_Force( reax_system *system, control_params *control,
#endif
for ( i = 0; i < system->N; ++i )
- fprintf(out_control->ftot, "%6d %23.15e %23.15e %23.15e\n",
+ {
+ fprintf( out_control->ftot, "%6d %23.15e %23.15e %23.15e\n",
//fprintf(out_control->ftot, "%6d %19.9e %19.9e %19.9e\n",
//fprintf(out_control->ftot, "%3d %12.6f %12.6f %12.6f\n",
workspace->orig_id[i],
- system->atoms[i].f[0], system->atoms[i].f[1], system->atoms[i].f[2]);
+ system->atoms[i].f[0], system->atoms[i].f[1], system->atoms[i].f[2] );
+ }
- fflush(out_control->ftot);
+ fflush( out_control->ftot );
#if !defined(TEST_FORCES)
- fclose(out_control->ftot);
+ fclose( out_control->ftot );
#endif
}
void Output_Results( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ list **lists, output_controls *out_control )
{
- int i, type_i, f_update;
- real q;
- real t_elapsed = 0;;
-
+ int i, type_i;
+ real e_pol, q, f_update;
+ real t_elapsed = 0;
/* Compute Polarization Energy */
- data->E_Pol = 0.0;
+ e_pol = 0.0;
+
+#ifdef _OPENMP
+ #pragma omp parallel for default(none) private(q, type_i,) shared(system) \
+ reduction(+: e_pol) schedule(static)
+#endif
for ( i = 0; i < system->N; i++ )
{
q = system->atoms[i].q;
type_i = system->atoms[i].type;
- data->E_Pol += ( system->reaxprm.sbp[ type_i ].chi * q +
- (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR(q) ) *
- KCALpMOL_to_EV;
- /* fprintf( stderr, "%6d%23.15e%23.15e%23.15e%23.15e\n",
- i, q, system->reaxprm.sbp[ type_i ].chi,
- system->reaxprm.sbp[ type_i ].eta, data->E_Pol ); */
+ e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
+ (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
+ KCALpMOL_to_EV;
}
+ data->E_Pol = e_pol;
+
data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
- data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
- data->E_Tor + data->E_Con +
- data->E_vdW + data->E_Ele + data->E_Pol;
+ data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
+ data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
data->E_Tot = data->E_Pot + E_CONV * data->E_Kin;
-
/* output energies if it is the time */
if ( out_control->energy_update_freq > 0 &&
data->step % out_control->energy_update_freq == 0 )
@@ -584,7 +614,7 @@ void Output_Results( reax_system *system, control_params *control,
"%-6d%24.15e%24.15e%24.15e%13.5f%13.5f%16.5f%13.5f%13.5f\n",
data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
data->therm.T, control->T, system->box.volume, data->iso_bar.P,
- (control->P[0] + control->P[1] + control->P[2]) / 3 );
+ (control->P[0] + control->P[1] + control->P[2]) / 3.0 );
fprintf( out_control->pot,
"%-6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
@@ -599,7 +629,7 @@ void Output_Results( reax_system *system, control_params *control,
"%-6d%16.2f%16.2f%16.2f%11.2f%11.2f%13.2f%13.5f%13.5f\n",
data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
data->therm.T, control->T, system->box.volume, data->iso_bar.P,
- (control->P[0] + control->P[1] + control->P[2]) / 3 );
+ (control->P[0] + control->P[1] + control->P[2]) / 3.0 );
fprintf( out_control->pot,
"%-6d%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f%13.2f\n",
@@ -613,24 +643,29 @@ void Output_Results( reax_system *system, control_params *control,
t_elapsed = Get_Timing_Info( data->timing.total );
if ( data->step == data->prev_steps )
- f_update = 1;
- else f_update = out_control->energy_update_freq;
+ {
+ f_update = 1.0;
+ }
+ else
+ {
+ f_update = 1.0 / out_control->energy_update_freq;
+ }
fprintf( out_control->log, "%6d %10.2f %10.2f %10.2f %10.2f %10.2f %10.4f %10.4f %10.2f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n",
- data->step, t_elapsed / f_update,
- data->timing.nbrs / f_update,
- data->timing.init_forces / f_update,
- data->timing.bonded / f_update,
- data->timing.nonb / f_update,
- data->timing.cm / f_update,
- data->timing.cm_sort_mat_rows / f_update,
- (double)data->timing.cm_solver_iters / f_update,
- data->timing.cm_solver_pre_comp / f_update,
- data->timing.cm_solver_pre_app / f_update,
- data->timing.cm_solver_spmv / f_update,
- data->timing.cm_solver_vector_ops / f_update,
- data->timing.cm_solver_orthog / f_update,
- data->timing.cm_solver_tri_solve / f_update );
+ data->step, t_elapsed * f_update,
+ data->timing.nbrs * f_update,
+ data->timing.init_forces * f_update,
+ data->timing.bonded * f_update,
+ data->timing.nonb * f_update,
+ data->timing.cm * f_update,
+ data->timing.cm_sort_mat_rows * f_update,
+ (double)data->timing.cm_solver_iters * f_update,
+ data->timing.cm_solver_pre_comp * f_update,
+ data->timing.cm_solver_pre_app * f_update,
+ data->timing.cm_solver_spmv * f_update,
+ data->timing.cm_solver_vector_ops * f_update,
+ data->timing.cm_solver_orthog * f_update,
+ data->timing.cm_solver_tri_solve * f_update );
data->timing.total = Get_Time( );
data->timing.nbrs = 0;
@@ -672,23 +707,21 @@ void Output_Results( reax_system *system, control_params *control,
system->box.box_norms[2],
data->tot_press[0], data->tot_press[1], data->tot_press[2],
control->P[0], control->P[1], control->P[2], system->box.volume );
- fflush( out_control->prs);
+ fflush( out_control->prs );
}
}
if ( out_control->write_steps > 0 &&
data->step % out_control->write_steps == 0 )
{
- // t_start = Get_Time( );
+ //t_start = Get_Time( );
out_control->append_traj_frame( system, control, data,
- workspace, lists, out_control );
+ workspace, lists, out_control );
//Write_PDB( system, *lists+BONDS, data, control, workspace, out_control );
- // t_elapsed = Get_Timing_Info( t_start );
- // fprintf(stdout, "append_frame took %.6f seconds\n", t_elapsed );
+ //t_elapsed = Get_Timing_Info( t_start );
+ //fprintf(stdout, "append_frame took %.6f seconds\n", t_elapsed );
}
-
- // fprintf( stderr, "output_results... done\n" );
}
@@ -780,9 +813,9 @@ void Print_Linear_System( reax_system *system, control_params *control,
void Print_Charges( reax_system *system, control_params *control,
- static_storage *workspace, int step )
+ static_storage *workspace, int step )
{
- int i;
+ int i;
char fname[100];
FILE *fout;
@@ -790,24 +823,28 @@ void Print_Charges( reax_system *system, control_params *control,
fout = fopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
+ {
fprintf( fout, "%6d%12.7f%12.7f%12.7f\n",
workspace->orig_id[i],
workspace->s[0][i], workspace->t[0][i], system->atoms[i].q );
+ }
fclose( fout );
}
void Print_Soln( static_storage *workspace,
- real *x, real *b_prm, real *b, int N )
+ real *x, real *b_prm, real *b, int N )
{
int i;
fprintf( stdout, "%6s%10s%10s%10s\n", "id", "x", "b_prm", "b" );
for ( i = 0; i < N; ++i )
+ {
fprintf( stdout, "%6d%10.4f%10.4f%10.4f\n",
workspace->orig_id[i], x[i], b_prm[i], b[i] );
+ }
fflush( stdout );
}
@@ -821,7 +858,9 @@ void Print_Sparse_Matrix( sparse_matrix *A )
{
fprintf( stderr, "i:%d j(val):", i );
for ( j = A->start[i]; j < A->start[i + 1]; ++j )
+ {
fprintf( stderr, "%d(%.4f) ", A->j[j], A->val[j] );
+ }
fprintf( stderr, "\n" );
}
}
@@ -843,15 +882,20 @@ void Print_Sparse_Matrix2( sparse_matrix *A, char *fname, char *mode )
for ( i = 0; i < A->n; ++i )
{
- for ( j = A->start[i]; j < A->start[i + 1]; ++j )
+ /* off-diagonals */
+ for ( j = A->start[i]; j < A->start[i + 1] - 1; ++j )
{
- //fprintf( f, "%d%d %.15e\n", A->entries[j].j, i, A->entries[j].val );
//Convert 0-based to 1-based (for Matlab)
- fprintf( f, "%6d %6d %24.15e\n", i+1, A->j[j]+1, A->val[j] );
+ fprintf( f, "%6d %6d %24.15e\n", i + 1, A->j[j] + 1, A->val[j] );
+ /* print symmetric entry */
+// fprintf( f, "%6d %6d %24.15e\n", A->j[j] + 1, i + 1, A->val[j] );
}
+
+ /* diagonal */
+ fprintf( f, "%6d %6d %24.15e\n", i + 1, A->j[A->start[i + 1] - 1] + 1, A->val[A->start[i + 1] - 1] );
}
- fclose(f);
+ fclose( f );
}
@@ -900,6 +944,7 @@ void Print_Bonds( reax_system *system, list *bonds, char *fname )
FILE *f = fopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
+ {
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
pbond = &(bonds->select.bond_list[pj]);
@@ -910,7 +955,9 @@ void Print_Bonds( reax_system *system, list *bonds, char *fname )
fprintf( f, "%6d%6d %9.5f %9.5f\n",
i + 1, pbond->nbr + 1, pbond->d, bo_ij->BO );
}
- fclose(f);
+ }
+
+ fclose( f );
}
@@ -931,35 +978,41 @@ void Print_Bond_List2( reax_system *system, list *bonds, char *fname )
nbr = bonds->select.bond_list[pj].nbr;
id_j = nbr + 1; //system->my_atoms[nbr].orig_id;
if ( id_i < id_j )
+ {
temp[num++] = id_j;
+ }
}
qsort(&temp, num, sizeof(int), fn_qsort_intcmp);
- for (j = 0; j < num; j++)
- fprintf(f, "%6d", temp[j] );
- fprintf(f, "\n");
+ for ( j = 0; j < num; j++ )
+ {
+ fprintf( f, "%6d", temp[j] );
+ }
+ fprintf( f, "\n" );
}
}
#ifdef LGJ
-Print_XYZ_Serial(reax_system* system, static_storage *workspace)
+Print_XYZ_Serial( reax_system* system, static_storage *workspace )
{
rvec p;
-
- char fname[100];
+ char fname[100];
FILE *fout;
- sprintf( fname, "READ_PDB.0" );
- fout = fopen( fname, "w" );
int i;
- for (i = 0; i < system->N; i++)
+
+ sprintf( fname, "READ_PDB.0" );
+ fout = fopen( fname, "w" );
+
+ for ( i = 0; i < system->N; i++ )
+ {
fprintf( fout, "%6d%24.15e%24.15e%24.15e\n",
workspace->orig_id[i],
p[0] = system->atoms[i].x[0],
p[1] = system->atoms[i].x[1],
- p[2] = system->atoms[i].x[2]);
-
+ p[2] = system->atoms[i].x[2] );
+ }
- fclose(fout);
+ fclose( fout );
}
#endif
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index 915a10db..0a364e69 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "reset_utils.h"
+
#include "list.h"
#include "vector.h"
@@ -29,36 +30,36 @@ void Reset_Atoms( reax_system* system )
int i;
for ( i = 0; i < system->N; ++i )
+ {
memset( system->atoms[i].f, 0.0, sizeof(rvec) );
+ }
}
void Reset_Pressures( simulation_data *data )
{
rtensor_MakeZero( data->flex_bar.P );
- data->iso_bar.P = 0;
+ data->iso_bar.P = 0.0;
rvec_MakeZero( data->int_press );
rvec_MakeZero( data->ext_press );
- /* fprintf( stderr, "reset: ext_press (%12.6f %12.6f %12.6f)\n",
- data->ext_press[0], data->ext_press[1], data->ext_press[2] ); */
}
void Reset_Simulation_Data( simulation_data* data )
{
- data->E_BE = 0;
- data->E_Ov = 0;
- data->E_Un = 0;
- data->E_Lp = 0;
- data->E_Ang = 0;
- data->E_Pen = 0;
- data->E_Coa = 0;
- data->E_HB = 0;
- data->E_Tor = 0;
- data->E_Con = 0;
- data->E_vdW = 0;
- data->E_Ele = 0;
- data->E_Kin = 0;
+ data->E_BE = 0.0;
+ data->E_Ov = 0.0;
+ data->E_Un = 0.0;
+ data->E_Lp = 0.0;
+ data->E_Ang = 0.0;
+ data->E_Pen = 0.0;
+ data->E_Coa = 0.0;
+ data->E_HB = 0.0;
+ data->E_Tor = 0.0;
+ data->E_Con = 0.0;
+ data->E_vdW = 0.0;
+ data->E_Ele = 0.0;
+ data->E_Kin = 0.0;
}
@@ -150,11 +151,7 @@ void Reset( reax_system *system, control_params *control,
Reset_Simulation_Data( data );
- if ( control->ensemble == NPT || control->ensemble == sNPT ||
- control->ensemble == iNPT )
- {
- Reset_Pressures( data );
- }
+ Reset_Pressures( data );
Reset_Workspace( system, workspace );
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 484104a8..9da3e21a 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -26,11 +26,8 @@
#include "vector.h"
-void LonePair_OverUnder_Coordination_Energy( reax_system *system,
- control_params *control,
- simulation_data *data,
- static_storage *workspace,
- list **lists,
+void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params *control,
+ simulation_data *data, static_storage *workspace, list **lists,
output_controls *out_control )
{
int i, j, pj, type_i, type_j;
@@ -44,8 +41,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system,
real e_un, CEunder1, CEunder2, CEunder3, CEunder4;
real p_lp1, p_lp2, p_lp3;
real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
-
- single_body_parameters *sbp_i, *sbp_j;
+ single_body_parameters *sbp_i;
two_body_parameters *twbp;
bond_data *pbond;
bond_order_data *bo_ij;
@@ -152,18 +148,11 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system,
j = bonds->select.bond_list[pj].nbr;
type_j = system->atoms[j].type;
bo_ij = &(bonds->select.bond_list[pj].bo_data);
- sbp_j = &(system->reaxprm.sbp[ type_j ]);
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]) *
( bo_ij->BO_pi + bo_ij->BO_pi2 );
-
- /*fprintf( stdout, "%4d%4d%23.15e%23.15e%23.15e\n",
- i+1, j+1,
- dfvl * workspace->Delta_lp_temp[j],
- sbp_j->nlp_opt,
- workspace->nlp_temp[j] );*/
}
exp_ovun1 = p_ovun3 * EXP( p_ovun4 * sum_ovun2 );
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index fc93a474..fad302b4 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -25,7 +25,7 @@
void Temperature_Control( control_params *control, simulation_data *data,
- output_controls *out_control )
+ output_controls *out_control )
{
real tmp;
@@ -34,17 +34,24 @@ void Temperature_Control( control_params *control, simulation_data *data,
if ( (data->step - data->prev_steps) %
((int)(control->T_freq / control->dt)) == 0 )
{
- if ( fabs( control->T - control->T_final ) >= fabs( control->T_rate ) )
+ if ( FABS( control->T - control->T_final ) >= FABS( control->T_rate ) )
+ {
control->T += control->T_rate;
- else control->T = control->T_final;
+ }
+ else
+ {
+ control->T = control->T_final;
+ }
}
}
else if ( control->T_mode == 2 ) // constant slope control
{
tmp = control->T_rate * control->dt / control->T_freq;
- if ( fabs( control->T - control->T_final ) >= fabs( tmp ) )
+ if ( FABS( control->T - control->T_final ) >= FABS( tmp ) )
+ {
control->T += tmp;
+ }
}
}
@@ -53,18 +60,19 @@ void Compute_Total_Mass( reax_system *system, simulation_data *data )
{
int i;
- data->M = 0;
+ data->M = 0.0;
for ( i = 0; i < system->N; i++ )
+ {
data->M += system->reaxprm.sbp[ system->atoms[i].type ].mass;
+ }
- //fprintf ( stderr, "Compute_total_Mass -->%f<-- \n", data->M );
- data->inv_M = 1. / data->M;
+ data->inv_M = 1.0 / data->M;
}
void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
- FILE *fout )
+ FILE *fout )
{
int i;
real m, xx, xy, xz, yy, yz, zz, det;
@@ -76,7 +84,6 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
rvec_MakeZero( data->amcm ); // angular momentum of CoM
rvec_MakeZero( data->avcm ); // angular velocity of CoM
-
/* Compute the position, velocity and angular momentum about the CoM */
for ( i = 0; i < system->N; ++i )
{
@@ -87,15 +94,6 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
rvec_Cross( tvec, system->atoms[i].x, system->atoms[i].v );
rvec_ScaledAdd( data->amcm, m, tvec );
-
- /*fprintf( fout,"%3d %g %g %g\n",
- i+1,
- system->atoms[i].v[0], system->atoms[i].v[1], system->atoms[i].v[2] );
- fprintf( fout, "vcm: %g %g %g\n",
- data->vcm[0], data->vcm[1], data->vcm[2] );
- */
- /* fprintf( stderr, "amcm: %12.6f %12.6f %12.6f\n",
- data->amcm[0], data->amcm[1], data->amcm[2] ); */
}
rvec_Scale( data->xcm, data->inv_M, data->xcm );
@@ -147,10 +145,14 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
inv[2][1] = mat[2][0] * mat[0][1] - mat[0][0] * mat[2][1];
inv[2][2] = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
- if ( fabs(det) > ALMOST_ZERO )
+ if ( FABS(det) > ALMOST_ZERO )
+ {
rtensor_Scale( inv, 1. / det, inv );
+ }
else
+ {
rtensor_MakeZero( inv );
+ }
/* Compute the angular velocity about the centre of mass */
rtensor_MatVec( data->avcm, inv, data->amcm );
@@ -186,7 +188,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
data->E_Kin = 0.0;
- for (i = 0; i < system->N; i++)
+ for ( i = 0; i < system->N; i++ )
{
m = system->reaxprm.sbp[system->atoms[i].type].mass;
@@ -200,8 +202,10 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
data->therm.T = (2. * data->E_Kin) / (data->N_f * K_B);
- if ( fabs(data->therm.T) < ALMOST_ZERO ) /* avoid T being an absolute zero! */
+ if ( FABS(data->therm.T) < ALMOST_ZERO ) /* avoid T being an absolute zero! */
+ {
data->therm.T = ALMOST_ZERO;
+ }
}
@@ -214,8 +218,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
* We may want to add that for more accuracy.
*/
void Compute_Pressure_Isotropic( reax_system* system, control_params *control,
- simulation_data* data,
- output_controls *out_control )
+ simulation_data* data, output_controls *out_control )
{
int i;
reax_atom *p_atom;
@@ -273,8 +276,7 @@ void Compute_Pressure_Isotropic( reax_system* system, control_params *control,
}
-void Compute_Pressure_Isotropic_Klein( reax_system* system,
- simulation_data* data )
+void Compute_Pressure_Isotropic_Klein( reax_system* system, simulation_data* data )
{
int i;
reax_atom *p_atom;
@@ -301,7 +303,7 @@ void Compute_Pressure_Isotropic_Klein( reax_system* system,
void Compute_Pressure( reax_system* system, simulation_data* data,
- static_storage *workspace )
+ static_storage *workspace )
{
int i;
reax_atom *p_atom;
diff --git a/sPuReMD/src/testmd.c b/sPuReMD/src/testmd.c
index 9ff007a9..d1159ef1 100644
--- a/sPuReMD/src/testmd.c
+++ b/sPuReMD/src/testmd.c
@@ -135,6 +135,9 @@ void static Read_System( char * const geo_file,
exit( INVALID_GEO );
}
+ fclose( ffield );
+ fclose( ctrl );
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "input files have been read...\n" );
Print_Box( &(system->box), stderr );
@@ -168,10 +171,10 @@ int main(int argc, char* argv[])
lists = (list*) malloc( sizeof(list) * LIST_N );
Read_System( argv[1], argv[2], argv[3], &system, &control,
- &data, &workspace, &out_control );
+ &data, &workspace, &out_control );
Initialize( &system, &control, &data, &workspace, &lists,
- &out_control, &Evolve );
+ &out_control, &Evolve );
/* compute f_0 */
//if( control.restart == 0 ) {
@@ -195,13 +198,15 @@ int main(int argc, char* argv[])
}
Evolve( &system, &control, &data, &workspace, &lists, &out_control );
Post_Evolve( &system, &control, &data, &workspace, &lists, &out_control );
- Output_Results(&system, &control, &data, &workspace, &lists, &out_control);
+ Output_Results( &system, &control, &data, &workspace, &lists, &out_control );
Analysis( &system, &control, &data, &workspace, &lists, &out_control );
steps = data.step - data.prev_steps;
if ( steps && out_control.restart_freq &&
steps % out_control.restart_freq == 0 )
+ {
Write_Restart( &system, &control, &data, &workspace, &out_control );
+ }
}
if ( out_control.write_steps > 0 )
@@ -214,5 +219,10 @@ int main(int argc, char* argv[])
data.timing.elapsed = Get_Timing_Info( data.timing.start );
fprintf( out_control.log, "total: %.2f secs\n", data.timing.elapsed );
+ Finalize( &system, &control, &data, &workspace, &lists,
+ &out_control );
+
+ free( lists );
+
return SUCCESS;
}
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 5a7eac49..c54def41 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -49,12 +49,14 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
rvec* dcos_theta_di, rvec* dcos_theta_dj, rvec* dcos_theta_dk )
{
int t;
- real sqr_d_ji = SQR(d_ji);
- real sqr_d_jk = SQR(d_jk);
- real inv_dists = 1.0 / (d_ji * d_jk);
- real inv_dists3 = POW( inv_dists, 3 );
- real dot_dvecs = rvec_Dot( dvec_ji, dvec_jk );
- real Cdot_inv3 = dot_dvecs * inv_dists3;
+ real sqr_d_ji, sqr_d_jk, inv_dists, inv_dists3, dot_dvecs, Cdot_inv3;
+
+ sqr_d_ji = SQR( d_ji );
+ sqr_d_jk = SQR( d_jk );
+ inv_dists = 1.0 / (d_ji * d_jk);
+ inv_dists3 = POW( inv_dists, 3 );
+ dot_dvecs = rvec_Dot( dvec_ji, dvec_jk );
+ Cdot_inv3 = dot_dvecs * inv_dists3;
for ( t = 0; t < 3; ++t )
{
@@ -80,494 +82,563 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
{
- int i, j, pi, k, pk, t;
- int type_i, type_j, type_k;
- int start_j, end_j, start_pk, end_pk;
- int flag, cnt, num_thb_intrs;
-
- real temp, temp_bo_jt, pBOjt7;
- real p_val1, p_val2, p_val3, p_val4, p_val5;
- real p_val6, p_val7, p_val8, p_val9, p_val10;
- real p_pen1, p_pen2, p_pen3, p_pen4;
- real p_coa1, p_coa2, p_coa3, p_coa4;
- real trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk;
- real exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2;
- real dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO;
- real CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8;
- real CEpen1, CEpen2, CEpen3;
- real e_ang, e_coa, e_pen;
- real CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5;
- real Cf7ij, Cf7jk, Cf8j, Cf9j;
- real f7_ij, f7_jk, f8_Dj, f9_Dj;
- real Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
- real r_ij, r_jk;
- real BOA_ij, BOA_jk;
- real vlpadj;
- rvec force, ext_press;
- // rtensor temp_rtensor, total_rtensor;
real *total_bo;
- three_body_header *thbh;
- three_body_parameters *thbp;
- three_body_interaction_data *p_ijk, *p_kji;
- bond_data *pbond_ij, *pbond_jk, *pbond_jt;
- bond_order_data *bo_ij, *bo_jk, *bo_jt;
list *bonds, *thb_intrs;
bond_data *bond_list;
three_body_interaction_data *thb_list;
+ real p_pen2, p_pen3, p_pen4;
+ real p_coa2, p_coa3, p_coa4;
+ real p_val6, p_val8, p_val9, p_val10;
+ int num_thb_intrs;
+ real e_ang_total, e_pen_total, e_coa_total;
total_bo = workspace->total_bond_order;
bonds = (*lists) + BONDS;
bond_list = bonds->select.bond_list;
thb_intrs = (*lists) + THREE_BODIES;
thb_list = thb_intrs->select.three_body_list;
-
/* global parameters used in these calculations */
+ p_pen2 = system->reaxprm.gp.l[19];
+ p_pen3 = system->reaxprm.gp.l[20];
+ p_pen4 = system->reaxprm.gp.l[21];
+ p_coa2 = system->reaxprm.gp.l[2];
+ p_coa3 = system->reaxprm.gp.l[38];
+ p_coa4 = system->reaxprm.gp.l[30];
p_val6 = system->reaxprm.gp.l[14];
p_val8 = system->reaxprm.gp.l[33];
p_val9 = system->reaxprm.gp.l[16];
p_val10 = system->reaxprm.gp.l[17];
num_thb_intrs = 0;
+ e_ang_total = 0.0;
+ e_pen_total = 0.0;
+ e_coa_total = 0.0;
- for ( j = 0; j < system->N; ++j )
+ //TODO: change interaction lists for parallelization
+#ifdef _OPENMP
+// #pragma omp parallel default(shared) reduction(+:total_Eang, total_Epen, total_Ecoa, num_thb_intrs)
+#endif
{
- // fprintf( out_control->eval, "j: %d\n", j );
- type_j = system->atoms[j].type;
- start_j = Start_Index(j, bonds);
- end_j = End_Index(j, bonds);
-
- p_val3 = system->reaxprm.sbp[ type_j ].p_val3;
- p_val5 = system->reaxprm.sbp[ type_j ].p_val5;
-
- SBOp = 0, prod_SBO = 1;
- for ( t = start_j; t < end_j; ++t )
- {
- 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;
- prod_SBO *= EXP( -temp );
- }
+ int i, j, pi, k, pk, t;
+ int type_i, type_j, type_k;
+ int start_j, end_j, start_pk, end_pk;
+ int cnt;
+ real temp, temp_bo_jt, pBOjt7;
+ real p_val1, p_val2, p_val3, p_val4, p_val5, p_val7;
+ real p_pen1;
+ real p_coa1;
+ real trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk;
+ real exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2;
+ real dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO;
+ real CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8;
+ real CEpen1, CEpen2, CEpen3;
+ real e_ang, e_coa, e_pen;
+ real CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5;
+ real Cf7ij, Cf7jk, Cf8j, Cf9j;
+ real f7_ij, f7_jk, f8_Dj, f9_Dj;
+ real Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
+ real BOA_ij, BOA_jk;
+ real vlpadj;
+ rvec force, ext_press;
+ //rtensor temp_rtensor, total_rtensor;
+ three_body_header *thbh;
+ three_body_parameters *thbp;
+ three_body_interaction_data *p_ijk, *p_kji;
+ bond_data *pbond_ij, *pbond_jk, *pbond_jt;
+ bond_order_data *bo_ij, *bo_jk, *bo_jt;
+ rvec *f_i, *f_j, *f_k;
+#ifdef _OPENMP
+// int tid = omp_get_thread_num( );
+#endif
- /* modifications to match Adri's code - 09/01/09 */
- if ( workspace->vlpex[j] >= 0 )
- {
- vlpadj = 0;
- dSBO2 = prod_SBO - 1;
- }
- else
+ for ( j = 0; j < system->N; ++j )
{
- vlpadj = workspace->nlp[j];
- dSBO2 = (prod_SBO - 1) * (1 - p_val8 * workspace->dDelta_lp[j]);
- }
-
- SBO = SBOp + (1 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
- dSBO1 = -8 * prod_SBO * ( workspace->Delta_boc[j] + p_val8 * vlpadj );
+ // fprintf( out_control->eval, "j: %d\n", j );
+ type_j = system->atoms[j].type;
+ start_j = Start_Index(j, bonds);
+ end_j = End_Index(j, bonds);
+//#ifdef _OPENMP
+// f_j = &(workspace->f_local[tid * system->N + j]);
+//#else
+ f_j = &(system->atoms[j].f);
+//#endif
+
+ p_val3 = system->reaxprm.sbp[ type_j ].p_val3;
+ p_val5 = system->reaxprm.sbp[ type_j ].p_val5;
+
+ 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);
+ temp = SQR( bo_jt->BO );
+ temp *= temp;
+ temp *= temp;
+ prod_SBO *= EXP( -temp );
+ }
- if ( SBO <= 0 )
- SBO2 = 0, CSBO2 = 0;
- else if ( SBO > 0 && SBO <= 1 )
- {
- SBO2 = POW( SBO, p_val9 );
- CSBO2 = p_val9 * POW( SBO, p_val9 - 1 );
- }
- else if ( SBO > 1 && SBO < 2 )
- {
- SBO2 = 2 - POW( 2 - SBO, p_val9 );
- CSBO2 = p_val9 * POW( 2 - SBO, p_val9 - 1 );
- }
- else
- SBO2 = 2, CSBO2 = 0;
+ /* modifications to match Adri's code - 09/01/09 */
+ if ( workspace->vlpex[j] >= 0.0 )
+ {
+ vlpadj = 0.0;
+ dSBO2 = prod_SBO - 1.0;
+ }
+ else
+ {
+ vlpadj = workspace->nlp[j];
+ dSBO2 = (prod_SBO - 1.0) * (1.0 - p_val8 * workspace->dDelta_lp[j]);
+ }
- expval6 = EXP( p_val6 * workspace->Delta_boc[j] );
+ SBO = SBOp + (1.0 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
+ dSBO1 = -8.0 * prod_SBO * ( workspace->Delta_boc[j] + p_val8 * vlpadj );
- /* unlike 2-body intrs where we enforce i<j, we cannot put any such
- restrictions here. such a restriction would prevent us from producing
- all 4-body intrs correctly */
- 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);
- BOA_ij = bo_ij->BO - control->thb_cut;
+ if ( SBO <= 0.0 )
+ {
+ SBO2 = 0.0;
+ CSBO2 = 0.0;
+ }
+ else if ( SBO > 0.0 && SBO <= 1.0 )
+ {
+ SBO2 = POW( SBO, p_val9 );
+ CSBO2 = p_val9 * POW( SBO, p_val9 - 1.0 );
+ }
+ else if ( SBO > 1.0 && SBO < 2.0 )
+ {
+ SBO2 = 2.0 - POW( 2.0 - SBO, p_val9 );
+ CSBO2 = p_val9 * POW( 2.0 - SBO, p_val9 - 1.0 );
+ }
+ else
+ {
+ SBO2 = 2.0;
+ CSBO2 = 0.0;
+ }
+ expval6 = EXP( p_val6 * workspace->Delta_boc[j] );
- if ( BOA_ij/*bo_ij->BO*/ > 0.0 )
+ /* unlike 2-body intrs where we enforce i<j, we cannot put any such
+ restrictions here. such a restriction would prevent us from producing
+ all 4-body intrs correctly */
+ for ( pi = start_j; pi < end_j; ++pi )
{
- i = pbond_ij->nbr;
- r_ij = pbond_ij->d;
- type_i = system->atoms[i].type;
- // fprintf( out_control->eval, "i: %d\n", i );
-
-
- /* first copy 3-body intrs from previously computed ones where i>k.
- IMPORTANT: if it is less costly to compute theta and its
- derivative, we should definitely re-compute them,
- instead of copying!
- in the second for-loop below, we compute only new 3-body intrs
- where i < k */
- for ( pk = start_j; pk < pi; ++pk )
+ Set_Start_Index( pi, num_thb_intrs, thb_intrs );
+ pbond_ij = &(bond_list[pi]);
+ bo_ij = &(pbond_ij->bo_data);
+ BOA_ij = bo_ij->BO - control->thb_cut;
+
+ if ( BOA_ij > 0.0 )
{
- // fprintf( out_control->eval, "pk: %d\n", pk );
- start_pk = Start_Index( pk, thb_intrs );
- end_pk = End_Index( pk, thb_intrs );
+ i = pbond_ij->nbr;
+ type_i = system->atoms[i].type;
+//#ifdef _OPENMP
+// f_i = &(workspace->f_local[tid * system->N + i]);
+//#else
+ f_i = &(system->atoms[i].f);
+//#endif
+
+ /* first copy 3-body intrs from previously computed ones where i>k.
+ IMPORTANT: if it is less costly to compute theta and its
+ derivative, we should definitely re-compute them,
+ instead of copying!
+ in the second for-loop below, we compute only new 3-body intrs
+ where i < k */
+ for ( pk = start_j; pk < pi; ++pk )
+ {
+ // fprintf( out_control->eval, "pk: %d\n", pk );
+ start_pk = Start_Index( pk, thb_intrs );
+ end_pk = End_Index( pk, thb_intrs );
- for ( t = start_pk; t < end_pk; ++t )
- if ( thb_list[t].thb == i )
+ for ( t = start_pk; t < end_pk; ++t )
{
- p_ijk = &(thb_list[num_thb_intrs]);
- p_kji = &(thb_list[t]);
-
- p_ijk->thb = bond_list[pk].nbr;
- p_ijk->pthb = pk;
- p_ijk->theta = p_kji->theta;
- rvec_Copy( p_ijk->dcos_di, p_kji->dcos_dk );
- rvec_Copy( p_ijk->dcos_dj, p_kji->dcos_dj );
- rvec_Copy( p_ijk->dcos_dk, p_kji->dcos_di );
-
- ++num_thb_intrs;
- break;
+ if ( thb_list[t].thb == i )
+ {
+ p_ijk = &(thb_list[num_thb_intrs]);
+ p_kji = &(thb_list[t]);
+
+ p_ijk->thb = bond_list[pk].nbr;
+ p_ijk->pthb = pk;
+ p_ijk->theta = p_kji->theta;
+ rvec_Copy( p_ijk->dcos_di, p_kji->dcos_dk );
+ rvec_Copy( p_ijk->dcos_dj, p_kji->dcos_dj );
+ rvec_Copy( p_ijk->dcos_dk, p_kji->dcos_di );
+
+ ++num_thb_intrs;
+ break;
+ }
}
- }
-
-
- /* 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);
- 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] );
-
- //CHANGE ORIGINAL
- if (BOA_jk <= 0) continue;
- //CHANGE ORIGINAL
-
-
- Calculate_Theta( pbond_ij->dvec, pbond_ij->d,
- pbond_jk->dvec, pbond_jk->d,
- &theta, &cos_theta );
-
- Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d,
- pbond_jk->dvec, pbond_jk->d,
- &(p_ijk->dcos_di), &(p_ijk->dcos_dj),
- &(p_ijk->dcos_dk) );
+ }
- p_ijk->thb = k;
- p_ijk->pthb = pk;
- p_ijk->theta = theta;
- sin_theta = SIN( theta );
- if ( sin_theta < 1.0e-5 )
- sin_theta = 1.0e-5;
+ /* 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);
+ 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] );
+//#ifdef _OPENMP
+// f_k = &(workspace->f_local[tid * system->N + k]);
+//#else
+ f_k = &(system->atoms[k].f);
+//#endif
+
+ //CHANGE ORIGINAL
+ if ( BOA_jk <= 0 )
+ {
+ continue;
+ }
+ //CHANGE ORIGINAL
- ++num_thb_intrs;
+ Calculate_Theta( pbond_ij->dvec, pbond_ij->d,
+ pbond_jk->dvec, pbond_jk->d,
+ &theta, &cos_theta );
+ Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d,
+ pbond_jk->dvec, pbond_jk->d,
+ &(p_ijk->dcos_di), &(p_ijk->dcos_dj),
+ &(p_ijk->dcos_dk) );
- if ( BOA_jk > 0.0 &&
- (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut)/*0*/)
- {
- r_jk = pbond_jk->d;
- thbh = &( system->reaxprm.thbp[type_i][type_j][type_k] );
- flag = 0;
+ p_ijk->thb = k;
+ p_ijk->pthb = pk;
+ p_ijk->theta = theta;
- /* if( workspace->orig_id[i] < workspace->orig_id[k] )
- fprintf( stdout, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
- workspace->orig_id[i], workspace->orig_id[j],
- workspace->orig_id[k], bo_ij->BO, bo_jk->BO, p_ijk->theta );
- else
- fprintf( stdout, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
- workspace->orig_id[k], workspace->orig_id[j],
- workspace->orig_id[i], bo_jk->BO, bo_ij->BO, p_ijk->theta ); */
+ sin_theta = SIN( theta );
+ if ( sin_theta < 1.0e-5 )
+ {
+ sin_theta = 1.0e-5;
+ }
+ ++num_thb_intrs;
- for ( cnt = 0; cnt < thbh->cnt; ++cnt )
+ if ( BOA_jk > 0.0 &&
+ (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut)/*0*/)
{
- // fprintf( out_control->eval,
- // "%6d%6d%6d -- exists in thbp\n", i+1, j+1, k+1 );
+ 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",
+ workspace->orig_id[i], workspace->orig_id[j],
+ workspace->orig_id[k], bo_ij->BO, bo_jk->BO, p_ijk->theta );
+ else
+ fprintf( stdout, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
+ workspace->orig_id[k], workspace->orig_id[j],
+ workspace->orig_id[i], bo_jk->BO, bo_ij->BO, p_ijk->theta ); */
- if ( fabs(thbh->prm[cnt].p_val1) > 0.001 )
+ for ( cnt = 0; cnt < thbh->cnt; ++cnt )
{
- thbp = &( thbh->prm[cnt] );
-
- /* ANGLE ENERGY */
- p_val1 = thbp->p_val1;
- p_val2 = thbp->p_val2;
- p_val4 = thbp->p_val4;
- p_val7 = thbp->p_val7;
- theta_00 = thbp->theta_00;
-
- exp3ij = EXP( -p_val3 * POW( BOA_ij, p_val4 ) );
- f7_ij = 1.0 - exp3ij;
- Cf7ij = p_val3 * p_val4 *
+ if ( FABS(thbh->prm[cnt].p_val1) > 0.001 )
+ {
+ thbp = &( thbh->prm[cnt] );
+
+ /* ANGLE ENERGY */
+ p_val1 = thbp->p_val1;
+ p_val2 = thbp->p_val2;
+ p_val4 = thbp->p_val4;
+ p_val7 = thbp->p_val7;
+ theta_00 = thbp->theta_00;
+
+ exp3ij = EXP( -p_val3 * POW( BOA_ij, p_val4 ) );
+ f7_ij = 1.0 - exp3ij;
+ Cf7ij = p_val3 * p_val4 *
POW( BOA_ij, p_val4 - 1.0 ) * exp3ij;
- exp3jk = EXP( -p_val3 * POW( BOA_jk, p_val4 ) );
- f7_jk = 1.0 - exp3jk;
- Cf7jk = p_val3 * p_val4 *
+ exp3jk = EXP( -p_val3 * POW( BOA_jk, p_val4 ) );
+ f7_jk = 1.0 - exp3jk;
+ Cf7jk = p_val3 * p_val4 *
POW( BOA_jk, p_val4 - 1.0 ) * exp3jk;
- expval7 = EXP( -p_val7 * workspace->Delta_boc[j] );
- trm8 = 1.0 + expval6 + expval7;
- f8_Dj = p_val5 - ( (p_val5 - 1.0) * (2.0 + expval6) / trm8 );
- Cf8j = ( (1.0 - p_val5) / SQR(trm8) ) *
- (p_val6 * expval6 * trm8 -
- (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
-
- theta_0 = 180.0 -
- theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
- theta_0 = DEG2RAD( theta_0 );
-
- expval2theta = EXP(-p_val2 * SQR(theta_0 - theta));
- if ( p_val1 >= 0 )
- expval12theta = p_val1 * (1.0 - expval2theta);
- else // To avoid linear Me-H-Me angles (6/6/06)
- expval12theta = p_val1 * -expval2theta;
-
- 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);
-
- Ctheta_0 = p_val10 * DEG2RAD(theta_00) *
- exp( -p_val10 * (2.0 - SBO2) );
-
- CEval5 = -CEval4 * Ctheta_0 * CSBO2;
- CEval6 = CEval5 * dSBO1;
- CEval7 = CEval5 * dSBO2;
- CEval8 = -CEval4 / sin_theta;
-
- data->E_Ang += e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
- /* END ANGLE ENERGY*/
-
-
- /* PENALTY ENERGY */
- p_pen1 = thbp->p_pen1;
- p_pen2 = system->reaxprm.gp.l[19];
- p_pen3 = system->reaxprm.gp.l[20];
- p_pen4 = system->reaxprm.gp.l[21];
-
- exp_pen2ij = EXP( -p_pen2 * SQR( BOA_ij - 2.0 ) );
- exp_pen2jk = EXP( -p_pen2 * SQR( BOA_jk - 2.0 ) );
- exp_pen3 = EXP( -p_pen3 * workspace->Delta[j] );
- exp_pen4 = EXP( p_pen4 * workspace->Delta[j] );
- trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
- f9_Dj = ( 2.0 + exp_pen3 ) / trm_pen34;
- Cf9j = (-p_pen3 * exp_pen3 * trm_pen34 -
- (2.0 + exp_pen3) * ( -p_pen3 * exp_pen3 +
- p_pen4 * exp_pen4 )) /
- SQR( trm_pen34 );
-
- data->E_Pen += e_pen =
- p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
-
- CEpen1 = e_pen * Cf9j / f9_Dj;
- temp = -2.0 * p_pen2 * e_pen;
- CEpen2 = temp * (BOA_ij - 2.0);
- CEpen3 = temp * (BOA_jk - 2.0);
- /* END PENALTY ENERGY */
-
-
- /* COALITION ENERGY */
- p_coa1 = thbp->p_coa1;
- p_coa2 = system->reaxprm.gp.l[2];
- p_coa3 = system->reaxprm.gp.l[38];
- p_coa4 = system->reaxprm.gp.l[30];
-
- exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
- data->E_Coa += e_coa =
- p_coa1 / (1. + exp_coa2) *
- EXP( -p_coa3 * SQR(total_bo[i] - BOA_ij) ) *
- EXP( -p_coa3 * SQR(total_bo[k] - BOA_jk) ) *
- EXP( -p_coa4 * SQR(BOA_ij - 1.5) ) *
- EXP( -p_coa4 * SQR(BOA_jk - 1.5) );
-
- CEcoa1 = -2 * p_coa4 * (BOA_ij - 1.5) * e_coa;
- CEcoa2 = -2 * p_coa4 * (BOA_jk - 1.5) * e_coa;
- CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1 + exp_coa2);
- CEcoa4 = -2 * p_coa3 * (total_bo[i] - BOA_ij) * e_coa;
- CEcoa5 = -2 * p_coa3 * (total_bo[k] - BOA_jk) * e_coa;
- /* END COALITION ENERGY */
-
- /* FORCES */
- bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
- bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
- workspace->CdDelta[j] += ((CEval3 + CEval7) +
- CEpen1 + CEcoa3);
- workspace->CdDelta[i] += CEcoa4;
- workspace->CdDelta[k] += CEcoa5;
-
- for ( t = start_j; t < end_j; ++t )
- {
- 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;
-
- // fprintf( out_control->eval, "%6d%12.8f\n",
- // workspace->orig_id[ bond_list[t].nbr ],
- // (CEval6 * pBOjt7) );
-
- bo_jt->Cdbo += (CEval6 * pBOjt7);
- bo_jt->Cdbopi += CEval5;
- bo_jt->Cdbopi2 += CEval5;
- }
+ expval7 = EXP( -p_val7 * workspace->Delta_boc[j] );
+ trm8 = 1.0 + expval6 + expval7;
+ f8_Dj = p_val5 - ( (p_val5 - 1.0) * (2.0 + expval6) / trm8 );
+ Cf8j = ( (1.0 - p_val5) / SQR(trm8) ) *
+ (p_val6 * expval6 * trm8 -
+ (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
+
+ theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
+ theta_0 = DEG2RAD( theta_0 );
+
+ expval2theta = EXP(-p_val2 * SQR(theta_0 - theta));
+ if ( p_val1 >= 0 )
+ {
+ expval12theta = p_val1 * (1.0 - expval2theta);
+ }
+ /* To avoid linear Me-H-Me angles (6/6/06) */
+ else
+ {
+ expval12theta = p_val1 * -expval2theta;
+ }
+
+ 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);
+
+ Ctheta_0 = p_val10 * DEG2RAD(theta_00) *
+ EXP( -p_val10 * (2.0 - SBO2) );
+
+ CEval5 = -CEval4 * Ctheta_0 * CSBO2;
+ CEval6 = CEval5 * dSBO1;
+ CEval7 = CEval5 * dSBO2;
+ CEval8 = -CEval4 / sin_theta;
+
+ e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
+ e_ang_total += e_ang;
+ /* END ANGLE ENERGY*/
+
+ /* PENALTY ENERGY */
+ p_pen1 = thbp->p_pen1;
+
+ exp_pen2ij = EXP( -p_pen2 * SQR( BOA_ij - 2.0 ) );
+ exp_pen2jk = EXP( -p_pen2 * SQR( BOA_jk - 2.0 ) );
+ exp_pen3 = EXP( -p_pen3 * workspace->Delta[j] );
+ exp_pen4 = EXP( p_pen4 * workspace->Delta[j] );
+ trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
+ f9_Dj = ( 2.0 + exp_pen3 ) / trm_pen34;
+ Cf9j = (-p_pen3 * exp_pen3 * trm_pen34 -
+ (2.0 + exp_pen3) * ( -p_pen3 * exp_pen3 +
+ p_pen4 * exp_pen4 )) / SQR( trm_pen34 );
+
+ e_pen = p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
+ e_pen_total += e_pen;
+
+ CEpen1 = e_pen * Cf9j / f9_Dj;
+ temp = -2.0 * p_pen2 * e_pen;
+ CEpen2 = temp * (BOA_ij - 2.0);
+ CEpen3 = temp * (BOA_jk - 2.0);
+ /* END PENALTY ENERGY */
+
+ /* COALITION ENERGY */
+ p_coa1 = thbp->p_coa1;
+
+ exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
+ e_coa = p_coa1 / (1. + exp_coa2) *
+ EXP( -p_coa3 * SQR(total_bo[i] - BOA_ij) ) *
+ EXP( -p_coa3 * SQR(total_bo[k] - BOA_jk) ) *
+ EXP( -p_coa4 * SQR(BOA_ij - 1.5) ) *
+ EXP( -p_coa4 * SQR(BOA_jk - 1.5) );
+ e_coa_total += e_coa;
+
+ CEcoa1 = -2 * p_coa4 * (BOA_ij - 1.5) * e_coa;
+ CEcoa2 = -2 * p_coa4 * (BOA_jk - 1.5) * e_coa;
+ CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1 + exp_coa2);
+ CEcoa4 = -2 * p_coa3 * (total_bo[i] - BOA_ij) * e_coa;
+ CEcoa5 = -2 * p_coa3 * (total_bo[k] - BOA_jk) * e_coa;
+ /* END COALITION ENERGY */
+
+ /* FORCES */
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ workspace->CdDelta[i] += CEcoa4;
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ workspace->CdDelta[k] += CEcoa5;
+
+ for ( t = start_j; t < end_j; ++t )
+ {
+ 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;
+
+ // fprintf( out_control->eval, "%6d%12.8f\n",
+ // workspace->orig_id[ bond_list[t].nbr ],
+ // (CEval6 * pBOjt7) );
+
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ bo_jt->Cdbo += (CEval6 * pBOjt7);
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ bo_jt->Cdbopi += CEval5;
+#ifdef _OPENMP
+// #pragma omp atomic
+#endif
+ bo_jt->Cdbopi2 += CEval5;
+ }
+
+
+ if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT)
+ {
+ rvec_ScaledAdd( *f_i, CEval8, p_ijk->dcos_di );
+ rvec_ScaledAdd( *f_j, CEval8, p_ijk->dcos_dj );
+ rvec_ScaledAdd( *f_k, CEval8, p_ijk->dcos_dk );
+ }
+ else
+ {
+ /* terms not related to bond order derivatives
+ are added directly into
+ forces and pressure vector/tensor */
+ rvec_Scale( force, CEval8, p_ijk->dcos_di );
+ rvec_Add( *f_i, force );
+ rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
+#ifdef _OPENMP
+// #pragma omp critical (Three_Body_Interactions_ext_press)
+#endif
+ {
+ rvec_Add( data->ext_press, ext_press );
+ }
+ rvec_ScaledAdd( *f_j, CEval8, p_ijk->dcos_dj );
- if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT)
- {
- rvec_ScaledAdd( system->atoms[i].f, CEval8, p_ijk->dcos_di );
- rvec_ScaledAdd( system->atoms[j].f, CEval8, p_ijk->dcos_dj );
- rvec_ScaledAdd( system->atoms[k].f, CEval8, p_ijk->dcos_dk );
- }
- else
- {
- /* terms not related to bond order derivatives
- are added directly into
- forces and pressure vector/tensor */
- rvec_Scale( force, CEval8, p_ijk->dcos_di );
- rvec_Add( system->atoms[i].f, force );
- rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
- rvec_Add( data->ext_press, ext_press );
-
- rvec_ScaledAdd( system->atoms[j].f,
- CEval8, p_ijk->dcos_dj );
-
- rvec_Scale( force, CEval8, p_ijk->dcos_dk );
- rvec_Add( system->atoms[k].f, force );
- rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
- rvec_Add( data->ext_press, ext_press );
-
-
- /* This part is for a fully-flexible box */
- /* rvec_OuterProduct( temp_rtensor,
- p_ijk->dcos_di, system->atoms[i].x );
- rtensor_Scale( total_rtensor, +CEval8, temp_rtensor );
-
- rvec_OuterProduct( temp_rtensor,
- p_ijk->dcos_dj, system->atoms[j].x );
- rtensor_ScaledAdd(total_rtensor, CEval8, temp_rtensor);
-
- rvec_OuterProduct( temp_rtensor,
- p_ijk->dcos_dk, system->atoms[k].x );
- rtensor_ScaledAdd(total_rtensor, CEval8, temp_rtensor);
-
- if( pbond_ij->imaginary || pbond_jk->imaginary )
- rtensor_ScaledAdd( data->flex_bar.P,
- -1.0, total_rtensor );
- else
- rtensor_Add( data->flex_bar.P, total_rtensor ); */
- }
+ rvec_Scale( force, CEval8, p_ijk->dcos_dk );
+ rvec_Add( *f_k, force );
+ rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
+#ifdef _OPENMP
+// #pragma omp critical (Three_Body_Interactions_ext_press)
+#endif
+ {
+ rvec_Add( data->ext_press, ext_press );
+ }
+
+ /* This part is for a fully-flexible box */
+ /* rvec_OuterProduct( temp_rtensor,
+ p_ijk->dcos_di, system->atoms[i].x );
+ rtensor_Scale( total_rtensor, +CEval8, temp_rtensor );
+
+ rvec_OuterProduct( temp_rtensor,
+ p_ijk->dcos_dj, system->atoms[j].x );
+ rtensor_ScaledAdd(total_rtensor, CEval8, temp_rtensor);
+
+ rvec_OuterProduct( temp_rtensor,
+ p_ijk->dcos_dk, system->atoms[k].x );
+ rtensor_ScaledAdd(total_rtensor, CEval8, temp_rtensor);
+
+ if( pbond_ij->imaginary || pbond_jk->imaginary )
+ rtensor_ScaledAdd( data->flex_bar.P,
+ -1.0, total_rtensor );
+ else
+ rtensor_Add( data->flex_bar.P, total_rtensor ); */
+ }
#ifdef TEST_ENERGY
- fprintf( out_control->eval,
- //"%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e",
- "%6d%6d%6d%23.15e%23.15e%23.15e\n",
- i + 1, j + 1, k + 1,
- //workspace->orig_id[i]+1,
- //workspace->orig_id[j]+1,
- //workspace->orig_id[k]+1,
- //workspace->Delta_boc[j],
- RAD2DEG(theta), /*BOA_ij, BOA_jk, */
- e_ang, data->E_Ang );
-
- /*fprintf( out_control->eval,
- "%23.15e%23.15e%23.15e%23.15e",
- p_val3, p_val4, BOA_ij, BOA_jk );
- fprintf( out_control->eval,
- "%23.15e%23.15e%23.15e%23.15e",
- f7_ij, f7_jk, f8_Dj, expval12theta );
- fprintf( out_control->eval,
- "%23.15e%23.15e%23.15e%23.15e%23.15e\n",
- CEval1, CEval2, CEval3, CEval4, CEval5
- //CEval6, CEval7, CEval8 );*/
-
- /*fprintf( out_control->eval,
- "%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e\n",
- -p_ijk->dcos_di[0]/sin_theta,
- -p_ijk->dcos_di[1]/sin_theta,
- -p_ijk->dcos_di[2]/sin_theta,
- -p_ijk->dcos_dj[0]/sin_theta,
- -p_ijk->dcos_dj[1]/sin_theta,
- -p_ijk->dcos_dj[2]/sin_theta,
- -p_ijk->dcos_dk[0]/sin_theta,
- -p_ijk->dcos_dk[1]/sin_theta,
- -p_ijk->dcos_dk[2]/sin_theta );*/
-
- /* fprintf( out_control->epen,
- "%23.15e%23.15e%23.15e\n",
- CEpen1, CEpen2, CEpen3 );
- fprintf( out_control->epen,
- "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
- workspace->orig_id[i], workspace->orig_id[j],
- workspace->orig_id[k], RAD2DEG(theta),
- BOA_ij, BOA_jk, e_pen, data->E_Pen ); */
-
- fprintf( out_control->ecoa,
- "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
- workspace->orig_id[i],
- workspace->orig_id[j],
- workspace->orig_id[k],
- RAD2DEG(theta), BOA_ij, BOA_jk,
- e_coa, data->E_Coa );
+ fprintf( out_control->eval,
+ //"%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e",
+ "%6d%6d%6d%23.15e%23.15e%23.15e\n",
+ i + 1, j + 1, k + 1,
+ //workspace->orig_id[i]+1,
+ //workspace->orig_id[j]+1,
+ //workspace->orig_id[k]+1,
+ //workspace->Delta_boc[j],
+ RAD2DEG(theta), /*BOA_ij, BOA_jk, */
+ e_ang, data->E_Ang );
+
+ /*fprintf( out_control->eval,
+ "%23.15e%23.15e%23.15e%23.15e",
+ p_val3, p_val4, BOA_ij, BOA_jk );
+ fprintf( out_control->eval,
+ "%23.15e%23.15e%23.15e%23.15e",
+ f7_ij, f7_jk, f8_Dj, expval12theta );
+ fprintf( out_control->eval,
+ "%23.15e%23.15e%23.15e%23.15e%23.15e\n",
+ CEval1, CEval2, CEval3, CEval4, CEval5
+ //CEval6, CEval7, CEval8 );*/
+
+ /*fprintf( out_control->eval,
+ "%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e%23.15e\n",
+ -p_ijk->dcos_di[0]/sin_theta,
+ -p_ijk->dcos_di[1]/sin_theta,
+ -p_ijk->dcos_di[2]/sin_theta,
+ -p_ijk->dcos_dj[0]/sin_theta,
+ -p_ijk->dcos_dj[1]/sin_theta,
+ -p_ijk->dcos_dj[2]/sin_theta,
+ -p_ijk->dcos_dk[0]/sin_theta,
+ -p_ijk->dcos_dk[1]/sin_theta,
+ -p_ijk->dcos_dk[2]/sin_theta );*/
+
+ /* fprintf( out_control->epen,
+ "%23.15e%23.15e%23.15e\n",
+ CEpen1, CEpen2, CEpen3 );
+ fprintf( out_control->epen,
+ "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
+ workspace->orig_id[i], workspace->orig_id[j],
+ workspace->orig_id[k], RAD2DEG(theta),
+ BOA_ij, BOA_jk, e_pen, data->E_Pen ); */
+
+ fprintf( out_control->ecoa,
+ "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
+ workspace->orig_id[i],
+ workspace->orig_id[j],
+ workspace->orig_id[k],
+ RAD2DEG(theta), BOA_ij, BOA_jk,
+ e_coa, data->E_Coa );
#endif
-#ifdef TEST_FORCES /* angle forces */
- Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
- Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
- Add_dDelta( system, lists,
- j, CEval3 + CEval7, workspace->f_ang );
-
- for ( t = start_j; t < end_j; ++t )
- {
- 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;
-
- Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
- workspace->f_ang );
- Add_dBOpinpi2( system, lists, j, t,
- CEval5, CEval5,
- workspace->f_ang, workspace->f_ang );
- }
-
- rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
- rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
- rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
- /* end angle forces */
-
- /* penalty forces */
- Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
- Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
- Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
- /* end penalty forces */
-
- /* coalition forces */
- Add_dBO( system, lists,
- j, pi, CEcoa1 - CEcoa4, workspace->f_coa );
- Add_dBO( system, lists,
- j, pk, CEcoa2 - CEcoa5, workspace->f_coa );
- Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
- Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
- Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
- /* end coalition forces */
+#ifdef TEST_FORCES
+ /* angle forces */
+ Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
+ Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
+ Add_dDelta( system, lists,
+ j, CEval3 + CEval7, workspace->f_ang );
+
+ for ( t = start_j; t < end_j; ++t )
+ {
+ 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;
+
+ Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
+ workspace->f_ang );
+ Add_dBOpinpi2( system, lists, j, t,
+ CEval5, CEval5,
+ workspace->f_ang, workspace->f_ang );
+ }
+
+ rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
+ rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
+ rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
+ /* end angle forces */
+
+ /* penalty forces */
+ Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
+ Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
+ Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
+ /* end penalty forces */
+
+ /* coalition forces */
+ Add_dBO( system, lists,
+ j, pi, CEcoa1 - CEcoa4, workspace->f_coa );
+ Add_dBO( system, lists,
+ j, pk, CEcoa2 - CEcoa5, workspace->f_coa );
+ Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
+ Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
+ Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
+ /* end coalition forces */
#endif
+ }
}
}
}
}
- }
- Set_End_Index(pi, num_thb_intrs, thb_intrs );
+ Set_End_Index( pi, num_thb_intrs, thb_intrs );
+ }
}
}
+ data->E_Ang += e_ang_total;
+ data->E_Pen += e_pen_total;
+ data->E_Coa += e_coa_total;
if ( num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE )
{
@@ -594,217 +665,272 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
}
-
void Hydrogen_Bonds( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
{
- int i, j, k, pi, pk, itr, top;
- int type_i, type_j, type_k;
- int start_j, end_j, hb_start_j, hb_end_j;
- int hblist[MAX_BONDS];
- int num_hb_intrs = 0;
- real r_ij, r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
- real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
- rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
- rvec dvec_jk, force, ext_press;
- ivec rel_jk;
- // rtensor temp_rtensor, total_rtensor;
- hbond_parameters *hbp;
- bond_order_data *bo_ij;
- bond_data *pbond_ij;
- far_neighbor_data *nbr_jk;
- list *bonds, *hbonds;
- bond_data *bond_list;
- hbond_data *hbond_list;
+ real e_hb_total;
- bonds = (*lists) + BONDS;
- bond_list = bonds->select.bond_list;
+ e_hb_total = 0.0;
- hbonds = (*lists) + HBONDS;
- hbond_list = hbonds->select.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.
- Hydrogen bond is between j and k.
- so in this function i->X, j->H, k->Z when we map
- variables onto the ones in the handout.*/
- for ( j = 0; j < system->N; ++j )
+#ifdef _OPENMP
+ #pragma omp parallel default(shared) reduction(+: e_hb_total)
+#endif
{
- /* j must be H */
- if ( system->reaxprm.sbp[system->atoms[j].type].p_hbond == 1 )
- {
- /*set j's variables */
- type_j = system->atoms[j].type;
- start_j = Start_Index( j, bonds );
- end_j = End_Index( j, bonds );
- hb_start_j = Start_Index( workspace->hbond_index[j], hbonds );
- hb_end_j = End_Index( workspace->hbond_index[j], hbonds );
+ int i, j, k, pi, pk, itr, top;
+ int type_i, type_j, type_k;
+ int start_j, end_j, hb_start_j, hb_end_j;
+ int hblist[MAX_BONDS];
+#ifdef TEST_FORCES
+ int num_hb_intrs;
+#endif
+ real r_ij, r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
+ real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
+ rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
+ rvec dvec_jk, force, ext_press;
+ ivec rel_jk;
+ //rtensor temp_rtensor, total_rtensor;
+ hbond_parameters *hbp;
+ bond_order_data *bo_ij;
+ bond_data *pbond_ij;
+ far_neighbor_data *nbr_jk;
+ list *bonds, *hbonds;
+ bond_data *bond_list;
+ hbond_data *hbond_list;
+ rvec *f_i, *f_j, *f_k;
+#ifdef _OPENMP
+ int tid = omp_get_thread_num( );
+#endif
- top = 0;
- for ( pi = start_j; pi < end_j; ++pi )
+#ifdef TEST_FORCES
+ num_hb_intrs = 0;
+#endif
+ bonds = (*lists) + BONDS;
+ bond_list = bonds->select.bond_list;
+ hbonds = (*lists) + HBONDS;
+ hbond_list = hbonds->select.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.
+ Hydrogen bond is between j and k.
+ so in this function i->X, j->H, k->Z when we map
+ variables onto the ones in the handout.*/
+#ifdef _OPENMP
+ #pragma omp for schedule(guided)
+#endif
+ for ( j = 0; j < system->N; ++j )
+ {
+ /* j must be H */
+ if ( system->reaxprm.sbp[system->atoms[j].type].p_hbond == 1 )
{
- pbond_ij = &( bond_list[pi] );
- i = pbond_ij->nbr;
- bo_ij = &(pbond_ij->bo_data);
- type_i = system->atoms[i].type;
+ /* set j's variables */
+ type_j = system->atoms[j].type;
+ start_j = Start_Index( j, bonds );
+ end_j = End_Index( j, bonds );
+ 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]);
+#else
+ f_j = &(system->atoms[j].f);
+#endif
- if ( system->reaxprm.sbp[type_i].p_hbond == 2 &&
- bo_ij->BO >= HB_THRESHOLD )
+ top = 0;
+ for ( pi = start_j; pi < end_j; ++pi )
{
- hblist[top++] = pi;
- }
- }
+ pbond_ij = &( bond_list[pi] );
+ i = pbond_ij->nbr;
+ bo_ij = &(pbond_ij->bo_data);
+ type_i = system->atoms[i].type;
- // fprintf( stderr, "j: %d, top: %d, hb_start_j: %d, hb_end_j:%d\n",
- // j, top, hb_start_j, hb_end_j );
+ if ( system->reaxprm.sbp[type_i].p_hbond == 2 &&
+ bo_ij->BO >= HB_THRESHOLD )
+ {
+ hblist[top++] = pi;
+ }
+ }
- for ( pk = hb_start_j; pk < hb_end_j; ++pk )
- {
- /* set k's varibles */
- k = hbond_list[pk].nbr;
- type_k = system->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 );
-
- for ( itr = 0; itr < top; ++itr )
+ for ( pk = hb_start_j; pk < hb_end_j; ++pk )
{
- pi = hblist[itr];
- pbond_ij = &( bond_list[pi] );
- i = pbond_ij->nbr;
+ /* set k's varibles */
+ k = hbond_list[pk].nbr;
+ type_k = system->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 );
+#ifdef _OPENMP
+ f_k = &(workspace->f_local[tid * system->N + k]);
+#else
+ f_k = &(system->atoms[k].f);
+#endif
- if ( i != k )
+ for ( itr = 0; itr < top; ++itr )
{
- 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 ]);
- ++num_hb_intrs;
+ pi = hblist[itr];
+ pbond_ij = &( bond_list[pi] );
+ i = pbond_ij->nbr;
- Calculate_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
- &theta, &cos_theta );
- /* the derivative of cos(theta) */
- Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
- &dcos_theta_di, &dcos_theta_dj, &dcos_theta_dk );
-
- /* hydrogen bond energy*/
- sin_theta2 = SIN( theta / 2.0 );
- sin_xhz4 = SQR(sin_theta2);
- sin_xhz4 *= sin_xhz4;
- cos_xhz1 = ( 1.0 - cos_theta );
- exp_hb2 = EXP( -hbp->p_hb2 * bo_ij->BO );
- exp_hb3 = EXP( -hbp->p_hb3 * ( hbp->r0_hb / r_jk +
- r_jk / hbp->r0_hb - 2.0 ) );
-
- e_hb = hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
- data->E_HB += e_hb;
-
- CEhb1 = hbp->p_hb1 * hbp->p_hb2 * exp_hb2 * exp_hb3 * sin_xhz4;
- CEhb2 = -hbp->p_hb1 / 2.0 * (1.0 - exp_hb2) * exp_hb3 * cos_xhz1;
- CEhb3 = -hbp->p_hb3 * e_hb * (-hbp->r0_hb / SQR(r_jk) +
- 1.0 / hbp->r0_hb);
-
- /* hydrogen bond forces */
- bo_ij->Cdbo += CEhb1; // dbo term
-
- if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT)
+ if ( i != k )
{
- rvec_ScaledAdd( system->atoms[i].f,
- +CEhb2, dcos_theta_di ); //dcos terms
- rvec_ScaledAdd( system->atoms[j].f,
- +CEhb2, dcos_theta_dj );
- rvec_ScaledAdd( system->atoms[k].f,
- +CEhb2, dcos_theta_dk );
- //dr terms
- rvec_ScaledAdd( system->atoms[j].f, -CEhb3 / r_jk, dvec_jk );
- rvec_ScaledAdd( system->atoms[k].f, +CEhb3 / r_jk, dvec_jk );
- }
- else
- {
- /* for pressure coupling, terms that are not related
- to bond order derivatives are added directly into
- pressure vector/tensor */
- rvec_Scale( force, +CEhb2, dcos_theta_di ); // dcos terms
- rvec_Add( system->atoms[i].f, force );
- rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
- rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
-
- rvec_ScaledAdd( system->atoms[j].f, +CEhb2, dcos_theta_dj );
-
- ivec_Scale( rel_jk, hbond_list[pk].scl, nbr_jk->rel_box );
- rvec_Scale( force, +CEhb2, dcos_theta_dk );
- rvec_Add( system->atoms[k].f, force );
- rvec_iMultiply( ext_press, rel_jk, force );
- rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
-
- //dr terms
- rvec_ScaledAdd( system->atoms[j].f, -CEhb3 / r_jk, dvec_jk );
-
- rvec_Scale( force, CEhb3 / r_jk, dvec_jk );
- rvec_Add( system->atoms[k].f, force );
- rvec_iMultiply( ext_press, rel_jk, force );
- rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
-
- /* This part is intended for a fully-flexible box */
- /* rvec_OuterProduct( temp_rtensor,
- dcos_theta_di, system->atoms[i].x );
- rtensor_Scale( total_rtensor, -CEhb2, temp_rtensor );
-
- rvec_ScaledSum( temp_rvec, -CEhb2, dcos_theta_dj,
- -CEhb3/r_jk, pbond_jk->dvec );
- rvec_OuterProduct( temp_rtensor,
- temp_rvec, system->atoms[j].x );
- rtensor_Add( total_rtensor, temp_rtensor );
-
- rvec_ScaledSum( temp_rvec, -CEhb2, dcos_theta_dk,
- +CEhb3/r_jk, pbond_jk->dvec );
- rvec_OuterProduct( temp_rtensor,
- temp_rvec, system->atoms[k].x );
- rtensor_Add( total_rtensor, temp_rtensor );
-
- if( pbond_ij->imaginary || pbond_jk->imaginary )
- rtensor_ScaledAdd( data->flex_bar.P, -1.0, total_rtensor );
- else
- rtensor_Add( data->flex_bar.P, total_rtensor ); */
- }
+ 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 ]);
+#ifdef _OPENMP
+ f_i = &(workspace->f_local[tid * system->N + i]);
+#else
+ f_i = &(system->atoms[i].f);
+#endif
-#ifdef TEST_ENERGY
- /*fprintf( out_control->ehb,
- "%23.15e%23.15e%23.15e\n%23.15e%23.15e%23.15e\n%23.15e%23.15e%23.15e\n",
- dcos_theta_di[0], dcos_theta_di[1], dcos_theta_di[2],
- dcos_theta_dj[0], dcos_theta_dj[1], dcos_theta_dj[2],
- dcos_theta_dk[0], dcos_theta_dk[1], dcos_theta_dk[2]);
- fprintf( out_control->ehb, "%23.15e%23.15e%23.15e\n",
- CEhb1, CEhb2, CEhb3 ); */
- fprintf( stderr, //out_control->ehb,
- "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
- workspace->orig_id[i],
- workspace->orig_id[j],
- workspace->orig_id[k],
- r_jk, theta, bo_ij->BO, e_hb, data->E_HB );
+#ifdef TEST_FORCES
+ ++num_hb_intrs;
+#endif
+ Calculate_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
+ &theta, &cos_theta );
+ /* the derivative of cos(theta) */
+ Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
+ &dcos_theta_di, &dcos_theta_dj, &dcos_theta_dk );
+
+ /* hydrogen bond energy */
+ sin_theta2 = SIN( theta / 2.0 );
+ sin_xhz4 = SQR( sin_theta2 );
+ sin_xhz4 *= sin_xhz4;
+ cos_xhz1 = ( 1.0 - cos_theta );
+ exp_hb2 = EXP( -hbp->p_hb2 * bo_ij->BO );
+ exp_hb3 = EXP( -hbp->p_hb3 * ( hbp->r0_hb / r_jk +
+ r_jk / hbp->r0_hb - 2.0 ) );
+
+ e_hb = hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
+ e_hb_total += e_hb;
+
+ CEhb1 = hbp->p_hb1 * hbp->p_hb2 * exp_hb2 * exp_hb3 * sin_xhz4;
+ CEhb2 = -hbp->p_hb1 / 2.0 * (1.0 - exp_hb2) * exp_hb3 * cos_xhz1;
+ CEhb3 = -hbp->p_hb3 * e_hb * (-hbp->r0_hb / SQR( r_jk ) +
+ 1.0 / hbp->r0_hb);
+
+ /* hydrogen bond forces */
+ /* dbo term,
+ * note: safe to update across threads as this points
+ * to the bond_order_data struct inside atom j's list,
+ * and threads are partitioned across all j's */
+#ifdef _OPENMP
+ #pragma omp atomic
+#endif
+ bo_ij->Cdbo += CEhb1;
+
+ if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT)
+ {
+ /* dcos terms */
+ rvec_ScaledAdd( *f_i, +CEhb2, dcos_theta_di );
+ rvec_ScaledAdd( *f_j, +CEhb2, dcos_theta_dj );
+ rvec_ScaledAdd( *f_k, +CEhb2, dcos_theta_dk );
+
+ /* dr terms */
+ rvec_ScaledAdd( *f_j, -CEhb3 / r_jk, dvec_jk );
+ rvec_ScaledAdd( *f_k, +CEhb3 / r_jk, dvec_jk );
+ }
+ else
+ {
+ /* for pressure coupling, terms that are not related
+ to bond order derivatives are added directly into
+ pressure vector/tensor */
+
+ /* dcos terms */
+ rvec_Scale( force, +CEhb2, dcos_theta_di );
+ rvec_Add( *f_i, force );
+ rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
+#ifdef _OPENMP
+ #pragma omp critical (Hydrogen_Bonds_ext_press)
+#endif
+ {
+ rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
+ }
+
+ rvec_ScaledAdd( *f_j, +CEhb2, dcos_theta_dj );
+
+ ivec_Scale( rel_jk, hbond_list[pk].scl, nbr_jk->rel_box );
+ rvec_Scale( force, +CEhb2, dcos_theta_dk );
+ rvec_Add( *f_k, force );
+ rvec_iMultiply( ext_press, rel_jk, force );
+#ifdef _OPENMP
+ #pragma omp critical (Hydrogen_Bonds_ext_press)
+#endif
+ {
+ rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
+ }
+
+ /* dr terms */
+ rvec_ScaledAdd( *f_j, -CEhb3 / r_jk, dvec_jk );
+
+ rvec_Scale( force, CEhb3 / r_jk, dvec_jk );
+ rvec_Add( *f_k, force );
+ rvec_iMultiply( ext_press, rel_jk, force );
+#ifdef _OPENMP
+ #pragma omp critical (Hydrogen_Bonds_ext_press)
+#endif
+ {
+ rvec_ScaledAdd( data->ext_press, 1.0, ext_press );
+ }
+
+ /* This part is intended for a fully-flexible box */
+ /* rvec_OuterProduct( temp_rtensor,
+ dcos_theta_di, system->atoms[i].x );
+ rtensor_Scale( total_rtensor, -CEhb2, temp_rtensor );
+
+ rvec_ScaledSum( temp_rvec, -CEhb2, dcos_theta_dj,
+ -CEhb3/r_jk, pbond_jk->dvec );
+ rvec_OuterProduct( temp_rtensor,
+ temp_rvec, system->atoms[j].x );
+ rtensor_Add( total_rtensor, temp_rtensor );
+
+ rvec_ScaledSum( temp_rvec, -CEhb2, dcos_theta_dk,
+ +CEhb3/r_jk, pbond_jk->dvec );
+ rvec_OuterProduct( temp_rtensor,
+ temp_rvec, system->atoms[k].x );
+ rtensor_Add( total_rtensor, temp_rtensor );
+
+ if( pbond_ij->imaginary || pbond_jk->imaginary )
+ rtensor_ScaledAdd( data->flex_bar.P, -1.0, total_rtensor );
+ else
+ rtensor_Add( data->flex_bar.P, total_rtensor ); */
+ }
+
+#ifdef TEST_ENERGY
+ /*fprintf( out_control->ehb,
+ "%23.15e%23.15e%23.15e\n%23.15e%23.15e%23.15e\n%23.15e%23.15e%23.15e\n",
+ dcos_theta_di[0], dcos_theta_di[1], dcos_theta_di[2],
+ dcos_theta_dj[0], dcos_theta_dj[1], dcos_theta_dj[2],
+ dcos_theta_dk[0], dcos_theta_dk[1], dcos_theta_dk[2]);
+ fprintf( out_control->ehb, "%23.15e%23.15e%23.15e\n",
+ CEhb1, CEhb2, CEhb3 ); */
+ fprintf( stderr, //out_control->ehb,
+ "%6d%6d%6d%23.15e%23.15e%23.15e%23.15e%23.15e\n",
+ workspace->orig_id[i],
+ workspace->orig_id[j],
+ workspace->orig_id[k],
+ r_jk, theta, bo_ij->BO, e_hb, data->E_HB );
#endif
#ifdef TEST_FORCES
- // dbo term
- Add_dBO( system, lists, j, pi, +CEhb1, workspace->f_hb );
- // dcos terms
- rvec_ScaledAdd( workspace->f_hb[i], +CEhb2, dcos_theta_di );
- rvec_ScaledAdd( workspace->f_hb[j], +CEhb2, dcos_theta_dj );
- rvec_ScaledAdd( workspace->f_hb[k], +CEhb2, dcos_theta_dk );
- // dr terms
- rvec_ScaledAdd( workspace->f_hb[j], -CEhb3 / r_jk, dvec_jk );
- rvec_ScaledAdd( workspace->f_hb[k], +CEhb3 / r_jk, dvec_jk );
+ /* dbo term */
+ Add_dBO( system, lists, j, pi, +CEhb1, workspace->f_hb );
+ /* dcos terms */
+ rvec_ScaledAdd( workspace->f_hb[i], +CEhb2, dcos_theta_di );
+ rvec_ScaledAdd( workspace->f_hb[j], +CEhb2, dcos_theta_dj );
+ rvec_ScaledAdd( workspace->f_hb[k], +CEhb2, dcos_theta_dk );
+ /* dr terms */
+ rvec_ScaledAdd( workspace->f_hb[j], -CEhb3 / r_jk, dvec_jk );
+ rvec_ScaledAdd( workspace->f_hb[k], +CEhb3 / r_jk, dvec_jk );
#endif
+ }
}
}
}
}
}
- /* fprintf( stderr, "hydbonds: ext_press (%23.15e %23.15e %23.15e)\n",
- data->ext_press[0], data->ext_press[1], data->ext_press[2] ); */
+ data->E_HB += e_hb_total;
#ifdef TEST_FORCES
fprintf( stderr, "Number of hydrogen bonds: %d\n", num_hb_intrs );
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index cc1da28a..245a4b76 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -383,6 +383,21 @@ int Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
}
+void Deallocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
+{
+ int i;
+
+ for ( i = 0; i < MAX_TOKENS; i++ )
+ {
+ free( (*tokens)[i] );
+ }
+
+ free( *line );
+ free( *backup );
+ free( *tokens );
+}
+
+
int Tokenize( char* s, char*** tok )
{
char test[MAX_LINE];
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index aec68895..38e6ba0c 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -57,6 +57,7 @@ int Get_Atom_Type( reax_interaction*, char* );
char *Get_Element( reax_system*, int );
char *Get_Atom_Name( reax_system*, int );
int Allocate_Tokenizer_Space( char**, char**, char*** );
+void Deallocate_Tokenizer_Space( char **, char **, char *** );
int Tokenize( char*, char*** );
/* from lammps */
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index 73f25e7c..d5379777 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "two_body_interactions.h"
+
#include "bond_orders.h"
#include "list.h"
#include "lookup.h"
@@ -30,6 +31,7 @@ void Bond_Energy( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
{
+ int i;
real gp3, gp4, gp7, gp10, gp37, ebond_total;
list *bonds;
@@ -41,9 +43,11 @@ void Bond_Energy( reax_system *system, control_params *control,
gp37 = (int) system->reaxprm.gp.l[37];
ebond_total = 0.0;
- #pragma omp parallel default(shared) reduction(+: ebond_total)
+#ifdef _OPENMP
+// #pragma omp parallel default(shared) reduction(+: ebond_total)
+#endif
{
- int i, j, pj;
+ int j, pj;
int start_i, end_i;
int type_i, type_j;
real ebond, pow_BOs_be2, exp_be12, CEbo;
@@ -53,13 +57,16 @@ void Bond_Energy( reax_system *system, control_params *control,
two_body_parameters *twbp;
bond_order_data *bo_ij;
- #pragma omp for schedule(guided)
+#ifdef _OPENMP
+// #pragma omp for schedule(guided)
+#endif
for ( i = 0; i < system->N; ++i )
{
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
- //fprintf( stderr, "i=%d start=%d end=%d\n", i, start_i, end_i );
+
for ( pj = start_i; pj < end_i; ++pj )
+ {
if ( i < bonds->select.bond_list[pj].nbr )
{
/* set the pointers */
@@ -81,7 +88,6 @@ void Bond_Energy( reax_system *system, control_params *control,
ebond = -twbp->De_s * bo_ij->BO_s * exp_be12
- twbp->De_p * bo_ij->BO_pi
- twbp->De_pp * bo_ij->BO_pi2;
-
ebond_total += ebond;
/* calculate derivatives of Bond Orders */
@@ -93,10 +99,7 @@ void Bond_Energy( reax_system *system, control_params *control,
fprintf( out_control->ebond, "%6d%6d%24.15e%24.15e\n",
workspace->orig_id[i], workspace->orig_id[j],
// i+1, j+1,
- bo_ij->BO, ebond/*, data->E_BE*/ );
- /* fprintf( out_control->ebond, "%6d%6d%12.6f%12.6f%12.6f\n",
- workspace->orig_id[i], workspace->orig_id[j],
- CEbo, -twbp->De_p, -twbp->De_pp );*/
+ bo_ij->BO, ebond );
#endif
#ifdef TEST_FORCES
@@ -113,7 +116,7 @@ void Bond_Energy( reax_system *system, control_params *control,
(sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990) ||
(sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) )
{
- // ba = SQR(bo_ij->BO - 2.50);
+ //ba = SQR(bo_ij->BO - 2.50);
exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.50) );
//oboa=abo(j1)-boa;
//obob=abo(j2)-boa;
@@ -129,11 +132,11 @@ void Bond_Energy( reax_system *system, control_params *control,
ebond_total += estriph;
decobdbo = gp10 * exphu * hulpov * (exphua1 + exphub1) *
- ( gp3 - 2.0 * gp7 * (bo_ij->BO - 2.50) );
+ ( gp3 - 2.0 * gp7 * (bo_ij->BO - 2.50) );
decobdboua = -gp10 * exphu * hulpov *
- (gp3 * exphua1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
+ (gp3 * exphua1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
decobdboub = -gp10 * exphu * hulpov *
- (gp3 * exphub1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
+ (gp3 * exphub1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
bo_ij->Cdbo += decobdbo;
workspace->CdDelta[i] += decobdboua;
@@ -155,6 +158,7 @@ void Bond_Energy( reax_system *system, control_params *control,
}
}
}
+ }
}
}
@@ -166,6 +170,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
{
+ int i;
real p_vdW1, p_vdW1i;
list *far_nbrs;
real e_vdW_total, e_ele_total;
@@ -176,9 +181,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_vdW_total = 0.0;
e_ele_total = 0.0;
+#ifdef _OPENMP
#pragma omp parallel default(shared) reduction(+: e_vdW_total, e_ele_total)
+#endif
{
- int i, j, pj;
+ int j, pj;
int start_i, end_i;
real self_coef;
real powr_vdW1, powgi_vdW1;
@@ -196,12 +203,14 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
tid = omp_get_thread_num( );
#endif
- e_ele = 0;
- e_vdW = 0;
- e_core = 0;
- de_core = 0;
+ e_ele = 0.0;
+ e_vdW = 0.0;
+ e_core = 0.0;
+ de_core = 0.0;
+#ifdef _OPENMP
#pragma omp for schedule(guided)
+#endif
for ( i = 0; i < system->N; ++i )
{
start_i = Start_Index( i, far_nbrs );
@@ -273,7 +282,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
if ( system->reaxprm.gp.vdw_type == 2 || system->reaxprm.gp.vdw_type == 3 )
{
/* innner wall */
- e_core = twbp->ecore * EXP(twbp->acore * (1.0 - (r_ij / twbp->rcore)));
+ e_core = twbp->ecore * EXP( twbp->acore * (1.0 - (r_ij / twbp->rcore)) );
e_vdW += self_coef * Tap * e_core;
e_vdW_total += self_coef * Tap * e_core;
@@ -322,7 +331,9 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
#endif
rvec_iMultiply( ext_press, nbr_pj->rel_box, temp );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (vdW_Coulomb_Energy_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
@@ -522,7 +533,9 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_vdW_total = 0.0;
e_ele_total = 0.0;
+#ifdef _OPENMP
#pragma omp parallel default(shared) reduction(+: e_vdW_total, e_ele_total)
+#endif
{
int i, j, pj, r;
int type_i, type_j, tmin, tmax;
@@ -537,9 +550,9 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
int tid;
tid = omp_get_thread_num( );
-#endif
#pragma omp for schedule(guided)
+#endif
for ( i = 0; i < system->N; ++i )
{
type_i = system->atoms[i].type;
@@ -617,7 +630,9 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
rvec_Add( workspace->f_local[tid * system->N + j], temp );
#endif
rvec_iMultiply( ext_press, nbr_pj->rel_box, temp );
- #pragma omp critical
+#ifdef _OPENMP
+ #pragma omp critical (Tabulated_vdW_Coulomb_Energy_ext_press)
+#endif
{
rvec_Add( data->ext_press, ext_press );
}
diff --git a/sPuReMD/src/vector.c b/sPuReMD/src/vector.c
index e4bd42ca..f320d857 100644
--- a/sPuReMD/src/vector.c
+++ b/sPuReMD/src/vector.c
@@ -32,12 +32,18 @@ inline int Vector_isZero( const real * const v, const unsigned int k )
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp single
+#endif
{
ret = TRUE;
}
+#ifdef _OPENMP
+ #pragma omp barrier
+
#pragma omp for reduction(&&: ret) schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
if ( FABS( v[i] ) > ALMOST_ZERO )
@@ -46,6 +52,10 @@ inline int Vector_isZero( const real * const v, const unsigned int k )
}
}
+#ifdef _OPENMP
+ #pragma omp barrier
+#endif
+
return ret;
}
@@ -54,7 +64,9 @@ inline void Vector_MakeZero( real * const v, const unsigned int k )
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
v[i] = ZERO;
@@ -66,7 +78,9 @@ inline void Vector_Copy( real * const dest, const real * const v, const unsigned
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
dest[i] = v[i];
@@ -78,7 +92,9 @@ inline void Vector_Scale( real * const dest, const real c, const real * const v,
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
dest[i] = c * v[i];
@@ -91,7 +107,9 @@ inline void Vector_Sum( real * const dest, const real c, const real * const v, c
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
dest[i] = c * v[i] + d * y[i];
@@ -103,7 +121,9 @@ inline void Vector_Add( real * const dest, const real c, const real * const v, c
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
dest[i] += c * v[i];
@@ -132,17 +152,27 @@ inline real Dot( const real * const v1, const real * const v2, const unsigned in
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp single
+#endif
{
ret2 = ZERO;
}
+#ifdef _OPENMP
+ #pragma omp barrier
+
#pragma omp for reduction(+: ret2) schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
ret2 += v1[i] * v2[i];
}
+#ifdef _OPENMP
+ #pragma omp barrier
+#endif
+
return ret2;
}
@@ -151,18 +181,37 @@ inline real Norm( const real * const v1, const unsigned int k )
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp single
+#endif
{
ret2 = ZERO;
}
+#ifdef _OPENMP
+ #pragma omp barrier
+
#pragma omp for reduction(+: ret2) schedule(static)
+#endif
for ( i = 0; i < k; ++i )
{
ret2 += SQR( v1[i] );
}
- return SQRT( ret2 );
+#ifdef _OPENMP
+ #pragma omp barrier
+
+ #pragma omp single
+#endif
+ {
+ ret2 = SQRT( ret2 );
+ }
+
+#ifdef _OPENMP
+ #pragma omp barrier
+#endif
+
+ return ret2;
}
@@ -303,9 +352,9 @@ inline real rvec_Norm( const rvec v )
inline int rvec_isZero( const rvec v )
{
- if ( fabs(v[0]) > ALMOST_ZERO ||
- fabs(v[1]) > ALMOST_ZERO ||
- fabs(v[2]) > ALMOST_ZERO )
+ if ( FABS(v[0]) > ALMOST_ZERO ||
+ FABS(v[1]) > ALMOST_ZERO ||
+ FABS(v[2]) > ALMOST_ZERO )
{
return FALSE;
}
diff --git a/sPuReMD/src/vector.h b/sPuReMD/src/vector.h
index 98ba7dd8..27ceb241 100644
--- a/sPuReMD/src/vector.h
+++ b/sPuReMD/src/vector.h
@@ -24,6 +24,7 @@
#include "mytypes.h"
+
int Vector_isZero( const real * const, const unsigned int );
void Vector_MakeZero( real * const, const unsigned int );
void Vector_Copy( real * const, const real * const, const unsigned int );
@@ -80,4 +81,5 @@ void ivec_Scale( ivec, const real, const ivec );
void ivec_rScale( ivec, const real, const rvec );
void ivec_Sum( ivec, const ivec, const ivec );
+
#endif
--
GitLab