/*
* Copyright 2009-2019 The VOTCA Development Team
* (http://www.votca.org)
*
* Licensed under the Apache License, Version 2.0 (the "License")
*
* You may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*/
#pragma once
#ifndef __VOTCA_XTP_POLARSITE_H
#define __VOTCA_XTP_POLARSITE_H
#include <votca/xtp/eigen.h>
#include <votca/xtp/staticsite.h>
namespace votca {
namespace xtp {
/**
\brief Class to represent Atom/Site in electrostatic+polarisation
The units are atomic units, e.g. Bohr, Hartree.
*/
class PolarSite : public StaticSite {
public:
// delete these two functions because we do not want to be able to read
// StaticSite::data but PolarSite::data
void WriteData(StaticSite::data& d) const = delete;
void ReadData(StaticSite::data& d) = delete;
PolarSite(Index id, std::string element, Eigen::Vector3d pos);
PolarSite(Index id, std::string element)
: PolarSite(id, element, Eigen::Vector3d::Zero()){};
~PolarSite() override = default;
void setPolarisation(const Eigen::Matrix3d& pol) override;
Eigen::Matrix3d getPolarisation() const { return _pinv.inverse(); }
const Eigen::Matrix3d& getPInv() const { return _pinv; }
// MULTIPOLES DEFINITION
Eigen::Vector3d getDipole() const override;
double getSqrtInvEigenDamp() const { return _eigendamp_invsqrt; }
void Rotate(const Eigen::Matrix3d& R,
const Eigen::Vector3d& ref_pos) override {
StaticSite::Rotate(R, ref_pos);
_pinv = R.transpose() * _pinv * R;
}
const Eigen::Vector3d& V() const { return _V; }
Eigen::Vector3d& V() { return _V; }
const Eigen::Vector3d& V_noE() const { return _V_noE; }
Eigen::Vector3d& V_noE() { return _V_noE; }
void Reset() {
_V.setZero();
_V_noE.setZero();
}
double deltaQ_V_ext() const { return _induced_dipole.dot(_V); }
double InternalEnergy() const {
return 0.5 * _induced_dipole.transpose() * _pinv * _induced_dipole;
}
const Eigen::Vector3d& Induced_Dipole() const { return _induced_dipole; }
void setInduced_Dipole(const Eigen::Vector3d& induced_dipole) {
_induced_dipole = induced_dipole;
}
struct data {
Index id;
char* element;
double posX;
double posY;
double posZ;
Index rank;
double Q00;
double Q11c;
double Q11s;
double Q10;
double Q20;
double Q21c;
double Q21s;
double Q22c;
double Q22s;
double Vx;
double Vy;
double Vz;
double Vx_noE;
double Vy_noE;
double Vz_noE;
double pxx;
double pxy;
double pxz;
double pyy;
double pyz;
double pzz;
double d_x_ind;
double d_y_ind;
double d_z_ind;
};
// do not move up has to be below data definition
PolarSite(const data& d);
double DipoleChange() const;
void SetupCptTable(CptTable& table) const override;
void WriteData(data& d) const;
void ReadData(const data& d);
std::string identify() const override { return "polarsite"; }
friend std::ostream& operator<<(std::ostream& out, const PolarSite& site) {
out << site.getId() << " " << site.getElement() << " " << site.getRank();
out << " " << site.getPos().transpose() << " "
<< site.Induced_Dipole().transpose() << "\n";
return out;
}
private:
std::string writePolarisation() const override;
// PolarSite has two external fields,
// the first is used for interaction with regions, which are further out, i.e.
// the interaction energy with it is included in the polar region energy
Eigen::Vector3d _V = Eigen::Vector3d::Zero();
// the second is used for interaction with regions, which are further inside,
// i.e. the interaction energy with it is included in the other region's
// energy
Eigen::Vector3d _V_noE = Eigen::Vector3d::Zero();
Eigen::Vector3d _induced_dipole = Eigen::Vector3d::Zero();
Eigen::Matrix3d _pinv = Eigen::Matrix3d::Zero();
double _eigendamp_invsqrt = 0.0;
};
} // namespace xtp
} // namespace votca
#endif