Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/md/1d/domain_md1d.cc
b05d5f2aba78features/reshape_comput…

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 "lm_common.hh"
#include "domain_md1d.hh"
#include "reference_manager.hh"
#include "reference_manager_md1d.hh"
#include "md_builder.hh"
#include "spatial_filter_atomic.hh"
/* -------------------------------------------------------------------------- */
//double nearbyint(double x);
__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,
const Real & domainSize){
Real rij = xi - xj;
// Real domainSize = this->getDomainSize();
if (this->periodic_flag)
rij -= domainSize*nearbyint(rij/domainSize);
rij = fabs(rij);
return rij;
}
/* -------------------------------------------------------------------------- */
Real DomainMD1D::getDomainSize(){
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->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);
dofs[i].a = 0.0;
dofs[i].f = 0.0;
}
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();
}
/* -------------------------------------------------------------------------- */
void DomainMD1D::performStep2() {
Real domainSize = this->getDomainSize();
if(elastic_flag){
quadraticForceComputation();
return;
}
LM_ASSERT(dofs[3].m != 0,"mass pb");
Real contribution = 0.0, rij ;
int nb_atoms = dofs.size();
int droite=-3*nb_atoms, j,i,gauche=3*nb_atoms;
for (i = 0 ;i < nb_atoms;++i){
droite = -3*nb_atoms;
for(j=i+1;j<nb_atoms;++j)
{
rij = this->computeDistance(dofs[j].p,dofs[i].p,domainSize);
if (rij <= rcut)
droite = j;
else if(droite != -3*nb_atoms)
break;
}
LM_ASSERT(dofs[3].m != 0,"mass pb");
if (droite == -3*nb_atoms && i!= nb_atoms-1)
{
// LM_FATAL("l'atome " << i << " a perdu tout ses voisins de droite");
droite = i;
}
else if (i==nb_atoms-1)
{droite = i;}
while(droite != i)
{
if(droite !=i ){
rij = computeDistance(dofs[i].p,dofs[droite].p,domainSize);
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[droite].alpha);
DUMP("contribution between " << i << " and " << droite << " = "
<< contribution << " [d=" << rij << "]",DBG_ALL);
dofs[i].a -= contribution ;
dofs[droite].a += contribution ;
dofs[i].f = dofs[i].a * dofs[i].m;
dofs[droite].f = dofs[droite].a * dofs[droite].m;
droite--;
}
}
//periodic boundary treatment
DUMP(i << " " << nb_atoms << " " << periodic_flag << " "
<< computeDistance(dofs[nb_atoms-1].p,dofs[i].p,domainSize) << " "
<< rcut,DBG_ALL);
if (periodic_flag && i > nb_atoms/2){
gauche = 3*nb_atoms;
for(j=0;j<2*((int)(rcut/r0));++j) {
rij = computeDistance(dofs[j].p,dofs[i].p,domainSize);
if (rij - domainSize*nearbyint(rij/domainSize) <=rcut) gauche = j;
else if(gauche != 3*nb_atoms) break;
}
if (gauche == 3*nb_atoms) continue;
DUMP("for atom " << i << " left is " << gauche,DBG_ALL);
while(gauche != -1) {
if(gauche != -1){
rij = computeDistance(dofs[i].p,dofs[gauche].p,domainSize);
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[gauche].alpha);
DUMP("contribution between " << i << " and " << gauche
<< " = " << contribution << " [d=" << rij << "]",DBG_ALL);
dofs[i].a -= contribution ;
dofs[gauche].a += contribution ;
dofs[i].f = dofs[i].a * dofs[i].m;
dofs[gauche].f = dofs[gauche].a * dofs[gauche].m;
gauche--;
}
}
}
}
DUMP("before mass weighting " << 3 << " : "
<< dofs[3].a << " " << dofs[3].m,DBG_DETAIL);
for (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() {
Real domainSize = this->getDomainSize();
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,domainSize);
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,domainSize);
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(0) << " 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