From 35394c7e03e00a0086960943be683625a861a5c4 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 10 Dec 2020 10:57:30 -0500
Subject: [PATCH] sPuReMD: fix memory leak with ICHOLT preconditioner (droptol
 for OpenMP only) for charge solver. Default to Jacobi preconditioner. Resolve
 several compiler warnings around C strings (snprintf, etc.). Other code
 clean-up.

---
 sPuReMD/src/analyze.c     | 34 ++++++++++++++++++++--------------
 sPuReMD/src/control.c     |  4 ++--
 sPuReMD/src/geo_tools.c   |  2 +-
 sPuReMD/src/lin_alg.c     |  7 +++++++
 sPuReMD/src/multi_body.c  | 10 +++++-----
 sPuReMD/src/reset_tools.c |  2 ++
 sPuReMD/src/restart.c     |  2 +-
 7 files changed, 38 insertions(+), 23 deletions(-)

diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index c02ac49c..fa831353 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -197,12 +197,12 @@ static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
 }
 
 
-static void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
-       size_t size )
+static void Print_Molecule( reax_system *system, molecule *m, int mode,
+        char *s, size_t s_size, char *s_out, size_t s_out_size )
 {
     int j, atom;
 
-    s[0] = 0;
+    s[0] = '\0';
 
     if ( mode == 1 )
     {
@@ -211,8 +211,9 @@ static void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
         {
             if ( m->mtypes[j] )
             {
-                snprintf( s, size, "%s%14s%3d", s, system->reax_param.sbp[j].name, m->mtypes[j] );
-                s[size - 1] = '\0';
+                snprintf( s_out, s_out_size, "%.*s%14s%3d", (int) MAX(s_out_size - 18, 0),
+                        s, system->reax_param.sbp[j].name, m->mtypes[j] );
+                s_out[s_out_size - 1] = '\0';
             }
         }
     }
@@ -222,8 +223,10 @@ static void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
         for ( j = 0; j < m->atom_count; ++j )
         {
             atom = m->atom_list[j];
-            snprintf( s, size, "%s%s(%d)",
-                     s, system->reax_param.sbp[ system->atoms[atom].type ].name, atom );
+            snprintf( s_out, s_out_size, "%.*s%.*s(%d)", (int) MAX((s_out_size - 5) / 2, 0), s,
+                    (int) MAX((s_out_size - 5) / 2, 0),
+                    system->reax_param.sbp[ system->atoms[atom].type ].name, atom );
+            s_out[s_out_size - 1] = '\0';
         }
     }
 }
@@ -238,6 +241,7 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
     int *old_mark = workspace->old_mark;
     int num_old, num_new;
     char s[MAX_MOLECULE_SIZE * 10];
+    char s_out[MAX_MOLECULE_SIZE * 10];
     reax_list *new_bonds = lists[BONDS];
     reax_list *old_bonds = lists[OLD_BONDS];
     molecule old_molecules[20], new_molecules[20];
@@ -306,8 +310,8 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
                 fprintf( fout, "old molecules: " );
                 for ( i = 0; i < num_old; ++i )
                 {
-                    Print_Molecule( system, &old_molecules[i], 1, &s[0],
-                            MAX_MOLECULE_SIZE * 10 );
+                    Print_Molecule( system, &old_molecules[i], 1, s,
+                            sizeof(s), s_out, sizeof(s_out) );
                     fprintf( fout, "%s\t", s );
                 }
                 fprintf( fout, "\n" );
@@ -315,8 +319,8 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
                 fprintf( fout, "new molecules: " );
                 for ( i = 0; i < num_new; ++i )
                 {
-                    Print_Molecule( system, &new_molecules[i], 1, &s[0],
-                            MAX_MOLECULE_SIZE * 10 );
+                    Print_Molecule( system, &new_molecules[i], 1, s,
+                            sizeof(s), s_out, sizeof(s_out) );
                     fprintf( fout, "%s\t", s );
                 }
                 fprintf( fout, "\n" );
@@ -595,6 +599,7 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
     int *mark = workspace->mark;
     int num_fragments, num_fragment_types;
     char fragment[MAX_ATOM_TYPES];
+    char fragment_out[MAX_ATOM_TYPES];
     char fragments[MAX_FRAGMENT_TYPES][MAX_ATOM_TYPES];
     int fragment_count[MAX_FRAGMENT_TYPES];
     molecule m;
@@ -618,13 +623,14 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
             memset( m.mtypes, 0, MAX_ATOM_TYPES * sizeof(int) );
             Visit_Bonds( atom, mark, m.mtypes, system, control, new_bonds, ignore );
             ++num_fragments;
-            Print_Molecule( system, &m, 1, fragment, MAX_FRAGMENT_TYPES );
+            Print_Molecule( system, &m, 1, fragment, sizeof(fragment),
+                    fragment_out, sizeof(fragment_out) );
 
             /* check if a similar fragment already exists */
             flag = 0;
             for ( i = 0; i < num_fragment_types; ++i )
             {
-                if ( !strncmp( fragments[i], fragment, MAX_ATOM_TYPES ) )
+                if ( !strncmp( fragments[i], fragment_out, MAX_ATOM_TYPES ) )
                 {
                     ++fragment_count[i];
                     flag = 1;
@@ -635,7 +641,7 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
             if ( flag == 0 )
             {
                 /* it is a new one, add to the fragments list */
-                strncpy( fragments[num_fragment_types], fragment,
+                strncpy( fragments[num_fragment_types], fragment_out,
                         sizeof(fragments[num_fragment_types]) - 1 );
                 fragments[num_fragment_types][sizeof(fragments[num_fragment_types]) - 1] = '\0';
                 fragment_count[num_fragment_types] = 1;
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index b8cceaa7..28978fa2 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -92,10 +92,10 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
     strncpy( control->cm_init_guess_gd_model, "frozen_model.pb", sizeof(control->cm_init_guess_gd_model) - 1 );
     control->cm_init_guess_gd_model[sizeof(control->cm_init_guess_gd_model) - 1] = '\0';
     control->cm_init_guess_win_size = 5;
-    control->cm_solver_pre_comp_type = ICHOLT_PC;
+    control->cm_solver_pre_comp_type = JACOBI_PC;
     control->cm_solver_pre_comp_sweeps = 3;
     control->cm_solver_pre_comp_sai_thres = 0.1;
-    control->cm_solver_pre_comp_refactor = 100;
+    control->cm_solver_pre_comp_refactor = 1;
     control->cm_solver_pre_comp_droptol = 0.01;
     control->cm_solver_pre_app_type = TRI_SOLVE_PA;
     control->cm_solver_pre_app_jacobi_iters = 50;
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 0fef6e03..94472a6d 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -275,7 +275,7 @@ static int Count_Atoms( reax_system *system, FILE *fp, int geo_format )
 
 #if defined(DEBUG)
     fprintf( stderr, "[INFO] Count_Atoms: " );
-    fprintf( stderr, "N = %d\n", system->N );
+    fprintf( stderr, "N = %d\n", n );
 #endif
 
     return n;
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 44ddce77..29f4eb8e 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -470,6 +470,13 @@ void Calculate_Droptol( const sparse_matrix * const A,
         {
             droptol[i] = SQRT( droptol[i] ) * dtol;
         }
+
+#if defined(_OPENMP)
+        #pragma omp master
+        {
+            sfree( droptol_local, "Calculate_Droptol::droptol_local" );
+        }
+#endif
     }
 }
 
diff --git a/sPuReMD/src/multi_body.c b/sPuReMD/src/multi_body.c
index f1216ac5..1dd047df 100644
--- a/sPuReMD/src/multi_body.c
+++ b/sPuReMD/src/multi_body.c
@@ -66,7 +66,7 @@ void Atom_Energy( reax_system *system, control_params *control,
         /* lone-pair Energy */
         p_lp2 = sbp_i->p_lp2;
         expvd2 = EXP( -75.0 * workspace->Delta_lp[i] );
-        inv_expvd2 = 1.0 / (1.0 + expvd2 );
+        inv_expvd2 = 1.0 / (1.0 + expvd2);
 
         /* calculate the energy */
         e_lp = p_lp2 * workspace->Delta_lp[i] * inv_expvd2;
@@ -178,16 +178,16 @@ void Atom_Energy( reax_system *system, control_params *control,
         exp_ovun2 = EXP( p_ovun2 * Delta_lpcorr );
         inv_exp_ovun2 = 1.0 / (1.0 + exp_ovun2);
 
-        DlpVi = 1.0 / (Delta_lpcorr + sbp_i->valency + 1e-8 );
+        DlpVi = 1.0 / (Delta_lpcorr + sbp_i->valency + 1.0e-8 );
         CEover1 = Delta_lpcorr * DlpVi * inv_exp_ovun2;
 
         e_ov = sum_ovun1 * CEover1;
         data->E_Ov += e_ov;
 
         CEover2 = sum_ovun1 * DlpVi * inv_exp_ovun2
-            * ( 1.0 - Delta_lpcorr * ( DlpVi + p_ovun2 * exp_ovun2 * inv_exp_ovun2 ) );
+            * (1.0 - Delta_lpcorr * (DlpVi + p_ovun2 * exp_ovun2 * inv_exp_ovun2));
 
-        CEover3 = CEover2 * (1.0 - dfvl * workspace->dDelta_lp[i] * inv_exp_ovun1 );
+        CEover3 = CEover2 * (1.0 - dfvl * workspace->dDelta_lp[i] * inv_exp_ovun1);
 
         CEover4 = CEover2 * (dfvl * workspace->Delta_lp_temp[i])
             * p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1);
@@ -205,7 +205,7 @@ void Atom_Energy( reax_system *system, control_params *control,
         e_un = -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
         data->E_Un += e_un;
 
-        CEunder1 = inv_exp_ovun2n * ( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8
+        CEunder1 = inv_exp_ovun2n * (p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8
                 + p_ovun2 * e_un * exp_ovun2n);
         CEunder2 = -e_un * p_ovun8 * exp_ovun8 * inv_exp_ovun8;
         CEunder3 = CEunder1 * (1.0 - dfvl * workspace->dDelta_lp[i] * inv_exp_ovun1);
diff --git a/sPuReMD/src/reset_tools.c b/sPuReMD/src/reset_tools.c
index ac2718f8..57ececef 100644
--- a/sPuReMD/src/reset_tools.c
+++ b/sPuReMD/src/reset_tools.c
@@ -27,7 +27,9 @@
 
 void Reset_Pressures( control_params *control, simulation_data *data )
 {
+#if defined(_OPENMP)
     int i;
+#endif
 
     rtensor_MakeZero( data->flex_bar.P );
     data->iso_bar.P = 0.0;
diff --git a/sPuReMD/src/restart.c b/sPuReMD/src/restart.c
index 32b3a622..4afaa367 100644
--- a/sPuReMD/src/restart.c
+++ b/sPuReMD/src/restart.c
@@ -37,7 +37,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
     restart_header res_header;
     restart_atom res_data;
 
-    snprintf( fname, MAX_STR, "%s.res%d", control->sim_name, data->step );
+    snprintf( fname, MAX_STR, "%.*s.res%d", MAX_STR - 12, control->sim_name, data->step );
     fres = sfopen( fname, "wb" );
 
     res_header.step = data->step;
-- 
GitLab