first version of basinhopping is running; results are incorrect

basinhopping.cc

@@ -7,11 +7,12 @@
 #include "potential.h"
 #include "iop.h"
 #include "parameter.h"
+#include "lina.h"
 
 using namespace std;
 
 
-int BasinHopping::run (unique_ptr<AcceptanceTest>& accept) 
+int BasinHopping::run () 
 {
 
     vector<double> p;
@@ -49,30 +50,34 @@ int BasinHopping::run (unique_ptr<AcceptanceTest>& accept)
         _currentStep.setEnergy(potential->calcEnergy(_currentStep)); 
 
         ofstream dummy;
+        dummy.open("out");
         _currentStep = potential->optimize(dummy, _currentStep, switches, opt); 
+        dummy.close();
+
+        vector< vector<double> >hessian = potential->calcHessian(_currentStep);
+        vector<double> eigenValues = diag(hessian);
+        _currentStep.setHessian(eigenValues);
 
-        if (this->checkConf())
+        if (_currentStep.isMinimum()) //check for true minimum
         {
-            double oldE, newE = _currentStep.getEnergy();
+            if (this->checkConf()) //check spherical container
+            {
+                double oldE, newE = _currentStep.getEnergy();
 
-            if (_iteration == 1) {oldE = newE;}
-            else {oldE = _previousStep.getEnergy();}
+                if (_iteration == 1) {oldE = newE;}
+                else {oldE = _previousStep.getEnergy();}
 
-            if (this->acceptStep(oldE, newE, accept)) 
-            {
-                accepted = true;
-            }
-            else
-            {
-                _currentStep = _previousStep;
+                if (this->acceptStep(oldE, newE)) //accept the step
+                {
+                    accepted = true;
+                    _uniqueStructures.addCluster(_currentStep);
+                }
+                else {_currentStep = _previousStep;} //rejected
             }
-        }
-        else
-        {
-            //cerr << "Cluster not within spherical container" << endl;
-            _currentStep = _previousStep;
+            else{_currentStep = _previousStep;} //rejected because spherical container
         }
 
+
         cout.precision(5);
         cout << fixed << left << setprecision(10);
         cout << setfill(' ')  << left <<
@@ -86,6 +91,15 @@ int BasinHopping::run (unique_ptr<AcceptanceTest>& accept)
 
     }
 
+    int n(0);
+    for (auto& i : _uniqueStructures)
+    {
+        cout << i.first << endl;
+        xyzout(i.second, to_string(n) + ".xyz");
+        n++;
+    }
+
+
     return 0;
 }
 
@@ -109,9 +123,9 @@ void BasinHopping::propagate()
     _currentStep.setCoordinates(coordinates);
 }
 
-bool BasinHopping::acceptStep (double oldE, double newE, unique_ptr<AcceptanceTest>& accept)
+bool BasinHopping::acceptStep (double oldE, double newE)
 {
-    return (*accept)(oldE, newE);
+    return (*_accept)(oldE, newE);
 }
 
 

basinhopping.h

@@ -4,25 +4,31 @@
 #include "structure.h"
 #include "potential.h"
 #include "acceptanceTest.h"
+#include "storage.h"
 
 class BasinHopping
 {
     public:
         BasinHopping(   structure initialCoordinates,
+                        std::shared_ptr<AcceptanceTest>& accept, 
                         int nsteps = 1000) :  
             _initialCoordinates(initialCoordinates),
             _previousStep(initialCoordinates),
             _currentStep(initialCoordinates),
             _size(initialCoordinates.nAtoms()),
             _iteration(0),
-            _nsteps(nsteps)
+            _accept(accept),
+            _nsteps(nsteps),
+            _uniqueStructures()
         {}
 
-        int run (std::unique_ptr<AcceptanceTest>& accept);
+        int nStructures() {return _uniqueStructures.getSize();}
+
+        int run ();
 
     private:
         bool checkConf();
-        bool acceptStep (double oldE, double newE, std::unique_ptr<AcceptanceTest>& accept);
+        bool acceptStep (double oldE, double newE);
         void propagate();
 
         const structure _initialCoordinates;
@@ -31,6 +37,8 @@ class BasinHopping
         const int _size;
         int _iteration;
         const int _nsteps;
+        std::shared_ptr<AcceptanceTest> _accept;
+        Storage _uniqueStructures;
 };
 
 #endif

global.cc

@@ -10,15 +10,19 @@ int main ()
 {
     timer T;
 
-    structure cluster(10); 
+    structure cluster(8); 
 
-    BasinHopping hop(cluster);
 
-    unique_ptr<AcceptanceTest> accept;
+    shared_ptr<AcceptanceTest> accept_shared;
     Metropolis* metro = new Metropolis(10);
-    accept.reset(metro);
+    accept_shared.reset(metro);
+
+    BasinHopping hop(cluster,accept_shared,1000);
+
+    hop.run();
+
+    cout << hop.nStructures() << endl;
 
-    hop.run(accept);
 
     cout << "Time: " << T.total_timing() << " s" << endl;
 

storage.cc

@@ -0,0 +1,20 @@
+#include <iostream>
+#include "storage.h"
+#include "structure.h"
+
+using namespace std;
+
+bool Storage::addCluster (structure &S)
+{
+    double energy = S.getEnergy();
+
+    energyStructureMap::iterator iter = _mapping.find(energy);
+
+    if (iter == _mapping.end())
+    {
+        _mapping[energy] = S;
+        return true;
+    }
+
+    return false;
+}

storage.h

@@ -0,0 +1,32 @@
+#ifndef STORAGE
+#define STORAGE
+
+#include <map>
+#include "structure.h"
+
+struct compareDouble
+{
+    compareDouble(double eps = 1e-2) : _eps(eps) {}
+    double _eps;
+    bool operator() (const double a, const double b) const
+    {
+       if (b - a > _eps) return true; 
+       return false;
+    }
+};
+
+class Storage
+{
+    private:
+        typedef std::map<double, structure, compareDouble> energyStructureMap;
+        energyStructureMap _mapping;
+
+    public:
+        Storage() : _mapping() {}
+        bool addCluster (structure &S);
+        int getSize() {return _mapping.size();}
+
+        energyStructureMap::iterator begin() {return _mapping.begin();}
+        energyStructureMap::iterator end() {return _mapping.end();}
+};
+#endif