MEAMPotential.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 "functions/Functions.h"

namespace atomicrex {

class MEAMPotential : public Potential
{
public:

    MEAMPotential(const FPString& id, FitJob* job) : Potential(id, job, "MEAM") {}

    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 parse(XML::Element potentialElement) override;

    virtual size_t perAtomDataSize() const override { return sizeof(double); }

    virtual size_t perBondDataSize() const override { return sizeof(MEAMBondData); }

    FunctionBase* pairPotential(int speciesA, int speciesB) const {
        return _pairPotentialMap[speciesA * job()->numAtomTypes() + speciesB];
    }

    FunctionBase* electronDensity(int speciesA, int speciesB) const {
        return _electronDensityMap[speciesA * job()->numAtomTypes() + speciesB];
    }

    FunctionBase* embeddingEnergy(int species) const {
        return _embeddingEnergyMap[species];
    }

    FunctionBase* ffunction(int speciesA, int speciesB) const {
        return _ffunctionMap[speciesA * job()->numAtomTypes() + speciesB];
    }

    FunctionBase* gfunction(int speciesA, int speciesB, int speciesC) const {
        return _gfunctionMap[speciesA * job()->numAtomTypes() * job()->numAtomTypes() + speciesB * job()->numAtomTypes() + speciesC];
    }

    FunctionBase* findFunctionById(const FPString& id) const {
        for(const auto& f : _functions)
            if(f && f->id() == id) return f.get();
        return nullptr;
    }

    virtual void outputResults() override;

    void writeTables(const FPString& basename) const;

private:

    std::vector<std::shared_ptr<FunctionBase>> _functions;

    std::vector<FunctionBase*> _pairPotentialMap;

    std::vector<FunctionBase*> _electronDensityMap;

    std::vector<FunctionBase*> _embeddingEnergyMap;

    std::vector<FunctionBase*> _ffunctionMap;

    std::vector<FunctionBase*> _gfunctionMap;

    double _cutoff;

    FPString _exportFunctionsFile;

    // Helper data structure used to store temporary per-bond values during energy/force calculation.
    struct MEAMBondData {
        double f, fprime;
    };
};

} // End of namespace