Fix hvap and prevent undefined behavior from using unitialized variables

src/BlockOutput.cpp

@@ -217,7 +217,7 @@ void BlockAverages::InitWatchSingle(config_setup::TrackedVars const& tracked)
     blocks[out::PRESSURE_IDX].SetRef(&var->pressure[b], b);
     blocks[out::MOL_NUM_IDX].SetRef(&var->numByBox[b], b);
     blocks[out::DENSITY_IDX].SetRef(&var->densityTot[b], b);
-    blocks[out::COMPRESSIBILITY_IDX].SetRef(&var->compressability[b], b);
+    blocks[out::COMPRESSIBILITY_IDX].SetRef(&var->compressibility[b], b);
     blocks[out::ENTHALPY_IDX].SetRef(&var->enthalpy[b], b);
     blocks[out::SURF_TENSION_IDX].SetRef(&var->surfaceTens[b], b);
 #if ENSEMBLE == GEMC

src/ConsoleOutput.cpp

@@ -229,7 +229,11 @@ void ConsoleOutput::PrintStatistic(const uint box, const ulong step) const
 
   if(enablePressure) {
     printElement(var->pressure[box], elementWidth);
-    printElement((var->pressure[box]) * (var->volumeRef[box]) / (var->numByBox[box]) / (var->T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR), elementWidth);
+    printElement(var->compressibility[box], elementWidth); 
+    printElement(var->enthalpy[box], elementWidth); 
+    #if ENSEMBLE == GEMC
+    printElement(var->heatOfVap, elementWidth);
+    #endif
   }
   if(enableMol) {
     printElement(var->numByBox[box], elementWidth);
@@ -340,6 +344,10 @@ void ConsoleOutput::PrintStatisticTitle()
   if(enablePressure) {
     printElement("PRESSURE", elementWidth);
     printElement("COMPRESSIBILITY", elementWidth);
+    printElement("ENTHALPY", elementWidth);
+    #if ENSEMBLE == GEMC
+    printElement("HEAT_VAP", elementWidth);
+    #endif
   }
 
   if(enableMol) {

src/OutputVars.cpp

@@ -18,6 +18,19 @@ OutputVars::OutputVars(System & sys, StaticVals const& statV, const std::vector<
   T_in_K(statV.forcefield.T_in_K), calc(sys.calcEnergy), molKindNames(molKindNames)
 {
   InitRef(sys, statV);
+  for (int b = 0; b < BOXES_WITH_U_NB; ++b){
+    compressibility[b] = 0.0;
+    enthalpy[b] = 0.0;
+  }
+  #if ENSEMBLE == GEMC
+  liqBox = 0;
+  vapBox = 0;
+  heatOfVap = 0.0;
+  for (int b = 0; b < BOXES_WITH_U_NB; ++b){
+    heatOfVap_energy_term_box[b] = 0.0;
+    heatOfVap_density_term_box[b] = 0.0;
+  }
+  #endif
 }
 
 void OutputVars::InitRef(System & sys, StaticVals const& statV)
@@ -101,13 +114,19 @@ void OutputVars::CalcAndConvert(ulong step)
   molLookupRef->TotalAndDensity(numByBox,  numByKindBox, molFractionByKindBox,
                                 densityByKindBox, volInvRef);
 
-#if ENSEMBLE == GEMC
-  //Determine which box is liquid for purposes of heat of vap.
-  if (densityByKindBox[numKinds] > densityByKindBox[0]) {
-    vapBox = mv::BOX0;
-    liqBox = mv::BOX1;
+  for (uint b = 0; b < BOX_TOTAL; b++) {
+    densityTot[b] = 0.0;
+    for (uint k = 0; k < numKinds; k++) {
+      double density = densityByKindBox[k + numKinds * b];
+
+      // Convert density to g/ml (which is equivalent to g/cm3)
+      // To get kg/m3, multiply output densities by 1000.
+      density *= unit::MOLECULES_PER_A3_TO_MOL_PER_CM3 *
+                 kindsRef[k].molMass;
+      densityTot[b] += density;
+    }
+    densityTot[b] *= 1000;
   }
-#endif
 
   for (uint b = 0; b < BOXES_WITH_U_NB; b++) {
     //Account for dimensionality of virial (raw "virial" is actually a
@@ -160,52 +179,50 @@ void OutputVars::CalcAndConvert(ulong step)
         pressure[b] = rawPressure[b];
         pressure[b] *= unit::K_MOLECULE_PER_A3_TO_BAR;
         if(numByBox[b] != 0){
-          compressability[b] = (pressure[b]) * (volumeRef[b]) / numByBox[b] / (T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR);
+          compressibility[b] = (pressure[b]) * (volumeRef[b]) / numByBox[b] / (T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR);
           enthalpy[b] = (energyRef[b].total / numByBox[b] + rawPressure[b] * volumeRef[b] / numByBox[b]) * UNIT_CONST_H::unit::K_TO_KJ_PER_MOL;
         } else {
-          compressability[b] = 0;
-          enthalpy[b] = 0;
-        }
-#if ENSEMBLE == GEMC
-        // delta Hv = (Uv-Ul) + P(Vv-Vl)
-        if (numByBox[vapBox] != 0){
-          heatOfVap_energy_term_box[vapBox] = energyRef[vapBox].total / numByBox[vapBox];
-          heatOfVap_density_term_box[vapBox] = volumeRef[vapBox] / numByBox[vapBox];
-        } else {
-          heatOfVap_energy_term_box[vapBox] = 0.0;
-          heatOfVap_density_term_box[vapBox] = 0.0;
+          compressibility[b] = 0.0;
+          enthalpy[b] = 0.0;
         }
-
-        if (numByBox[liqBox] != 0){
-          heatOfVap_energy_term_box[liqBox] = energyRef[liqBox].total / numByBox[liqBox];
-          heatOfVap_density_term_box[liqBox] = volumeRef[liqBox] / numByBox[liqBox];
-        } else {
-          heatOfVap_energy_term_box[liqBox] = 0.0;
-          heatOfVap_density_term_box[liqBox] = 0.0;
-        }     
-
-        heatOfVap = heatOfVap_energy_term_box[vapBox] -
-                    heatOfVap_energy_term_box[liqBox];
-        heatOfVap += rawPressure[vapBox] *
-                     (heatOfVap_density_term_box[vapBox] -
-                      heatOfVap_density_term_box[liqBox]);
-        heatOfVap *= unit::K_TO_KJ_PER_MOL;
-#endif
       }
     }
   }
+#if ENSEMBLE == GEMC
+  if (pressureCalc) {
+    if((step + 1) % pCalcFreq == 0 || step == 0) {
+      //Determine which box is liquid for purposes of heat of vap.
+      if (densityTot[mv::BOX1] >= densityTot[mv::BOX0]) {
+        vapBox = mv::BOX0;
+        liqBox = mv::BOX1;
+      } else {
+        vapBox = mv::BOX1;
+        liqBox = mv::BOX0;
+      }
+      // delta Hv = (Uv-Ul) + P(Vv-Vl)
+      if (numByBox[vapBox] != 0){
+        heatOfVap_energy_term_box[vapBox] = energyRef[vapBox].total / numByBox[vapBox];
+        heatOfVap_density_term_box[vapBox] = volumeRef[vapBox] / numByBox[vapBox];
+      } else {
+        heatOfVap_energy_term_box[vapBox] = 0.0;
+        heatOfVap_density_term_box[vapBox] = 0.0;
+      }
 
-  for (uint b = 0; b < BOX_TOTAL; b++) {
-    densityTot[b] = 0.0;
-    for (uint k = 0; k < numKinds; k++) {
-      double density = densityByKindBox[k + numKinds * b];
-
-      // Convert density to g/ml (which is equivalent to g/cm3)
-      // To get kg/m3, multiply output densities by 1000.
-      density *= unit::MOLECULES_PER_A3_TO_MOL_PER_CM3 *
-                 kindsRef[k].molMass;
-      densityTot[b] += density;
+      if (numByBox[liqBox] != 0){
+        heatOfVap_energy_term_box[liqBox] = energyRef[liqBox].total / numByBox[liqBox];
+        heatOfVap_density_term_box[liqBox] = volumeRef[liqBox] / numByBox[liqBox];
+      } else {
+        heatOfVap_energy_term_box[liqBox] = 0.0;
+        heatOfVap_density_term_box[liqBox] = 0.0;
+      }     
+
+      heatOfVap = heatOfVap_energy_term_box[vapBox] -
+                  heatOfVap_energy_term_box[liqBox];
+      heatOfVap += rawPressure[vapBox] *
+                    (heatOfVap_density_term_box[vapBox] -
+                    heatOfVap_density_term_box[liqBox]);
+      heatOfVap *= unit::K_TO_KJ_PER_MOL;
     }
-    densityTot[b] *= 1000;
   }
+#endif
 }

src/OutputVars.h

@@ -40,7 +40,7 @@ class OutputVars
   //Intermediate vars.
   uint * numByBox, * numByKindBox;
   double *molFractionByKindBox, *densityByKindBox,
-      pressure[BOXES_WITH_U_NB], densityTot[BOX_TOTAL], compressability[BOXES_WITH_U_NB], enthalpy[BOXES_WITH_U_NB];
+      pressure[BOXES_WITH_U_NB], densityTot[BOX_TOTAL], compressibility[BOXES_WITH_U_NB], enthalpy[BOXES_WITH_U_NB];
   double pressureTens[BOXES_WITH_U_NB][3][3];
   double surfaceTens[BOXES_WITH_U_NB];
   ulong pCalcFreq;

test/Run_Examples.py

@@ -0,0 +1,214 @@
+#!/usr/bin/python
+
+HEADER = '\033[95m'
+OKBLUE = '\033[94m'
+OKGREEN = '\033[92m'
+WARNING = '\033[93m'
+FAIL = '\033[91m'
+ENDC = '\033[0m'
+BOLD = '\033[1m'
+UNDERLINE = '\033[4m'
+
+def substring_after(s, delim):
+    return s.partition(delim)[2]
+
+import os
+import glob
+import pandas as pd
+import time
+import datetime
+import subprocess
+from subprocess import STDOUT
+from filecmp import cmp
+binaries_dict = {}
+GPU_binaries_dict = {}
+
+os.chdir("new_binaries")
+pathToDir = os.getcwd()
+binaries_cpu_new = sorted(glob.glob('GOMC_CPU_*'), key=os.path.getmtime)
+binaries_gpu_new = sorted(glob.glob('GOMC_GPU_*'), key=os.path.getmtime)
+for binary in binaries_cpu_new:
+    binaries_dict[substring_after(binary, "GOMC_")+"_new"] = os.path.join(pathToDir, binary)
+for binary in binaries_gpu_new:
+    binaries_dict[substring_after(binary, "GOMC_")+"_new"] = os.path.join(pathToDir, binary)
+
+os.chdir("../ref_binaries")
+pathToDir = os.getcwd()
+binaries_cpu_ref = sorted(glob.glob('GOMC_CPU_*'), key=os.path.getmtime)
+binaries_gpu_ref = sorted(glob.glob('GOMC_GPU_*'), key=os.path.getmtime)
+for binary in binaries_cpu_ref:
+    binaries_dict[substring_after(binary, "GOMC_")+"_ref"] = os.path.join(pathToDir, binary)
+for binary in binaries_gpu_ref:
+    binaries_dict[substring_after(binary, "GOMC_")+"_ref"] = os.path.join(pathToDir, binary)
+
+print("Binaries Dict", binaries_dict)
+
+os.chdir("../integration")
+
+confFileName = "in.conf"
+
+Log_Template_file = open("IntegrationTest.log", 'w')
+print("opened Log")
+
+listOfTests = []
+
+# traverse root directory, and list directories as dirs and files as files
+# traverse root directory, and list directories as dirs and files as files
+for root, dirs, files in os.walk("."):
+    path = root.split(os.sep)
+#    print((len(path) - 1) * '---', os.path.basename(root))
+    for file in files:
+#        print(len(path) * '---', file)
+        if file==confFileName:
+            newOrRef = ""
+            cpuOrGpu = ""
+            if "new_cpu" in path:
+                newOrRef = "_new"
+                cpuOrGpu = "CPU_"
+            elif "ref_cpu" in path:
+                newOrRef = "_ref"
+                cpuOrGpu = "CPU_"
+            elif "new_gpu" in path:
+                newOrRef = "_new"
+                cpuOrGpu = "GPU_"
+            elif "ref_gpu" in path:
+                newOrRef = "_ref"
+                cpuOrGpu = "GPU_"
+
+            if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'NVT_GEMC' not in path:
+                print("Call GOMC")
+                command = binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root)
+                print(binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root))
+                listOfTests.append(command)
+            elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'NPT_GEMC' not in path:
+                print("Call GOMC")
+                print(binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root))
+                command = binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root)
+                listOfTests.append(command)
+            elif cpuOrGpu+"GCMC"+newOrRef in binaries_dict and 'GCMC' in path:
+                print("Call GOMC")
+                print(binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,cpuOrGpu,newOrRef,"GCMC_"+os.path.basename(root))
+                command = binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,cpuOrGpu,newOrRef,"GCMC_"+os.path.basename(root)
+                listOfTests.append(command)
+            elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'NVT_GEMC' in path:
+                print("Call GOMC")
+                print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root))
+                command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NVT_GEMC_"+os.path.basename(root)
+                listOfTests.append(command)
+            elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'NPT_GEMC' in path:
+                print("Call GOMC")
+                print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root))
+                command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root)
+                listOfTests.append(command)
+print(listOfTests)
+# Create the pandas DataFrame
+colNames = ['PathToBinary', 'PathToExample', 'Binary', 'CPU_or_GPU','New_or_Ref','Example']
+df = pd.DataFrame(listOfTests, columns = colNames)
+df = df.sort_values(by=['Example'])
+print(df)
+#for index, row in df.iterrows():
+#    print(index)
+#    print(row)
+grouped = df.groupby(['Example'])
+all_examples = df['Example'].unique()
+CPU_v_GPU_global = True
+New_v_Ref_global = True
+cross_bool_global = True
+CPU_v_GPU_exists_global = False
+New_v_Ref_exists_global = False
+cross_exists_global = False
+for example in all_examples:
+    ex_df = grouped.get_group(example)
+    for index, row in ex_df.iterrows():
+        #print("run file {}".format(row['PathToExample']+" conf"))
+        os.chdir(row['PathToExample'])
+        write_log_data = "Changing directory to {}\n".format(row['PathToExample'])
+        Log_Template_file.write(str(write_log_data))
+        command = (row['PathToBinary'] + " in.conf > out.log")
+        write_log_data = "Issuing command: {}\n".format(command)
+        Log_Template_file.write(str(write_log_data))
+        start = time.time()
+        exec_GOMC_run_command = subprocess.Popen(command, stderr=STDOUT)
+        write_log_data = "Waiting for GOMC Example {} {} to finish.\n".format(row['Binary'],row['Example'])
+        print(str(write_log_data))
+        Log_Template_file.write(str(write_log_data))
+        GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED)  # pauses python until box 0 sim done
+        end = time.time()
+        write_log_data = "Elapsed time: {}.\n".format(datetime.timedelta(seconds=end-start))
+        Log_Template_file.write(str(write_log_data))
+        print(str(write_log_data))
+        write_log_data = "The GOMC Example {} {} has finished.\n".format(row['Binary'],row['Example'])
+        print(str(write_log_data))
+        Log_Template_file.write(str(write_log_data))
+        Log_Template_file.flush()
+
+    # Create a list of the PDB files in this example
+    full_path_pdb_files = sorted(glob.glob(os.path.join(row['PathToExample'],'*.pdb')), key=os.path.getmtime)
+    just_file_names = []
+    for path in full_path_pdb_files:
+        just_file_names.append(os.path.basename(path))
+    print(just_file_names)
+    cross = ex_df.merge(ex_df, on=['Example'],how='outer')
+    print('---', ex_df['Example'].iloc[0])
+    CPU_v_GPU = True
+    New_v_Ref = True
+    cross_bool = True
+    CPU_v_GPU_exists = False
+    New_v_Ref_exists = False
+    cross_exists = False
+    for pdb_file in just_file_names:
+        print(2 * '---', pdb_file)
+        my_tuples = []
+        for index, row in cross.iterrows():
+            f1 = os.path.join(row['PathToExample_x'],pdb_file)
+            f2 = os.path.join(row['PathToExample_y'],pdb_file)
+            if ((row['CPU_or_GPU_x'] != row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] == row['New_or_Ref_y'])):
+                CPU_v_GPU_exists = True
+                CPU_v_GPU_exists_global = True
+                result = cmp(f1, f2, shallow=False)
+                CPU_v_GPU = CPU_v_GPU and result
+                CPU_v_GPU_global = CPU_v_GPU_global and result
+            elif ((row['CPU_or_GPU_x'] == row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] != row['New_or_Ref_y'])):
+                New_v_Ref_exists = True
+                New_v_Ref_exists_global = True
+                result = cmp(f1, f2, shallow=False)
+                New_v_Ref = New_v_Ref and result
+                New_v_Ref_global = New_v_Ref_global and result
+            elif ((row['CPU_or_GPU_x'] != row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] != row['New_or_Ref_y'])):
+                cross_exists = True
+                cross_exists_global = True
+                result = cmp(f1, f2, shallow=False)
+                cross_bool = cross_bool and result
+                cross_bool_global = cross_bool_global and result
+    if(CPU_v_GPU_exists):
+        if(CPU_v_GPU):
+            print((3 * '---')+"CPU_v_GPU: "+ OKGREEN + "PASS" + ENDC)
+        else:
+            print((3 * '---')+"CPU_v_GPU: "+ FAIL + "FAIL" + ENDC)
+    if(New_v_Ref_exists):
+        if(New_v_Ref):
+            print((3 * '---')+"New vs Ref: "+ OKGREEN + "PASS" + ENDC)
+        else:
+            print((3 * '---')+"New vs Ref: "+ FAIL + "FAIL" + ENDC)
+    if(cross_exists):
+        if(cross_bool):
+            print((3 * '---')+"CPU vs GPU X New vs Ref: "+ OKGREEN + "PASS" + ENDC)
+        else:
+            print((3 * '---')+"CPU vs GPU X New vs Ref: "+ FAIL + "FAIL" + ENDC)
+
+
+if(CPU_v_GPU_exists_global):
+    if(CPU_v_GPU_global):
+        print("CPU_v_GPU Global: "+ OKGREEN + "PASS" + ENDC)
+    else:
+        print("CPU_v_GPU Global: "+ FAIL + "FAIL" + ENDC)
+if(New_v_Ref_exists_global):
+    if(New_v_Ref_global):
+        print("New vs Ref Global: "+ OKGREEN + "PASS" + ENDC)
+    else:
+        print("New vs Ref Global: "+ FAIL + "FAIL" + ENDC)
+if(cross_exists_global):
+    if(cross_bool_global):
+        print("CPU vs GPU X New vs Ref Global: "+ OKGREEN + "PASS" + ENDC)
+    else:
+        print("CPU vs GPU X New vs Ref Global: "+ FAIL + "FAIL" + ENDC)