Site.h

#ifndef SITE_H_
#define SITE_H_
 
#include "utils/mardyn_assert.h"
 
#include "utils/Logger.h"
#include "utils/xmlfileUnits.h"
#include <array>
#include <cmath>
#include <cstdint>
 
using Log::global_log;
 
class Site {
public:
    double rx() const { return _r[0]; }  
    double ry() const { return _r[1]; }  
    double rz() const { return _r[2]; }  
    std::array<double, 3> r() const { return _r; }  
    double m() const { return _m; }         
    void setR(int d, double r) {
        mardyn_assert(d < 3);
        _r[d] = r;
    }
 
    void setM(double m) { _m = m; }
 
    virtual void readXML(XMLfileUnits& xmlconfig) {
        xmlconfig.getNodeValueReduced("coords/x", _r[0]);
        xmlconfig.getNodeValueReduced("coords/y", _r[1]);
        xmlconfig.getNodeValueReduced("coords/z", _r[2]);
        xmlconfig.getNodeValueReduced("mass", _m);
    }
 
    virtual ~~Site() {}
 
    virtual void write(std::ostream& ostrm) const {
        ostrm << _r[0] << " " << _r[1] << " " << _r[2];
    }
 
protected:
 
    Site(double x = 0., double y = 0., double z = 0., double m = 0.) : _r{{x, y, z}}, _m(m)  {}
 
    std::array<double, 3> _r;  
    double _m;  
};
 
 
class LJcenter : public Site {
public:
    LJcenter(): Site(0., 0., 0., 0.), _epsilon(0.), _sigma(0.), _uLJshift6(0.) {}
    LJcenter(double x, double y, double z, double m, double epsilon, double sigma, double shift)
        : Site(x, y, z, m), _epsilon(epsilon), _sigma(sigma), _uLJshift6(shift) {}
 
    void readXML(XMLfileUnits& xmlconfig) {
        Site::readXML(xmlconfig);
        xmlconfig.getNodeValueReduced("epsilon", _epsilon);
        xmlconfig.getNodeValueReduced("sigma", _sigma);
        xmlconfig.getNodeValueReduced("shifted", _uLJshift6);
    }
    
    void write(std::ostream& ostrm) const {
        Site::write(ostrm);
        ostrm << "\t" << eps() << " " << sigma() << " " << shift6();
    }
    /* @todo TODO: rename to epsilon */
    double eps() const { return _epsilon; }  
    double sigma() const { return _sigma; }  
    double shift6() const { return _uLJshift6; }  
    void setEps(double epsilon) { _epsilon = epsilon; }
    void setSigma(double sigma) { _sigma = sigma; }
    void setULJShift6(double uLJshift6) { _uLJshift6 = uLJshift6; }
 
 
private:
    double _epsilon;  
    double _sigma;  
    double _uLJshift6; 
};
 
class Charge : public Site {
public:
    Charge() : Site(0., 0., 0., 0.), _q(0.) {}
    Charge(double x, double y, double z, double m, double q) : Site(x, y, z, m), _q(q) {}
 
 
    void readXML(XMLfileUnits& xmlconfig) {
        Site::readXML(xmlconfig);
        xmlconfig.getNodeValueReduced("charge", _q);
    }
 
    void write(std::ostream& ostrm) const {
        ostrm << _r[0] << " " << _r[1] << " " << _r[2] << "\t" << _m << " " << _q;
    }
 
    double q() const { return _q; }  
    void setQ(double q) { _q = q; }
 
private:
    double _q;  
};
 
 
class OrientedSite : public Site {
public:
    double ex() const { return _e[0]; }  
    double ey() const { return _e[1]; }  
    double ez() const { return _e[2]; }  
    std::array<double, 3> e() const { return _e; }  
    double abs() const { return _abs; }  
    void write(std::ostream& ostrm) const {
        Site::write(ostrm);
        ostrm << "\t" << ex() << " " << ey() << " " << ez() << "\t" << abs();
    }
 
    void setOrientationVectorByPolarAngles(const double& theta_deg, const double& phi_deg) {
        const double fac = M_PI / 180.;  // translate: degrees --> rad
        double theta_rad = theta_deg * fac;
        double phi_rad   = phi_deg   * fac;
        _e[0] = sin(theta_rad) * cos(phi_rad);
        _e[1] = sin(theta_rad) * sin(phi_rad);
        _e[2] = cos(theta_rad);
    }
    void setAbs(double abs) {_abs = abs; }
 
    void setE(int d, double e) {
        mardyn_assert(d < 3);
        _e[d] = e;
    }
    void normalize_e() {
        _abs = sqrt(ex()*ex() + ey()*ey() + ez()*ez());
        for(int d = 0; d < 3; ++d) { _e[d] /= _abs; }
    }
 
protected:
    OrientedSite(double x = 0., double y = 0., double z = 0., double m = 0., double ex = 0., double ey = 0., double ez = 0., double abs = 0.)
        : Site(x, y, z, m), _e{{ex, ey, ez}}, _abs(abs) {}
 
    std::array<double, 3> _e;  
    double _abs;  
};
 
 
class Dipole : public OrientedSite {
public:
    Dipole() : OrientedSite(0., 0., 0., 0., 0., 0., 0., 0.) {}
    Dipole(double x, double y, double z, double eMyx, double eMyy, double eMyz, double absMy)
        : OrientedSite(x, y, z, 0., eMyx, eMyy, eMyz, absMy) {}
 
    void readXML(XMLfileUnits& xmlconfig) {
        Site::readXML(xmlconfig);
        bool bAngleInput = true;
        double theta, phi;
        bAngleInput = bAngleInput && xmlconfig.getNodeValueReduced("dipolemoment/theta", theta);
        bAngleInput = bAngleInput && xmlconfig.getNodeValueReduced("dipolemoment/phi",   phi);
        if(true == bAngleInput) {
            this->setOrientationVectorByPolarAngles(theta, phi);
        } else {
            xmlconfig.getNodeValueReduced("dipolemoment/x", _e[0]);
            xmlconfig.getNodeValueReduced("dipolemoment/y", _e[1]);
            xmlconfig.getNodeValueReduced("dipolemoment/z", _e[2]);
            normalize_e();
        }
        xmlconfig.getNodeValueReduced("dipolemoment/abs", _abs);
    }
 
    double absMy() const { return abs(); }  
    void setAbyMy(double my) { setAbs(my); }
};
 
class Quadrupole : public OrientedSite {
public:
    Quadrupole() : OrientedSite(0., 0., 0., 0., 0., 0., 0., 0.) {}
    Quadrupole(double x, double y, double z, double eQx, double eQy, double eQz, double absQ)
        : OrientedSite(x, y, z, 0., eQx, eQy, eQz, absQ) {}
 
    void readXML(XMLfileUnits& xmlconfig) {
        Site::readXML(xmlconfig);
        bool bAngleInput = true;
        double theta, phi;
        bAngleInput = bAngleInput && xmlconfig.getNodeValueReduced("quadrupolemoment/theta", theta);
        bAngleInput = bAngleInput && xmlconfig.getNodeValueReduced("quadrupolemoment/phi",   phi);
        if(true == bAngleInput) {
            setOrientationVectorByPolarAngles(theta, phi);
        } else {
            xmlconfig.getNodeValueReduced("quadrupolemoment/x", _e[0]);
            xmlconfig.getNodeValueReduced("quadrupolemoment/y", _e[1]);
            xmlconfig.getNodeValueReduced("quadrupolemoment/z", _e[2]);
            normalize_e();
        }
        xmlconfig.getNodeValueReduced("quadrupolemoment/abs", _abs);
    }
 
    double absQ() const { return abs(); }  
    void setAbsQ(double q) { setAbs(q); }
};
 
#endif  /* SITE_H_ */