Merge pull request #67 from Zhaoli2042/push-branch

Push branch: Widom insertion adsorption energy

Author: Xiaoyi-ZHANG23

src_clean/axpy.cu

@@ -170,13 +170,18 @@ void RunMoves(Variables& Vars, size_t box_index, int Cycle)
     double2& newScale = SystemComponents.TempVal.Scale; 
     newScale = SystemComponents.Lambda[comp].SET_SCALE(1.0); //Set scale for full molecule (lambda = 1.0)//
     double Rosenbluth = MOVES.INSERTION.WidomMove(Vars, box_index);
+    //Also record move energy (delta energy)
+    //MoveEnergy widom_e = MOVES.INSERTION.energy;
+    //printf("Widom E: "); widom_e.print();
 
     if(Vars.SimulationMode == PRODUCTION)
     {
       size_t blockID = Cycle/Vars.BlockAverageSize;
       if(blockID >= SystemComponents.Nblock) blockID --;
       SystemComponents.Moves[comp].BlockID = blockID;
       SystemComponents.Moves[comp].RecordRosen(Rosenbluth, WIDOM);
+      //weight with Rosenbluth (heavy weight on low (attractive) energy, lower on high energy)
+      SystemComponents.Moves[comp].Rosen[blockID].widom_energy += MOVES.INSERTION.energy * Rosenbluth;
     }
   }
   else if(RANDOMNUMBER < SystemComponents.Moves[comp].ReinsertionProb)

src_clean/data_struct.h

@@ -35,7 +35,7 @@ enum INTERACTION_TYPES {HH = 0, HG, GG};
 enum RESTART_FILE_TYPES {RASPA_RESTART = 0, LAMMPS_DATA};
 
 //Zhao's note: For the stage of evaluating total energy of the system//
-enum ENERGYEVALSTAGE {INITIAL = 0, CREATEMOL, FINAL, CREATEMOL_DELTA, DELTA, CREATEMOL_DELTA_CHECK, DELTA_CHECK, DRIFT, GPU_DRIFT, AVERAGE, AVERAGE_ERR};
+enum ENERGYEVALSTAGE {INITIAL = 0, CREATEMOL, FINAL, CREATEMOL_DELTA, DELTA, CREATEMOL_DELTA_CHECK, DELTA_CHECK, DRIFT, GPU_DRIFT, AVERAGE, AVERAGE_ERR, WIDOM_AVG, WIDOM_ERR};
 
 struct EnergyComplex
 {
@@ -409,6 +409,82 @@ struct TMMC
   }
 };
 
+//Temporary energies for a move//
+//Zhao's note: consider making this the default return variable for moves, like RASPA-3?//
+struct MoveEnergy
+{
+  double storedHGVDW=0.0;
+  double storedHGReal=0.0;
+  double storedHGEwaldE=0.0;
+  // van der Waals //
+  double HHVDW=0.0;
+  double HGVDW=0.0;
+  double GGVDW=0.0;
+  // Real Part of Coulomb //
+  double HHReal=0.0;
+  double HGReal=0.0;
+  double GGReal=0.0;
+  // Long-Range Ewald Energy //
+  double HHEwaldE=0.0;
+  double HGEwaldE=0.0;
+  double GGEwaldE=0.0;
+  // Other Energies //
+  double TailE =0.0;
+  double DNN_E =0.0;
+  double total()
+  {
+    return HHVDW + HGVDW + GGVDW + 
+           HHReal + HGReal + GGReal + 
+           HHEwaldE + HGEwaldE + GGEwaldE + 
+           TailE + DNN_E;
+  };
+  void take_negative()
+  {
+    storedHGVDW *= -1.0;
+    storedHGReal*= -1.0;
+    storedHGEwaldE  *= -1.0;
+    HHVDW *= -1.0; HHReal *= -1.0;
+    HGVDW *= -1.0; HGReal *= -1.0;
+    GGVDW *= -1.0; GGReal *= -1.0;
+    HHEwaldE *= -1.0; HGEwaldE *= -1.0; GGEwaldE  *= -1.0;
+    TailE     *= -1.0;
+    DNN_E     *= -1.0;
+  };
+  void zero()
+  {
+    storedHGVDW =0.0;
+    storedHGReal=0.0;
+    storedHGEwaldE=0.0;
+    HHVDW=0.0; HHReal=0.0;
+    HGVDW=0.0; HGReal=0.0;
+    GGVDW=0.0; GGReal=0.0;
+    HHEwaldE =0.0;
+    HGEwaldE =0.0;
+    GGEwaldE =0.0;
+    TailE    =0.0;
+    DNN_E    =0.0;
+  };
+  void print()
+  {
+    printf("HHVDW: %.5f, HHReal: %.5f, HGVDW: %.5f, HGReal: %.5f, GGVDW: %.5f, GGReal: %.5f, HHEwaldE: %.5f,\n HGEwaldE: %.5f,\n GGEwaldE: %.5f, TailE: %.5f, DNN_E: %.5f\n", HHVDW, HHReal, HGVDW, HGReal, GGVDW, GGReal, HHEwaldE, HGEwaldE, GGEwaldE, TailE, DNN_E);
+    printf("Stored HGVDW: %.5f, Stored HGReal: %.5f, Stored HGEwaldE: %.5f\n", storedHGVDW, storedHGReal, storedHGEwaldE);
+  };
+  void DNN_Replace_Energy()
+  {
+    storedHGVDW = HGVDW;
+    storedHGReal= HGReal;
+    storedHGEwaldE  = HGEwaldE;
+    HGVDW = 0.0;
+    HGReal= 0.0;
+    HGEwaldE = 0.0;
+  }
+  double DNN_Correction() //Using DNN energy to replace HGVDW, HGReal and HGEwaldE//
+  {
+    double correction = DNN_E - storedHGVDW - storedHGReal - storedHGEwaldE;
+    return correction;
+  };
+};
+
 //Zhao's note: keep track of the Rosenbluth weights during the simulation//
 struct RosenbluthWeight
 {
@@ -417,6 +493,7 @@ struct RosenbluthWeight
   double3 Insertion      = {0.0, 0.0, 0.0};
   double3 Deletion       = {0.0, 0.0, 0.0};
   //NOTE: DO NOT RECORD FOR REINSERTIONS, SINCE THE DELETION PART OF REINSERTION IS MODIFIED//
+  MoveEnergy widom_energy;
 };
 
 struct Move_Statistics
@@ -483,6 +560,10 @@ struct Move_Statistics
   std::vector<std::vector<double>>MolSQPerComponent;
   //x: average; y: average^2; z: Number of Widom insertion performed//
   std::vector<RosenbluthWeight>Rosen; //vector over Nblocks//
+
+  MoveEnergy WidomEnergy;
+  MoveEnergy WidomEnergy_ERR;
+
   void NormalizeProbabilities()
   {
     //Zhao's note: the probabilities here are what we defined in simulation.input, raw values//
@@ -640,81 +721,7 @@ struct Tail
   double Energy  = {0.0};
 };
 
-//Temporary energies for a move//
-//Zhao's note: consider making this the default return variable for moves, like RASPA-3?//
-struct MoveEnergy
-{
-  double storedHGVDW=0.0;
-  double storedHGReal=0.0;
-  double storedHGEwaldE=0.0;
-  // van der Waals //
-  double HHVDW=0.0;
-  double HGVDW=0.0;
-  double GGVDW=0.0;
-  // Real Part of Coulomb //
-  double HHReal=0.0;
-  double HGReal=0.0;
-  double GGReal=0.0;
-  // Long-Range Ewald Energy //
-  double HHEwaldE=0.0;
-  double HGEwaldE=0.0;
-  double GGEwaldE=0.0;
-  // Other Energies //
-  double TailE =0.0;
-  double DNN_E =0.0;
-  double total()
-  {
-    return HHVDW + HGVDW + GGVDW + 
-           HHReal + HGReal + GGReal + 
-           HHEwaldE + HGEwaldE + GGEwaldE + 
-           TailE + DNN_E;
-  };
-  void take_negative()
-  {
-    storedHGVDW *= -1.0;
-    storedHGReal*= -1.0;
-    storedHGEwaldE  *= -1.0;
-    HHVDW *= -1.0; HHReal *= -1.0;
-    HGVDW *= -1.0; HGReal *= -1.0;
-    GGVDW *= -1.0; GGReal *= -1.0;
-    HHEwaldE *= -1.0; HGEwaldE *= -1.0; GGEwaldE  *= -1.0;
-    TailE     *= -1.0;
-    DNN_E     *= -1.0;
-  };
-  void zero()
-  {
-    storedHGVDW =0.0;
-    storedHGReal=0.0;
-    storedHGEwaldE=0.0;
-    HHVDW=0.0; HHReal=0.0;
-    HGVDW=0.0; HGReal=0.0;
-    GGVDW=0.0; GGReal=0.0;
-    HHEwaldE =0.0;
-    HGEwaldE =0.0;
-    GGEwaldE =0.0;
-    TailE    =0.0;
-    DNN_E    =0.0;
-  };
-  void print()
-  {
-    printf("HHVDW: %.5f, HHReal: %.5f, HGVDW: %.5f, HGReal: %.5f, GGVDW: %.5f, GGReal: %.5f, HHEwaldE: %.5f,\n HGEwaldE: %.5f,\n GGEwaldE: %.5f, TailE: %.5f, DNN_E: %.5f\n", HHVDW, HHReal, HGVDW, HGReal, GGVDW, GGReal, HHEwaldE, HGEwaldE, GGEwaldE, TailE, DNN_E);
-    printf("Stored HGVDW: %.5f, Stored HGReal: %.5f, Stored HGEwaldE: %.5f\n", storedHGVDW, storedHGReal, storedHGEwaldE);
-  };
-  void DNN_Replace_Energy()
-  {
-    storedHGVDW = HGVDW;
-    storedHGReal= HGReal;
-    storedHGEwaldE  = HGEwaldE;
-    HGVDW = 0.0;
-    HGReal= 0.0;
-    HGEwaldE = 0.0;
-  }
-  double DNN_Correction() //Using DNN energy to replace HGVDW, HGReal and HGEwaldE//
-  {
-    double correction = DNN_E - storedHGVDW - storedHGReal - storedHGEwaldE;
-    return correction;
-  };
-};
+
 /*
 __host__ MoveEnergy MoveEnergy_Multiply(MoveEnergy A, MoveEnergy B)
 {

src_clean/fxn_main.h

@@ -442,7 +442,7 @@ inline void Copy_AtomData_from_Device(Atoms* System, Atoms* Host_System, Atoms*
   }
 }
 
-inline void PRINT_ENERGY_AT_STAGE(Components& SystemComponents, int stage, Units& Constants)
+inline void PRINT_ENERGY_AT_STAGE(Components& SystemComponents, int stage, Units& Constants, size_t comp)
 {
   std::string stage_name;
   MoveEnergy  E;
@@ -459,6 +459,8 @@ inline void PRINT_ENERGY_AT_STAGE(Components& SystemComponents, int stage, Units
     case GPU_DRIFT: {stage_name = "GPU DRIFT (GPU FINAL - CPU FINAL)"; E = SystemComponents.Final_Energy - SystemComponents.GPU_Energy; break;}
     case AVERAGE: {stage_name = "PRODUCTION PHASE AVERAGE ENERGY"; E = SystemComponents.AverageEnergy + SystemComponents.Initial_Energy;break;}
     case AVERAGE_ERR: {stage_name = "PRODUCTION PHASE AVERAGE ENERGY ERRORBAR"; E = SystemComponents.AverageEnergy_Errorbar; break;}
+    case WIDOM_AVG: {stage_name = "AVERAGE WIDOM ADSORPTION ENERGY for " + SystemComponents.MoleculeName[comp]; E = SystemComponents.Moves[comp].WidomEnergy;break;}
+    case WIDOM_ERR: {stage_name = "WIDOM ADSORPTION ENERGY ERRORBAR for " + SystemComponents.MoleculeName[comp]; E = SystemComponents.Moves[comp].WidomEnergy_ERR;break;}
   }
   fprintf(SystemComponents.OUTPUT, " *** %s *** \n", stage_name.c_str());
   fprintf(SystemComponents.OUTPUT, "========================================================================\n");
@@ -501,19 +503,28 @@ void ENERGY_SUMMARY(Variables& Vars)
   for(size_t i = 0; i < NumberOfSimulations; i++)
   {
     fprintf(SystemComponents[i].OUTPUT, "======================== ENERGY SUMMARY (Simulation %zu) =========================\n", i);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], INITIAL, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA_CHECK, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], FINAL, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA_CHECK, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DRIFT, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], GPU_DRIFT, Constants);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], INITIAL, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA_CHECK, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], FINAL, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA_CHECK, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], DRIFT, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], GPU_DRIFT, Constants, 0);
     fprintf(SystemComponents[i].OUTPUT, "================================================================================\n");
     fprintf(SystemComponents[i].OUTPUT, "======================== PRODUCTION PHASE AVERAGE ENERGIES (Simulation %zu) =========================\n", i);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE, Constants);
-    PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE_ERR, Constants);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE, Constants, 0);
+    PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE_ERR, Constants, 0);
+    fprintf(SystemComponents[i].OUTPUT, "================================================================================\n");
+    
+    fprintf(SystemComponents[i].OUTPUT, "======================== WIDOM ADSORPTION ENERGIES (Simulation %zu) =========================\n", i);
+    for(size_t j = 1; j < SystemComponents[i].NComponents.x; j++)
+    {
+      if(SystemComponents[i].Moves[j].WidomProb < 1e-10) continue;
+      PRINT_ENERGY_AT_STAGE(SystemComponents[i], WIDOM_AVG, Constants, j);
+      PRINT_ENERGY_AT_STAGE(SystemComponents[i], WIDOM_ERR, Constants, j);
+    }
     fprintf(SystemComponents[i].OUTPUT, "================================================================================\n");
 
     fprintf(SystemComponents[i].OUTPUT, "DNN Rejection Summary:\nTranslation+Rotation: %zu\nReinsertion: %zu\nInsertion: %zu\nDeletion: %zu\nSingleSwap: %zu\n", SystemComponents[i].TranslationRotationDNNReject, SystemComponents[i].ReinsertionDNNReject, SystemComponents[i].InsertionDNNReject, SystemComponents[i].DeletionDNNReject, SystemComponents[i].SingleSwapDNNReject);

src_clean/print_statistics.cuh

@@ -89,17 +89,21 @@ static inline void ComputeFugacity(double& AverageWr, double& AverageMu, double&
 
 static inline void Print_Widom_Statistics(Components& SystemComponents, Boxsize Box, Units& Constants, size_t comp)
 {
-  Move_Statistics MoveStats = SystemComponents.Moves[comp];
+  Move_Statistics& MoveStats = SystemComponents.Moves[comp];
   double2 totRosen = {0.0, 0.0};
   double2 totMu    = {0.0, 0.0};
   double2 totHenry = {0.0, 0.0};
-  double2 totFuga  = {0.0, 0.0};
+  double2 totFuga  = {0.0, 0.0};    
+  MoveEnergy AverageEnergy;
+  MoveEnergy AverageEnergy_SQ;
+
   fprintf(SystemComponents.OUTPUT, "=====================Rosenbluth Summary For Component [%zu] (%s)=====================\n", comp, SystemComponents.MoleculeName[comp].c_str());
   fprintf(SystemComponents.OUTPUT, "There are %zu blocks\n", SystemComponents.Nblock);
   for(size_t i = 0; i < SystemComponents.Nblock; i++)
   {
     fprintf(SystemComponents.OUTPUT, "=====BLOCK %zu=====\n", i);
     fprintf(SystemComponents.OUTPUT, "Widom Performed: %.1f\n", MoveStats.Rosen[i].Total.z);
+    double OverallAvgWr = 0.0;
     if(MoveStats.Rosen[i].Total.z > 0)
     {
       double AverageWr = 0.0; double AverageMu = 0.0; double AverageHenry = 0.0; double Fugacity = 0.0;
@@ -112,6 +116,8 @@ static inline void Print_Widom_Statistics(Components& SystemComponents, Boxsize
       totMu.x    += AverageMu; totMu.y    += AverageMu * AverageMu;
       totHenry.x += AverageHenry; totHenry.y += AverageHenry * AverageHenry;
       totFuga.x  += Fugacity;  totFuga.y  += Fugacity * Fugacity;
+
+      OverallAvgWr = AverageWr;
     }
     if(MoveStats.Rosen[i].Insertion.z > 0)
     {
@@ -129,7 +135,21 @@ static inline void Print_Widom_Statistics(Components& SystemComponents, Boxsize
       fprintf(SystemComponents.OUTPUT, "(Deletion) Averaged Excess Mu: %.10f\n", AverageMu);
       fprintf(SystemComponents.OUTPUT, "(Deletion) Converted to Fugacity: %.10f [Pascal], Temp: %.5f [K]\n", Fugacity, SystemComponents.Temperature);
     }
+
+    //Add Widom Delta Energy
+    //Calculate just the overall, now//
+    if(MoveStats.Rosen[i].Total.z > 0)
+    {
+      printf("Raw WIDOM %zu", i); MoveStats.Rosen[i].widom_energy.print();
+      MoveEnergy Average = MoveStats.Rosen[i].widom_energy / static_cast<double>(MoveStats.Rosen[i].Total.z) / OverallAvgWr;
+      printf("AVG WIDOM %zu", i); Average.print();
+      AverageEnergy += Average / static_cast<double>(SystemComponents.Nblock);
+      AverageEnergy_SQ += Average * Average / static_cast<double>(SystemComponents.Nblock);
+    }
   }
+  MoveEnergy AverageEnergy_Errorbar = sqrt_MoveEnergy(AverageEnergy_SQ - AverageEnergy * AverageEnergy) * 2.0;
+  SystemComponents.Moves[comp].WidomEnergy = AverageEnergy;
+  SystemComponents.Moves[comp].WidomEnergy_ERR = AverageEnergy_Errorbar;
 
   double2 AvgBlockRosen = {0.0, 0.0};
   AvgBlockRosen.x = totRosen.x / (double) SystemComponents.Nblock;