CDIPotential.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 for pair-wise interactions that are concentration
 * dependent.
 *
 * Each interaction (between two species A and B) consists of
 * a potential function V(r) and an interpolation function h(x).
 * The concentration x is computed from the local density
 * of the matrix species S.
 ************************************************************/
class CDIPairInteraction : public FitObject
{
public:

    CDIPairInteraction(const FPString& id, FitJob* job, int atomTypeA, int atomTypeB, int matrix, const std::shared_ptr<FunctionBase>& V, const std::shared_ptr<FunctionBase>& h)
        : FitObject(id, job), _atomTypeA(atomTypeA), _atomTypeB(atomTypeB), _matrix(matrix), _V(V), _sigmaInverseArgument(new FunctionInverseArgument("CDIPHFunction", job, h)) {
        registerSubObject(V.get());
        registerSubObject(h.get());
    }

    int atomTypeA() const { return _atomTypeA; }
    int atomTypeB() const { return _atomTypeB; }

    int matrix() const { return _matrix; }

    bool interacting(int type_i, int type_j) const {
        if(type_i == _atomTypeA && type_j == _atomTypeB) return true;
        if(type_i == _atomTypeB && type_j == _atomTypeA) return true;
        return false;
    }

    const std::shared_ptr<FunctionBase>& V() const { return _V; }

    const std::shared_ptr<FunctionInverseArgument>& h() const { return _sigmaInverseArgument; }

private:

    int _atomTypeA, _atomTypeB, _matrix;

    std::shared_ptr<FunctionBase> _V;

    std::shared_ptr<FunctionInverseArgument> _sigmaInverseArgument;
};


/************************************************************
 * Composition-dependent interatomic potential (CDIP).
 ************************************************************/
class CDIPotential : public Potential
{
public:

    CDIPotential(const FPString& id, FitJob* job) : Potential(id, job, "CDIP"),
        _linearizedMode(false), _legacyMode(false) {}

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

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

    virtual void outputResults() override;

    void writePotential(const FPString& filename) const;

    void writeTables(const FPString& basename) const;

public:

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

    virtual XML::OElement generateXMLDefinition() override;

private:

    void parseSigmas(XML::Element parentElement);

    void parseInteractions(XML::Element parentElement);

    double _cutoff;

    int _linearizedMode;

    int _legacyMode;
  
    std::vector<std::shared_ptr<FunctionBase>> _sigmas;

    std::vector<std::shared_ptr<CDIPairInteraction>> _pairInteractions;

    FPString _exportPotentialFile;

    FPString _exportFunctionsFile;
};

} // End of namespace