Changeset 644ba1 for molecuilder/src/builder.cpp

Timestamp:

Aug 18, 2008, 8:57:58 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

04c868

Parents:

e6971b

Message:

Adaptivity fixes, MD by VerletForceIntegration introduced, MD molecule::Trajectories, atom Max::Order, no more recursive going down the fragmentation level

MD
==
molecule::Trajectories is now a map to a struct trajectory list of all the MD steps.
struct Trajectory contains STL vectors of coordinates, velocities and forces. Both are needed for the new VerletForceIntegration.
Parsing of coordinates, velocities and forces from the config file was completely rewritten:

• in the FastParsing case, we just scan for IonType1_1 to find the last step, set the file pointer there and scan all remaining ones
• in the other case, we create the atoms in the first step and put them in a hash for lookup on later steps and read in sequentially (with file pointer moving on).
• This is a lot faster than the old variant.

VerletForceIntegration() implemented in a working manner with force smoothing (mean of actual and last one).
OutputTrajectoriesXYZ() now concatenates the single MD steps into one xyz file, so that the animation can be viewed with e.g. jmol or vmd
config:Deltat is now public (lazy me) and set to 1 instead of 0 initially, also Deltat is parsed accordingly (if not present, defaults to 1)
MatrixContainer::ParseMatrix from parser.cpp is used during VerletForceIntegration() to parse the forces file. Consequently, we have included parser.* in the Makefile.am.
Fix: MoleculeListClass::OutputConfigForListOfFragments() stores config file under config::configpath, before it backup'd the path twice to PathBackup.

Adaptivity
==========
Adaptivity (CheckOrderAtSite()) had compared not absolute, but real values which caused sign problems and faulty behaviour.
Adapatvity (CheckOrderAtSite()) had included atoms by Order (desired one) not by FragOrder (current one) into the list of candidates, which caused faulty behaviour.
CheckOrderAtSite() did single stepping wrong as the mask bit (last in AtomMask) was checked for true not false! Also bit was not set to false initially in FragmentMolecule().
Adaptivity: FragmentMolecule now returns 1 - continue, 2 - don't ... to tell whether we still have to continue with the adaptive cycle (is given as return value from molecuilder)
introduced atom::MaxOrder
StoreOrderAtSiteFile() now also stores the MaxOrder and ParseOrderAtSiteFromFile() parses it back into atom class

Removed Fragmentation Recursion
===============================
As we switched in the prelude of the adaptivity to a monotonous increase from order 1 to the desired one, we don't need to recursively go down each level from a given current bond order, as all these fragments have already been created on the lower orders. Consequently, FragmentBOSSANOVA does not contain NumLevels or alike anymore. This simplified the code a bit, but probably is not yet completely done. (but is working, checked on heptan).
SPFragmentGenerator() does not need Labels anymore, global ones are used for checks. Consequently, PowerSetGenerator() and FragmentSearch structure don't initialise/contain them anymore. We always compare against ...->GetTrueFather()->nr.

File:

• molecuilder/src/builder.cpp (modified) (18 diffs)

molecuilder/src/builder.cpp

@@ -703 +703 @@
     output.open(filename, ios::trunc);
   } 
-  if (mol->OutputXYZ(&output))
-    cout << Verbose(0) << "Saving of XYZ file successful." << endl;
-  else
-    cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  if (mol->MDSteps <= 1) {
+    if (mol->OutputXYZ(&output))
+      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+    else
+      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  } else {
+    if (mol->OutputTrajectoriesXYZ(&output))
+      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+    else
+      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  }
   output.close();
   output.clear();
@@ -733 +740 @@
   string line;
   atom *first;
-  int ExitFlag = 0;
+  bool SaveFlag = false;
+  int ExitFlag = 1;
   int j;
   double volume = 0.;
@@ -758 +766 @@
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
             cout << "\t-O\tCenter atoms in origin." << endl;
             cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
-            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl;
+            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
@@ -767 +776 @@
             cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
-            cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl;
+            cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
@@ -860 +869 @@
           switch(argv[argptr-1][1]) {
             case 'p':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
               if (!mol->AddXYZFile(argv[argptr]))
@@ -875 +884 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
+            case 'D': 
+              {
+                cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
+                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
+                int *MinimumRingSize = NULL;
+                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
+                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
+                Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
+                delete[](MinimumRingSize);
+                if (Subgraphs != NULL) {
+                  while (Subgraphs->next != NULL) {
+                    Subgraphs = Subgraphs->next;
+                    delete(Subgraphs->previous);
+                  }
+                  delete(Subgraphs);
+                }
+              }
+              argptr+=1;
+              break;
             case 'P':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
-              if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1])))
+              if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
                 cout << Verbose(2) << "File not found." << endl;
               else
                 cout << Verbose(2) << "File found and parsed." << endl;
-              argptr+=2;
+              argptr+=1;
               break;
             case 't':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Translating all ions to new origin." << endl;
               for (int i=NDIM;i--;)
@@ -893 +921 @@
               break;
             case 'a':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
               first = new atom;
@@ -908 +936 @@
               break;
             case 's':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
@@ -929 +957 @@
               break;
             case 'b':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
@@ -944 +972 @@
               break;
             case 'c':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
@@ -961 +989 @@
               break;
             case 'O':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Centering atoms in origin." << endl;
               mol->CenterOrigin((ofstream *)&cout, &x);
@@ -967 +995 @@
               break;
             case 'r':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
               break;
             case 'F':
             case 'f':
-              if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
               cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
               if (argc >= argptr+2) {
@@ -980 +1007 @@
                 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
                 if (mol->first->next != mol->last) {
-                  mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
+                  ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
                 }
                 end = clock();
@@ -990 +1017 @@
               break;
             case 'm':
-              ExitFlag = 1;
               j = atoi(argv[argptr++]);
               if ((j<0) || (j>1)) {
@@ -996 +1022 @@
                 j = 0;
               }
-              if (j)
+              if (j) {
+                SaveFlag = true;
                 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
-              else
+              } else
                 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
               mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
               break;
             case 'o':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(0) << "Evaluating volume of the convex envelope.";
               VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
@@ -1013 +1040 @@
               {
                 double density;
-                ExitFlag = 1;
+                SaveFlag = true;
                 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
                 density = atof(argv[argptr++]);
@@ -1030 +1057 @@
               break;
             case 'd':
+              SaveFlag = true;
               for (int axis = 1; axis <= NDIM; axis++) {
                 int faktor = atoi(argv[argptr++]);
@@ -1089 +1117 @@
       } else argptr++;
     } while (argptr < argc);
-    if (ExitFlag == 1) // 1 means save and exit
+    if (SaveFlag)
       SaveConfig(ConfigFileName, &configuration, periode, mol);
-    if (ExitFlag >= 1) { // 2 means just exit
+    if ((ExitFlag >= 1)) {
       delete(mol);
       delete(periode);
-      return (1);
+      return (ExitFlag);
     }
   } else {  // no arguments, hence scan the elements db