ElasticBand.cpp

#include "ElasticBand.h"
#include "CVs/CVManager.h"
#include "Snapshot.h" 
#include <math.h>
#include <iostream>

namespace SSAGES
{
    // Post-integration hook.
    void ElasticBand::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
    {
        auto& forces = snapshot->GetForces();
        auto cvs = cvmanager.GetCVs(cvmask_);

        // Apply umbrella to cvs
        for(size_t i = 0; i < cvs.size(); ++i)
        {
            // Get current CV and gradient.
            auto& cv = cvs[i];
            auto& grad = cv->GetGradient();

            // Compute dV/dCV.
            auto diff = cv->GetDifference(centers_[i]);
            auto D = cvspring_[i] * diff;

            // Update forces.
            for(size_t j = 0; j < forces.size(); ++j)
                for(int k = 0; k < forces[j].size(); ++k)
                    forces[j][k] -= D*grad[j][k];

            // If not equilibrating and has evolved enough steps,
            // generate the gradient
            if(iterator_ >= equilibrate_ && iterator_ % evolution_ == 0)
            {
                newcenters_[i] += diff;
                nsampled_++;
            }
        }

        // Restart iteration and zero gradients when moving onto
        // next elastic band iteration
        if(iterator_ > (equilibrate_ + evolution_ * nsamples_))
        {   
            StringUpdate();
            CheckEnd(cvs);
            UpdateWorldString(cvs); 
            PrintString(cvs); 

            iterator_ = 0;

            for(size_t ii = 0; ii < newcenters_.size(); ii++)
                newcenters_[ii] = 0;

            nsampled_ = 0;
            iteration_++;
        }

        iterator_++;
    }

    void ElasticBand::StringUpdate()
    {
        double dot=0, norm=0;
        std::vector<double> lcv0, ucv0, tngt;
        lcv0.resize(centers_.size(), 0);
        ucv0.resize(centers_.size(), 0);

        GatherNeighbors(&lcv0, &ucv0);

        tngt.resize(centers_.size());

        for(size_t ii = 0; ii<centers_.size(); ii++)
        {
            tngt[ii] = ucv0[ii] - lcv0[ii];
            norm+=tngt[ii]*tngt[ii];
        }

        norm=sqrt(norm);

        for(size_t ii = 0; ii < centers_.size(); ii++)  {
            tngt[ii] /= norm;
            newcenters_[ii] /= ((double) nsampled_);
            dot+=newcenters_[ii]*tngt[ii];
        }
        
        // Evolution of the images and reparametrirized of the string
        for(size_t ii = 0; ii < centers_.size(); ii++)
        {
            if((mpiid_ != 0) && (mpiid_ != world_.size()-1))
            {
                newcenters_[ii] -= dot*tngt[ii];
                centers_[ii] = centers_[ii] + tau_ * 
                (newcenters_[ii] + stringspring_ * (ucv0[ii] + lcv0[ii] - 2 * centers_[ii]));
            }
            else
            {
                centers_[ii] = centers_[ii] + tau_ * newcenters_[ii];
            }
        }
    }
}