RDF.h

#ifndef RDF_H
#define RDF_H
 
#include <cmath>
#include <string>
#include <utility>
#include <vector>
 
#include "plugins/PluginBase.h"
#include "molecules/Molecule.h"
#include "utils/CommVar.h"
 
class Component;
class RDFCellProcessor;
 
class RDF : public PluginBase {
 
    friend class RDFTest;
 
public:
 
    RDF();
    virtual ~~RDF();
 
    void readXML(XMLfileUnits& xmlconfig);
 
    void afterForces(ParticleContainer* particleContainer, DomainDecompBase* domainDecomp, unsigned long simstep);
 
    void init(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, Domain *domain);
 
    void finish(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, Domain *domain);
 
    std::string getPluginName() {
        return std::string("RDF");
    }
    static PluginBase* createInstance() { return new RDF(); }
 
    void setOutputTimestep(unsigned int timestep) { _writeFrequency = timestep; }
 
    void setOutputPrefix(std::string prefix) { _outputPrefix = prefix; }
 
    void endStep(ParticleContainer *particleContainer, DomainDecompBase *domainDecomposition, Domain *domain,
                 unsigned long simStep);
 
    void tickRDF() {
        _numberOfRDFTimesteps++;
    }
 
    void accumulateNumberOfMolecules(std::vector<Component>& components);
    
    void observeARDFMolecule(double dd, double cosPhi, double cosPhiReverse, unsigned cid1, unsigned cid2) {
        
        if(dd > _maxDistanceSquare) { return; }
        size_t distanceBinID = floor( sqrt(dd) / binwidth() );
        size_t angularBinID = floor( (-cosPhi + 1.)/ angularbinwidth() );
        size_t angularBinIDReverse = floor((-cosPhiReverse + 1.)/ angularbinwidth() );
        size_t binID = distanceBinID * _angularBins + angularBinID;
        size_t binIDReverse = distanceBinID * _angularBins + angularBinIDReverse;
        #if defined _OPENMP
        #pragma omp atomic
        #endif
        _ARDFdistribution.local[cid1][cid2][binID]++;
        #if defined _OPENMP
        #pragma omp atomic
        #endif
        _ARDFdistribution.local[cid2][cid1][binIDReverse]++;
    }
 
    void observeRDF(Molecule const& mi, Molecule const& mj, double dd) {
        observeRDFMolecule(dd, mi.componentid(), mj.componentid());
 
        if(isEnabledSiteRDF()) {
            double drs[3];
            double dr2;
            unsigned si = mi.numSites();
            unsigned sj = mj.numSites();
            if(si+sj > 2) {
                for(unsigned m = 0; m < si; m++) {
                    // when interacting two molecules with the same component id, the interaction of site m with site n is calculated twice if m!=n.
                    // this duplication is corrected by correctionFactor in the writeToFile(...) function.
                    // fixing this here, by starting n at m will break the unit tests and some symmetry.
                    for(unsigned n = 0; n < sj; n++) {
                        const std::array<double,3> dii = mi.site_d_abs(m);
                        const std::array<double,3> djj = mj.site_d_abs(n);
                        SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
                        observeRDFSite(dr2, mi.componentid(), mj.componentid(), m, n);
                    }
                }
            }
        }
    }
 
    void observeRDFMolecule(double dd, unsigned i, unsigned j) {
        if(dd > _maxDistanceSquare) { return; }
        if(i > j) { std::swap(j, i); }
        size_t binId = floor( sqrt(dd) / binwidth() );
        #if defined _OPENMP
        #pragma omp atomic
        #endif
        _distribution.local[i][j-i][binId]++;
    }
 
    inline void observeRDFSite(double dd, unsigned i, unsigned j, unsigned m, unsigned n) {
        if(dd > _maxDistanceSquare) { return; }
        if (i > j) {
            std::swap(j, i);
            std::swap(m, n);
        }
 
        unsigned int binId = floor( sqrt(dd) / binwidth() );
        #if defined _OPENMP
        #pragma omp atomic
        #endif
        _siteDistribution.local[i][j-i][m][n][binId] ++;
        if((i == j) && (m != n)){
            #if defined _OPENMP
            #pragma omp atomic
            #endif
            _siteDistribution.local[i][j-i][n][m][binId] ++;
        }
    }
 
    bool isEnabledSiteRDF() const { return _doCollectSiteRDF; }
    bool doARDF() const { return _doARDF; }
 
    void reset();  
 
private:
 
    template<typename T>
    void resizeExactly(std::vector<T>& v, unsigned int numElements) const {
        v.reserve(numElements);
        v.resize(numElements);
    }
 
    void init();
    unsigned int numBins() const { return _bins; }
    unsigned int numARDFBins() const { return _ARDFBins; }
    double binwidth() const { return _intervalLength; }
    double angularbinwidth() const { return _angularIntervalLength; }
    void collectRDF(DomainDecompBase* domainDecomp);  
 
    void accumulateRDF();
 
    void writeToFile(const Domain* domain, const std::string& filename, unsigned int i, unsigned int j) const;
    void writeToFileARDF(const Domain* domain, const std::string& filename, unsigned int i, unsigned int j) const;
    double _intervalLength;
    
    double _angularIntervalLength;
 
    unsigned long _bins;
    
    unsigned long _angularBins {1};
    
    unsigned long _ARDFBins;
 
    unsigned int _numberOfComponents;
 
    std::vector<Component>* _components;
 
    int _samplingFrequency;
 
    int _numberOfRDFTimesteps;
 
    int _accumulatedNumberOfRDFTimesteps;
 
    double _maxDistanceSquare;
 
    std::vector<unsigned long> _globalCtr;
 
    std::vector<unsigned long> _globalAccumulatedCtr;
 
    CommVar<std::vector<std::vector<std::vector<unsigned long>>>> _distribution;
    CommVar<std::vector<std::vector<std::vector<unsigned long>>>> _ARDFdistribution;
    std::vector<std::vector<std::vector<unsigned long>>> _globalAccumulatedDistribution;
    std::vector<std::vector<std::vector<unsigned long>>> _globalAccumulatedARDFDistribution;
    
    bool _doCollectSiteRDF;
    bool _doARDF;
    // vector indices:
    // first: component i
    // second: component j, but shifted by -i, so that we only have to save it once for each component interaction (i<->j is the same as j<->i)
    // third: site m of first component
    // fourth: site n of second component
    // fifth: bin to sort the particles into to get the actual RDF.
    CommVar<std::vector<std::vector<std::vector<std::vector<std::vector<unsigned long>>>>>> _siteDistribution;
 
    std::vector<std::vector<std::vector<std::vector<std::vector<unsigned long>>>>> _globalAccumulatedSiteDistribution;
 
    unsigned int _writeFrequency;  
    std::string _outputPrefix;  
 
    bool _initialized;
    bool _readConfig;
 
    RDFCellProcessor * _cellProcessor;
};
 
#endif /* RDF_H */