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] += D;
                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, lnorm = 0, unorm = 0;
        std::vector<double> lcv0, ucv0, tngt, ltngt, utngt;
        lcv0.resize(centers_.size(), 0);
        ucv0.resize(centers_.size(), 0);
 
        GatherNeighbors(&lcv0, &ucv0);
 
        tngt.resize(centers_.size());
        ltngt.resize(centers_.size());
        utngt.resize(centers_.size());
 
        for(size_t ii = 0; ii < centers_.size(); ii++)
        {
            utngt[ii] = ucv0[ii] - centers_[ii];
            ltngt[ii] = centers_[ii] - lcv0[ii];
 
            unorm += utngt[ii]*utngt[ii];
            lnorm += ltngt[ii]*ltngt[ii];
        }
 
        unorm = sqrt(unorm);
        lnorm = sqrt(lnorm);
 
        for(size_t ii = 0; ii<centers_.size(); ii++)
        {
            tngt[ii] = ltngt[ii]/lnorm + utngt[ii]/unorm;
            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 within the band
        // Endpoints evolve due to gradient alone.
        for(size_t ii = 0; ii < centers_.size(); ii++)
        {
            if((mpiid_ != 0) && (mpiid_ != world_.size()-1))
            {
                //subtract out tangent from "real" force
                newcenters_[ii] -= dot*tngt[ii];
 
                //centers evolve according to perpendicular "real" force and parallel "spring" force; reuses tngt def.
                centers_[ii] += tau_ * (newcenters_[ii] + stringspring_ * (unorm - lnorm) * tngt[ii]);
            }
            else
            {
                centers_[ii] += tau_ * newcenters_[ii];
            }
        }
    }
}