Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82358301
domain_md1d.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 10, 23:20
Size
17 KB
Mime Type
text/x-c
Expires
Thu, Sep 12, 23:20 (2 d)
Engine
blob
Format
Raw Data
Handle
20684640
Attached To
rLIBMULTISCALE LibMultiScale
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__
Event Timeline
Log In to Comment