diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 783c73852c03a4a0aa7790830d71009c881a7fff..bec04123391c1a165a1252d8d7ccb66a8deb37c0 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -22,6 +22,7 @@
#include "control.h"
#include <ctype.h>
+#include <zlib.h>
#include "traj.h"
#include "tool_box.h"
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index e0ea4d5f09c50d7a7f331332d4261e730460483c..bcf785139b285daa63854785c97860e3fd83584b 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -29,10 +29,10 @@
/* Intel MKL */
#if defined(HAVE_LAPACKE_MKL)
-#include "mkl.h"
+ #include "mkl.h"
/* reference LAPACK */
#elif defined(HAVE_LAPACKE)
-#include "lapacke.h"
+ #include "lapacke.h"
#endif
@@ -3233,7 +3233,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
Vector_Copy( u[0], z[0], N );
u[0][0] += ( u[0][0] < 0.0 ? -1 : 1 ) * w[0];
- Vector_Scale( u[0], 1 / Norm( u[0], N ), u[0], N );
+ Vector_Scale( u[0], 1.0 / Norm( u[0], N ), u[0], N );
w[0] *= ( u[0][0] < 0.0 ? 1 : -1 );
t_vops += Get_Timing_Info( t_start );
@@ -3243,12 +3243,13 @@ int GMRES_HouseHolder( const static_storage * const workspace,
{
/* compute v_j */
t_start = Get_Time( );
- Vector_Scale( z[j], -2 * u[j][j], u[j], N );
+ Vector_Scale( z[j], -2.0 * u[j][j], u[j], N );
z[j][j] += 1.; /* due to e_j */
for ( i = j - 1; i >= 0; --i )
{
- Vector_Add( z[j] + i, -2 * Dot( u[i] + i, z[j] + i, N - i ), u[i] + i, N - i );
+ Vector_Add( z[j] + i, -2.0 * Dot( u[i] + i, z[j] + i, N - i ),
+ u[i] + i, N - i );
}
t_vops += Get_Timing_Info( t_start );
@@ -3267,7 +3268,8 @@ int GMRES_HouseHolder( const static_storage * const workspace,
t_start = Get_Time( );
for ( i = 0; i <= j; ++i )
{
- Vector_Add( v + i, -2 * Dot( u[i] + i, v + i, N - i ), u[i] + i, N - i );
+ Vector_Add( v + i, -2.0 * Dot( u[i] + i, v + i, N - i ),
+ u[i] + i, N - i );
}
if ( !Vector_isZero( v + (j + 1), N - (j + 1) ) )
@@ -3279,16 +3281,17 @@ int GMRES_HouseHolder( const static_storage * const workspace,
#ifdef _OPENMP
#pragma omp single
#endif
- u[j + 1][j + 1] += ( v[j + 1] < 0.0 ? -1 : 1 ) * temp;
+ u[j + 1][j + 1] += ( v[j + 1] < 0.0 ? -1.0 : 1.0 ) * temp;
#ifdef _OPENMP
#pragma omp barrier
#endif
- Vector_Scale( u[j + 1], 1 / Norm( u[j + 1], N ), u[j + 1], N );
+ Vector_Scale( u[j + 1], 1.0 / Norm( u[j + 1], N ), u[j + 1], N );
/* overwrite v with P_m+1 * v */
- temp = 2 * Dot( u[j + 1] + (j + 1), v + (j + 1), N - (j + 1) ) * u[j + 1][j + 1];
+ temp = 2.0 * Dot( u[j + 1] + (j + 1), v + (j + 1), N - (j + 1) )
+ * u[j + 1][j + 1];
#ifdef _OPENMP
#pragma omp single
#endif
@@ -3299,7 +3302,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
#endif
Vector_MakeZero( v + (j + 2), N - (j + 2) );
-// Vector_Add( v, -2 * Dot( u[j+1], v, N ), u[j+1], N );
+// Vector_Add( v, -2.0 * Dot( u[j+1], v, N ), u[j+1], N );
}
t_vops += Get_Timing_Info( t_start );
@@ -3314,7 +3317,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
tmp1 = workspace->hc[i] * v[i] + workspace->hs[i] * v[i + 1];
tmp2 = -workspace->hs[i] * v[i] + workspace->hc[i] * v[i + 1];
- v[i] = tmp1;
+ v[i] = tmp1;
v[i + 1] = tmp2;
}
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index a173ecca352c8460181947f859c9d17491f1c236..dfadda825b596470163f7b396abb01688954af1d 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -28,13 +28,9 @@
#endif
#include <math.h>
-#include <random.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <sys/time.h>
-#include <time.h>
-#include <zlib.h>
#ifdef _OPENMP
#include <omp.h>
@@ -43,13 +39,13 @@
//#define DEBUG_FOCUS
//#define TEST_FORCES
//#define TEST_ENERGY
-#define REORDER_ATOMS // turns on nbrgen opt by re-ordering atoms
+#define REORDER_ATOMS /* turns on nbrgen opt by re-ordering atoms */
//#define LGJ
-#define SUCCESS 1
-#define FAILURE 0
-#define TRUE 1
-#define FALSE 0
+#define SUCCESS (1)
+#define FAILURE (0)
+#define TRUE (1)
+#define FALSE (0)
#define LOG log
#define EXP exp
@@ -81,56 +77,56 @@
#endif
#if !defined(PI)
- #define PI 3.14159265
+ #define PI (3.14159265)
#endif
-#define C_ele 332.06371
-//#define K_B 503.398008 // kcal/mol/K
-#define K_B 0.831687 // amu A^2 / ps^2 / K
+#define C_ele (332.06371)
+//#define K_B (503.398008) // kcal/mol/K
+#define K_B (0.831687) // amu A^2 / ps^2 / K
#define F_CONV (1e6 / 48.88821291 / 48.88821291) // --> amu A / ps^2
-#define E_CONV 0.002391 // amu A^2 / ps^2 --> kcal/mol
-#define EV_to_KCALpMOL 14.400000 // ElectronVolt --> KCAL per MOLe
-#define KCALpMOL_to_EV 23.060549 // 23.020000//KCAL per MOLe --> ElectronVolt
-#define ECxA_to_DEBYE 4.803204 // elem. charge * angstrom -> debye conv
-#define CAL_to_JOULES 4.184000 // CALories --> JOULES
-#define JOULES_to_CAL (1/4.184000) // JOULES --> CALories
-#define AMU_to_GRAM 1.6605e-24
-#define ANG_to_CM 1.0e-8
-#define AVOGNR 6.0221367e23
+#define E_CONV (0.002391) // amu A^2 / ps^2 --> kcal/mol
+#define EV_to_KCALpMOL (14.400000) // ElectronVolt --> KCAL per MOLe
+#define KCALpMOL_to_EV (23.060549) // 23.020000//KCAL per MOLe --> ElectronVolt
+#define ECxA_to_DEBYE (4.803204) // elem. charge * angstrom -> debye conv
+#define CAL_to_JOULES (4.184000) // CALories --> JOULES
+#define JOULES_to_CAL (1.0 / 4.184000) // JOULES --> CALories
+#define AMU_to_GRAM (1.6605e-24)
+#define ANG_to_CM (1.0e-8)
+#define AVOGNR (6.0221367e23)
#define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
-#define MAX_STR 1024
-#define MAX_LINE 1024
-#define MAX_TOKENS 1024
-#define MAX_TOKEN_LEN 1024
-
-#define MAX_ATOM_ID 100000
-#define MAX_RESTRICT 15
-#define MAX_MOLECULE_SIZE 20
-#define MAX_ATOM_TYPES 25
-
-#define MAX_GRID 50
-#define MAX_3BODY_PARAM 5
-#define MAX_4BODY_PARAM 5
-#define NO_OF_INTERACTIONS 10
-
-#define MAX_dV 1.01
-#define MIN_dV 0.99
-#define MAX_dT 4.00
-#define MIN_dT 0.00
-
-#define ZERO 0.000000000000000e+00
-#define ALMOST_ZERO 1e-10
-#define NEG_INF -1e10
-#define NO_BOND 1e-3
-#define HB_THRESHOLD 1e-2
-#define MAX_BONDS 40
-#define MIN_BONDS 15
-#define MIN_HBONDS 50
-#define SAFE_HBONDS 1.4
-#define MIN_GCELL_POPL 50
-#define SAFE_ZONE 1.2
-#define DANGER_ZONE 0.95
-#define LOOSE_ZONE 0.75
+#define MAX_STR (1024)
+#define MAX_LINE (1024)
+#define MAX_TOKENS (1024)
+#define MAX_TOKEN_LEN (1024)
+
+#define MAX_ATOM_ID (100000)
+#define MAX_RESTRICT (15)
+#define MAX_MOLECULE_SIZE (20)
+#define MAX_ATOM_TYPES (25)
+
+#define MAX_GRID (50)
+#define MAX_3BODY_PARAM (5)
+#define MAX_4BODY_PARAM (5)
+#define NO_OF_INTERACTIONS (10)
+
+#define MAX_dV (1.01)
+#define MIN_dV (0.99)
+#define MAX_dT (4.00)
+#define MIN_dT (0.00)
+
+#define ZERO (0.000000000000000e+00)
+#define ALMOST_ZERO (1e-10)
+#define NEG_INF (-1e10)
+#define NO_BOND (1e-3)
+#define HB_THRESHOLD (1e-2)
+#define MAX_BONDS (40)
+#define MIN_BONDS (15)
+#define MIN_HBONDS (50)
+#define SAFE_HBONDS (1.4)
+#define MIN_GCELL_POPL (50)
+#define SAFE_ZONE (1.2)
+#define DANGER_ZONE (0.95)
+#define LOOSE_ZONE (0.75)
/* config params */
@@ -251,7 +247,8 @@ typedef int ivec[3];
typedef real rtensor[3][3];
-/* Force field global params mapping:
+/* Force field global parameters mapping
+ * (contained in section 1 of file):
*
* l[0] = p_boc1
* l[1] = p_boc2
@@ -287,7 +284,7 @@ typedef real rtensor[3][3];
* l[31] = p_ovun4
* l[32] = p_ovun3
* l[33] = p_val8
- * l[34] = ACKS2 bond softness
+ * l[34] = b_s_acks2 (ACKS2 bond softness)
* l[35] = N/A
* l[36] = N/A
* l[37] = version number
@@ -303,11 +300,14 @@ typedef struct
typedef struct
{
/* Line one in field file */
- char name[15]; /* Two character atom name */
+ /* Two character atom name */
+ char name[15];
real r_s;
- real valency; /* Valency of the atom */
- real mass; /* Mass of atom */
+ /* Valency of the atom */
+ real valency;
+ /* Mass of atom */
+ real mass;
real r_vdw;
real epsilon;
real gamma;
@@ -322,8 +322,8 @@ typedef struct
real p_ovun5;
real chi;
real eta;
- /* 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:
+ * 1 for H donor, 2 for acceptors (O,S,N), 0 for others */
int p_hbond;
/* Line three in field file */
@@ -358,7 +358,8 @@ typedef struct
real p_bo6;
real r_s;
real r_p;
- real r_pp; /* r_o distances in BO formula */
+ /* r_o distances in BO formula */
+ real r_pp;
real p_boc3;
real p_boc4;
real p_boc5;
@@ -383,9 +384,11 @@ typedef struct
real acore;
/* electrostatic parameters */
- real gamma; // note: this parameter is gamma^-3 and not gamma.
+ /* note: this parameter is gamma^-3 and not gamma */
+ real gamma;
- real v13cor, ovc;
+ real v13cor;
+ real ovc;
} two_body_parameters;
@@ -483,8 +486,11 @@ typedef struct
rvec side_prop;
rvec nbr_box_press[27];
- rtensor box, box_inv, old_box;
- rtensor trans, trans_inv;
+ rtensor box;
+ rtensor box_inv;
+ rtensor old_box;
+ rtensor trans;
+ rtensor trans_inv;
rtensor g;
} simulation_box;
@@ -501,13 +507,13 @@ typedef struct
rvec len;
rvec inv_len;
- int**** atoms;
- int*** top;
- int*** mark;
- int*** start;
- int*** end;
- ivec**** nbrs;
- rvec**** nbrs_cp;
+ int **** atoms;
+ int *** top;
+ int *** mark;
+ int *** start;
+ int *** end;
+ ivec **** nbrs;
+ rvec **** nbrs_cp;
} grid;
@@ -554,9 +560,10 @@ typedef struct
int reneighbor;
real vlist_cut;
real nbr_cut;
+ /* upper and lower taper */
real r_cut;
real r_sp_cut;
- real r_low; // upper and lower taper
+ real r_low;
real bo_cut;
real thb_cut;
real hb_cut;
@@ -637,7 +644,6 @@ typedef struct
real v_eps;
real v_eps_old;
real a_eps;
-
} isotropic_barostat;
@@ -688,8 +694,10 @@ typedef struct
int prev_steps;
real time;
- real M; /* Total Mass */
- real inv_M; /* 1 / Total Mass */
+ /* Total Mass */
+ real M;
+ /* 1.0 / Total Mass */
+ real inv_M;
rvec xcm; /* Center of mass */
rvec vcm; /* Center of mass velocity */
@@ -720,7 +728,7 @@ typedef struct
real E_Ele; /* Total electrostatics energy */
real E_Pol; /* Polarization energy */
- real N_f; /*Number of degrees of freedom */
+ real N_f; /* Number of degrees of freedom */
rvec t_scale;
rtensor p_scale;
thermostat therm; /* Used in Nose_Hoover method */
@@ -1038,6 +1046,8 @@ typedef struct
/* local forces per thread */
rvec *f_local;
#endif
+ unsigned int temp_int_omp;
+ real temp_real_omp;
reallocate_data realloc;
@@ -1127,7 +1137,7 @@ typedef struct
/* trajectory output function pointer definitions */
int (* write_header)( reax_system*, control_params*, static_storage*, void* );
- int (* append_traj_frame)(reax_system*, control_params*,
+ int (* append_traj_frame)( reax_system*, control_params*,
simulation_data*, static_storage*, reax_list **, void* );
int (* write)( FILE *, const char *, ... );
@@ -1146,7 +1156,6 @@ typedef struct
FILE *ecou;
#endif
- FILE *ftot;
#ifdef TEST_FORCES
FILE *fbo;
FILE *fdbo;
@@ -1157,6 +1166,7 @@ typedef struct
FILE *fhb;
FILE *f4body;
FILE *fnonb;
+ FILE *ftot;
FILE *ftot2;
#endif
} output_controls;
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 447e7cc326a7b10e4f0c60aca94bfe16829c30be..b019dfa0fb7514636d7213c7859be800ec790951 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -546,16 +546,12 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
}
+#if defined(TEST_FORCES)
void Print_Total_Force( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
int i;
-#if !defined(TEST_FORCES)
- char temp[1000];
- snprintf( temp, 1000, "%.*s.ftot", 994, control->sim_name );
- out_control->ftot = fopen( temp, "w" );
-#endif
for ( i = 0; i < system->N; ++i )
{
@@ -567,10 +563,8 @@ void Print_Total_Force( reax_system *system, control_params *control,
}
fflush( out_control->ftot );
-#if !defined(TEST_FORCES)
- fclose( out_control->ftot );
-#endif
}
+#endif
void Output_Results( reax_system *system, control_params *control,
diff --git a/sPuReMD/src/random.c b/sPuReMD/src/random.c
index 4afe6aeac249822b2303a553bbf4938b08892f80..b4ce0ee6dccd19465a905ad991658f8f34f12b05 100644
--- a/sPuReMD/src/random.c
+++ b/sPuReMD/src/random.c
@@ -21,40 +21,13 @@
#include "random.h"
-
-/* System random number generator used linear congruance method with
- * large periodicity for generation of pseudo random number. function
- * Random returns this random number appropriately scaled so that
- * 0 <= Random(range) < range */
-double Random(double range)
-{
- return (random() * range) / 2147483647L;
-}
+#include <time.h>
/* This function seeds the system pseudo random number generator with
* current time. Use this function once in the begining to initialize
- * the system */
-void Randomize()
+ * the system. */
+void Randomize( )
{
- srandom(time(NULL));
-}
-
-
-/* GRandom return random number with gaussian distribution with mean
- * and standard deviation "sigma" */
-double GRandom(double mean, double sigma)
-{
- double v1 = Random(2.0) - 1.0;
- double v2 = Random(2.0) - 1.0;
- double rsq = v1 * v1 + v2 * v2;
-
- while (rsq >= 1.0 || rsq == 0.0)
- {
- v1 = Random(2.0) - 1.0;
- v2 = Random(2.0) - 1.0;
- rsq = v1 * v1 + v2 * v2;
- }
-
- return mean + v1 * sigma * SQRT(-2.0 * LOG(rsq) / rsq);
+ srandom( time(NULL) );
}
diff --git a/sPuReMD/src/random.h b/sPuReMD/src/random.h
index 1fd1642f4e5a7044f898fcc003ded9d428e4d07b..cf8517ab151c3b8628a0c9715476f06db3cfc428 100644
--- a/sPuReMD/src/random.h
+++ b/sPuReMD/src/random.h
@@ -22,23 +22,43 @@
#ifndef __RANDOM_H_
#define __RANDOM_H_
+/* Includes <stdlib.h> for random( ) */
#include "mytypes.h"
-/* System random number generator used linear congruance method with
- large periodicity for generation of pseudo random number. function
- Random returns this random number appropriately scaled so that
- 0 <= Random(range) < range */
-double Random( double );
-
/* This function seeds the system pseudo random number generator with
current time. Use this function once in the begining to initialize
the system */
void Randomize( );
+
+/* System random number generator used linear congruance method with
+ * large periodicity for generation of pseudo random number. Function
+ * Random returns this random number appropriately scaled so that
+ * 0 <= Random(range) < range. */
+static inline double Random( double range )
+{
+ return (random( ) * range) / 2147483647L;
+}
+
+
/* GRandom return random number with gaussian distribution with mean
- and standard deviation "sigma" */
-double GRandom( double, double );
+ * and standard deviation "sigma." */
+static inline double GRandom( double mean, double sigma )
+{
+ double v1 = Random(2.0) - 1.0;
+ double v2 = Random(2.0) - 1.0;
+ double rsq = v1 * v1 + v2 * v2;
+
+ while (rsq >= 1.0 || rsq == 0.0)
+ {
+ v1 = Random(2.0) - 1.0;
+ v2 = Random(2.0) - 1.0;
+ rsq = v1 * v1 + v2 * v2;
+ }
+
+ return mean + v1 * sigma * SQRT(-2.0 * LOG(rsq) / rsq);
+}
#endif
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index c0a61f51d633b63a5c3d8d360d08659cda3e9a80..1cc6ba72f6fa14fe5d501194932f5f1a8434e633 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -22,6 +22,7 @@
#include "tool_box.h"
#include <ctype.h>
+#include <sys/time.h>
/************** taken from box.c **************/
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index 03cc5d74e1a644fdc10d84e8c6b57d9331c1a644..35e893584f524d5bab66ffc88a120aa61c87fd3b 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -23,6 +23,8 @@
#include "list.h"
+#include <zlib.h>
+
/************************************************/
/* CUSTOM FORMAT ROUTINES */
diff --git a/sPuReMD/src/traj.h b/sPuReMD/src/traj.h
index 738bc19333b244a586e8a4bb63d29c275f7915cf..f4d074c8fdd8242737bb38693adad6c03c7638f5 100644
--- a/sPuReMD/src/traj.h
+++ b/sPuReMD/src/traj.h
@@ -24,8 +24,6 @@
#include "mytypes.h"
-#include <zlib.h>
-
#define BLOCK_MARK "REAX_BLOCK_MARK "
#define BLOCK_MARK_LEN 16
diff --git a/sPuReMD/src/vector.h b/sPuReMD/src/vector.h
index f150b84cc933bb4f1f23d400473546794ce4cd00..5550a2d7354362ffcfe24bbdbd94bc77af0afc13 100644
--- a/sPuReMD/src/vector.h
+++ b/sPuReMD/src/vector.h
@@ -24,6 +24,8 @@
#include "mytypes.h"
+#include "random.h"
+
/* file scope to make OpenMP shared (Vector_isZero) */
static unsigned int ret_omp;
@@ -200,10 +202,6 @@ static inline void rvec_Copy( rvec dest, const rvec src )
{
dest[i] = src[i];
}
-
-// dest[0] = src[0];
-// dest[1] = src[1];
-// dest[2] = src[2];
}
static inline void rvec_Scale( rvec ret, const real c, const rvec v )
@@ -217,10 +215,6 @@ static inline void rvec_Scale( rvec ret, const real c, const rvec v )
{
ret[i] = c * v[i];
}
-
-// ret[0] = c * v[0];
-// ret[1] = c * v[1];
-// ret[2] = c * v[2];
}
@@ -235,10 +229,6 @@ static inline void rvec_Add( rvec ret, const rvec v )
{
ret[i] += v[i];
}
-
-// ret[0] += v[0];
-// ret[1] += v[1];
-// ret[2] += v[2];
}
@@ -253,10 +243,6 @@ static inline void rvec_ScaledAdd( rvec ret, const real c, const rvec v )
{
ret[i] += c * v[i];
}
-
-// ret[0] += c * v[0];
-// ret[1] += c * v[1];
-// ret[2] += c * v[2];
}
@@ -271,10 +257,6 @@ static inline void rvec_Sum( rvec ret, const rvec v1 , const rvec v2 )
{
ret[i] = v1[i] + v2[i];
}
-
-// ret[0] = v1[0] + v2[0];
-// ret[1] = v1[1] + v2[1];
-// ret[2] = v1[2] + v2[2];
}
@@ -290,10 +272,6 @@ static inline void rvec_ScaledSum( rvec ret, const real c1, const rvec v1,
{
ret[i] = c1 * v1[i] + c2 * v2[i];
}
-
-// ret[0] = c1 * v1[0] + c2 * v2[0];
-// ret[1] = c1 * v1[1] + c2 * v2[1];
-// ret[2] = c1 * v1[2] + c2 * v2[2];
}
@@ -313,8 +291,6 @@ static inline real rvec_Dot( const rvec v1, const rvec v2 )
}
return ret;
-
-// return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
}
@@ -335,8 +311,6 @@ static inline real rvec_ScaledDot( const real c1, const rvec v1,
}
return c1 * c2 * ret;
-
-// return (c1 * c2) * (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
}
@@ -351,10 +325,6 @@ static inline void rvec_Multiply( rvec r, const rvec v1, const rvec v2 )
{
r[i] = v1[i] * v2[i];
}
-
-// r[0] = v1[0] * v2[0];
-// r[1] = v1[1] * v2[1];
-// r[2] = v1[2] * v2[2];
}
@@ -369,10 +339,6 @@ static inline void rvec_iMultiply( rvec r, const ivec v1, const rvec v2 )
{
r[i] = v1[i] * v2[i];
}
-
-// r[0] = v1[0] * v2[0];
-// r[1] = v1[1] * v2[1];
-// r[2] = v1[2] * v2[2];
}
@@ -387,10 +353,6 @@ static inline void rvec_Divide( rvec r, const rvec v1, const rvec v2 )
{
r[i] = v1[i] / v2[i];
}
-
-// r[0] = v1[0] / v2[0];
-// r[1] = v1[1] / v2[1];
-// r[2] = v1[2] / v2[2];
}
@@ -405,10 +367,6 @@ static inline void rvec_iDivide( rvec r, const rvec v1, const ivec v2 )
{
r[i] = v1[i] / v2[i];
}
-
-// r[0] = v1[0] / v2[0];
-// r[1] = v1[1] / v2[1];
-// r[2] = v1[2] / v2[2];
}
@@ -423,10 +381,6 @@ static inline void rvec_Invert( rvec r, const rvec v )
{
r[i] = 1.0 / v[i];
}
-
-// r[0] = 1.0 / v[0];
-// r[1] = 1.0 / v[1];
-// r[2] = 1.0 / v[2];
}
@@ -442,10 +396,6 @@ static inline void rvec_Cross( rvec ret, const rvec v1, const rvec v2 )
ret[i] = v1[(i + 1) % 3] * v2[(i + 2) % 3]
- v1[(i + 2) % 3] * v2[(i + 1) % 3];
}
-
-// ret[0] = v1[1] * v2[2] - v1[2] * v2[1];
-// ret[1] = v1[2] * v2[0] - v1[0] * v2[2];
-// ret[2] = v1[0] * v2[1] - v1[1] * v2[0];
}
@@ -479,8 +429,6 @@ static inline real rvec_Norm_Sqr( const rvec v )
}
return ret;
-
-// return SQR(v[0]) + SQR(v[1]) + SQR(v[2]);
}
@@ -500,8 +448,6 @@ static inline real rvec_Norm( const rvec v )
}
return SQRT( ret );
-
-// return SQRT( SQR(v[0]) + SQR(v[1]) + SQR(v[2]) );
}
@@ -520,14 +466,6 @@ static inline int rvec_isZero( const rvec v )
}
return ret;
-
-// if ( FABS(v[0]) > ALMOST_ZERO ||
-// FABS(v[1]) > ALMOST_ZERO ||
-// FABS(v[2]) > ALMOST_ZERO )
-// {
-// return FALSE;
-// }
-// return TRUE;
}
@@ -542,13 +480,11 @@ static inline void rvec_MakeZero( rvec v )
{
v[i] = ZERO;
}
-
-// v[0] = ZERO;
-// v[1] = ZERO;
-// v[2] = ZERO;
}
+/* Note: currently not thread-safe since uses random( ) underneath,
+ * change to reentrant version later if needed */
static inline void rvec_Random( rvec v )
{
v[0] = Random(2.0) - 1.0;
@@ -788,10 +724,6 @@ static inline void ivec_MakeZero( ivec v )
{
v[i] = 0;
}
-
-// v[0] = 0;
-// v[1] = 0;
-// v[2] = 0;
}
@@ -806,10 +738,6 @@ static inline void ivec_Copy( ivec dest, const ivec src )
{
dest[i] = src[i];
}
-
-// dest[0] = src[0];
-// dest[1] = src[1];
-// dest[2] = src[2];
}
@@ -824,10 +752,6 @@ static inline void ivec_Scale( ivec dest, const real C, const ivec src )
{
dest[i] = C * src[i];
}
-
-// dest[0] = C * src[0];
-// dest[1] = C * src[1];
-// dest[2] = C * src[2];
}
@@ -842,10 +766,6 @@ static inline void ivec_rScale( ivec dest, const real C, const rvec src )
{
dest[i] = (int)(C * src[i]);
}
-
-// dest[0] = (int)(C * src[0]);
-// dest[1] = (int)(C * src[1]);
-// dest[2] = (int)(C * src[2]);
}
@@ -864,12 +784,6 @@ static inline int ivec_isZero( const ivec v )
}
return ret;
-
-// if ( v[0] == 0 && v[1] == 0 && v[2] == 0 )
-// {
-// return TRUE;
-// }
-// return FALSE;
}
@@ -888,13 +802,6 @@ static inline int ivec_isEqual( const ivec v1, const ivec v2 )
}
return ret;
-
-// if ( v1[0] == v2[0] && v1[1] == v2[1] && v1[2] == v2[2] )
-// {
-// return TRUE;
-// }
-//
-// return FALSE;
}
@@ -909,10 +816,6 @@ static inline void ivec_Sum( ivec dest, const ivec v1, const ivec v2 )
{
dest[i] = v1[i] * v2[i];
}
-
-// dest[0] = v1[0] + v2[0];
-// dest[1] = v1[1] + v2[1];
-// dest[2] = v1[2] + v2[2];
}
diff --git a/tools/driver.py b/tools/driver.py
index 59a8c4a8ba7b6eda19fe08b8dcb9e61fd18e38fc..93163deaec3d723284d3d61d73d57d3d23df3fb3 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -1,6 +1,9 @@
#!/bin/python3
-from ctypes import *
+from ctypes import c_int, c_double, c_char, c_char_p, c_void_p, \
+ Structure, Union, POINTER, CFUNCTYPE, cdll
+import sqlite3 as sq3
+from os import path
class BondOrderData(Structure):
@@ -220,7 +223,6 @@ class SimulationData(Structure):
("ext_press", c_double * 3),
("kin_press", c_double),
("tot_press", c_double * 3),
- ("tot_press", c_double * 3),
("timing", ReaxTiming),
]
@@ -236,6 +238,141 @@ class ReaxAtom(Structure):
]
+def create_db(name='spuremd.db'):
+ conn = sq3.connect(name)
+
+ conn.executescript("""
+ CREATE TABLE simulation(
+ id integer,
+ date text,
+ name text,
+ ensemble_type integer,
+ steps integer,
+ time_step integer,
+ restart_format integer,
+ random_velocity integer,
+ reposition_atoms integer,
+ peroidic_boundary integer,
+ geo_format integer,
+ restrict_bonds integer,
+ tabulate_long_range integer,
+ reneighbor integer,
+ vlist_cutoff real,
+ neighbor_cutoff real,
+ three_body_cutoff real,
+ hydrogen_bond_cutoff real,
+ bond_graph_cutoff real,
+ charge_method integer,
+ cm_q_net real,
+ cm_solver_type integer,
+ cm_solver_max_iters integer,
+ cm_solver_restart integer,
+ cm_solver_q_err real,
+ cm_domain_sparsity real,
+ cm_solver_pre_comp_type integer,
+ cm_solver_pre_comp_refactor integer,
+ cm_solver_pre_comp_droptol real,
+ cm_solver_pre_comp_sweeps integer,
+ cm_solver_pre_comp_sai_thres real,
+ cm_solver_pre_app_type integer,
+ cm_solver_pre_app_jacobi_iters integer,
+ temp_init real,
+ temp_final real,
+ temp_mass real,
+ temp_mode integer,
+ temp_rate real,
+ temp_freq integer,
+ pressure real,
+ pressure_mass real,
+ compress integer,
+ pressure_mode integer,
+ remove_center_of_mass integer,
+ debug_level integer,
+ write_freq integer,
+ traj_compress integer,
+ traj_format integer,
+ traj_title text,
+ atom_info integer,
+ atom_velocities integer,
+ atom_forces integer,
+ bond_info integer,
+ angle_info integer,
+ test_forces integer,
+ molecule_analysis integer,
+ freq_molecule_analysis integer,
+ ignore integer,
+ dipole_analysis integer,
+ freq_dipole_analysis integer,
+ diffusion_coefficient integer,
+ freq_diffusion_coefficient integer,
+ restrict_type integer,
+ PRIMARY KEY (id)
+ );
+
+ CREATE TABLE system_properties(
+ id integer,
+ step integer,
+ total_energy real,
+ potential_energy real,
+ kinetic_energy real,
+ temperature real,
+ volume real,
+ pressure real,
+ PRIMARY KEY (id, step)
+ );
+
+ CREATE TABLE potential(
+ id integer,
+ step integer,
+ bond_energy real,
+ atom_energy real,
+ lone_pair_energy real,
+ angle_energy real,
+ coa_energy real,
+ hydrogen_bond_energy real,
+ torsion_energy real,
+ conjugation_energy real,
+ van_der_waals_energy real,
+ coulombic_energy real,
+ polarization_energy real,
+ PRIMARY KEY (id, step)
+ );
+
+ CREATE TABLE trajectory(
+ id integer,
+ step integer,
+ atom_id integer,
+ position_x real,
+ position_y real,
+ position_z real,
+ charge real,
+ PRIMARY KEY (id, step, atom_id)
+ );
+
+ CREATE TABLE performance(
+ id integer,
+ step integer,
+ time_total real,
+ time_nbrs real,
+ time_init real,
+ time_bonded real,
+ time_nonbonded real,
+ time_cm real,
+ time_cm_sort real,
+ cm_solver_iters integer,
+ time_cm_pre_comp real,
+ time_cm_pre_app real,
+ time_cm_solver_spmv real,
+ time_cm_solver_vec_ops real,
+ time_cm_solver_orthog real,
+ time_cm_solver_tri_solve real,
+ PRIMARY KEY (id, step)
+ );
+ """)
+
+ conn.close()
+
+
if __name__ == '__main__':
lib = cdll.LoadLibrary("libspuremd.so.1")
@@ -258,10 +395,6 @@ if __name__ == '__main__':
CALLBACKFUNC = CFUNCTYPE(None, POINTER(ReaxAtom),
POINTER(SimulationData), POINTER(ReaxList))
- def get_simulation_step_results(atoms, data, lists):
- print("{0:24.15f} {1:24.15f} {2:24.15f}".format(
- data[0].E_Tot, data[0].E_Kin, data[0].E_Pot))
-
setup_callback = lib.setup_callback
setup_callback.restype = c_int
@@ -269,6 +402,57 @@ if __name__ == '__main__':
set_output_enabled.argtypes = [c_void_p, c_int]
set_output_enabled.restype = c_int
+ db_file = "spuremd.db"
+
+ if not path.isfile(db_file):
+ create_db(db_file)
+
+ conn = sq3.connect(db_file)
+ record_potential = True
+ record_trajectory = False
+ record_performance = True
+
+ def get_simulation_step_results(atoms, data, lists):
+ print("{0:24.15f} {1:24.15f} {2:24.15f}".format(
+ data[0].E_Tot, data[0].E_Kin, data[0].E_Pot))
+
+ if data[0].step == 0:
+ #TODO: insert data into simulation table
+ pass
+
+ with conn:
+ conn.execute("INSERT INTO system_properties VALUES (?,?,?,?,?,?,?,?)",
+ (0, data[0].step, data[0].E_Tot, data[0].E_Pot, data[0].E_Kin,
+ data[0].therm.T, 0.0, data[0].iso_bar.P))
+ # MISSING: ID, system->box.volume
+
+ if record_potential:
+ with conn:
+ conn.execute("INSERT INTO potential VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?)",
+ (0, data[0].step, data[0].E_BE, data[0].E_Ov + data[0].E_Un,
+ data[0].E_Lp, data[0].E_Ang + data[0].E_Pen, data[0].E_Coa,
+ data[0].E_HB, data[0].E_Tor, data[0].E_Con, data[0].E_vdW,
+ data[0].E_Ele, data[0].E_Pol))
+
+ if record_trajectory:
+ with conn:
+ for i in range(6540):
+ conn.execute("INSERT INTO trajectory VALUES (?,?,?,?,?,?,?)",
+ (0, data[0].step, i, atoms[i].x[0], atoms[i].x[1],
+ atoms[i].x[2], atoms[i].q))
+
+ if record_performance:
+ with conn:
+ conn.execute("INSERT INTO performance VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)",
+ (0, data[0].step, data[0].timing.total, data[0].timing.nbrs, data[0].timing.init_forces,
+ data[0].timing.bonded, data[0].timing.nonb, data[0].timing.cm,
+ data[0].timing.cm_sort_mat_rows, data[0].timing.cm_solver_iters,
+ data[0].timing.cm_solver_pre_comp,
+ data[0].timing.cm_solver_pre_app, data[0].timing.cm_solver_spmv,
+ data[0].timing.cm_solver_vector_ops, data[0].timing.cm_solver_orthog,
+ data[0].timing.cm_solver_tri_solve))
+
+
handle = setup(b"data/benchmarks/water/water_6540.pdb",
b"data/benchmarks/water/ffield.water",
b"environ/param.gpu.water")
@@ -288,4 +472,5 @@ if __name__ == '__main__':
print("{0:9d} {1:24.15f} {2:24.15f} {3:24.15f} {4:24.15f}".format(
i + 1, atoms[i].x[0], atoms[i].x[1], atoms[i].x[2], atoms[i].q))
+ conn.close()
cleanup(handle)