Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/md/lammps/ref_atom_lammps.hh
1d1b2315bf45master

Last Change
Raw File

ref_atom_lammps.hh
View Options

/**
* @file ref_atom_lammps.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief Reference over LAMMPS atom structure
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#ifndef __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__
#define __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__
/* -------------------------------------------------------------------------- */
#include "epot_hook.hh"
#include "ref_atom.hh"
#include "stress_hook.hh"
#include <Eigen/Dense>
#include <cmath>
/* -------------------------------------------------------------------------- */
#include "atom.h"
/* -------------------------------------------------------------------------- */
extern ::LAMMPS_NS::LAMMPS *lammps_main_object;
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <UInt Dim> class ReferenceManagerLammps;
template <UInt Dim> class ContainerLammps;
/* -------------------------------------------------------------------------- */
class ComparatorRefLammps;
template <UInt Dim> class DomainLammps;
/**
* Class RefLammps
*
*/
template <UInt _Dim> class RefLammps : public RefAtom<_Dim, RefLammps<_Dim>> {
public:
static constexpr UInt Dim = _Dim;
using RefComparator = ComparatorRefLammps;
using Domain = DomainLammps<_Dim>;
using fields =
field_list<_position0, _position, _displacement, _velocity, _force,
_stress, _mass, _angular_velocity, _torque, _radius>;
RefLammps(UInt id)
: index_atom(id), _weight(1), stress_hook(NULL), epot_hook(NULL){};
RefLammps() { RefLammps(0); }
virtual ~RefLammps(){};
public:
inline VectorView<_Dim> position0();
inline VectorView<_Dim> position();
inline AccessorAtomDof<RefLammps<_Dim>, _displacement> displacement();
inline VectorView<_Dim> velocity();
inline VectorView<_Dim> force();
inline VectorView<_Dim> angular_velocity();
inline VectorView<_Dim> angular_acceleration();
inline VectorView<_Dim> torque();
inline Real &radius();
inline MatrixView<_Dim> stress();
inline MatrixView<_Dim> kinetic_stress();
inline Real electronic_density();
inline Real getEPot();
inline UInt typeID() { return lammps_main_object->atom->type[index_atom]; };
inline UInt id();
inline Real &mass();
inline void setMass(Real mass);
inline Real &charge();
inline Real &alpha();
inline UInt tag();
inline int getIndex() const { return index_atom; };
inline bool operator==(const RefLammps &ref) const {
return (ref.index_atom == index_atom);
};
inline void setIndex(UInt n) { index_atom = n; };
friend class ComparatorRefLammps;
friend class ContainerLammps<Dim>;
friend class ReferenceManagerLammps<Dim>;
private:
int index_atom;
Real _weight;
StressHookLammps *stress_hook;
EpotHookLammps *epot_hook;
};
/* -------------------------------------------------------------------------- */
template <UInt _Dim> inline Real RefLammps<_Dim>::electronic_density() {
#ifdef PATCHED_LAMMPS
CHECK_INDEX;
LM_ASSERT(dynamic_cast<LAMMPS_NS::PairEAM *>(lammps_main_object->force->pair),
"the potential has to be EAM to get electronic density");
if (current_step == 0)
return 0;
LAMMPS_NS::PairEAM *ptr =
dynamic_cast<LAMMPS_NS::PairEAM *>(lammps_main_object->force->pair);
return ptr->rho[index_atom];
#else
LM_FATAL("to use this option you have to patch lammps to let 'rho'"
<< " array within EAM strucuture to be public");
#endif
return 0;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::position0() {
return VectorView<Dim>(lammps_main_object->atom->x0[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::angular_velocity() {
return VectorView<Dim>(lammps_main_object->atom->omega[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::torque() {
return VectorView<Dim>(lammps_main_object->atom->torque[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real &RefLammps<Dim>::radius() {
return lammps_main_object->atom->radius[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline AccessorAtomDof<RefLammps<Dim>, _displacement>
RefLammps<Dim>::displacement() {
return AccessorAtomDof<RefLammps<Dim>, _displacement>(*this);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::position() {
return VectorView<Dim>(lammps_main_object->atom->x[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline UInt RefLammps<Dim>::id() {
return lammps_main_object->atom->type[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::velocity() {
return VectorView<Dim>(lammps_main_object->atom->v[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline VectorView<Dim> RefLammps<Dim>::force() {
return VectorView<Dim>(lammps_main_object->atom->f[index_atom]);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline UInt RefLammps<Dim>::tag() {
return lammps_main_object->atom->tag[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline MatrixView<Dim> RefLammps<Dim>::stress() {
LM_TOIMPLEMENT;
// LM_ASSERT(stress_hook,
// "stress hook not initialized : stress cannot be computed");
// if (!stress_hook->Computed())
// stress_hook->ComputeStress();
// stress_hook->stress(index_atom, i);
// return stress_hook->stress(index_atom, i);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline MatrixView<Dim> RefLammps<Dim>::kinetic_stress() {
LM_TOIMPLEMENT;
// LM_ASSERT(stressKin_hook, "stress hook(Kin) not initialized :"
// << " stress cannot be computed");
// if (!stressKin_hook->Computed())
// stressKin_hook->ComputeStress();
// stressKin_hook->stress(index_atom, i);
// return stressKin_hook->stress(index_atom, i);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline Real RefLammps<Dim>::getEPot() {
if (!epot_hook) {
LM_ASSERT(epot_hook, "epot hook not initialized : epot cannot be computed");
}
if (!epot_hook->isComputed())
epot_hook->computeEpot();
return epot_hook->epot(index_atom);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline Real &RefLammps<Dim>::charge() {
LM_TOIMPLEMENT;
// return lammps_main_object->atom
// ->dipole[lammps_main_object->atom->type[index_atom]];
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline Real &RefLammps<Dim>::alpha() { return _weight; }
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline Real &RefLammps<Dim>::mass() {
if (lammps_main_object->atom->mass)
return lammps_main_object->atom
->mass[lammps_main_object->atom->type[index_atom]];
return lammps_main_object->atom->rmass[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> inline void RefLammps<Dim>::setMass(Real mass) {
lammps_main_object->atom->mass[lammps_main_object->atom->type[index_atom]] =
mass;
}
/* -------------------------------------------------------------------------- */
class ComparatorRefLammps {
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <UInt Dim>
bool operator()(RefLammps<Dim> &r1, RefLammps<Dim> &r2) const {
return r1.index_atom < r2.index_atom;
}
template <UInt Dim>
bool operator()(const RefLammps<Dim> &r1, const RefLammps<Dim> &r2) const {
return r1.index_atom < r2.index_atom;
}
};
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline std::ostream &operator<<(std::ostream &os, RefLammps<Dim> &ref) {
UInt n = ref.getIndex();
Eigen::IOFormat fmt(Eigen::FullPrecision, 0, ", ", "", "", "", "[", "]");
os << "Lammps reference : " << n << " P=" << ref.position().format(fmt)
<< " P0=" << ref.position0().format(fmt)
<< " V=" << ref.velocity().format(fmt) << " F=" << ref.force().format(fmt);
return os;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline std::ostream &operator<<(std::ostream &os,
const RefLammps<Dim> &ref [[gnu::unused]]) {
LM_TOIMPLEMENT;
return os;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__ */
c4science · Help