Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/md/lammps/ref_atom_lammps.hh
b05d5f2aba78features/coupler_compon…

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 "ref_atom.hh"
#include "stress_hook.hh"
#include "epot_hook.hh"
#include "pair_eam.h"
#include <cmath>
/* -------------------------------------------------------------------------- */
#include "atom.h"
/* -------------------------------------------------------------------------- */
extern ::LAMMPS_NS::LAMMPS * lammps_main_object;
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename T > class ReferenceManager;
/* -------------------------------------------------------------------------- */
#ifndef LM_OPTIMIZED
#define CHECK_INDEX check_index()
#else
#define CHECK_INDEX
#endif
class ComparatorRefLammps;
template <UInt Dim> class DomainLammps;
/**
* Class RefLammps
*
*/
template <UInt _Dim>
class RefLammps : public RefAtom<_Dim> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
RefLammps(UInt id):
index_atom(id),_weight(1),
stress_hook(NULL),stressKin_hook(NULL),epot_hook(NULL){};
RefLammps(){
RefLammps(0);
}
virtual ~RefLammps(){};
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
typedef ComparatorRefLammps RefComparator;
typedef DomainLammps<_Dim> Domain;
static const UInt Dim = _Dim;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
inline Real & position0(UInt );
inline Real & position(UInt );
inline AccessorAtomDof<_Dim,_displacement> displacement(UInt );
inline Real & velocity(UInt );
inline Real & force(UInt );
inline Real & stress(UInt i);
inline Real & kinetic_stress(UInt i);
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(){return index_atom;};
inline bool operator == (RefLammps & ref){
return (ref.index_atom == index_atom);
};
inline void setIndex(UInt n){
index_atom = n;
};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
private:
#ifndef LM_OPTIMIZED
inline void check_index(){
LM_ASSERT(index_atom >= 0 &&
index_atom < lammps_main_object->atom->nlocal
+ lammps_main_object->atom->nghost,
"attempt to access to atom " << index_atom
<< " out of range 0 - "
<< lammps_main_object->atom->nlocal+lammps_main_object->atom->nghost-1);
}
#endif
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
public:
int index_atom;
Real _weight;
StressHookLammps * stress_hook;
StressHookLammps * stressKin_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 Real & RefLammps<_Dim>::position0(UInt i){
CHECK_INDEX;
return lammps_main_object->atom->x0[index_atom][i];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline AccessorAtomDof<_Dim,_displacement> RefLammps<_Dim>::displacement(UInt i){
CHECK_INDEX;
return AccessorAtomDof<_Dim,_displacement>(*this,i);
// Real disp = position(i)-position0(i);
// UInt pbc_flag = false;
// switch (i){
// case 0 : pbc_flag = lammps_main_object->domain->xperiodic;break;
// case 1 : pbc_flag = lammps_main_object->domain->yperiodic;break;
// case 2 : pbc_flag = lammps_main_object->domain->zperiodic;break;
// default: LM_FATAL("internal error");
// }
// if (pbc_flag)
// {
// Real size = lammps_main_object->domain->boxhi[i] -
// lammps_main_object->domain->boxlo[i];
// if (disp >= size/2) disp -= size;
// if (disp < -size/2) disp += size;
// LM_ASSERT(fabs(disp) <= size/2,"this should not happend");
// }
// return disp;
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline Real & RefLammps<_Dim>::position(UInt i){
CHECK_INDEX;
return lammps_main_object->atom->x[index_atom][i];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline UInt RefLammps<_Dim>::id(){
CHECK_INDEX;
return lammps_main_object->atom->type[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline Real & RefLammps<_Dim>::velocity(UInt i){
CHECK_INDEX;
return lammps_main_object->atom->v[index_atom][i];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline Real & RefLammps<_Dim>::force(UInt i){
CHECK_INDEX;
return lammps_main_object->atom->f[index_atom][i];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline UInt RefLammps<_Dim>::tag(){
return lammps_main_object->atom->type[index_atom];
}
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
inline Real & RefLammps<_Dim>::stress(UInt i){
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 Real & RefLammps<_Dim>::kinetic_stress(UInt i){
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(){
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(){
return lammps_main_object->atom->mass[lammps_main_object->atom->type[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();
os << "Lammps reference : " << n
<< " [P] " << ref.position(0) << " " << ref.position(1) << " " << ref.position(2)
<< " [P0] " << ref.position0(0) << " " << ref.position0(1) << " " << ref.position0(2)
<< " [V] " << ref.velocity(0) << " " << ref.velocity(1) << " " << ref.velocity(2)
<< " [F] " << ref.force(0) << " " << ref.force(1) << " " << ref.force(2);
return os;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__ */
c4science · Help