Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/math/dispersion_relation.hh
8ae23a91b711lammps_2020

Last Change
Raw File

dispersion_relation.hh
View Options

/**
* @file dispersion_relation.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief This function implements the dispersion relation of various crytals
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_DISPERSION_RELATION_HH__
#define __LIBMULTISCALE_DISPERSION_RELATION_HH__
/* -------------------------------------------------------------------------- */
#include "function_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/** Class DispersionRelation
*
* This class implements access to the dispersion relation
*
*/
class DispersionRelationInterface : public FunctionInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DispersionRelationInterface() {
this->elastic_stiffness = 1.0;
this->mass = 1.0;
this->r0 = 1.0;
};
DispersionRelationInterface(Real stiffness, const Quantity<Mass> &mass,
const Quantity<Length> &r0) {
this->elastic_stiffness = stiffness;
DUMP(" STIFFNESS " << this->elastic_stiffness, DBG_MESSAGE);
this->mass = mass;
DUMP(" ATOM MASS " << this->mass, DBG_MESSAGE);
this->r0 = r0;
DUMP(" R0 " << this->r0, DBG_MESSAGE);
};
virtual ~DispersionRelationInterface(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
// inline Real compute(Real * k, UInt branch);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
Real elastic_stiffness;
Quantity<Length> r0;
Quantity<Mass> mass;
};
/* -------------------------------------------------------------------------- */
class DispersionRelation : public DispersionRelationInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DispersionRelation(){};
DispersionRelation(Real stiffness, const Quantity<Mass> &mass,
const Quantity<Length> &r0)
: DispersionRelationInterface(stiffness, mass, r0){};
virtual ~DispersionRelation(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <LatticeType lattice_type> inline Real compute(Real *k, UInt branch);
inline std::vector<Real> compute(std::vector<Real> &wave_numbers,
const LatticeType &type, UInt branch);
template <LatticeType lattice_type>
inline std::vector<Real> compute(std::vector<Real> &wave_numbers,
UInt branch);
template <LatticeType lattice_type>
inline void computePolarization(Real *k, Real *pol, UInt branch);
inline std::vector<Real> computePolarization(std::vector<Real> &wave_numbers,
const LatticeType &type,
UInt branch);
template <LatticeType lattice_type>
inline std::vector<Real> computePolarization(std::vector<Real> &wave_numbers,
UInt branch);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
};
/* -------------------------------------------------------------------------- */
template <LatticeType lattice_type>
inline Real DispersionRelation::compute(Real *k, UInt branch) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <LatticeType lattice_type>
inline void DispersionRelation::computePolarization(Real *k, Real *pol,
UInt branch) {
constexpr UInt Dim = LatticeDim<lattice_type>::value;
pol[0] = 1.;
for (UInt i = 1; i < Dim; ++i)
pol[i] = 0.;
}
/* -------------------------------------------------------------------------- */
template <>
inline void DispersionRelation::computePolarization<hex_2d>(Real *k, Real *pol,
UInt branch) {
Real A = cos(k[0] * r0 / 2) * cos(k[1] * sqrt(3) * r0 / 2);
Real B = sin(k[0] * r0 / 2) * sin(sqrt(3) * k[1] * r0 / 2);
Real C = cos(k[0] * r0);
Real D = C - A;
Real D2 = D * D;
Real B2 = B * B;
Real sign = pow(-1., branch);
DUMP("A = " << A, DBG_DETAIL);
DUMP("B = " << B, DBG_DETAIL);
DUMP("C = " << C, DBG_DETAIL);
DUMP("D = " << D, DBG_DETAIL);
DUMP("sign = " << sign, DBG_DETAIL);
if (B != 0) {
Real Delta = D2 + 3. * B2;
DUMP("Delta = " << Delta, DBG_DETAIL);
pol[0] = A - C - sqrt(Delta);
pol[1] = sqrt(3.) * B;
Real n = sqrt(pol[0] * pol[0] + pol[1] * pol[1]);
DUMP("n = " << n, DBG_DETAIL);
pol[0] /= n;
pol[1] /= n;
} else {
if (sign < 0) {
pol[0] = 0;
pol[1] = 1;
} else {
pol[0] = 1;
pol[1] = 0;
}
}
DUMP("k = " << k[0] << " " << k[1] << " polarization = " << pol[0] << " "
<< pol[1],
DBG_MESSAGE);
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<mono_atom_1d>(Real *k, UInt branch) {
Real omega = 0.0;
omega = 2.0 * sqrt(elastic_stiffness * code_unit_system.e_dd2m_tt / mass) *
(std::abs(sin(k[0] / 2.0 * r0)));
DUMP("k = " << k[0] << " omega = " << omega, DBG_MESSAGE);
return omega;
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<hex_2d>(Real *k, UInt branch) {
Real A = cos(k[0] * r0 / 2) * cos(k[1] * sqrt(3) * r0 / 2);
Real B = sin(k[0] * r0 / 2) * sin(sqrt(3) * k[1] * r0 / 2);
Real C = cos(k[0] * r0);
Real D = C - A;
Real D2 = D * D;
Real B2 = B * B;
Real Delta = D2 + 3. * B2;
Real sign = pow(-1., branch);
Real omega = sqrt((elastic_stiffness * code_unit_system.e_dd2m_tt / mass) *
(3.0 - 2.0 * A - C + sign * sqrt(Delta)));
// Real omega1D = compute<mono_atom_1d>(k,branch);
DUMP("k = " << k[0] << " " << k[1] << " omega = " << omega, DBG_MESSAGE);
// DUMP("k = " << k[0] << " " << k[1] << " omega1D = " <<
// omega1D,DBG_MESSAGE);
return omega;
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<hex_2d_radial>(Real *k, UInt branch) {
Real temp0 = (cos(2 * k[0]) - cos(k[0]) * cos(k[1]));
Real temp1 = 3.0 * sin(k[0]) * sin(k[0]) * sin(k[1]) * sin(k[1]);
Real omega = sqrt((elastic_stiffness / mass) *
(3.0 - cos(2 * k[0]) - 2.0 * cos(k[0]) * cos(k[1]) -
sqrt(temp0 * temp0 + temp1)));
return omega;
}
/* -------------------------------------------------------------------------- */
inline std::vector<Real>
DispersionRelation::compute(std::vector<Real> &wave_numbers,
const LatticeType &type, UInt branch) {
switch (type) {
case mono_atom_1d:
return compute<mono_atom_1d>(wave_numbers, branch);
break;
case hex_2d:
return compute<hex_2d>(wave_numbers, branch);
break;
case hex_2d_radial:
return compute<hex_2d_radial>(wave_numbers, branch);
break;
case fcc_3d:
return compute<fcc_3d>(wave_numbers, branch);
break;
default:
LM_FATAL("internal error: unknown lattice");
}
}
/* -------------------------------------------------------------------------- */
inline std::vector<Real>
DispersionRelation::computePolarization(std::vector<Real> &wave_numbers,
const LatticeType &type, UInt branch) {
switch (type) {
case mono_atom_1d:
return computePolarization<mono_atom_1d>(wave_numbers, branch);
break;
case hex_2d:
return computePolarization<hex_2d>(wave_numbers, branch);
break;
case hex_2d_radial:
return computePolarization<hex_2d_radial>(wave_numbers, branch);
break;
case fcc_3d:
return computePolarization<fcc_3d>(wave_numbers, branch);
break;
default:
LM_FATAL("internal error: unknown lattice");
}
}
/* -------------------------------------------------------------------------- */
template <LatticeType type>
inline std::vector<Real>
DispersionRelation::compute(std::vector<Real> &wave_numbers, UInt branch) {
constexpr UInt Dim = LatticeDim<type>::value;
UInt total_nb_modes = wave_numbers.size() / Dim;
std::vector<Real> res;
for (UInt nb = 0; nb < total_nb_modes; ++nb) {
Real *k = &wave_numbers[nb * Dim];
res.push_back(this->compute<type>(k, branch));
}
return res;
}
/* -------------------------------------------------------------------------- */
template <LatticeType type>
inline std::vector<Real>
DispersionRelation::computePolarization(std::vector<Real> &wave_numbers,
UInt branch) {
constexpr UInt Dim = LatticeDim<type>::value;
Real pol[Dim];
UInt total_nb_modes = wave_numbers.size() / Dim;
std::vector<Real> res;
for (UInt nb = 0; nb < total_nb_modes; ++nb) {
Real *k = &wave_numbers[nb * Dim];
this->computePolarization<type>(k, pol, branch);
for (UInt i = 0; i < Dim; ++i)
res.push_back(pol[i]);
}
return res;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_DISPERSION_RELATION_HH__ */
c4science · Help