Merge pull request #304 from GregorySchwing/rebasedpsfparsertest

PSF Parser can handle MosDef format Porous Materials
GOMC-WSU · Jan 7, 2021 · 3bd6a3f · 3bd6a3f
2 parents e6984f3 + d381921
commit 3bd6a3f
Show file tree

Hide file tree

Showing 23 changed files with 3,450 additions and 91 deletions.
diff --git a/src/BondAdjacencyList.cpp b/src/BondAdjacencyList.cpp
@@ -22,7 +22,7 @@ adjNode* BondAdjacencyList::getAdjListNode(int value, int weight, adjNode* head)
     return newNode;
 }
 
-BondAdjacencyList::BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::vector< std::vector<uint> > & moleculeXAtomIDY){
+BondAdjacencyList::BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::vector< std::vector<int> > & moleculeXAtomIDY){
 
     this->nAtoms = nAtoms;
     this->nBonds = nBonds;
@@ -75,18 +75,19 @@ BondAdjacencyList::BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::v
 
     connectedComponents(moleculeXAtomIDY);
     /* Sort Atom Indices in N connected components, then sort N connected components by first atom index*/
-    for (std::vector< std::vector<uint> >::iterator it = moleculeXAtomIDY.begin();
+    for (std::vector< std::vector<int> >::iterator it = moleculeXAtomIDY.begin();
         it != moleculeXAtomIDY.end(); it++){
         std::sort(it->begin(), it->end());
     }
     std::sort(moleculeXAtomIDY.begin(), moleculeXAtomIDY.end());
+/*
 #ifndef NDEBUG
     std::cout << "Adjacency List" << std::endl;
     for (uint i = 0; i < nAtoms; i++)
         display_AdjList(head[i], i);
 
     std::cout << "Connected Components" << std::endl;
-    for (std::vector< std::vector<uint> >::iterator it = moleculeXAtomIDY.begin();
+    for (std::vector< std::vector<int> >::iterator it = moleculeXAtomIDY.begin();
         it != moleculeXAtomIDY.end(); it++){
         for (std::vector<uint>::iterator it2 = it->begin();
             it2 != it->end(); it2++){
@@ -95,7 +96,7 @@ BondAdjacencyList::BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::v
         std::cout << std::endl;
     }
 #endif    
-  
+*/
 }
 
 // Destructor
@@ -124,7 +125,7 @@ void BondAdjacencyList::display_AdjList(adjNode* ptr, int i)
 
 // Method to print connected components in an 
 // undirected graph 
-void BondAdjacencyList::connectedComponents(std::vector< std::vector<uint> > & moleculeXAtomIDY) 
+void BondAdjacencyList::connectedComponents(std::vector< std::vector<int> > & moleculeXAtomIDY) 
 { 
     // Mark all the vertices as not visited 
     bool *visited = new bool[nAtoms]; 
@@ -139,7 +140,7 @@ void BondAdjacencyList::connectedComponents(std::vector< std::vector<uint> > & m
             // from v 
             /* For debugging
             std::cout << "Calling DFSUtil" << std::endl; */
-            std::vector<uint> moleculeX;
+            std::vector<int> moleculeX;
             DFSUtil(v, this->head[v], this->head, visited, moleculeX); 
             moleculeXAtomIDY.push_back(moleculeX);
             /* For debugging std::cout << "\n"; */
@@ -148,7 +149,7 @@ void BondAdjacencyList::connectedComponents(std::vector< std::vector<uint> > & m
     delete[] visited; 
 } 
 
-void BondAdjacencyList::DFSUtil(int v, adjNode * node, adjNode ** head, bool * visited, std::vector<uint> & moleculeX) 
+void BondAdjacencyList::DFSUtil(int v, adjNode * node, adjNode ** head, bool * visited, std::vector<int> & moleculeX) 
 { 
     // Mark the current node as visited and print it 
     visited[v] = true; 

diff --git a/src/BondAdjacencyList.h b/src/BondAdjacencyList.h
@@ -37,12 +37,12 @@ uint nAtoms;  // number of nodes in the graph
 uint nBonds; // number of edges in the graph
 
 
-BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::vector< std::vector<uint> > & moleculeXAtomIDY);
+BondAdjacencyList(FILE* psf, uint nAtoms, uint nBonds, std::vector< std::vector<int> > & moleculeXAtomIDY);
 ~BondAdjacencyList();
 adjNode* getAdjListNode(int value, int weight, adjNode* head);
 void display_AdjList(adjNode* ptr, int i);
-void connectedComponents(std::vector< std::vector<uint> > & moleculeXAtomIDY);
-void DFSUtil(int v, adjNode * node, adjNode ** head, bool * visited, std::vector<uint> & moleculeX);
+void connectedComponents(std::vector< std::vector<int> > & moleculeXAtomIDY);
+void DFSUtil(int v, adjNode * node, adjNode ** head, bool * visited, std::vector<int> & moleculeX);
 
 graphEdge * edges;
 

diff --git a/src/MolSetup.cpp b/src/MolSetup.cpp
diff --git a/src/MolSetup.h b/src/MolSetup.h
@@ -27,6 +27,14 @@ class FFSetup;
 
 namespace mol_setup
 {
+
+struct MoleculeVariables {
+  std::vector<uint> startIdxMolecules, moleculeKinds;
+  std::vector<std::string> moleculeNames, moleculeKindNames;
+  uint lastAtomIndexInBox0 = 0;
+  uint lastResKindIndex = 0;
+};
+
 //!structure to contain an atom's data during initialization
 class Atom
 {
@@ -146,7 +154,7 @@ typedef std::map<std::size_t, std::vector<std::string> > SizeMap;
 *\param numFiles number of files to read
 *\return -1 if failed, 0 if successful
 */
-int ReadCombinePSF(MolMap& kindMap, SizeMap& sizeMap, const std::string* psfFilename,
+int ReadCombinePSF(MoleculeVariables & molVars, MolMap& kindMap, SizeMap& sizeMap, const std::string* psfFilename,
                    const bool* psfDefined, pdb_setup::Atoms& pdbAtoms);
 
 void PrintMolMapVerbose(const MolMap& kindMap);
@@ -159,14 +167,15 @@ class MolSetup
 public:
   class Atom;
   int read_atoms(FILE *, unsigned int nAtoms, std::vector<mol_setup::Atom> & allAtoms);
-  void createMapAndModifyPDBAtomDataStructure(const BondAdjacencyList & bondAdjList,
-                                              const std::vector< std::vector<uint> > & moleculeXAtomIDY, 
-                                              std::vector<mol_setup::Atom> & allAtoms,
-                                              mol_setup::MolMap & kindMap,
-                                              mol_setup::SizeMap & sizeMap,
-                                              pdb_setup::Atoms& pdbAtoms,
-                                              mol_setup::MolMap * kindMapFromBox1,
-                                              mol_setup::SizeMap * sizeMapFromBox1);
+  void createKindMap (pdb_setup::Atoms& pdbAtoms,
+                      mol_setup::MoleculeVariables & molVars,
+                      const BondAdjacencyList & bondAdjList,
+                      std::vector< std::vector<int> > & moleculeXAtomIDY, 
+                      std::vector<mol_setup::Atom> & allAtoms,
+                      mol_setup::MolMap & kindMap,
+                      mol_setup::SizeMap & sizeMap,
+                      mol_setup::MolMap * kindMapFromBox1,
+                      mol_setup::SizeMap * sizeMapFromBox1);
 
   static void copyBondInfoIntoMapEntry(const BondAdjacencyList & bondAdjList, mol_setup::MolMap & kindMap, std::string fragName);
 
@@ -178,10 +187,11 @@ class MolSetup
            const bool* psfDefined, 
            pdb_setup::Atoms& pdbAtoms);
 
-  void AssignKinds(const pdb_setup::Atoms& pdbAtoms, const FFSetup& ffData);
+  void AssignKinds(const mol_setup::MoleculeVariables& molVars, const FFSetup& ffData);
 
 //private:
   mol_setup::MolMap kindMap;
   mol_setup::SizeMap sizeMap;
+  mol_setup::MoleculeVariables molVars;
 };
 #endif
diff --git a/src/Molecules.cpp b/src/Molecules.cpp
@@ -43,23 +43,25 @@ void Molecules::Init(Setup & setup, Forcefield & forcefield,
   kinds = new MoleculeKind[kindsCount];
 
   //Molecule instance arrays/data
-  count = atoms.startIdxRes.size();
+  count = setup.mol.molVars.startIdxMolecules.size();
   if (count == 0) {
     std::cerr << "Error: No Molecule was found in the PDB file(s)!" << std::endl;
     exit(EXIT_FAILURE);
   }
 
   start = new uint [count + 1];
-  start = vect::TransferInto<uint>(start, atoms.startIdxRes);
+  start = vect::TransferInto<uint>(start, setup.mol.molVars.startIdxMolecules);
   start[count] = atoms.x.size();
-  kIndex = vect::transfer<uint>(atoms.resKinds);
-  kIndexCount = atoms.resKinds.size();
+
+  kIndex = vect::transfer<uint>(setup.mol.molVars.moleculeKinds);
+  kIndexCount = setup.mol.molVars.moleculeKinds.size();
+
   chain = vect::transfer<char>(atoms.chainLetter);
   for (uint mk = 0 ; mk < kindsCount; mk++) {
     countByKind[mk] =
-      std::count(atoms.resNames.begin(), atoms.resNames.end(),
-                 atoms.resKindNames[mk]);
-    kinds[mk].Init(atoms.resKindNames[mk], setup, forcefield, sys);
+      std::count(setup.mol.molVars.moleculeNames.begin(), setup.mol.molVars.moleculeNames.end(),
+                 setup.mol.molVars.moleculeKindNames[mk]);
+    kinds[mk].Init(setup.mol.molVars.moleculeKindNames[mk], setup, forcefield, sys);
   }
 
 #if ENSEMBLE == GCMC

diff --git a/src/PDBSetup.h b/src/PDBSetup.h
@@ -132,8 +132,6 @@ class Atoms : public FWReadableBase
   uint numAtomsInBox[BOX_TOTAL]; // number of atom in each box
   std::string currResname;
 
-  uint lastAtomIndexInBox0 = 0;
-  uint lastResKindIndex = 0;
 };
 
 }

diff --git a/src/PSFOutput.cpp b/src/PSFOutput.cpp
@@ -37,12 +37,12 @@ const int dihPerLine = 2;
 
 PSFOutput::PSFOutput(const Molecules& molecules, const System &sys,
                      Setup & set) :
-  molecules(&molecules), molNames(set.pdb.atoms.resKindNames),
+  molecules(&molecules), molNames(set.mol.molVars.moleculeKindNames),
   molLookRef(sys.molLookup)
 {
   molKinds.resize(set.mol.kindMap.size());
-  for(uint i = 0; i < set.pdb.atoms.resKindNames.size(); ++i) {
-    molKinds[i] = set.mol.kindMap[set.pdb.atoms.resKindNames[i]];
+  for(uint i = 0; i < set.mol.molVars.moleculeKindNames.size(); ++i) {
+    molKinds[i] = set.mol.kindMap[set.mol.molVars.moleculeKindNames[i]];
   }
   CountMolecules();
   PrintPSF(set.config.out.state.files.psf.name);

diff --git a/src/Setup.h b/src/Setup.h
@@ -55,7 +55,7 @@ class Setup
     if(mol.Init(config.in.restart, config.in.files.psf.name, config.in.files.psf.defined, pdb.atoms) != 0) {
       exit(EXIT_FAILURE);
     }
-    mol.AssignKinds(pdb.atoms, ff);
+    mol.AssignKinds(mol.molVars, ff);
 
   }
 };

diff --git a/test/FileList.cmake b/test/FileList.cmake
@@ -1,5 +1,6 @@
 set(TestSources
     test/src/BasicTypesTest.cpp
     test/src/BitLibTest.cpp
+    test/src/PSFParserTest.cpp
     test/src/EndianTest.cpp
 )
diff --git a/test/GoogleTest.cmake b/test/GoogleTest.cmake
@@ -36,4 +36,5 @@ include(test/FileList.cmake)
 # Now simply link against gtest or gtest_main as needed. Eg
 add_executable(GOMC_Test ${TestSources})
 target_link_libraries(GOMC_Test gtest_main)
-add_test(NAME BasicTypesTest COMMAND BasicTypesTest)
+add_test(NAME BasicTypesTest COMMAND BasicTypesTest)
+add_test(NAME PSFParserTest COMMAND PSFParserTest)
diff --git a/test/input/GOMC_GEMC.conf b/test/input/GOMC_GEMC.conf
@@ -0,0 +1,157 @@
+########################
+#Control file for SPCE water at 350K
+########################
+
+############################################################################
+#  ========-------------------- INPUT --------------------------===========
+############################################################################
+
+#########################
+# enable, step
+#########################
+Restart	 	false  
+
+
+####################################
+# kind {RESTART, RANDOM, INTSEED}
+####################################
+PRNG		RANDOM
+
+####################################
+# FORCE FIELD
+####################################
+ParaTypeCHARMM	 true
+
+Parameters     	 ./water.par
+Parameters     	 ./GOMC_water_FF.inp
+####################################
+# INPUT PDB FILES
+####################################
+Coordinates 0    ./PSFParser/BOX_0.pdb
+Coordinates 1	 ./PSFParser/BOX_1.pdb
+
+####################################
+# INPUT PSF FILES
+####################################
+Structure 0  	 ./PSFParser/BOX_0.psf
+Structure 1	 ./PSFParser/BOX_1.psf
+
+############################################################################
+#  =======--------------------- SYSTEM --------------------------===========
+############################################################################
+##################################
+# GEMC TYPE (DEFULT IS NVT_GEMC)  
+##################################
+GEMC 	NVT
+
+#############################
+# SIMULATION CONDITION   
+#############################
+Temperature     350
+Potential       SWITCH 
+Rswitch		10.0
+LRC		true
+Rcut		12
+RcutLow		0.1   # was set to 1.2
+Exclude 	1-4 
+
+#############################
+# ELECTROSTATIC   
+#############################
+Ewald		true
+ElectroStatic   true
+CachedFourier	false
+Tolerance      	0.00001
+1-4scaling     false
+
+
+###############################                                                
+# PRESSURE CALCULATION                                                         
+################################                                               
+PressureCalc  true 5000
+
+################################
+# STEPS 
+################################
+RunSteps           	10
+EqSteps 		5
+AdjSteps	   	2 
+
+
+
+################################
+# MOVE FREQUENCY              
+################################
+
+VolFreq		     0.01
+MultiParticleFreq     0.01
+DisFreq              0.18
+RotFreq              0.10
+RegrowthFreq         0.10
+#IntraMEMC-2Freq      0.10
+IntraSwapFreq        0.10
+#MEMC-2Freq           0.25
+SwapFreq             0.50	
+
+################################
+# BOX DIMENSION #, X, Y, Z
+################################
+CellBasisVector1  0   40.0  	0.00  		0.00
+CellBasisVector2  0   0.00  		40.0  	0.00
+CellBasisVector3  0   0.00  		0.00  		40.0
+
+CellBasisVector1  1   200.00  	0.00  		0.00
+CellBasisVector2  1   0.00  		200.00  	0.00
+CellBasisVector3  1   0.00  		0.00  		200.00
+
+
+
+##############################
+# CBMC TRIALS
+##############################
+CBMC_First   16    
+CBMC_Nth     8
+CBMC_Ang     50
+CBMC_Dih     50
+
+
+############################################################################
+#  =======-------------------- OUTPUT --------------------------===========
+############################################################################
+
+##########################
+# statistics filename add
+##########################
+OutputName Output_data
+
+#####################################
+# enable, frequency           
+#####################################
+RestartFreq 	   	true	10
+CheckpointFreq		true	10
+CoordinatesFreq   	true  	10
+ConsoleFreq		true	10
+BlockAverageFreq	true	10
+HistogramFreq		true	10
+
+DCDRestartFreq		true 	10
+
+################################
+# OutHistSettings
+################################
+DistName	 dis
+HistName	 his
+RunNumber	 1
+RunLetter	 a
+SampleFreq	 100
+
+##################################
+# enable: blk avg., fluct.
+##################################
+OutEnergy         true    true   
+OutPressure       true    true  
+OutMolNum         true    true   
+OutDensity        true    true
+OutVolume         true    true
+OutSurfaceTention false false
+