Bug-fix: computation of average excess loading at small pressures (di…

…fference is at most 1 molecule per uc)
iRASPA · May 19, 2021 · beb278b · beb278b
1 parent 405beff
commit beb278b
Show file tree

Hide file tree

Showing 5 changed files with 14 additions and 8 deletions.
diff --git a/README.md b/README.md
@@ -49,6 +49,10 @@ control volume grand canonical molecular dynamics study",
 * Helium potential for computing void-fractions.\
   J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1954, p. 1114.
   [reference in article](https://www.sciencedirect.com/science/article/abs/pii/S0927775701006288)
+* Linear/Branched alkane model.\
+  "United atom force field for alkanes in nanoporous materials",
+  D. Dubbeldam, S. Calero, T.J.H. Vlugt, R. Krishna, T.L.M. Maesen, B. Smit, J. Phys. Chem. B 2004, 108(33), 12301-12313
+  [link to article](https://pubs.acs.org/doi/abs/10.1021/jp0376727)
 
 rigid MOFs
 ------------

diff --git a/src/framework.c b/src/framework.c
@@ -4138,7 +4138,7 @@ void ExpandAsymmetricFrameworkToFullFramework(void)
 
               dr=ConvertFromABCtoXYZ(ds); 
               rr=SQR(dr.x)+SQR(dr.y)+SQR(dr.z);
-              if(rr<1e-3) index=k;
+              if(rr<1e-2) index=k;
             }
 
             // the symmetry position is new -> insert it as a framework atom

diff --git a/src/output.c b/src/output.c
@@ -260,7 +260,7 @@ void PrintPreSimulationStatusCurrentSystem(int system)
   fprintf(FilePtr,"Compiler and run-time data\n");
   fprintf(FilePtr,"===========================================================================\n");
 
-  fprintf(FilePtr,"%s\n","RASPA 2.0.44");
+  fprintf(FilePtr,"%s\n","RASPA 2.0.45");
 
   #if defined (__LP64__) || defined (__64BIT__) || defined (_LP64) || (__WORDSIZE == 64)
     fprintf(FilePtr,"Compiled as a 64-bits application\n");

diff --git a/src/potentials.c b/src/potentials.c
@@ -11957,11 +11957,13 @@ REAL TailMolecularEnergyDifferenceAddRemove(int ComponentToAdd,int ComponentToRe
     {
       if(TailCorrection[i][j])
       {
-        energy_new+=2.0*M_PI*NumberOfPseudoAtomsTypeNew[i]*NumberOfPseudoAtomsTypeNew[j]*
-                    PotentialCorrection(i,j,CutOffVDW);
+        energy_new+=2.0*M_PI*(NumberOfPseudoAtomsTypeNew[i]-NumberOfFractionalPseudoAtomsType[CurrentSystem][i])*
+                             (NumberOfPseudoAtomsTypeNew[j]-NumberOfFractionalPseudoAtomsType[CurrentSystem][j])*
+                             PotentialCorrection(i,j,CutOffVDW);
 
-        energy_old+=2.0*M_PI*NumberOfPseudoAtomsType[CurrentSystem][i]*NumberOfPseudoAtomsType[CurrentSystem][j]*
-                    PotentialCorrection(i,j,CutOffVDW);
+        energy_old+=2.0*M_PI*(NumberOfPseudoAtomsType[CurrentSystem][i]-NumberOfFractionalPseudoAtomsType[CurrentSystem][i])*
+                             (NumberOfPseudoAtomsType[CurrentSystem][j]-NumberOfFractionalPseudoAtomsType[CurrentSystem][j])*
+                             PotentialCorrection(i,j,CutOffVDW);
       }
     }
   return (energy_new-energy_old)/Volume[CurrentSystem];

diff --git a/src/statistics.c b/src/statistics.c
@@ -1120,7 +1120,7 @@ void UpdateEnergyAveragesCurrentSystem(void)
   REAL UCFMCAdsorbate,UCFMCCation;
   REAL weight,lambda;
   int numberOfAbsoluteIntegerMoleculesForComponent;
-  int numberOfExcessIntegerMoleculesForComponent;
+  double numberOfExcessIntegerMoleculesForComponent;
 
   // check for new block
   if(CurrentCycle==BlockCycle[Block])
@@ -1754,7 +1754,7 @@ void PrintIntervalStatusInit(long long CurrentCycle,long long NumberOfCycles,FIL
          (double)(Components[i].MOLEC_PER_UC_TO_CC_STP_G[CurrentSystem]*loading),
          (double)(Components[i].MOLEC_PER_UC_TO_CC_STP_CC[CurrentSystem]*loading));
 
-      loading=((REAL)Components[i].NumberOfMolecules[CurrentSystem]
+      loading=(REAL)(Components[i].NumberOfMolecules[CurrentSystem]
                     -Components[i].AmountOfExcessMolecules[CurrentSystem]
                     -(Components[i].FractionalMolecule[CurrentSystem]>=0?1:0)
                     -Components[i].NumberOfRXMCMoleculesPresent[CurrentSystem])/(REAL)number_of_unit_cells;