AtomicStructure.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 "../Atomicrex.h"
#include "NeighborList.h"
#include "../job/FitObject.h"
#include "../properties/FitProperty.h"
#include "../properties/AtomVectorProperty.h"
#include "../dof/AtomCoordinatesDOF.h"

namespace atomicrex {

class FitJob;                   // Declared in FitJob.h

class AtomicStructure : public FitObject
{
public:

    enum DirtyFlags {
        NEIGHBOR_LISTS = (1<<0),
        ATOM_POSITIONS = (1<<1),
        STRUCTURE_DOF = (1<<4),
    };

    enum OutputFormat {
        POSCAR,
        LAMMPS
    };


public:

    /*************************************** Constructor / destructor ***********************************/

    AtomicStructure(const FPString& id, FitJob* job);

    /*************************************** Simulation cell geometry ************************************/

    void setupSimulationCell(const Matrix3& cellVectors, const Point3& cellOrigin = Point3::Origin(), const std::array<bool,3>& pbc = std::array<bool,3>{{true,true,true}});

    void deformSimulationCell(const Matrix3& deformation);

    const Matrix3& simulationCell() const { return _simulationCell; }

    const Point3& simulationCellOrigin() const { return _simulationCellOrigin; }

    const Matrix3& reciprocalSimulationCell() const { return _reciprocalSimulationCell; }

    bool hasPBC(size_t dimension) const { return _pbc[dimension]; }

    const std::array<bool,3>& pbc() const { return _pbc; }

    /*************************************** Coordinate utility functions ************************************/

    Vector3I periodicImage(const Point3& p) const;

    Point3 wrapPoint(Point3 p) const;

    Point3 wrapReducedPoint(Point3 p) const;

    Point3 reducedToAbsolute(const Point3& reducedPoint) const { return simulationCellOrigin() + (simulationCell() * (reducedPoint - Point3::Origin())); }

    Point3 absoluteToReduced(const Point3& worldPoint) const { return Point3::Origin() + (reciprocalSimulationCell() * (worldPoint - simulationCellOrigin())); }

    Vector3 reducedToAbsolute(const Vector3& reducedVec) const { return simulationCell() * reducedVec; }

    Vector3 absoluteToReduced(const Vector3& worldVec) const { return reciprocalSimulationCell() * worldVec; }

    /************************************************ Atoms ***********************************************/

    void setAtomCount(int numLocalAtoms);

    int numLocalAtoms() const { return _numLocalAtoms; }

    int numGhostAtoms() const { return _numGhostAtoms; }

    int numAtoms() const { return _numAtoms; }

    const std::vector<Point3>& atomPositions() const { return _atomPositions; }

    std::vector<Point3>& atomPositions() { return _atomPositions; }

    void setAtomPositions(const std::vector<Point3>& newPositions) {
        BOOST_ASSERT(newPositions.size() == numLocalAtoms());
        std::copy(newPositions.begin(), newPositions.end(), _atomPositions.begin());
        setDirty(ATOM_POSITIONS);
    }

    const std::vector<Point3>& atomDisplacements() const { return _atomDisplacements; }

    std::vector<Point3>& atomDisplacements() { return _atomDisplacements; }

    const std::vector<Point3>& atomInitialPositions() const { return _atomInitialPositions; }

    int atomType(int i) const { return _atomTypes[i]; }

    const std::vector<int>& atomTypes() const { return _atomTypes; }

    std::vector<int>& atomTypes() { return _atomTypes; }

    const std::vector<int>& atomTags() const { return _atomTags; }

    std::vector<int>& atomTags() { return _atomTags; }

    const std::vector<Vector3>& atomForces() const { return _atomForcesProperty.values(); }

    std::vector<Vector3>& atomForces() { return _atomForcesProperty.values(); }

    AtomVectorProperty& atomForcesProperty() { return _atomForcesProperty; }

    const std::vector< std::pair<int,int> >& forwardMapping() const { return _forwardMapping; }

    const std::vector<int>& reverseMapping() const { return _reverseMapping; }

    template<typename T>
    T* perAtomData() { return _perAtomPotentialData.empty() ? nullptr : reinterpret_cast<T*>(&_perAtomPotentialData.front()); }

    // Returns the neighbor list built for this structure and the given potential.
    NeighborList& neighborList(const Potential* potential) {
        auto nl = _neighborLists.find(potential);
        BOOST_ASSERT(nl != _neighborLists.end());
        return nl->second;
    }

    /************************************ Updating and modifying the structure ***************************************/

    void setDirty(DirtyFlags flags) { _dirtyFlags |= flags; }

    bool isDirty(DirtyFlags parts) const { return (_dirtyFlags & parts) != 0; }

    void clearDirty(DirtyFlags parts) { _dirtyFlags &= ~parts; }

    virtual void dofValueChanged(DegreeOfFreedom& dof) override {
        FitObject::dofValueChanged(dof);
        setDirty(STRUCTURE_DOF);
    }

    virtual void updateStructure();

    /******************************************** Property calculation ***********************************************/

    virtual bool computeProperties(bool isFitting) override;

    double computeEnergy(bool computeForces, bool isFitting, bool suppressRelaxation = false);

    bool relax(bool isFitting);

    /******************************************** Input/output ***********************************************/

    void writeToDumpFile(const FPString& filename, bool includeGhostAtoms = false) const;

    void writeToPoscarFile(const FPString& filename, bool includeGhostAtoms = false) const;

    void writeToFile();

    /************************************* Access to calculated quantities *****************************************/

    double totalEnergy() const { return _totalEnergy; }

    const std::array<double,6>& virial() const { return _virial; }

    std::array<double,6>& virial() { return _virial; }

    double pressureTensor(int voigtIndex) const { BOOST_ASSERT(voigtIndex >= 0 && voigtIndex < 6); return _pressureTensorProperty[voigtIndex]; }

    double pressure() const { return _pressureProperty; }

public:

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

private:

    void computeElasticConstants(bool isFitting);

    int _dirtyFlags;

    Matrix3 _simulationCell;
    Point3 _simulationCellOrigin;
    Matrix3 _reciprocalSimulationCell;
    std::array<bool,3> _pbc;

    int _numLocalAtoms;
    int _numGhostAtoms;
    int _numAtoms;
    double _skinThickness;
    std::vector<Point3> _atomPositions;
    std::vector<Point3> _lastNeighborUpdatePositions;
    Matrix3 _lastSimulationCell;
    Vector3I _ghostAtomImages;
    std::vector<Point3> _atomDisplacements;
    std::vector<Point3> _atomInitialPositions;
    std::vector<int> _atomTypes;
    std::vector<int> _atomTags;
    std::vector<std::pair<int,int>> _forwardMapping;
    std::vector<int> _reverseMapping;
    std::vector<unsigned char> _perAtomPotentialData;
    std::vector<Vector3> _ghostAtomPBCImages;

    std::map<const Potential*, NeighborList> _neighborLists;

    double _totalEnergy;

    int _lastUpdate = 0;

    FPString _outputFile = "";
    OutputFormat _outputFormat;

    std::array<double,6> _virial;

protected:

    ScalarFitProperty _totalEnergyProperty;
    ScalarFitProperty _atomicEnergyProperty;

    ScalarFitProperty _totalVolumeProperty;
    ScalarFitProperty _atomicVolumeProperty;

    ScalarFitProperty _pressureTensorProperty[6];
    ScalarFitProperty _pressureProperty;

    ScalarFitProperty _bulkModulusProperty;
    ScalarFitProperty _elasticConstantProperties[21];

    AtomVectorProperty _atomForcesProperty;

    AtomCoordinatesDOF _atomCoordinatesDOF;

};

} // End of namespace