Changeset d067d45 for src/molecules.cpp

Timestamp:

Jul 23, 2009, 1:45:24 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

51c910

Parents:

ce5ac3 (diff), 437922 (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.

git-author:

Frederik Heber <heber@…> (07/23/09 12:34:47)

git-committer:

Frederik Heber <heber@…> (07/23/09 13:45:24)

Message:

Merge branch 'MultipleMolecules'

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/atom.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/helpers.cpp
molecuilder/src/joiner.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp
molecuilder/src/vector.cpp
molecuilder/src/verbose.cpp

merges:

• analyzer.cpp: all from HEAD, Hessian stuff
• atom.cpp: ostream & operator << (ostream &ost, const atom &a), 2nd param is const
• boundary.cpp: DoSingleStepOutput, VRMLSUffix, HULLEPSILON, Get_center_of_sphere(), Find_next_suitable_point_via_Angle_of_Sphere() and all of tesselation stuff from MultipleMolecules; DoRaster3DOutput, ~BoundaryLineSet(), ~Tesselation() from HEAD; write_raster3d_file with some changes from MultipleMolecules
• boundary.hpp: all from MultipleMolecules
• builder.cpp: all in ParseCommandLineOptions() from MultipleMolecules
• config.cpp: all from HEAD, is stuff from ThermoStat and ConstrainedMD
• datacreator.hpp: all from HEAD, Hessian stuff
• helpers.cpp: all from HEAD, ReAlloc() differed just by indentation
• moleculelist.cpp: all from MultipleMolecules, has to do with introduced molecule centers
• molecules.cpp: all from HEAD, was ConstrainedMD stuff
• molecules.hpp: ConstrainedMD and VerletForceIntegration from HEAD, moleculelist stuff from MultipleMolecules
• parser.cpp: all from HEAD, ColumnCounter is *ColumnCounter there
• parser.hpp: all from HEAD
• vector.cpp: completely from MultipleMolecules (more lines)
• verbose.cpp: all from HEAD

compilation fixes:

• CheckPresenceOfTriangle() from MultipleMolecules
• MoleculeListClass::ListOfMolecules has changed to STL list, hence use insert(), ReturnIndex() and MoleculeListClass::OutputConfigForListOfFragments
• MoleculeListClass::OutputConfigForListOfFragments() had fragmentprefix instead of FRAGMENTPREFIX
• parser.hpp: EnergyMatrix::ParseIndices() declaration was missing

File:

• src/molecules.cpp (modified) (27 diffs)

src/molecules.cpp

@@ -62 +62 @@
   cell_size[0] = cell_size[2] = cell_size[5]= 20.;
   cell_size[1] = cell_size[3] = cell_size[4]= 0.;
+  strcpy(name,"none");
+  IndexNr  = -1;
+  ActiveFlag = false;
 };
 
@@ -598 +601 @@
 };
 
+/** Set molecule::name from the basename without suffix in the given \a *filename.
+ * \param *filename filename
+ */
+void molecule::SetNameFromFilename(const char *filename)
+{
+  int length = 0;
+  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
+  char *endname = strchr(molname, '.');
+  if ((endname == NULL) || (endname < molname))
+    length = strlen(molname);
+  else
+    length = strlen(molname) - strlen(endname);
+  strncpy(name, molname, length);
+  name[length]='\0';
+};
+
 /** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
  * \param *dim vector class
@@ -652 +671 @@
     max->Scale(0.5);
     Translate(max);
+    Center.Zero();
   }
 
@@ -691 +711 @@
     max->AddVector(min);
     Translate(min);
+    Center.Zero();
   }
   delete(min);
@@ -700 +721 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterOrigin(ofstream *out, Vector *center)
+void molecule::CenterOrigin(ofstream *out)
 {
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
 
-  for(int i=NDIM;i--;) // zero center vector
-    center->x[i] = 0.;
+  Center.Zero();
 
   if (ptr != end) {   //list not empty?
@@ -712 +732 @@
       ptr = ptr->next;
       Num++;
-      center->AddVector(&ptr->x);
-    }
-    center->Scale(-1./Num); // divide through total number (and sign for direction)
-    Translate(center);
+      Center.AddVector(&ptr->x);
+    }
+    Center.Scale(-1./Num); // divide through total number (and sign for direction)
+    Translate(&Center);
+    Center.Zero();
   }
 };
@@ -780 +801 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterGravity(ofstream *out, Vector *center)
-{
-  if (center == NULL) {
-    DetermineCenter(*center);
-    Translate(center);
-    delete(center);
-  } else {
-    Translate(center);
-  }
+void molecule::CenterPeriodic(ofstream *out)
+{
+  DeterminePeriodicCenter(Center);
+};
+
+/** Centers the center of gravity of the atoms at (0,0,0).
+ * \param *out output stream for debugging
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterAtVector(ofstream *out, Vector *newcenter)
+{
+  Center.CopyVector(newcenter);
 };
 
@@ -837 +861 @@
 
 /** Determines center of molecule (yet not considering atom masses).
- * \param Center reference to return vector
- */
-void molecule::DetermineCenter(Vector &Center)
+ * \param center reference to return vector
+ */
+void molecule::DeterminePeriodicCenter(Vector &center)
 {
   atom *Walker = start;
@@ -849 +873 @@
 
   do {
-    Center.Zero();
+    center.Zero();
     flag = true;
     while (Walker->next != end) {
@@ -876 +900 @@
         Testvector.AddVector(&Translationvector);
         Testvector.MatrixMultiplication(matrix);
-        Center.AddVector(&Testvector);
+        center.AddVector(&Testvector);
         cout << Verbose(1) << "vector is: ";
         Testvector.Output((ofstream *)&cout);
@@ -889 +913 @@
             Testvector.AddVector(&Translationvector);
             Testvector.MatrixMultiplication(matrix);
-            Center.AddVector(&Testvector);
+            center.AddVector(&Testvector);
             cout << Verbose(1) << "Hydrogen vector is: ";
             Testvector.Output((ofstream *)&cout);
@@ -900 +924 @@
   } while (!flag);
   Free((void **)&matrix, "molecule::DetermineCenter: *matrix");
-  Center.Scale(1./(double)AtomCount);
+  center.Scale(1./(double)AtomCount);
 };
 
@@ -913 +937 @@
   Vector *CenterOfGravity = DetermineCenterOfGravity(out);
 
-  CenterGravity(out, CenterOfGravity);
+  CenterAtVector(out, CenterOfGravity);
 
   // reset inertia tensor
@@ -1409 +1433 @@
 bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
 {
+  molecule *mol = NULL;
   bool status = true;
   int MaxSteps = configuration.MaxOuterStep; 
-  MoleculeListClass *MoleculePerStep = new MoleculeListClass(MaxSteps+1, AtomCount);
+  MoleculeListClass *MoleculePerStep = new MoleculeListClass();
   // Get the Permutation Map by MinimiseConstrainedPotential
   atom **PermutationMap = NULL;
@@ -1443 +1468 @@
   *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
   for (int step = 0; step <= MaxSteps; step++) {
-    MoleculePerStep->ListOfMolecules[step] = new molecule(elemente);
+    mol = new molecule(elemente);
+    MoleculePerStep->insert(mol);
     Walker = start;
     while (Walker->next != end) {
       Walker = Walker->next;
       // add to molecule list
-      Sprinter = MoleculePerStep->ListOfMolecules[step]->AddCopyAtom(Walker);
+      Sprinter = mol->AddCopyAtom(Walker);
       for (int n=NDIM;n--;) {
         Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
@@ -1467 +1493 @@
   for (int i=AtomCount; i--; )
     SortIndex[i] = i;
-  status = MoleculePerStep->OutputConfigForListOfFragments(out, "ConstrainedStep", &configuration, SortIndex, false, false);
+  status = MoleculePerStep->OutputConfigForListOfFragments(out, &configuration, SortIndex);
   
   // free and return
@@ -1828 +1854 @@
 };
 
-/** Removes atom from molecule list.
+/** Removes atom from molecule list and deletes it.
  * \param *pointer atom to be removed
  * \return true - succeeded, false - atom not found in list
@@ -1842 +1868 @@
   Trajectories.erase(pointer);
   return remove(pointer, start, end);
+};
+
+/** Removes atom from molecule list, but does not delete it.
+ * \param *pointer atom to be removed
+ * \return true - succeeded, false - atom not found in list
+ */
+bool molecule::UnlinkAtom(atom *pointer)
+{
+  if (pointer == NULL)
+    return false;
+  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
+    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
+  else
+    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
+    ElementCount--;
+  Trajectories.erase(pointer);
+  unlink(pointer);
+  return true;
 };
 
@@ -2368 +2413 @@
             };
             AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-            
+
           }
 
@@ -2484 +2529 @@
                         //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
                         /// \todo periodic check is missing here!
-                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
+                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                         MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
                         MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
                         MaxDistance = MinDistance + BONDTHRESHOLD;
                         MinDistance -= BONDTHRESHOLD;
-                        distance =   OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size);
+                        distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
                         if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
-                          //*out << Verbose(0) << "Adding Bond between " << *Walker << " and " << *OtherWalker << "." << endl;
+                          //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
                           AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
-                          BondCount++;
                         } else {
                           //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
@@ -2564 +2608 @@
       *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
     *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-    
+
     // output bonds for debugging (if bond chain list was correctly installed)
     *out << Verbose(1) << endl << "From contents of bond chain list:";
@@ -3171 +3215 @@
           cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
         }
-        //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
       }
     }
@@ -3731 +3774 @@
   //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
     // allocate memory for the pointer array and transmorph graphs into full molecular fragments
-    BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
+    BondFragments = new MoleculeListClass();
     int k=0;
     for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
       KeySet test = (*runner).first;
       *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
-      BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
+      BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
       k++;
     }
-    *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
+    *out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;
 
     // ===== 9. Save fragments' configuration and keyset files et al to disk ===
-    if (BondFragments->NumberOfMolecules != 0) {
+    if (BondFragments->ListOfMolecules.size() != 0) {
       // create the SortIndex from BFS labels to order in the config file
       CreateMappingLabelsToConfigSequence(out, SortIndex);
 
-      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-      if (BondFragments->OutputConfigForListOfFragments(out, FRAGMENTPREFIX, configuration, SortIndex, true, true))
+      *out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
+      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
         *out << Verbose(1) << "All configs written." << endl;
       else
@@ -5318 +5361 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenter(CenterOfGravity);
-    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
+    DeterminePeriodicCenter(CenterOfGravity);
+    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
     *out << Verbose(5) << "Center of Gravity: ";
     CenterOfGravity.Output(out);
@@ -5325 +5368 @@
     OtherCenterOfGravity.Output(out);
     *out << endl;
-    if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
+    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
       *out << Verbose(4) << "Centers of gravity don't match." << endl;
       result = false;
@@ -5339 +5382 @@
     while (Walker->next != end) {
       Walker = Walker->next;
-      Distances[Walker->nr] = CenterOfGravity.Distance(&Walker->x);
+      Distances[Walker->nr] = CenterOfGravity.DistanceSquared(&Walker->x);
     }
     Walker = OtherMolecule->start;
     while (Walker->next != OtherMolecule->end) {
       Walker = Walker->next;
-      OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
+      OtherDistances[Walker->nr] = OtherCenterOfGravity.DistanceSquared(&Walker->x);
     }
 
@@ -5362 +5405 @@
     flag = 0;
     for (int i=0;i<AtomCount;i++) {
-      *out << Verbose(5) << "Distances: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
-      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold)
+      *out << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
+      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
         flag = 1;
     }