Page MenuHomec4science
Files F65193730

lattice.cc
No OneTemporary
Actions

File Metadata

Created
Sat, Jun 1, 17:39

lattice.cc
View Options

/**
* @file lattice.cc
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Fri Apr 13 15:55:56 2012
*
* @brief
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#include "lattice.hh"
#include <gsl/gsl_blas.h>
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include <string>
#include <limits>
#include <sstream>
#include <cmath>
template<Uint DIM>
Real determinant(gsl_matrix * rotation) throw(std::string){
Real determinant = 0;
Uint num_terms;
if (DIM == 2) {
num_terms = 1;
} else if (DIM == 3) {
num_terms = 3;
} else {
std::stringstream error_stream ;
error_stream << "determinant for matrices other than 2x2 or 3x3 are not "
<< "implemented" << std::endl;
std::cout << error_stream.str();
throw(error_stream.str());
}
for (Uint i = 0; i < num_terms ; ++i) {
Real increment = 1;
Real decrement = 1;
for (Uint j = 0; j < DIM ; ++j) {
increment *= gsl_matrix_get(rotation, j, (j+i)%DIM);
decrement *= gsl_matrix_get(rotation, DIM-1-j, (j+i)%DIM);
}
determinant += increment;
determinant -= decrement;
}
return determinant;
}
template<Uint DIM>
Lattice<DIM>::Lattice(Real * _constants, Real * miller) throw(std::string)
:atoms("atoms")
{
this->constants[DIM] = 1;
for (Uint i = 0; i < DIM; ++i) {
this->constants[i] = _constants[i];
this->constants[DIM] *= _constants[i];
}
this->constants[DIM] = pow(this->constants[DIM], 1/static_cast<Real>(DIM));
this->millers = gsl_matrix_calloc(DIM, DIM);
this->rotation_matrix = gsl_matrix_calloc(DIM, DIM);
if (miller == NULL) {
gsl_vector_view diag = gsl_matrix_diagonal(this->millers);
gsl_vector_set_all(&diag.vector, 1);
diag = gsl_matrix_diagonal(this->rotation_matrix);
gsl_vector_set_all(&diag.vector, 1);
} else {
for (Uint i = 0; i < DIM; ++i) {
for (Uint j = 0; j < DIM; ++j) {
gsl_matrix_set(this->millers, i, j, miller[i*DIM + j]);
}
gsl_vector_const_view miller_row =
gsl_matrix_const_row(this->millers, i);
Real norm = gsl_blas_dnrm2(&miller_row.vector);
gsl_vector_view rotation_column =
gsl_matrix_column(this->rotation_matrix, i);
gsl_vector_memcpy(&rotation_column.vector, &miller_row.vector);
gsl_vector_scale(&rotation_column.vector, 1./norm);
}
}
Real det = determinant<DIM>(this->rotation_matrix);
if (fabs(det - 1) > std::numeric_limits<Real>::epsilon()){
std::stringstream error_stream;
error_stream << "Error: the miller indices you specified yield an invalid "
<< "rotation matrix (The determinant is " << det << " instead "
<< "of 1)." << std::endl;
for (Uint i = 0; i < DIM ; ++i) {
error_stream << "[";
for (Uint j = 0; j < DIM-1 ; ++j) {
error_stream << gsl_matrix_get(this->rotation_matrix, i, j) << ", ";
}
error_stream << gsl_matrix_get(this->rotation_matrix, i, DIM-1) << "]"
<< std::endl;
}
error_stream <<" Make sure they for a normal basis and comply with "
<< "the right hand rule. Lattice details: "
<< std::endl;
this->printself(error_stream, indent_incr);
std::cout << error_stream.str();
throw error_stream.str();
}
// compute spatial constants
this->spatial_constants[DIM] = 1;
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->spatial_constants[dir] = 0;
for (Uint j = 0 ; j < DIM ; ++j) {
this->spatial_constants[dir] +=
square(this->constants[j] * gsl_matrix_get(this->millers, dir, j));
}
this->spatial_constants[dir] = sqrt(spatial_constants[dir]);
this->spatial_constants[DIM] *= this->spatial_constants[dir];
}
this->spatial_constants[DIM] = pow(this->spatial_constants[DIM],
1/static_cast<Real>(DIM));
}
template<Uint DIM>
Lattice<DIM>::Lattice(const Lattice<DIM> & other) {
for (Uint i = 0 ; i < DIM+1 ; ++i) {
this->constants[i] = other.constants[i];
this->spatial_constants[i] = other.spatial_constants[i];
}
this->millers = gsl_matrix_calloc(DIM, DIM);
gsl_matrix_memcpy (this->millers, other.millers);
this->rotation_matrix = gsl_matrix_calloc(DIM, DIM);
gsl_matrix_memcpy(this->rotation_matrix, other.rotation_matrix);
this->atoms = other.atoms;
}
template<Uint DIM>
Lattice<DIM>::~Lattice<DIM>()
{
gsl_matrix_free(this->millers);
gsl_matrix_free(this->rotation_matrix);
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Lattice<DIM>::fillAtoms(int * lattice_coords, PointContainer<DIM> & container)
{
gsl_vector * current_point = gsl_vector_alloc(DIM);
gsl_vector * rotated_point = gsl_vector_alloc(DIM);
for (Uint i = 0; i < this->atoms.getSize(); ++i) {
PointRef<DIM> current_atom = this->atoms.getPoint(i);
for (Uint j = 0; j < DIM; ++j) {
gsl_vector_set(current_point, j, current_atom.getComponent(j) +
lattice_coords[j] * this->constants[j]);
}
gsl_blas_dgemv (CblasNoTrans,
1.0, this->rotation_matrix, current_point,
0.0, rotated_point);
container.addPoint(rotated_point->data);
}
gsl_vector_free(current_point);
gsl_vector_free(rotated_point);
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Lattice<DIM>::fillAtoms(const std::vector<Real> sw_coords,
PointContainer<DIM> & container) {
if (sw_coords.size() != DIM) {
std::stringstream error;
error << "expected a vector of dimension " << DIM << " but got "
<< sw_coords.size() << ".";
std::runtime_error(error.str());
}
Real coords[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
coords[i] = sw_coords[i];
}
this->fillAtoms(coords, container);
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Lattice<DIM>::fillAtoms(Real * sw_coords, PointContainer<DIM> & container)
{
gsl_vector * current_point = gsl_vector_alloc(DIM);
gsl_vector * rotated_point = gsl_vector_alloc(DIM);
for (Uint i = 0; i < this->atoms.getSize(); ++i) {
PointRef<DIM> current_atom = this->atoms.getPoint(i);
for (Uint j = 0; j < DIM; ++j) {
gsl_vector_set(current_point, j, current_atom.getComponent(j));
}
gsl_blas_dgemv (CblasTrans,
1.0, this->rotation_matrix, current_point,
0.0, rotated_point);
for (Uint i = 0; i < DIM ; ++i) {
rotated_point->data[i] += sw_coords[i];
}
container.addPoint(rotated_point->data);
}
gsl_vector_free(current_point);
gsl_vector_free(rotated_point);
}
template <Uint DIM>
const Real & Lattice<DIM>::getMinConstant() const {
Real tmp = std::numeric_limits<Real>::max();
Uint index = 0;
for (Uint dir = 0 ; dir < DIM ; ++dir) {
if (tmp < this->constants[dir]) {
tmp = this->constants[dir];
index = dir;
}
}
return this->constants[index];
}
template <Uint DIM>
void Lattice<DIM>::printself(std::ostream & stream, int indent) const
{
indent_space(stream, indent);
stream << this->lattice_type << DIM << "D Lattice, a = (";
for (Uint i = 0; i < DIM-1; ++i) {
stream << this->constants[i] << ", ";
}
stream << this->constants[DIM-1] << "):" << std::endl;
this->atoms.printself(stream, indent+indent_incr);
indent_space(stream, indent+indent_incr);
stream << "Miller orientations :" << std::endl;
for (Uint i = 0; i < DIM; ++i) {
indent_space(stream, indent + 2*indent_incr);
stream << "[";
for (Uint j = 0; j < DIM-1; ++j) {
stream << gsl_matrix_get(this->millers, i, j) << ", ";
}
stream << gsl_matrix_get(this->millers, i, DIM-1) << "]" << std::endl;
}
stream << "Transposed Rotation matrix :" << std::endl;
for (Uint i = 0; i < DIM; ++i) {
indent_space(stream, indent + 2*indent_incr);
stream << "[";
for (Uint j = 0; j < DIM-1; ++j) {
stream << gsl_matrix_get(this->rotation_matrix, i, j) << ", ";
}
stream << gsl_matrix_get(this->rotation_matrix, i, DIM-1) << "]" << std::endl;
}
}
template <Uint DIM>
std::string Lattice<DIM>::lattice_type = "";
template class Lattice<2>;
template class Lattice<3>;

Event Timeline

c4science · Help