Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85130430
create_atoms.cpp
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
Thu, Sep 26, 23:49
Size
7 KB
Mime Type
text/x-c
Expires
Sat, Sep 28, 23:49 (2 d)
Engine
blob
Format
Raw Data
Handle
21132769
Attached To
rLAMMPS lammps
create_atoms.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "create_atoms.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_hybrid.h"
#include "comm.h"
#include "domain.h"
#include "lattice.h"
#include "region.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
#define MAXATOMS 0x7FFFFFFF
#define BIG 1.0e30
#define EPSILON 1.0e-6
/* ---------------------------------------------------------------------- */
CreateAtoms
::
CreateAtoms
(
LAMMPS
*
lmp
)
:
Pointers
(
lmp
)
{}
/* ---------------------------------------------------------------------- */
void
CreateAtoms
::
command
(
int
narg
,
char
**
arg
)
{
if
(
domain
->
box_exist
==
0
)
error
->
all
(
"Create_atoms command before simulation box is defined"
);
if
(
domain
->
lattice
==
NULL
)
error
->
all
(
"Cannot create atoms with undefined lattice"
);
// parse arguments
int
nbasis
=
domain
->
lattice
->
nbasis
;
int
*
basistype
=
new
int
[
nbasis
];
if
(
narg
<
1
)
error
->
all
(
"Illegal create_atoms command"
);
int
itype
=
atoi
(
arg
[
0
]);
if
(
itype
<=
0
||
itype
>
atom
->
ntypes
)
error
->
all
(
"Invalid atom type in create_atoms command"
);
for
(
int
i
=
0
;
i
<
nbasis
;
i
++
)
basistype
[
i
]
=
itype
;
regionflag
=
-
1
;
nhybrid
=
0
;
int
iarg
=
1
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"region"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal create_atoms command"
);
int
iregion
;
for
(
iregion
=
0
;
iregion
<
domain
->
nregion
;
iregion
++
)
if
(
strcmp
(
arg
[
iarg
+
1
],
domain
->
regions
[
iregion
]
->
id
)
==
0
)
break
;
if
(
iregion
==
domain
->
nregion
)
error
->
all
(
"Create_atoms region ID does not exist"
);
regionflag
=
iregion
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"basis"
)
==
0
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
"Illegal create_atoms command"
);
int
ibasis
=
atoi
(
arg
[
iarg
+
1
]);
itype
=
atoi
(
arg
[
iarg
+
2
]);
if
(
ibasis
<=
0
||
ibasis
>
nbasis
||
itype
<=
0
||
itype
>
atom
->
ntypes
)
error
->
all
(
"Illegal create_atoms command"
);
basistype
[
ibasis
-
1
]
=
itype
;
iarg
+=
3
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"hybrid"
)
==
0
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
"Illegal create_atoms command"
);
AtomVecHybrid
*
avec_hybrid
=
(
AtomVecHybrid
*
)
atom
->
avec
;
int
ihybrid
;
for
(
ihybrid
=
0
;
ihybrid
<
avec_hybrid
->
nstyles
;
ihybrid
++
)
if
(
strcmp
(
avec_hybrid
->
keywords
[
ihybrid
],
arg
[
iarg
+
1
])
==
0
)
break
;
if
(
ihybrid
==
avec_hybrid
->
nstyles
)
error
->
all
(
"Create atoms hybrid sub-style does not exist"
);
nhybrid
=
ihybrid
;
iarg
+=
3
;
}
else
error
->
all
(
"Illegal create_atoms command"
);
}
// convert 8 corners of my sub-box from box coords to lattice coords
// min to max = bounding box around the pts in lattice space
subxlo
=
domain
->
subxlo
;
subxhi
=
domain
->
subxhi
;
subylo
=
domain
->
subylo
;
subyhi
=
domain
->
subyhi
;
subzlo
=
domain
->
subzlo
;
subzhi
=
domain
->
subzhi
;
double
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
;
xmin
=
ymin
=
zmin
=
BIG
;
xmax
=
ymax
=
zmax
=
-
BIG
;
domain
->
lattice
->
bbox
(
1
,
subxlo
,
subylo
,
subzlo
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxhi
,
subylo
,
subzlo
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxlo
,
subyhi
,
subzlo
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxhi
,
subyhi
,
subzlo
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxlo
,
subylo
,
subzhi
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxhi
,
subylo
,
subzhi
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxlo
,
subyhi
,
subzhi
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
domain
->
lattice
->
bbox
(
1
,
subxhi
,
subyhi
,
subzhi
,
xmin
,
ymin
,
zmin
,
xmax
,
ymax
,
zmax
);
// ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my sub-box
// overlap = any part of a unit cell (face,edge,pt) in common with my sub-box
// in lattice space, sub-box is a tilted box
// but bbox of sub-box is aligned with lattice axes
// so ilo:khi unit cells should completely tile bounding box
// decrement lo values if min < 0, since static_cast(-1.5) = -1
int
ilo
,
ihi
,
jlo
,
jhi
,
klo
,
khi
;
ilo
=
static_cast
<
int
>
(
xmin
);
jlo
=
static_cast
<
int
>
(
ymin
);
klo
=
static_cast
<
int
>
(
zmin
);
ihi
=
static_cast
<
int
>
(
xmax
);
jhi
=
static_cast
<
int
>
(
ymax
);
khi
=
static_cast
<
int
>
(
zmax
);
if
(
xmin
<
0.0
)
ilo
--
;
if
(
ymin
<
0.0
)
jlo
--
;
if
(
zmin
<
0.0
)
klo
--
;
// iterate on 3d periodic lattice using loop bounds
// invoke add_atom for nbasis atoms in each unit cell
// add_atom converts lattice coords to box coords, checks if in my sub-box
boxxhi
=
domain
->
boxxhi
;
boxyhi
=
domain
->
boxyhi
;
boxzhi
=
domain
->
boxzhi
;
double
natoms_previous
=
atom
->
natoms
;
int
nlocal_previous
=
atom
->
nlocal
;
double
**
basis
=
domain
->
lattice
->
basis
;
int
i
,
j
,
k
,
m
;
for
(
k
=
klo
;
k
<=
khi
;
k
++
)
for
(
j
=
jlo
;
j
<=
jhi
;
j
++
)
for
(
i
=
ilo
;
i
<=
ihi
;
i
++
)
for
(
m
=
0
;
m
<
nbasis
;
m
++
)
add_atom
(
basistype
[
m
],
i
+
basis
[
m
][
0
],
j
+
basis
[
m
][
1
],
k
+
basis
[
m
][
2
]);
delete
[]
basistype
;
// new total # of atoms
double
rlocal
=
atom
->
nlocal
;
MPI_Allreduce
(
&
rlocal
,
&
atom
->
natoms
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
// print status
if
(
comm
->
me
==
0
)
{
if
(
screen
)
fprintf
(
screen
,
"Created %.15g atoms
\n
"
,
atom
->
natoms
-
natoms_previous
);
if
(
logfile
)
fprintf
(
logfile
,
"Created %.15g atoms
\n
"
,
atom
->
natoms
-
natoms_previous
);
}
// reset simulation now that more atoms are defined
// add tags for newly created atoms if possible
// if global map exists, reset it
// if a molecular system, set nspecial to 0 for new atoms
if
(
atom
->
natoms
>
MAXATOMS
)
atom
->
tag_enable
=
0
;
if
(
atom
->
natoms
<=
MAXATOMS
)
atom
->
tag_extend
();
if
(
atom
->
map_style
)
{
atom
->
map_init
();
atom
->
map_set
();
}
if
(
atom
->
molecular
)
{
int
**
nspecial
=
atom
->
nspecial
;
for
(
i
=
nlocal_previous
;
i
<
atom
->
nlocal
;
i
++
)
{
nspecial
[
i
][
0
]
=
0
;
nspecial
[
i
][
1
]
=
0
;
nspecial
[
i
][
2
]
=
0
;
}
}
}
/* ----------------------------------------------------------------------
add an atom of type at lattice coords x,y,z if it meets all criteria
------------------------------------------------------------------------- */
void
CreateAtoms
::
add_atom
(
int
ntype
,
double
x
,
double
y
,
double
z
)
{
// convert from lattice coords to box coords
domain
->
lattice
->
lattice2box
(
x
,
y
,
z
);
// if a region was specified, test if atom is in it
if
(
regionflag
>=
0
)
if
(
!
domain
->
regions
[
regionflag
]
->
match
(
x
,
y
,
z
))
return
;
// test if atom is in my subbox
if
(
x
<
subxlo
||
x
>=
subxhi
||
y
<
subylo
||
y
>=
subyhi
||
z
<
subzlo
||
z
>=
subzhi
)
return
;
// don't put atoms within EPSILON of upper periodic boundary
// else may overlap image atom at lower boundary
if
(
domain
->
xperiodic
&&
fabs
(
x
-
boxxhi
)
<
EPSILON
)
return
;
if
(
domain
->
yperiodic
&&
fabs
(
y
-
boxyhi
)
<
EPSILON
)
return
;
if
(
domain
->
zperiodic
&&
fabs
(
z
-
boxzhi
)
<
EPSILON
)
return
;
// add the atom to my list of atoms
atom
->
avec
->
create_atom
(
ntype
,
x
,
y
,
z
,
nhybrid
);
}
Event Timeline
Log In to Comment