very first version of eigenvector following transition state search

eigenvectorFollowing.cc

@@ -0,0 +1,134 @@
+#include <cmath>
+#include "eigenvectorFollowing.h"
+#include "iop.h"
+
+using namespace std;
+
+void EigenvectorFollowing::updateStructure (structure &S)
+{
+    _gradientIsCurrent = false;
+    _currentStructure = S;
+}
+
+column_vector EigenvectorFollowing::getTrueGradient ()
+{
+    if (_gradientIsCurrent) 
+        return _trueGradient;
+    else
+    {
+        _gradientIsCurrent = true;
+        _trueGradient = _tPot.getTrueGradient(_currentStructure);
+        return _trueGradient;
+    }
+}
+
+structure EigenvectorFollowing::stepUphill (structure &S)
+{
+    double h;
+    column_vector v = this->_tPot.getVector();
+        double eval = fabs(_tPot.get_eval());
+    if (_currentIter == 0) h = 0.3;
+    else
+    {
+        column_vector g = getTrueGradient();
+        double overlap = dot(g,v);
+        h = 2 * overlap / ( eval * ( 1 + sqrt( 1 + 4 * pow(2, overlap / eval) ) ) );
+    }
+    cout << "h: " << h << endl << "eval: " << eval << endl;
+
+    vector<coord3d> newCoord(S.nAtoms());
+    for (int i = 0; i < S.nAtoms(); i++)
+        for (int j = 0; j < 3; j++)
+            newCoord[i][j] = S[i][j] + h * v(3 * i + j);
+
+    structure newS(0, newCoord, false);
+    xyzout(newS, "uphill.xyz");
+
+    return newS;
+}
+
+void EigenvectorFollowing::run (unsigned int n) 
+{
+        for (unsigned int i = 0; i < n; i++)
+        {
+            _currentIter = i;
+            _tPot.calcTransverseDirection(_currentStructure);
+            structure trialS = stepUphill(_currentStructure);
+            structure newS = _tPot.optimize(trialS);
+
+            _previousStructure = _currentStructure;
+            updateStructure(newS);
+        }
+}
+
+//------------------------------------------------------------------------------//
+//                           main for testing                                   //
+//------------------------------------------------------------------------------//
+
+
+int main ()
+{
+
+    LJ *potential = new LJ();
+    TransversePotential T(*potential);
+
+    vector<coord3d> coords;
+    /*
+    coords.push_back(coord3d(0.25877219650832,-0.51611072802221,0));
+    coords.push_back(coord3d(0.31757890372974,0.48215865980443,0));
+    coords.push_back(coord3d(-0.57635110023806,0.033952068217783,0));
+    */
+    coords.push_back(coord3d(1,0,0));
+    coords.push_back(coord3d(-1,0,0));
+    coords.push_back(coord3d(0,1,0));
+    coords.push_back(coord3d(0,-1,0));
+    coords.push_back(coord3d(0,0,1));
+    coords.push_back(coord3d(0,0,-1));
+    structure S1(1,coords);
+
+    ofstream dummy;
+    structure S2 = potential->optimize(dummy,S1);
+
+    vector< vector<double> > hessian = potential->calcHessian(S2);
+    vector<double> eval = diag(hessian);
+    for (auto& i : eval) cout << "e " << i << endl;
+    xyzout(S2, "S2.xyz");
+
+
+
+
+    //std::cout << T.getEnergy(S) << std::endl;
+    //std::cout << "g2: " << T.getTrueGradient(S2) << std::endl;
+
+    EigenvectorFollowing E(T, S2);
+    E.run(1);
+
+    structure newS = E.getCurrentStructure();
+    xyzout(newS, "newS.xyz");
+
+    /*
+    vector<coord3d> coordinates = newS.getCoordinates();
+    for (auto& i : coordinates)
+        cout << i << endl;
+    */
+
+    /*
+    column_vector gradient = potential->calcGradient(newS);
+    column_vector gradientT = T.getTransverseGradient(newS);
+    vector< vector<double> > hessian = potential->calcHessian(newS);
+    vector<double> eval = diag(hessian);
+    */
+
+    //xyzout(newS, "newS.xyz");
+
+    /*
+    for (auto& i : gradientT)
+        cout << "gT: " << i << endl;
+    for (auto& i : gradient)
+        cout << "g: " << i << endl;
+    for (auto& i : eval)
+        cout << "H2: " << i << endl;
+    */
+
+    return 0;
+}

eigenvectorFollowing.h

@@ -0,0 +1,34 @@
+#include "transversePotential.h"
+#include "structure.h"
+
+class EigenvectorFollowing
+{
+    private:
+        TransversePotential _tPot;
+        structure _previousStructure;
+        structure _currentStructure; 
+        double _max_stepsize;
+        column_vector _trueGradient;
+        bool _gradientIsCurrent; 
+
+        unsigned int _currentIter;
+
+    public:
+        EigenvectorFollowing (TransversePotential &tPot, structure &S) :
+            _tPot(tPot),
+            _previousStructure(S),
+            _currentStructure(S),
+            _max_stepsize(1),
+            _trueGradient(0),
+            _gradientIsCurrent(false),
+            _currentIter(0)
+        {}
+
+        void updateStructure (structure &S);
+        column_vector getTrueGradient();
+        structure getCurrentStructure() { return _currentStructure; }
+        structure stepUphill (structure &S);
+
+
+        void run(unsigned int n);
+};

potential.cc

@@ -157,21 +157,17 @@ vector< vector<double> > pairPotential::calcHessian (structure &S)
     return hessianMatrix;
 }
 
-vector<double> pairPotential::getLowestEvec (vector< pair< double,vector<double> > > V)
+pair<double,vector<double> > pairPotential::getLowestEvec (vector< pair< double,vector<double> > > &V)
 {
-    vector<double> evec;
     sort(V.begin(),V.end());
     for (auto& i : V)
     {
-        cout << i.first << endl;
-        if (i.first > 0.001) 
+        if (i.first > 0.001 || i.first < -0.001) 
         {
-            evec = i.second;
+            return i;
             break;
         }
     }
-
-    return evec;
 }
 
 /*----------------------------------Optimization----------------------------------------*/
@@ -180,9 +176,6 @@ structure pairPotential::optimize (ostream &min, structure &S)
 {
     min.precision(16);
 
-    const size_t nsteps = static_cast<const size_t>(_nsteps);
-
-
     column_vector x((S.nAtoms()) * 3);
     for (int i = 0; i < S.nAtoms(); i++)
     {
@@ -204,7 +197,7 @@ structure pairPotential::optimize (ostream &min, structure &S)
                 try
                 {
                     dlib::find_min( dlib::bfgs_search_strategy(),
-                                    dlib::stop_strategy(_stop_crit, nsteps).be_verbose(min),
+                                    dlib::stop_strategy(_stop_crit, _nsteps).be_verbose(min),
                                     f, df, x, -(S.nAtoms()) * 1000);
                 }
                 catch (std::exception &e)
@@ -220,7 +213,7 @@ structure pairPotential::optimize (ostream &min, structure &S)
                 try
                 {
                     dlib::find_min( dlib::cg_search_strategy(),
-                                    dlib::stop_strategy(_stop_crit, nsteps).be_verbose(min),
+                                    dlib::stop_strategy(_stop_crit, _nsteps).be_verbose(min),
                                     f, df, x, -(S.nAtoms()) * 1000);
                 }
                 catch (std::exception &e)
@@ -257,7 +250,6 @@ structure pairPotential::optimize (ostream &min, structure &S)
     newS.setEnergy(finalEnergy);
     newS.setConverged();
 
-
     return newS;
 }
 

potential.h

@@ -2,6 +2,7 @@
 #define POTENTIAL
 
 #include <dlib/optimization.h>
+#include "stop_strategy.h"
 #include "structure.h"
 
 
@@ -21,7 +22,7 @@ class pairPotential
     protected:
         pairPotential (
                         const int algo_switch = 1, 
-                        const double stop_crit = 1e-3,
+                        const double stop_crit = 1e-15,
                         const int nsteps = 1000         ) 
                         : 
                         _algo_switch(algo_switch), 
@@ -38,9 +39,13 @@ class pairPotential
         const column_vector calcGradient (const column_vector &v);
         const column_vector calcGradient (structure &S);
         std::vector< std::vector<double> > calcHessian (structure &S);
+
+        const int getAlgoSwitch() const {return _algo_switch;}
+        const double getStopCrit() const {return _stop_crit;}
+        const int getNsteps() const {return _nsteps;}
 
         structure optimize (std::ostream &min, structure &S);
-        std::vector<double> getLowestEvec (std::vector< std::pair< double,std::vector<double> > > V);
+        std::pair<double,std::vector<double> > getLowestEvec (std::vector< std::pair< double,std::vector<double> > > &V);
 };
 
 

transversePotential.cc

@@ -0,0 +1,130 @@
+#include <dlib/optimization.h>
+#include <cassert>
+#include "stop_strategy.h"
+#include "potential.h"
+#include "transversePotential.h"
+#include "lina.h"
+#include "iop.h"
+
+using namespace std;
+
+void TransversePotential::setVector (double &eval, column_vector &v) 
+{ 
+    _eval = eval;
+    _vector = v / sqrt(dot(v, v)); 
+}
+
+void TransversePotential::setVector (double &eval, vector<double> &v)
+{
+    _eval = eval;
+    column_vector vector(v.size());
+    for (int i = 0; i < v.size() ; i++) 
+    {
+        vector(i) = v[i];
+    }
+
+    _vector = vector / sqrt(dot(vector, vector));
+}
+
+column_vector TransversePotential::getTransverseGradient (const column_vector &v)
+{
+    assert(_vector.size() == v.size());
+    const column_vector g = getTrueGradient(v);
+    assert(_vector.size() == g.size());
+
+    return g - dot(g, _vector) * _vector; 
+
+}
+
+column_vector TransversePotential::getTransverseGradient (structure &S)
+{
+    assert(_vector.size()); 
+    column_vector g = getTrueGradient(S);
+    assert(_vector.size() == g.size());
+
+    return g - dot(g, _vector) * _vector;
+}
+
+void TransversePotential::calcTransverseDirection( structure &S )
+{
+    vector< vector<double> > hessian = _potential->calcHessian(S);
+
+    vector<pair<double, vector<double> > > eval_evec = diagv(hessian);
+
+    pair<double,vector<double> > lowest_eval_evec = _potential->getLowestEvec(eval_evec);
+    setVector(lowest_eval_evec.first, lowest_eval_evec.second);
+}
+
+
+structure TransversePotential::optimize (structure &S)
+{
+    //calcTransverseDirection(S);
+
+    column_vector x((S.nAtoms()) * 3);
+    for (int i = 0; i < S.nAtoms(); i++)
+    {
+        for (int j = 0; j < 3; j++)
+        {
+            x(3 * i + j) = S[i][j];
+        }
+    }
+
+
+    auto f = [this] (column_vector v) -> const double {return this->getEnergy(v);};
+    auto df = [this] (column_vector v) -> column_vector {return this->getTransverseGradient(v);};
+
+    try
+    {
+        switch (_potential->getAlgoSwitch())
+        {
+            case 1:
+            {
+                dlib::find_min( dlib::bfgs_search_strategy(),
+                                dlib::stop_strategy(_potential->getStopCrit(), _potential->getNsteps()).be_verbose(),
+                                f, df, x, -(S.nAtoms()) * 1000);
+                break;
+            }
+            case 2:
+            {
+                dlib::find_min( dlib::cg_search_strategy(),
+                                dlib::stop_strategy(_potential->getStopCrit(), _potential->getNsteps()),
+                                f, df, x, -(S.nAtoms()) * 1000);
+                break;
+            }
+            default:
+            {
+                cerr << "Bad input of algorithm name" << endl;
+                structure newS(S.getNumber(), S.getCoordinates());
+                return newS;
+            }
+        }
+    }
+    catch (std::exception &e)
+    {
+        structure newS(S.getNumber(), S.getCoordinates());
+        return newS;
+    }
+
+    vector<coord3d> newCoordinates;
+    for (long i = 0; i < x.size() / 3; i++)
+    {
+        coord3d sphere (x(3 * i), x(3 * i + 1), x(3 * i + 2));
+        newCoordinates.push_back(sphere);
+    }
+
+    structure newS(S.getNumber(), newCoordinates);
+    double finalEnergy = this->getEnergy(x);
+
+    newS.setEnergy(finalEnergy);
+    newS.setConverged();
+
+    cout << "gT: " << dlib::length(this->getTransverseGradient(x)) << endl;
+    cout << "g: " << dlib::length(this->getTrueGradient(x)) << endl;
+    vector< vector<double> > hessian = _potential->calcHessian(newS);
+    vector<double> eval = diag(hessian);
+    for (auto& i : eval)
+        cout << "H: " << i << endl;
+
+    return newS;
+}
+