Changeset 4bb871a

Timestamp:

Feb 2, 2010, 12:22:06 PM (16 years ago)

Author:

Tillmann Crueger <crueger@…>

Children:

45cd358

Parents:

adcdf8 (diff), 6d0fcaa (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge commit 'heber/Analysis_PairCorrelation' into MenuRefactoring

Files:

• molecuilder/src/analysis_correlation.cpp (modified) (4 diffs)
• molecuilder/src/analysis_correlation.hpp (modified) (2 diffs)
• molecuilder/src/atom_bondedparticle.cpp (modified) (2 diffs)
• molecuilder/src/atom_bondedparticle.hpp (modified) (1 diff)
• molecuilder/src/boundary.cpp (modified) (2 diffs)
• molecuilder/src/builder.cpp (modified) (3 diffs)
• molecuilder/src/molecule.hpp (modified) (1 diff)
• molecuilder/src/molecule_graph.cpp (modified) (1 diff)
• util/src/NanoCreator.c (modified) (11 diffs)
• util/src/NanoCreator.h (modified) (1 diff)

molecuilder/src/analysis_correlation.cpp

@@ -370 +370 @@
 {
   Info FunctionInfo(__func__);
-  *file << "# BinStart\tCount" << endl;
+  *file << "BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << runner->second << endl;
@@ -383 +383 @@
 {
   Info FunctionInfo(__func__);
-  *file << "# BinStart\tAtom1\tAtom2" << endl;
+  *file << "BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
@@ -396 +396 @@
 {
   Info FunctionInfo(__func__);
-  *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
+  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
   for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first;
@@ -412 +412 @@
 {
   Info FunctionInfo(__func__);
-  *file << "# BinStart\tTriangle" << endl;
+  *file << "BinStart\tTriangle" << endl;
   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;

molecuilder/src/analysis_correlation.hpp

@@ -90 +90 @@
 /** Puts given correlation data into bins of given size (histogramming).
  * Note that BinStart and BinEnd are desired to fill the complete range, even where Bins are zero. If this is
- * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater.
+ * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater. If a
+ * certain start value is desired, give BinStart and a BinEnd that is smaller than BinStart, then only BinEnd will be
+ * calculated automatically, i.e. BinStart = 0. and BinEnd = -1. .
  * Also note that the range is given inclusive, i.e. [ BinStart, BinEnd ].
  * \param *map map of doubles to count
@@ -114 +116 @@
   if (BinStart == BinEnd) { // if same, find range ourselves
     GetMinMax( map, start, end);
-  } else { // if not, initialise range to zero
+  } else if (BinEnd < BinStart) { // if BinEnd smaller, just look for new max
+    GetMinMax( map, start, end);
+    start = BinStart;
+  } else { // else take both values
     start = BinStart;
     end = BinEnd;

molecuilder/src/atom_bondedparticle.cpp

@@ -56 +56 @@
  * \param *AdjacencyFile output stream
  */
-void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
+void BondedParticle::OutputAdjacency(ofstream * const AdjacencyFile) const
 {
   *AdjacencyFile << nr << "\t";
@@ -62 +62 @@
     *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
   *AdjacencyFile << endl;
+};
+
+/** Output of atom::nr along each bond partner per line.
+ * Only bonds are printed where atom::nr is smaller than the one of the bond partner.
+ * \param *AdjacencyFile output stream
+ */
+void BondedParticle::OutputBonds(ofstream * const BondFile) const
+{
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    if (nr < (*Runner)->GetOtherAtom(this)->nr)
+      *BondFile << nr << "\t" << (*Runner)->GetOtherAtom(this)->nr << "\n";
 };
 

molecuilder/src/atom_bondedparticle.hpp

@@ -44 +44 @@
   int CorrectBondDegree();
   void OutputBondOfAtom() const;
-  void OutputAdjacency(ofstream *AdjacencyFile) const;
+  void OutputAdjacency(ofstream * const AdjacencyFile) const;
+  void OutputBonds(ofstream * const BondFile) const;
   void OutputOrder(ofstream *file) const;
 

molecuilder/src/boundary.cpp

@@ -19 +19 @@
 
 #include<gsl/gsl_poly.h>
+#include<time.h>
 
 // ========================================== F U N C T I O N S =================================
@@ -840 +841 @@
   Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
 
+  // initialize seed of random number generator to current time
+  srand ( time(NULL) );
+
   // go over [0,1]^3 filler grid
   for (n[0] = 0; n[0] < N[0]; n[0]++)

molecuilder/src/builder.cpp

@@ -1308 +1308 @@
             Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
             Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
+            Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl;
+            Log() << Verbose(0) << "\t-J\t<path> Store adjacency per atom to file." << endl;
             Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
             Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
@@ -1594 +1596 @@
                 int ranges[NDIM] = {1,1,1};
                 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
-                //OutputCorrelationToSurface(&output, surfacemap);
-                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
+                OutputCorrelationToSurface(&output, surfacemap);
+                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 20. );
                 OutputCorrelation ( &binoutput, binmap );
                 output.close();
@@ -1679 +1681 @@
               }
               break;
+
+            case 'J':
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag =255;
+                eLog() << Verbose(0) << "Missing path of adjacency file: -j <path>" << endl;
+                performCriticalExit();
+              } else {
+                Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl;
+                configuration.BG->ConstructBondGraph(mol);
+                mol->StoreAdjacencyToFile(argv[argptr]);
+                argptr+=1;
+              }
+              break;
+
+            case 'j':
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag =255;
+                eLog() << Verbose(0) << "Missing path of bonds file: -j <path>" << endl;
+                performCriticalExit();
+              } else {
+                Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl;
+                configuration.BG->ConstructBondGraph(mol);
+                mol->StoreBondsToFile(argv[argptr]);
+                argptr+=1;
+              }
+              break;
+
             case 'N':
               if (ExitFlag == 0) ExitFlag = 1;

molecuilder/src/molecule.hpp

@@ -277 +277 @@
   int FragmentMolecule(int Order, config *configuration);
   bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
+  bool StoreBondsToFile(char *path);
   bool StoreAdjacencyToFile(char *path);
   bool CheckAdjacencyFileAgainstMolecule(char *path, atom **ListOfAtoms);

molecuilder/src/molecule_graph.cpp

@@ -992 +992 @@
   Log() << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
+    AdjacencyFile << "m\tn" << endl;
     ActOnAllAtoms(&atom::OutputAdjacency, &AdjacencyFile);
     AdjacencyFile.close();
+    Log() << Verbose(1) << "done." << endl;
+  } else {
+    Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl;
+    status = false;
+  }
+
+  return status;
+}
+;
+
+/** Storing the bond structure of a molecule to file.
+ * Simply stores Atom::nr and then the Atom::nr of all bond partners, one per line.
+ * \param *out output stream for debugging
+ * \param *path path to file
+ * \return true - file written successfully, false - writing failed
+ */
+bool molecule::StoreBondsToFile(char *path)
+{
+  ofstream BondFile;
+  stringstream line;
+  bool status = true;
+
+  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  BondFile.open(line.str().c_str(), ios::out);
+  Log() << Verbose(1) << "Saving adjacency list ... ";
+  if (BondFile != NULL) {
+    BondFile << "m\tn" << endl;
+    ActOnAllAtoms(&atom::OutputBonds, &BondFile);
+    BondFile.close();
     Log() << Verbose(1) << "done." << endl;
   } else {

util/src/NanoCreator.c

@@ -520 +520 @@
 {
   int c;
+  if ((a == 0) && (b == 0))
+    return 1;
+  else if (a == 0)
+    return b;
+  else if (b == 0)
+    return a;
   do {
     c = a % b;        /* Rest of integer divison */
@@ -672 +678 @@
   double *x = NULL, *coord = NULL, *tempvector = NULL, *offset = NULL;
   //double *cog = NULL, *cog_projected = NULL;
-  char stage[6];
+  char *stage = NULL;
   int axis[NDIM] = {0,1,2}, chiral[2] = {1,1}, factors[NDIM] = {1,1,1};
   char name[255], line[255], input = 'n';
   char *CellBuffer = NULL, *SheetBuffer = NULL, *TubeBuffer = NULL, *bufptr = NULL;
   double *randomness = NULL, percentage; // array with percentages for presence in sheet and beyond
-  int i,j, ggT;
+  int i,j, ggT, ggTsmall;
   int length;
 
@@ -688 +694 @@
         } else {
                 // store variables
-                filename = argv[1];
-                strncpy(stage, argv[2], 5);
+                filename = Malloc(sizeof(char)*(strlen(argv[1])+1), "Main: filename");
+    stage = Malloc(sizeof(char)*(strlen(argv[2])+1), "Main: stage");
+                strcpy(filename, argv[1]);
+                strcpy(stage, argv[2]);
     fprintf(stdout, "I will begin at stage %s.\n", stage);
 
                 // build filenames
+    length = strlen(filename);
     char *ptr = strrchr(filename, '.');
     if (ptr == NULL) {
       ptr = filename;
     } else {
-      length = strlen(filename);
       if (ptr != NULL) {
         length = ((int)(ptr-filename) >= length-5) ? (int)(ptr-filename) : length;
@@ -841 +849 @@
 
   int numbersheet, biggestdiameter, sheetnr[NDIM], tmp, seed;
+  double dfactors[2];
   double x1,x2,x3, angle;
   char flag = 'n';
@@ -874 +883 @@
 
     do {
-      fprintf(stdout, "\nNow specify the two natural numbers (m n) defining the chiral angle, \nif the result is crap, try flipping to (m,n): ");
+      fprintf(stdout, "\nNow specify the two natural numbers (n m) defining the chiral angle, \nif the result is crap, try flipping to (n,m): ");
       fscanf(stdin, "%d %d", &chiral[0], &chiral[1]);
-      ggT = GCD(2*chiral[1]+chiral[0],2*chiral[0]+chiral[1]);
+      ggTsmall = GCD(chiral[0],chiral[1]);
+      ggT = GCD(2*chiral[0]+chiral[1],2*chiral[1]+chiral[0]);
+      fprintf(stdout, "Greatest Common Denominator of (n, m) is %d\n", ggTsmall);
       fprintf(stdout, "Greatest Common Denominator of (2n+m, 2m+n) is %d\n", ggT);
       fprintf(stdout, "chiral0: %d\tchiral1: %d\n", chiral[0], chiral[1]);
       for (i=0;i<NDIM;i++) {
-        Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
-        //Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
+//        Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
+        Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
         //Tubevector[axis[0]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
         //Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
-        Tubevector[axis[1]][i] = (double)chiral[0] * OrthoVector[axis[0]][i] - (double)chiral[1] * OrthoVector[axis[1]][i];
+//        Tubevector[axis[1]][i] = (double)chiral[0] * OrthoVector[axis[0]][i] - (double)chiral[1] * OrthoVector[axis[1]][i];
+        Tubevector[axis[1]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
         //Tubevector[axis[1]][i] = (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[0]][i] + (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[1]][i];
 //        fprintf(stderr, "Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i]\n = %lg * %lg + %lg * %lg = %lg + %lg = %lg\n\n",
@@ -924 +936 @@
         x1 = gsl_vector_get(u,0)*(double)i;
         x2 = gsl_vector_get(u,1)*(double)i;
-        x3 =
         fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
         if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-6) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-6 )) {
@@ -931 +942 @@
         }
       }
-      fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));
-
-      // get length
+      x1 = fabs(gsl_vector_get(u,0));
+      x2 = fabs(gsl_vector_get(u,1));
+      fprintf(stdout, "(c,d) = (%d,%d).\n",((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)));
+      j = GCD(((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)));
+      fprintf(stdout, "GCD(%d,%d) = %i.\n", ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), j);
+      for (i=0;i<2;i++)
+        gsl_vector_set(u,i, gsl_vector_get(u,i)/(double)j);
+        // get length
       double x[NDIM];
       for (i=0;i<NDIM;i++)
         x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
       angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
-      fprintf(stdout, "angle is %lg\n", angle);
+      fprintf(stdout, "Angle is %lg.\n", angle);
       for (i=0;i<NDIM;i++) {
         Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
@@ -994 +1010 @@
         );
       fprintf(stdout, "\nGive integer factors for length and radius of tube (multiple of %d suggested) : ", ggT);
-      fscanf(stdin, "%d %d", &factors[1], &factors[0]);
+      fscanf(stdin, "%lg %lg", &dfactors[1], &dfactors[0]);
+      for(i=0;i<2;i++)
+        factors[i] = (int)(dfactors[i]);
+      fprintf(stdout, "\nI got %d and %d.\n", factors[1], factors[0]);
       fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n",
         acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
@@ -1323 +1342 @@
               fprintf(SheetFile, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", atombuffer[nr].name, coord[0], coord[1], coord[2]);
               // rotate to align the sheet in xy plane
-              x1 = coord[0]*cos(-angle) + coord[1] * sin(-angle);
-              x2 = coord[0]*(-sin(-angle)) + coord[1] * cos(-angle);
+              x1 = coord[0]*cos(angle) + coord[1] * sin(angle);
+              x2 = coord[0]*(-sin(angle)) + coord[1] * cos(angle);
               x3 = coord[2];
               fprintf(SheetFileAligned, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", atombuffer[nr].name, x1, x2, x3);
@@ -1424 +1443 @@
       fprintf(TubeFile, "%s\t%lg\t%lg\t%lg\n", name, atom_transformed[0], atom_transformed[1], atom_transformed[2]);
       // rotate and flip to align tube in z-direction
-      x1 = atom_transformed[0]*cos(-angle) + atom_transformed[1] * sin(-angle);
-      x2 = atom_transformed[0]*(-sin(-angle)) + atom_transformed[1] * cos(-angle);
+      x1 = atom_transformed[0]*cos(angle) + atom_transformed[1] * sin(angle);
+      x2 = atom_transformed[0]*(-sin(angle)) + atom_transformed[1] * cos(angle);
       x3 = atom_transformed[2];
       fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x3, x2, x1);  // order so that symmetry is along z axis
@@ -1537 +1556 @@
   if (TubeBuffer != NULL) Free(TubeBuffer, "Main: end of stages - TubeBuffer");
 
-  Free(CellFilename, "Main: end of stafes - CellFilename");
-  Free(SheetFilename, "Main: end of stafes - CellFilename");
-  Free(TubeFilename, "Main: end of stafes - CellFilename");
-  Free(TorusFilename, "Main: end of stafes - CellFilename");
-  Free(SheetFilenameAligned, "Main: end of stafes - CellFilename");
-  Free(TubeFilenameAligned, "Main: end of stafes - CellFilename");
+  Free(CellFilename, "Main: end of stages - CellFilename");
+  Free(SheetFilename, "Main: end of stages - CellFilename");
+  Free(TubeFilename, "Main: end of stages - CellFilename");
+  Free(TorusFilename, "Main: end of stages - CellFilename");
+  Free(SheetFilenameAligned, "Main: end of stages - CellFilename");
+  Free(TubeFilenameAligned, "Main: end of stages - CellFilename");
+  Free(filename, "Main: end of stages - filename");
+  Free(stage, "Main: end of stages - stage");
 
   // exit

util/src/NanoCreator.h

@@ -2 +2 @@
 #define NANOCREATOR_H_
 
-#define MYEPSILON 1e-13
+#define MYEPSILON 2e-1
 
 struct Atoms