Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/md/1d/domain_md1d.cc
f7c9b4a7ad6amaster

Last Change
Raw File

domain_md1d.cc
View Options

/**
* @file domain_md1d.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Fri Jan 10 20:47:45 2014
*
* @brief This is the domain for 1d atomic chains
*
* @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.
*
*/
/* -------------------------------------------------------------------------- */
#include "domain_md1d.hh"
#include "lm_common.hh"
#include "md_builder.hh"
#include "reference_manager.hh"
#include "reference_manager_md1d.hh"
#include "spatial_filter_atomic.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
void DomainMD1D::correctPBCPositions() {
int nb_atoms = this->getNbDofs();
Real domainSize = this->getDomainSize();
for (int i = 0; i < nb_atoms; ++i) {
dofs[i].p -= domainSize * nearbyint(dofs[i].p / domainSize);
}
}
/* -------------------------------------------------------------------------- */
inline Real DomainMD1D::computeDistance(const Real &xi, const Real &xj) {
Real rij = xj - xi;
Real domainSize = this->getDomainSize();
if (this->periodic_flag)
rij -= domainSize * nearbyint(rij / domainSize);
return rij;
}
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getDomainSize() {
if (this->periodic_flag)
return this->getNbDofs() * r0;
return (this->getNbDofs() - 1) * r0;
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::init() {
if (!Communicator::getCommunicator().amIinGroup(this->getGroupID()))
return;
ZoneMDBuilder z;
Cube b = this->getGeom().getBoundingBox();
if (this->atom_file == "")
z.buildMD1D(b, this->r0, this->mass, this->shift, this->periodic_flag,
this->dofs);
else
z.readMD1D(this->atom_file, b, this->r0, this->mass, this->dofs);
DUMP("Size of domain " << this->getDomainSize(), DBG_MESSAGE);
DUMP("position offset " << -1.0 * this->getDomainSize() / 2 << " Center "
<< b.getCenter(0) << " Shift " << shift,
DBG_MESSAGE);
static_cast<ContainerMD1D &>(atoms).setDofs(dofs);
DUMP(" MD1D R0 " << r0 << " RCUT " << rcut << " Rcut/r0 "
<< (Real)(rcut / r0),
DBG_MESSAGE);
int nb_atoms = this->getNbDofs();
stress = new Real[nb_atoms];
stiffness_parameter = 36. * epsilon * pow(2., 2. / 3.) / sigma / sigma;
DUMP("Stiffness parameter " << stiffness_parameter << " epsilon " << epsilon
<< " sigma " << sigma,
DBG_MESSAGE);
this->computeNeighborLists();
this->atoms.setRelease(INITIAL_MODEL_RELEASE);
}
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getEcin() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getEpot() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
DomainMD1D::DomainMD1D(DomainID ID, CommGroup GID)
: DomainAtomic<ContainerMD1D, 1>(ID, GID), periodic_flag(false),
elastic_flag(false), shift(0.0) {}
/* -------------------------------------------------------------------------- */
void DomainMD1D::performStep1() {
DUMP("Position step avec dt = " << timeStep, DBG_INFO);
int nb_atoms = this->getNbDofs();
for (int i = 0; i < nb_atoms; ++i) {
/* mise ajour des positions */
dofs[i].v +=
code_unit_system.ft_m2v * .5 * timeStep * dofs[i].f / dofs[i].m;
dofs[i].p += timeStep * dofs[i].v;
LM_ASSERT(!isinf(dofs[i].p),
"one atomic position "
<< i << " is attributed an infinite position : " << dofs[i].v
<< " " << dofs[i].a << " " << timeStep);
DUMP(getID() << " dof " << i << " x0 = " << dofs[i].p0
<< " u = " << dofs[i].p - dofs[i].p0 << " v = " << dofs[i].v
<< " a = " << dofs[i].a,
DBG_DETAIL);
}
if (periodic_flag)
this->correctPBCPositions();
this->atoms.incRelease();
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::performStep3() {
DUMP("Velocity step avec dt = " << timeStep, DBG_INFO);
int nb_atoms = dofs.size();
for (int i = 0; i < nb_atoms; ++i) {
dofs[i].v +=
code_unit_system.ft_m2v * 0.5 * timeStep * dofs[i].f / dofs[i].m;
DUMP(getID() << " dof " << i << " v = " << dofs[i].v
<< " a = " << dofs[i].a,
DBG_DETAIL);
}
this->atoms.incRelease();
}
/* -------------------------------------------------------------------------- */
inline void DomainMD1D::computeNeighborLists() {
int nb_atoms = dofs.size();
neighbor_lists.resize(nb_atoms);
for (UInt i = 0; i < nb_atoms; ++i) {
std::vector<UInt> neighbors;
for (UInt j = i + 1; j < nb_atoms; ++j) {
Real rij = this->computeDistance(dofs[j].p, dofs[i].p);
LM_ASSERT(std::fabs(rij) / sigma > 1e-6,
"problem in neighbor list: " << i << " " << j);
if (rij <= rcut * 1.5)
neighbors.push_back(j);
}
neighbor_lists[i] = neighbors;
// for (auto j: neighbor_lists[i])
// DUMP("neighborlist " << i << ":" << j, DBG_MESSAGE);
}
}
/* -------------------------------------------------------------------------- */
inline void DomainMD1D::computeContribution(UInt i, UInt j) {
Real rij = computeDistance(dofs[i].p, dofs[j].p);
Real sign = rij < 0 ? -1. : 1.;
rij = std::fabs(rij);
LM_ASSERT(std::fabs(rij) / sigma > 1e-6, "problem in neighbor list");
Real contribution =
24.0 * epsilon * (2.0 * MathTools::ipow(sigma / rij, 12) / rij -
1.0 * MathTools::ipow(sigma / rij, 6) / rij);
if (arlequin_flag)
contribution *= 0.5 * (dofs[i].alpha + dofs[j].alpha);
DUMP("contribution between " << i << " and " << j << " = " << contribution
<< " [d=" << rij << "]",
DBG_ALL);
dofs[i].a -= contribution*sign;
dofs[j].a += contribution*sign;
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::performStep2() {
// if (elastic_flag) {
// quadraticForceComputation();
// return;
// }
LM_ASSERT(dofs[3].m != 0, "mass pb");
int nb_atoms = dofs.size();
for (UInt i = 0; i < nb_atoms; ++i) {
dofs[i].f = 0.;
dofs[i].a = 0.;
}
// loop over atoms to compute their force
for (UInt i = 0; i < nb_atoms; ++i) {
// loop over right neighbors
for (auto j : neighbor_lists[i]) {
Real rij = this->computeDistance(dofs[j].p, dofs[i].p);
if (rij <= rcut)
computeContribution(i, j);
}
}
for (UInt i = 0; i < nb_atoms; ++i) {
dofs[i].f = dofs[i].a * dofs[i].m;
}
for (UInt i = 0; i < nb_atoms; ++i) {
// if (arlequin_flag) {
// dofs[i].a /= dofs[i].m * dofs[i].alpha;
// dofs[i].f = dofs[i].a * dofs[i].m;
// } else {
// dofs[i].a /= dofs[i].m;
// dofs[i].f = dofs[i].a * dofs[i].m;
// }
LM_ASSERT(!isinf(dofs[i].a),
"one atomic "
<< i << " position is attributed an infinite acceleration : "
<< dofs[i].m);
}
this->atoms.incRelease();
}
/* -------------------------------------------------------------------------- */
// Added by Srinivasa Babu Ramisetti
void DomainMD1D::quadraticForceComputation() {
double contribution;
double rij;
double pos_right = 0, pos_left = 0;
int i;
int nb_atoms = dofs.size();
for (i = 0; i < nb_atoms; ++i) {
dofs[i].pe = 0.0;
double pos_i = dofs[i].p;
int left = (i - 1 + nb_atoms) % nb_atoms;
int right = (i + 1) % nb_atoms;
if (periodic_flag != 1) {
left = i - 1;
right = i + 1;
}
if (right < nb_atoms) {
pos_right = dofs[right].p;
rij = computeDistance(pos_i, pos_right);
contribution = (rij - r0) * stiffness_parameter;
dofs[i].f += contribution;
dofs[i].pe += (rij - r0) * contribution / 4.0;
}
if (left >= 0) {
pos_left = dofs[left].p;
rij = computeDistance(pos_i, pos_left);
contribution = (rij - r0) * stiffness_parameter;
dofs[i].f -= contribution;
dofs[i].pe += (rij - r0) * contribution / 4.0;
}
dofs[i].a = dofs[i].f / dofs[i].m;
}
}
/* -------------------------------------------------------------------------- */
// Added by Srinivasa Babu Ramisetti
// void DomainMD1D::getNeighborList(RefAtom & atom, Real distance,
// std::vector<RefAtom> & neigbhor_list)
void DomainMD1D::getNeighborList(RefAtom &atom, UInt count,
std::vector<RefAtom> &neigbhor_list) {
UInt index_count = 0;
neigbhor_list.clear();
RefMD1D r1d(dofs);
int current_id = atom.id();
int nb_atoms = dofs.size();
if (periodic_flag == 1) {
// DUMP(" Current atom id "<< atom.id() << " pos0 " << atom.position0(0) <<
// " nb_coefficients " << count << " nb_atoms " << nb_atoms, DBG_MESSAGE);
int id = 1;
int i_id1 = (current_id + id);
int i_id2 = (current_id - id);
// Real d = fabs(atom.position0(0) - dofs[i_id1].p0);
// if(current_id == (nb_atoms - 1))
// d = fabs(atom.position0(0) - dofs[i_id2].p0);
// while(d < distance){
while (index_count < count) {
i_id1 = (current_id + id);
if (i_id1 < 0 || i_id1 >= nb_atoms) {
i_id1 = (current_id + id) % nb_atoms;
// d = fabs(atom.position0(0) - dofs[current_id - id].p0);
}
r1d.setIndex(i_id1);
neigbhor_list.push_back(r1d);
i_id2 = (current_id - id);
if (i_id2 < 0 || i_id2 > nb_atoms) {
i_id2 = (current_id - id + nb_atoms) % nb_atoms;
// d = fabs(atom.position0(0) - dofs[current_id + id].p0);
}
r1d.setIndex(i_id2);
neigbhor_list.push_back(r1d);
id++;
index_count++;
// DUMP(" l_i "<< i_id2 << " r_i " << i_id1 << " vel[l] " << v2 << "
// vel[r] " << v1 , DBG_MESSAGE);
}
} else {
DUMP(" Current atom id " << atom.id() << " pos0 " << atom.position0()
<< " nb_atoms " << nb_atoms,
DBG_DETAIL);
int id = 1;
// Real d = fabs(atom.position0(0) - dofs[current_id + 1].p0);
// while(d < distance){
while (index_count < count) {
if ((current_id + id) < 0 || (current_id + id) > nb_atoms)
LM_FATAL("Neighbor id's are out of bound!");
r1d.setIndex(current_id + id);
neigbhor_list.push_back(r1d);
if ((current_id - id) < 0 || (current_id - id) > nb_atoms)
LM_FATAL("Neighbor id's are out of bound!");
r1d.setIndex(current_id - id);
neigbhor_list.push_back(r1d);
DUMP(" neigbhor id on left : " << current_id - id
<< " neigbhor id on right : "
<< current_id + id,
DBG_DETAIL);
id++;
// d = fabs(atom.position0(0) - dofs[current_id + id].p0);
// count = count + 2;
index_count++;
}
}
DUMP("Number of neighbor's found " << index_count * 2, DBG_DETAIL);
}
Real DomainMD1D::getElasticConstant() { return this->stiffness_parameter; }
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getR0() { return this->r0; }
/* -------------------------------------------------------------------------- */
int DomainMD1D::getNbDofs() { return this->dofs.size(); }
/* -------------------------------------------------------------------------- */
UInt DomainMD1D::getDim() { return 1; }
/* -------------------------------------------------------------------------- */
void DomainMD1D::getDomainDimensions(Real *xmin, Real *xmax) {
xmin[0] = lm_real_max;
xmax[0] = -1. * lm_real_max;
int nb_atoms = dofs.size();
for (int i = 0; i < nb_atoms; ++i) {
if (xmin[0] > dofs[i].p0)
xmin[0] = dofs[i].p0;
if (xmax[0] < dofs[i].p0)
xmax[0] = dofs[i].p0;
}
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::getSubDomainDimensions(Real *xmin, Real *xmax) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::setCurrentStep(UInt ts) { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
UInt DomainMD1D::getCurrentStep() {
LM_TOIMPLEMENT;
return 0;
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::setTimeStep(Real ts) { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getTimeStep() { return timeStep; }
/* -------------------------------------------------------------------------- */
void DomainMD1D::setPeriodicFlag(bool flag) { periodic_flag = flag; }
/* -------------------------------------------------------------------------- */
void DomainMD1D::setReferenceManagerState(bool state) {}
/* -------------------------------------------------------------------------- */
bool DomainMD1D::isPeriodic(UInt dir) { return (periodic_flag && dir == 1); }
/* -------------------------------------------------------------------------- */
/* LMDESC MD1D
This plugin is the implementation of a very simple
one dimensional molecular dynamics class.
*/
/* LMHERITANCE domain_md */
/* LMEXAMPLE
Section MultiScale AtomsUnits \\
...\\
GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\
MODEL MD1D md\\
...\\
endSection\\ \\
Section MD1D:md AtomsUnits \\
RCUT rcut \\
R0 r0 \\
MASS 39.95e-3 \\
TIMESTEP 1 \\
EPSILON 1.6567944e-21/6.94769e-21 \\
SIGMA 1.1 \\
\\
DOMAIN_GEOMETRY 1 \\
endSection \\
*/
void DomainMD1D::declareParams() {
DomainAtomic<ContainerMD1D, 1>::declareParams();
/* LMKEYWORD RCUT
Sets the cutoff radius for force calculation
*/
this->parseKeyword("RCUT", rcut);
/* LMKEYWORD R0
Sets the interatomic distance
*/
this->parseKeyword("R0", r0);
/* LMKEYWORD SHIFT
used in geometry creation
*/
this->parseKeyword("SHIFT", shift, 0.);
/* LMKEYWORD MASS
Sets the per atom mass
*/
this->parseKeyword("MASS", mass);
/* LMKEYWORD TIMESTEP
Sets the timestep
*/
this->parseKeyword("TIMESTEP", timeStep, 1.);
/* LMKEYWORD EPSILON
Sets the LJ energy parameter
*/
this->parseKeyword("EPSILON", epsilon);
/* LMKEYWORD SIGMA
Sets the LJ distance parameter
*/
this->parseKeyword("SIGMA", sigma);
/* LMKEYWORD PERIODIC
Sets the periodic boundary condition
*/
this->parseTag("PERIODIC", periodic_flag, false);
/* LMKEYWORD ELASTIC
Sets the flag to switch elastic interactions
*/
this->parseTag("ELASTIC", elastic_flag, false);
/* LMKEYWORD ATOM_FILE
Sets the file to read the atomic positions
*/
this->parseKeyword("ATOM_FILE", atom_file, "");
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help