Merge pull request #61 from Zhaoli2042/main

PR: small change of code and initial pybind mix

Author: Xiaoyi-ZHANG23

README.md

@@ -28,6 +28,11 @@ A detailed installation note for gRASPA on CentOS/Ubuntu 24.04 is documented in
   * A40, A100, RTX 3080 Ti, RTX 3090, RTX 4090.
   * 🤯: RTX 3090/4090 is faster than A40/A100 for gRASPA
 * gRASPA has a SYCL version (experimental) that supports other devices, available in [Releases](https://github.com/snurr-group/gRASPA/releases)
+
+### Pybind Extension (testing)
+* Pybind extension for gRASPA allows user to interact with the **internal variables** of gRASPA, break down **MC moves**, add their **modifications**
+* Access the pybind-gRASPA extension [here](https://github.com/Zhaoli2042/gRASPA_pybind), as a patch to the original code
+
 ## Quick Start
 * Go to [```Examples/```](Examples/) folder and read more!
 
@@ -75,11 +80,11 @@ A detailed installation note for gRASPA on CentOS/Ubuntu 24.04 is documented in
 | Automatic Determination<br>of # unit cells | | | :heavy_check_mark: |
 
 ## Authors
-* Zhao Li (Northwestern University, currently at Purdue University)
-* Kaihang Shi (Northwestern University, currently at University at Buffalo)
-* David Dubbeldam (University of Amsterdam)
-* Mark Dewing (Argonne National Laboratory)
-* Christopher Knight (Argonne National Laboratory)
-* Alvaro Vazquez Mayagoitia (Argonne National Laboratory)
-* Randall Q. Snurr (Northwestern University)
+* Zhao Li (Northwestern University, currently at Purdue University/University of Notre Dame, [email protected])
+* Kaihang Shi (Northwestern University, currently at University at Buffalo, [email protected])
+* David Dubbeldam (University of Amsterdam, [email protected])
+* Mark Dewing (Argonne National Laboratory, [email protected])
+* Christopher Knight (Argonne National Laboratory, [email protected])
+* Alvaro Vazquez Mayagoitia (Argonne National Laboratory, [email protected])
+* Randall Q. Snurr (Northwestern University, [email protected])
 

src_clean/DNN_HostGuest_Energy_Functions.h

@@ -10,8 +10,15 @@
 #include "read_data.h"
 //###PATCH_LCLIN_INCLUDE_HEADER###//
 
-__global__ void Initialize_DNN_Positions(Atoms* d_a, Atoms New, Atoms Old, size_t Oldsize, size_t Newsize, size_t SelectedComponent, size_t Location, size_t chainsize, int MoveType, size_t CYCLE)
+__global__ void Initialize_DNN_Positions(Atoms* d_a, Atoms New, Atoms Old, double3* temp, size_t Oldsize, size_t Newsize, size_t SelectedComponent, size_t Location, size_t chainsize, int MoveType, size_t CYCLE)
 {
+  size_t ij = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if(ij < (Newsize + Oldsize))
+  {
+    Initialize_Copy_Positions_Together(d_a, New, Old, temp, Oldsize, Newsize, SelectedComponent, Location, chainsize, MoveType);
+  }
+  /*
   //Zhao's note: need to think about changing this boolean to switch//
   if(MoveType == TRANSLATION || MoveType == ROTATION || MoveType == SINGLE_INSERTION || MoveType == SINGLE_DELETION) // Translation/Rotation/single_insertion/single_deletion //
   {
@@ -53,6 +60,8 @@ __global__ void Initialize_DNN_Positions(Atoms* d_a, Atoms New, Atoms Old, size_
       Old.scaleCoul[i]     = d_a[SelectedComponent].scaleCoul[Location + i];
     }
   }
+  */
+
   /*
   if(CYCLE == 145) 
   {
@@ -62,7 +71,7 @@ __global__ void Initialize_DNN_Positions(Atoms* d_a, Atoms New, Atoms Old, size_
   */
 }
 
-void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, Components& SystemComponents, size_t SelectedComponent, int MoveType, size_t Location)
+void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, double3* temp, Components& SystemComponents, size_t SelectedComponent, int MoveType, size_t Location)
 {
   size_t Oldsize = 0; size_t Newsize = 0; size_t chainsize = 0;
   switch(MoveType)
@@ -90,8 +99,11 @@ void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, Component
     }
     case REINSERTION: // Reinsertion //
     {
-      throw std::runtime_error("Use the Special Function for Reinsertion");
-      //break;
+      Oldsize   = SystemComponents.Moleculesize[SelectedComponent];
+      Newsize   = SystemComponents.Moleculesize[SelectedComponent];
+      chainsize = SystemComponents.Moleculesize[SelectedComponent];
+      //throw std::runtime_error("Use the Special Function for Reinsertion");
+      break;
     }
     case IDENTITY_SWAP:
     {
@@ -120,7 +132,9 @@ void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, Component
       break;
     }
   }
-  Initialize_DNN_Positions<<<1,1>>>(d_a, New, Old, Oldsize, Newsize, SelectedComponent, Location, chainsize, MoveType, SystemComponents.CURRENTCYCLE);
+  //Initialize_DNN_Positions<<<1,1>>>(d_a, New, Old, Oldsize, Newsize, SelectedComponent, Location, chainsize, MoveType, SystemComponents.CURRENTCYCLE);
+  size_t Nblock = 0; size_t Nthread = 0; Setup_threadblock(Oldsize + Newsize, Nblock, Nthread);
+  Initialize_DNN_Positions<<<Nblock,Nthread>>>(d_a, New, Old, SystemComponents.tempMolStorage, Oldsize, Newsize, SelectedComponent, Location, chainsize, MoveType, SystemComponents.CURRENTCYCLE);
 }
 
 __global__ void Initialize_DNN_Positions_Reinsertion(double3* temp, Atoms* d_a, Atoms Old, size_t Oldsize, size_t Newsize, size_t realpos, size_t SelectedComponent)

src_clean/Ewald_Energy_Functions.h

@@ -102,46 +102,71 @@ __device__ void Initialize_Vectors_thread(Complex* eik, size_t numberOfAtoms, in
     eik[i + k * numberOfAtoms] = multiply(eik[i + (k - 1) * numberOfAtoms], eik[i + 1 * numberOfAtoms]);
   }
 }
-
-__global__ void Initialize_WaveVector_General(Boxsize Box, int3 kmax, Atoms* d_a, Atoms New, Atoms Old, size_t Oldsize, size_t Newsize, size_t SelectedComponent, size_t Location, size_t chainsize, size_t numberOfAtoms, int MoveType)
+//Copy Old positions and new positions together into Old//
+//including double3* temp, for reinsertion//
+__device__ void Initialize_Copy_Positions_Together(Atoms*& d_a, Atoms& New, Atoms& Old, double3* temp, size_t Oldsize, size_t Newsize, size_t SelectedComponent, size_t Location, size_t chainsize, int MoveType)
 {
-  //Zhao's note: need to think about changing this boolean to switch//
   size_t ij = blockIdx.x * blockDim.x + threadIdx.x;
-  if(ij < (Newsize + Oldsize))
+  if(MoveType == TRANSLATION || MoveType == ROTATION || MoveType == SPECIAL_ROTATION || MoveType == SINGLE_INSERTION || MoveType == SINGLE_DELETION) // Translation/Rotation/single_insertion/single_deletion //
   {
-    if(MoveType == TRANSLATION || MoveType == ROTATION || MoveType == SPECIAL_ROTATION || MoveType == SINGLE_INSERTION || MoveType == SINGLE_DELETION) // Translation/Rotation/single_insertion/single_deletion //
+    //For Translation/Rotation, the Old positions are already in the Old struct, just need to put the New positions into Old, after the Old positions//
+    if(ij >= Oldsize)
     {
-      //For Translation/Rotation, the Old positions are already in the Old struct, just need to put the New positions into Old, after the Old positions//
-      if(ij >= Oldsize)
-      {
-        Old.pos[ij]           = New.pos[ij - Oldsize];
-        Old.scale[ij]         = New.scale[ij - Oldsize];
-        Old.charge[ij]        = New.charge[ij - Oldsize];
-        Old.scaleCoul[ij]     = New.scaleCoul[ij - Oldsize];
-      }
+      Old.pos[ij]           = New.pos[ij - Oldsize];
+      Old.scale[ij]         = New.scale[ij - Oldsize];
+      Old.charge[ij]        = New.charge[ij - Oldsize];
+      Old.scaleCoul[ij]     = New.scaleCoul[ij - Oldsize];
     }
-    else if(MoveType == INSERTION || MoveType == CBCF_INSERTION) // Insertion & Fractional Insertion //
+  }
+  else if(MoveType == INSERTION || MoveType == CBCF_INSERTION) // Insertion & Fractional Insertion //
+  {
+    //Put the trial orientations in New to Old, right after the first bead position//
+    if(ij < chainsize)
     {
-      //Put the trial orientations in New to Old, right after the first bead position//
-      if(ij < chainsize)
-      {
-        Old.pos[ij + 1]       = New.pos[Location * chainsize + ij];
-        Old.scale[ij + 1]     = New.scale[Location * chainsize + ij];
-        Old.charge[ij + 1]    = New.charge[Location * chainsize + ij];
-        Old.scaleCoul[ij + 1] = New.scaleCoul[Location * chainsize + ij];
-      }
+      Old.pos[ij + 1]       = New.pos[Location * chainsize + ij];
+      Old.scale[ij + 1]     = New.scale[Location * chainsize + ij];
+      Old.charge[ij + 1]    = New.charge[Location * chainsize + ij];
+      Old.scaleCoul[ij + 1] = New.scaleCoul[Location * chainsize + ij];
     }
-    else if(MoveType == DELETION || MoveType == CBCF_DELETION) // Deletion //
+  }
+  else if(MoveType == DELETION || MoveType == CBCF_DELETION) // Deletion //
+  {
+    if(ij < Oldsize)
     {
-      if(ij < Oldsize)
-      {
-        // For deletion, Location = UpdateLocation, see Deletion Move //
-        Old.pos[ij]           = d_a[SelectedComponent].pos[Location + ij];
-        Old.scale[ij]         = d_a[SelectedComponent].scale[Location + ij];
-        Old.charge[ij]        = d_a[SelectedComponent].charge[Location + ij];
-        Old.scaleCoul[ij]     = d_a[SelectedComponent].scaleCoul[Location + ij];
-      }
+      // For deletion, Location = UpdateLocation, see Deletion Move //
+      Old.pos[ij]           = d_a[SelectedComponent].pos[Location + ij];
+      Old.scale[ij]         = d_a[SelectedComponent].scale[Location + ij];
+      Old.charge[ij]        = d_a[SelectedComponent].charge[Location + ij];
+      Old.scaleCoul[ij]     = d_a[SelectedComponent].scaleCoul[Location + ij];
     }
+  }
+  else if(MoveType == REINSERTION)
+  {
+    if(ij < Oldsize)
+    { 
+      Old.pos[ij]       = d_a[SelectedComponent].pos[Location + ij];
+      Old.scale[ij]     = d_a[SelectedComponent].scale[Location + ij];
+      Old.charge[ij]    = d_a[SelectedComponent].charge[Location + ij];
+      Old.scaleCoul[ij] = d_a[SelectedComponent].scaleCoul[Location + ij];
+    }
+    else
+    {
+      Old.pos[ij]       = temp[ij - Oldsize];
+      Old.scale[ij]     = d_a[SelectedComponent].scale[Location + ij - Oldsize];
+      Old.charge[ij]    = d_a[SelectedComponent].charge[Location + ij - Oldsize];
+      Old.scaleCoul[ij] = d_a[SelectedComponent].scaleCoul[Location + ij - Oldsize];
+    }
+  }
+}
+
+__global__ void Initialize_WaveVector_General(Boxsize Box, int3 kmax, Atoms* d_a, Atoms New, Atoms Old, double3* temp, size_t Oldsize, size_t Newsize, size_t SelectedComponent, size_t Location, size_t chainsize, size_t numberOfAtoms, int MoveType)
+{
+  //Zhao's note: need to think about changing this boolean to switch//
+  size_t ij = blockIdx.x * blockDim.x + threadIdx.x;
+  
+  if(ij < (Newsize + Oldsize))
+  {
+    Initialize_Copy_Positions_Together(d_a, New, Old, temp, Oldsize, Newsize, SelectedComponent, Location, chainsize, MoveType);
     //Old+New//
     Complex tempcomplex; tempcomplex.real = 1.0; tempcomplex.imag = 0.0;
     tempcomplex.real = 1.0; tempcomplex.imag = 0.0;
@@ -224,8 +249,6 @@ __global__ void Initialize_WaveVector_IdentitySwap(Boxsize Box, int3 kmax, doubl
       Old.charge[ij]    = d_a[NEWComponent].charge[ij - Oldsize];
       Old.scaleCoul[ij] = 1.0;
     }
-    
-    
     //Old+New//
     Complex tempcomplex; tempcomplex.real = 1.0; tempcomplex.imag = 0.0;
     tempcomplex.real = 1.0; tempcomplex.imag = 0.0;
@@ -467,8 +490,10 @@ double2 GPU_EwaldDifference_General(Simulations& Sim, ForceField& FF, Components
     }
     case REINSERTION: // Reinsertion //
     {
-      throw std::runtime_error("EWALD: Use the Special Function for Reinsertion");
-      //break;
+       Newsize   = SystemComponents.Moleculesize[SelectedComponent];
+       Oldsize   = SystemComponents.Moleculesize[SelectedComponent];
+       chainsize = SystemComponents.Moleculesize[SelectedComponent] - 1;
+       break;
     }
     case IDENTITY_SWAP:
     {
@@ -501,7 +526,7 @@ double2 GPU_EwaldDifference_General(Simulations& Sim, ForceField& FF, Components
   size_t numberOfAtoms = Oldsize + Newsize;
 
   size_t Nblock = 0; size_t Nthread = 0; Setup_threadblock(Oldsize + Newsize, Nblock, Nthread);
-  Initialize_WaveVector_General<<<Nblock,Nthread>>>(Box, Box.kmax, d_a, New, Old, Oldsize, Newsize, SelectedComponent, Location, chainsize, numberOfAtoms, MoveType); checkCUDAErrorEwald("error Initializing Ewald Vectors");
+  Initialize_WaveVector_General<<<Nblock,Nthread>>>(Box, Box.kmax, d_a, New, Old, SystemComponents.tempMolStorage, Oldsize, Newsize, SelectedComponent, Location, chainsize, numberOfAtoms, MoveType); checkCUDAErrorEwald("error Initializing Ewald Vectors");
 
   //Fourier Loop//
   size_t numberOfStructureFactors = (Box.kmax.x + 1) * (2 * Box.kmax.y + 1) * (2 * Box.kmax.z + 1);
@@ -552,39 +577,6 @@ double2 GPU_EwaldDifference_General(Simulations& Sim, ForceField& FF, Components
   return {SameSum, 2.0 * CrossSum};
 }
 
-//Zhao's note: THIS IS A SPECIAL FUNCTION JUST FOR REINSERTION//
-double2 GPU_EwaldDifference_Reinsertion(Boxsize& Box, Atoms*& d_a, Atoms& Old, double3* temp, ForceField& FF, double* Blocksum, Components& SystemComponents, size_t SelectedComponent, size_t UpdateLocation)
-{
-  if(FF.noCharges && !SystemComponents.hasPartialCharge[SelectedComponent]) return {0.0, 0.0};
-  double alpha = Box.Alpha; double alpha_squared = alpha * alpha;
-  double prefactor = Box.Prefactor * (2.0 * M_PI / Box.Volume);
-
-  size_t numberOfAtoms = SystemComponents.Moleculesize[SelectedComponent];
-  size_t Oldsize = 0; size_t Newsize = numberOfAtoms;
-  //Zhao's note: translation/rotation/reinsertion involves new + old states. Insertion/Deletion only has the new state.
-  Oldsize         = SystemComponents.Moleculesize[SelectedComponent];
-  numberOfAtoms  += Oldsize;
-
-  Complex* SameType = Box.AdsorbateEik; Complex* CrossType = Box.FrameworkEik;
-
-  // Construct exp(ik.r) for atoms and k-vectors kx, ky, kz = 0, 1 explicitly
-  size_t Nblock = 0; size_t Nthread = 0; Setup_threadblock(Oldsize + Newsize, Nblock, Nthread);
-  Initialize_WaveVector_Reinsertion<<<Nblock,Nthread>>>(Box, Box.kmax, temp, d_a, Old, Oldsize, Newsize, UpdateLocation, numberOfAtoms, SelectedComponent);
-
-  //Fourier Loop//
-  size_t numberOfStructureFactors = (Box.kmax.x + 1) * (2 * Box.kmax.y + 1) * (2 * Box.kmax.z + 1);
-  Nblock = 0; Nthread = 0; Setup_threadblock(numberOfStructureFactors, Nblock, Nthread);
-  Fourier_Ewald_Diff<<<Nblock * 2, Nthread, Nthread * sizeof(double)>>>(Box, SameType, CrossType, Old, alpha_squared, prefactor, Box.kmax, Oldsize, Newsize, Blocksum, false, Nblock);
-  //double sum[Nblock * 2]; 
-  double SameSum = 0.0;  double CrossSum = 0.0;
-  cudaMemcpy(SystemComponents.host_array, Blocksum, 2 * Nblock * sizeof(double), cudaMemcpyDeviceToHost);
-
-  for(size_t i = 0; i < Nblock; i++){SameSum += SystemComponents.host_array[i];}
-  for(size_t i = Nblock; i < 2 * Nblock; i++){CrossSum += SystemComponents.host_array[i];}
-
-  return {SameSum, 2.0 * CrossSum};
-}
-
 double2 GPU_EwaldDifference_IdentitySwap(Boxsize& Box, Atoms*& d_a, Atoms& Old, double3* temp, ForceField& FF, double* Blocksum, Components& SystemComponents, size_t OLDComponent, size_t NEWComponent, size_t UpdateLocation)
 {
   if(FF.noCharges && !SystemComponents.hasPartialCharge[NEWComponent] && !SystemComponents.hasPartialCharge[OLDComponent]) return {0.0, 0.0};

src_clean/VDW_Coulomb.cuh

@@ -98,8 +98,6 @@ __global__ void Energy_difference_LambdaChange(Boxsize Box, Atoms* System, Atoms
 //////////////////////
 double CPU_EwaldDifference(Boxsize& Box, Atoms& New, Atoms& Old, ForceField& FF, Components& SystemComponents, size_t SelectedComponent, bool Swap, size_t SelectedTrial);
 
-double2 GPU_EwaldDifference_Reinsertion(Boxsize& Box, Atoms*& d_a, Atoms& Old, double3* temp, ForceField& FF, double* Blocksum, Components& SystemComponents, size_t SelectedComponent, size_t UpdateLocation);
-
 double2 GPU_EwaldDifference_IdentitySwap(Boxsize& Box, Atoms*& d_a, Atoms& Old, double3* temp, ForceField& FF, double* Blocksum, Components& SystemComponents, size_t OLDComponent, size_t NEWComponent, size_t UpdateLocation);
 
 void Copy_Ewald_Vector(Simulations& Sim);
@@ -149,7 +147,7 @@ double TailCorrectionIdentitySwap(Components& SystemComponents, size_t NEWCompon
 // Deep Potential For Host-Guest //
 //             General Functions //
 ///////////////////////////////////
-void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, Components& SystemComponents, size_t SelectedComponent, int MoveType, size_t Location);
+void Prepare_DNN_InitialPositions(Atoms*& d_a, Atoms& New, Atoms& Old, double3* temp, Components& SystemComponents, size_t SelectedComponent, int MoveType, size_t Location);
 
 void Prepare_DNN_InitialPositions_Reinsertion(Atoms*& d_a, Atoms& Old, double3* temp, Components& SystemComponents, size_t SelectedComponent, size_t Location);
 

src_clean/axpy.cu

@@ -99,7 +99,7 @@ void Select_Box_Component_Molecule(Variables& Vars, size_t box_index)
 
   Vars.RandomNumber = Get_Uniform_Random();
 }
-inline void RunMoves(Variables& Vars, size_t box_index, int Cycle)
+void RunMoves(Variables& Vars, size_t box_index, int Cycle)
 {
   MC_MOVES MOVES;
 

src_clean/axpy.h

@@ -16,6 +16,8 @@ void Calculate_Exclusion_Energy_Rigid(Boxsize& Box, Atoms* System, ForceField FF
 void Check_WaveVector_CPUGPU(Boxsize& Box, Components& SystemComponents);
 */
 
+void RunMoves(Variables& Vars, size_t box_index, int Cycle);
+
 double CreateMolecule_InOneBox(Variables& Vars, size_t systemId, bool AlreadyHasFractionalMolecule);
 
 void Run_Simulation_MultipleBoxes(Variables& Vars);
@@ -49,7 +51,10 @@ __global__ void StoreNewLocation_Reinsertion(Atoms Mol, Atoms NewMol, double3* t
 __global__ void Update_Reinsertion_data(Atoms* d_a, double3* temp, size_t SelectedComponent, size_t UpdateLocation);
 
 double GetPrefactor(Components& SystemComponents, Simulations& Sims, size_t SelectedComponent, int MoveType);
-void AcceptInsertion(Variables& Vars, size_t systemId, int MoveType, CBMC_Variables& InsertionVariables);
+
+//void AcceptInsertion(Variables& Vars, size_t systemId, int MoveType, CBMC_Variables& InsertionVariables);
+void AcceptInsertion(Variables& Vars, CBMC_Variables& InsertionVariables, size_t systemId, int MoveType);
+
 void AcceptDeletion(Variables& Vars, size_t systemId, int MoveType);
 MoveEnergy Insertion_Body(Variables& Vars, size_t systemId, CBMC_Variables& CBMC);
 MoveEnergy Deletion_Body(Variables& Vars, size_t systemId, CBMC_Variables& CBMC);
@@ -72,14 +77,4 @@ struct SingleMove
   }
 };
 
-//struct InsertionMove;
-//struct InsertionMove {};
-//struct MC_MOVES {};
-/*
-struct MC_MOVES
-{
-  InsertionMove INSERTION;
-};
-*/
-
 #include "move_struct.h"

src_clean/data_struct.h

@@ -1229,6 +1229,7 @@ struct Simulations //For multiple simulations//
   Atoms*  d_a;                  // Pointer For Atom Data in the Simulation Box //
   Atoms   Old;                  // Temporary data storage for Old Configuration //
   Atoms   New;                  // Temporary data storage for New Configuration //
+  Atoms   Temp;                 // Temporary storage (for xfering data and accepting move) //
   int2*   ExcludeList;          // Atoms to exclude during energy calculations: x: component, y: molecule-ID (may need to add z and make it int3, z: atom-ID)
   double* Blocksum;             // Block sums for partial reduction //
   bool*   device_flag;          // flags for overlaps on the device //

src_clean/main.cpp

@@ -272,8 +272,9 @@ Variables Initialize(void) //for pybind
     Prepare_TempSystem_On_Host(Vars.SystemComponents[a].TempSystem);
     cudaMemcpy(Vars.Sims[a].d_a, device_System, sizeof(Atoms)*NComponents.x, cudaMemcpyHostToDevice);
     // SET UP TEMPORARY ARRAYS //
-    Setup_Temporary_Atoms_Structure(Vars.Sims[a].Old, Vars.SystemComponents[a].HostSystem);
-    Setup_Temporary_Atoms_Structure(Vars.Sims[a].New, Vars.SystemComponents[a].HostSystem);
+    Setup_Temporary_Atoms_Structure(Vars.Sims[a].Old,  Vars.SystemComponents[a].HostSystem);
+    Setup_Temporary_Atoms_Structure(Vars.Sims[a].New,  Vars.SystemComponents[a].HostSystem);
+    Setup_Temporary_Atoms_Structure(Vars.Sims[a].Temp, Vars.SystemComponents[a].HostSystem);
 
     if(Vars.SystemComponents[a].UseDNNforHostGuest)
     {

src_clean/main.h

@@ -0,0 +1,6 @@
+Variables Initialize(void);
+void RunSimulation(Variables& Vars);
+void EndOfSimulationWrapUp(Variables& Vars);
+
+MoveEnergy check_energy_wrapper(Variables& Var, size_t i);
+void ENERGY_SUMMARY(Variables& Vars);