ABOPotential.h

//
//  Copyright (C) 2017, Alexander Stukowski and Paul Erhart
//
//  This file is part of atomicrex.
//
//  Atomicrex is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 3 of the License, or
//  (at your option) any later version.
//
//  Atomicrex is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
//

#pragma once

#include "Potential.h"
#include <cmath>
#include "../dof/DegreeOfFreedom.h"

namespace atomicrex {

class ABOPotential : public Potential
{
public:

    ABOPotential(const FPString& id, FitJob* job, const FPString& tag = FPString("Tersoff"));

    virtual double cutoff() const override { return _cutoff; }

    virtual double computeEnergyAndForces(AtomicStructure& structure, NeighborList& neighborList) override;

    virtual double computeEnergy(AtomicStructure& structure, NeighborList& neighborList) override;

    virtual void outputResults() override;

    void writePotential(const FPString& filename) const;

public:

    virtual void parse(XML::Element potentialElement) override;

protected:

    struct ABOParamSet
    {
        ScalarDOF r0, D0, beta, S;      
        ScalarDOF gamma, c, d, h, twomu;    
        ScalarDOF beta2;            
        ScalarDOF powern;           
        double powerm;              
        ScalarDOF bigr;             
        ScalarDOF bigd;             

        double c1,c2,c3,c4;
        int powermint;

        ABOParamSet() :
            r0("r0", 0.0, 0.0),
            D0("D0", 0.0, 0.0),
            beta("beta"),
            S("S", 0.0, 0.0),
            gamma("gamma", 0.0, 0.0),
            c("c"),
            d("d", 0.0, 0.0),
            h("h", 0.0, 0.0),
            twomu("twomu", 0.0, 0.0),
            beta2("beta2", 0.0, 0.0),
            powern("powern", 0.0, 0.0),
            bigr("bigr", 0.0, 0.0),
            bigd("bigd", 0.0, 0.0)
            {}

        double ters_fc(double r) const {
            if(r < bigr-bigd) return 1.0;
            if(r > bigr+bigd) return 0.0;
            return 0.5*(1.0 - sin(M_PI/2*(r - bigr)/bigd));
        }
        double ters_fc_d(double r) const {
            if(r < bigr-bigd) return 0.0;
            if(r > bigr+bigd) return 0.0;
            return -(M_PI/4/bigd) * cos(M_PI/2*(r - bigr)/bigd);
        }
        double cut() const {
            return bigr + bigd;
        }
        double ters_fa(double r) const {
            if(r > bigr + bigd) return 0.0;
            return -D0*S/(S-1.0) * exp(-beta*sqrt(2/S)*(r-r0)) * ters_fc(r);
        }
        double ters_fa_d(double r) const {
            if(r > bigr + bigd) return 0.0;
            return D0*S/(S-1.0) * exp(-beta*sqrt(2/S)*(r-r0)) * (beta*sqrt(2/S) * ters_fc(r) - ters_fc_d(r));
        }
        double ters_gijk(double costheta) const {
            double ters_c = c * c;
            double ters_d = d * d;
            double hcth = h + costheta;
            return gamma * (1.0 + ters_c/ters_d - ters_c / (ters_d + hcth * hcth));
        }
        double ters_gijk_d(const double costheta) const {
            double ters_c = c * c;
            double ters_d = d * d;
            double hcth = h + costheta;
            double numerator = 2.0 * ters_c * hcth;
            double denominator = 1.0/(ters_d + hcth*hcth);
            return gamma * numerator * denominator * denominator;
        }
        double ters_bij(double zeta) const {
            double tmp = beta2 * zeta;
            if(tmp > c1) return 1.0 / sqrt(tmp);
            if(tmp > c2) return (1.0 - pow(tmp,-powern) / (2.0*powern)) / sqrt(tmp);
            if(tmp < c4) return 1.0;
            if(tmp < c3) return 1.0 - pow(tmp,powern)/(2.0*powern);
            return pow(1.0 + pow(tmp,powern), -1.0/(2.0*powern));
        }
        double ters_bij_d(double zeta) const {
            double tmp = beta2 * zeta;
            if(tmp > c1) return beta2 * -0.5*pow(tmp,-1.5);
            if(tmp > c2) return beta2 * (-0.5*pow(tmp,-1.5) * (1.0 - 0.5*(1.0 +  1.0/(2.0*powern)) * pow(tmp,-powern)));
            if(tmp < c4) return 0.0;
            if(tmp < c3) return -0.5*beta2 * pow(tmp,powern-1.0);
            double tmp_n = pow(tmp,powern);
            return -0.5 * pow(1.0+tmp_n, -1.0-(1.0/(2.0*powern)))*tmp_n / zeta;
        }
        void setCorrectedInitialValue(ScalarDOF& dof, double value) {
            // Due to the conversion between tersoff and abop format a parameter
            // that is exactly on the lower or upper boundary may be a tiny bit
            // outside of the boundary. This conversion errors are fixed here.
            // #include <cmath> to import the floating point version of abs. 
            if(dof.hasLowerBound() && std::abs(value - dof.lowerBound()) < 1e-9) {
                value = dof.lowerBound();
            }
            if(dof.hasUpperBound() && std::abs(value - dof.upperBound()) < 1e-9) {
                value = dof.upperBound();
            }
            dof.setInitialValue(value);
        }
        void setInitialParams(double r0,double D0, double beta, double S,
                double gamma, double c,double d,double h, double twomu,
                double beta2, double powern, double powerm, double bigd, double bigr) {
            // Store parameters.
            setCorrectedInitialValue(this->r0, r0);
            setCorrectedInitialValue(this->D0, D0);
            setCorrectedInitialValue(this->beta, beta);
            setCorrectedInitialValue(this->S, S);
            this->gamma.setInitialValue(gamma);
            this->c.setInitialValue(c);
            this->d.setInitialValue(d);
            this->h.setInitialValue(h);
            this->twomu.setInitialValue(twomu);
            this->beta2.setInitialValue(beta2);
            this->powern.setInitialValue(powern);
            this->bigd.setInitialValue(bigd);
            this->bigr.setInitialValue(bigr);
            this->powerm = powerm;

            // Compute dependent parameters.
            this->powermint = (int)powerm;
            this->c1 = pow(2.0*powern*1.0e-16,-1.0/powern);
            this->c2 = pow(2.0*powern*1.0e-8,-1.0/powern);
            this->c3 = 1.0/this->c2;
            this->c4 = 1.0/this->c1;
        }
    };

    const ABOParamSet& getParamSet(int itype, int jtype, int ktype) const {
        BOOST_ASSERT(itype >= 0 && itype < _numAtomTypes);
        BOOST_ASSERT(jtype >= 0 && jtype < _numAtomTypes);
        BOOST_ASSERT(ktype >= 0 && ktype < _numAtomTypes);
        return _params[itype * _numAtomTypes * _numAtomTypes + jtype * _numAtomTypes + ktype];
    }

    // The maximum cutoff radius of the potential.
    double _cutoff;

    // Parameter sets for I-J-K interactions.
    std::vector<ABOParamSet> _params;

    int _numAtomTypes;

    void parseTersoffFile(const FPString& filename);

    void linkDOF();

    FPString _exportPotentialFile;

    int _dofMode;
};

} // End of namespace