Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83057246
set.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
Sat, Sep 14, 23:40
Size
17 KB
Mime Type
text/x-c
Expires
Mon, Sep 16, 23:40 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20790284
Attached To
rLAMMPS lammps
set.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 "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "set.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "region.h"
#include "group.h"
#include "comm.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "random_park.h"
#include "math_extra.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
enum
{
ATOM
,
GROUP
,
REGION
};
enum
{
TYPE
,
TYPE_FRACTION
,
MOLECULE
,
X
,
Y
,
Z
,
VX
,
VY
,
VZ
,
CHARGE
,
DIPOLE
,
DIPOLE_RANDOM
,
QUAT
,
QUAT_RANDOM
,
DIAMETER
,
DENSITY
,
VOLUME
,
BOND
,
ANGLE
,
DIHEDRAL
,
IMPROPER
};
/* ---------------------------------------------------------------------- */
Set
::
Set
(
LAMMPS
*
lmp
)
:
Pointers
(
lmp
)
{
PI
=
4.0
*
atan
(
1.0
);
}
/* ---------------------------------------------------------------------- */
void
Set
::
command
(
int
narg
,
char
**
arg
)
{
if
(
domain
->
box_exist
==
0
)
error
->
all
(
"Set command before simulation box is defined"
);
if
(
atom
->
natoms
==
0
)
error
->
all
(
"Set command with no atoms existing"
);
if
(
narg
<
3
)
error
->
all
(
"Illegal set command"
);
// style and ID
if
(
strcmp
(
arg
[
0
],
"atom"
)
==
0
)
style
=
ATOM
;
else
if
(
strcmp
(
arg
[
0
],
"group"
)
==
0
)
style
=
GROUP
;
else
if
(
strcmp
(
arg
[
0
],
"region"
)
==
0
)
style
=
REGION
;
else
error
->
all
(
"Illegal set command"
);
int
n
=
strlen
(
arg
[
1
])
+
1
;
id
=
new
char
[
n
];
strcpy
(
id
,
arg
[
1
]);
select
=
NULL
;
// loop over keyword/value pairs
// call appropriate routine to reset attributes
if
(
comm
->
me
==
0
&&
screen
)
fprintf
(
screen
,
"Setting atom values ...
\n
"
);
int
allcount
,
origarg
;
int
iarg
=
2
;
while
(
iarg
<
narg
)
{
count
=
0
;
origarg
=
iarg
;
if
(
strcmp
(
arg
[
iarg
],
"type"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
ivalue
<=
0
||
ivalue
>
atom
->
ntypes
)
error
->
all
(
"Invalid value in set command"
);
set
(
TYPE
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"type/fraction"
)
==
0
)
{
if
(
iarg
+
4
>
narg
)
error
->
all
(
"Illegal set command"
);
newtype
=
atoi
(
arg
[
iarg
+
1
]);
fraction
=
atof
(
arg
[
iarg
+
2
]);
ivalue
=
atoi
(
arg
[
iarg
+
3
]);
if
(
fraction
<
0.0
||
fraction
>
1.0
)
error
->
all
(
"Invalid value in set command"
);
if
(
ivalue
<=
0
)
error
->
all
(
"Invalid random number seed in set command"
);
setrandom
(
TYPE_FRACTION
);
iarg
+=
4
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"mol"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
molecule_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
MOLECULE
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"x"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
X
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"y"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
Y
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"z"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
Z
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"vx"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
VX
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"vy"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
VY
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"vz"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
set
(
VZ
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"charge"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
q_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
CHARGE
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"dipole"
)
==
0
)
{
if
(
iarg
+
4
>
narg
)
error
->
all
(
"Illegal set command"
);
xvalue
=
atof
(
arg
[
iarg
+
1
]);
yvalue
=
atof
(
arg
[
iarg
+
2
]);
zvalue
=
atof
(
arg
[
iarg
+
3
]);
if
(
!
atom
->
mu_flag
||
atom
->
dipole
==
NULL
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
DIPOLE
);
iarg
+=
4
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"dipole/random"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
mu_flag
||
atom
->
shape
==
NULL
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
)
error
->
all
(
"Invalid random number seed in set command"
);
setrandom
(
DIPOLE_RANDOM
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"quat"
)
==
0
)
{
if
(
iarg
+
5
>
narg
)
error
->
all
(
"Illegal set command"
);
xvalue
=
atof
(
arg
[
iarg
+
1
]);
yvalue
=
atof
(
arg
[
iarg
+
2
]);
zvalue
=
atof
(
arg
[
iarg
+
3
]);
wvalue
=
atof
(
arg
[
iarg
+
4
]);
if
(
!
atom
->
quat_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
QUAT
);
iarg
+=
5
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"quat/random"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
quat_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
)
error
->
all
(
"Invalid random number seed in set command"
);
setrandom
(
QUAT_RANDOM
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"diameter"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
radius_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
DIAMETER
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"density"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
/*
if (!atom->density_flag)
error->all("Cannot set this attribute for this atom style");
*/
set
(
DENSITY
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"volume"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
dvalue
=
atof
(
arg
[
iarg
+
1
]);
if
(
!
atom
->
vfrac_flag
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
set
(
VOLUME
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"bond"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
atom
->
avec
->
bonds_allow
==
0
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
||
ivalue
>
atom
->
nbondtypes
)
error
->
all
(
"Invalid value in set command"
);
topology
(
BOND
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"angle"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
atom
->
avec
->
angles_allow
==
0
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
||
ivalue
>
atom
->
nangletypes
)
error
->
all
(
"Invalid value in set command"
);
topology
(
ANGLE
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"dihedral"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
atom
->
avec
->
dihedrals_allow
==
0
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
||
ivalue
>
atom
->
ndihedraltypes
)
error
->
all
(
"Invalid value in set command"
);
topology
(
DIHEDRAL
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"improper"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
"Illegal set command"
);
ivalue
=
atoi
(
arg
[
iarg
+
1
]);
if
(
atom
->
avec
->
impropers_allow
==
0
)
error
->
all
(
"Cannot set this attribute for this atom style"
);
if
(
ivalue
<=
0
||
ivalue
>
atom
->
nimpropertypes
)
error
->
all
(
"Invalid value in set command"
);
topology
(
IMPROPER
);
iarg
+=
2
;
}
else
error
->
all
(
"Illegal set command"
);
// statistics
MPI_Allreduce
(
&
count
,
&
allcount
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
comm
->
me
==
0
)
{
if
(
screen
)
fprintf
(
screen
,
" %d settings made for %s
\n
"
,
allcount
,
arg
[
origarg
]);
if
(
logfile
)
fprintf
(
logfile
,
" %d settings made for %s
\n
"
,
allcount
,
arg
[
origarg
]);
}
}
// free local memory
delete
[]
id
;
delete
[]
select
;
}
/* ----------------------------------------------------------------------
select atoms according to ATOM, GROUP, REGION style
n = nlocal or nlocal+nghost depending on keyword
------------------------------------------------------------------------- */
void
Set
::
selection
(
int
n
)
{
delete
[]
select
;
select
=
new
int
[
n
];
if
(
style
==
ATOM
)
{
if
(
atom
->
tag_enable
==
0
)
error
->
all
(
"Cannot use set atom with no atom IDs defined"
);
int
idatom
=
atoi
(
id
);
int
*
tag
=
atom
->
tag
;
for
(
int
i
=
0
;
i
<
n
;
i
++
)
if
(
idatom
==
tag
[
i
])
select
[
i
]
=
1
;
else
select
[
i
]
=
0
;
}
else
if
(
style
==
GROUP
)
{
int
igroup
=
group
->
find
(
id
);
if
(
igroup
==
-
1
)
error
->
all
(
"Could not find set group ID"
);
int
groupbit
=
group
->
bitmask
[
igroup
];
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
n
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
select
[
i
]
=
1
;
else
select
[
i
]
=
0
;
}
else
{
int
iregion
=
domain
->
find_region
(
id
);
if
(
iregion
==
-
1
)
error
->
all
(
"Set region ID does not exist"
);
double
**
x
=
atom
->
x
;
for
(
int
i
=
0
;
i
<
n
;
i
++
)
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
select
[
i
]
=
1
;
else
select
[
i
]
=
0
;
}
}
/* ----------------------------------------------------------------------
set an owned atom property directly
------------------------------------------------------------------------- */
void
Set
::
set
(
int
keyword
)
{
if
(
keyword
==
DIPOLE
)
atom
->
check_dipole
();
selection
(
atom
->
nlocal
);
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
if
(
keyword
==
TYPE
)
atom
->
type
[
i
]
=
ivalue
;
else
if
(
keyword
==
MOLECULE
)
atom
->
molecule
[
i
]
=
ivalue
;
else
if
(
keyword
==
X
)
atom
->
x
[
i
][
0
]
=
dvalue
;
else
if
(
keyword
==
Y
)
atom
->
x
[
i
][
1
]
=
dvalue
;
else
if
(
keyword
==
Z
)
atom
->
x
[
i
][
2
]
=
dvalue
;
else
if
(
keyword
==
VX
)
atom
->
v
[
i
][
0
]
=
dvalue
;
else
if
(
keyword
==
VY
)
atom
->
v
[
i
][
1
]
=
dvalue
;
else
if
(
keyword
==
VZ
)
atom
->
v
[
i
][
2
]
=
dvalue
;
else
if
(
keyword
==
CHARGE
)
atom
->
q
[
i
]
=
dvalue
;
else
if
(
keyword
==
DIAMETER
)
{
atom
->
radius
[
i
]
=
0.5
*
dvalue
;
if
(
domain
->
dimension
==
3
)
atom
->
rmass
[
i
]
=
4.0
*
PI
/
3.0
*
atom
->
radius
[
i
]
*
atom
->
radius
[
i
]
*
atom
->
radius
[
i
]
*
atom
->
density
[
i
];
else
atom
->
rmass
[
i
]
=
PI
*
atom
->
radius
[
i
]
*
atom
->
radius
[
i
]
*
atom
->
density
[
i
];
}
else
if
(
keyword
==
DENSITY
)
{
atom
->
rmass
[
i
]
=
dvalue
;
/*
atom->density[i] = dvalue;
if (domain->dimension == 3)
atom->rmass[i] = 4.0*PI/3.0 *
atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i];
else
atom->rmass[i] = PI *
atom->radius[i]*atom->radius[i] * atom->density[i];
*/
}
else
if
(
keyword
==
VOLUME
)
atom
->
vfrac
[
i
]
=
dvalue
;
else
if
(
keyword
==
DIPOLE
)
{
if
(
atom
->
dipole
[
atom
->
type
[
i
]]
>
0.0
)
{
double
**
mu
=
atom
->
mu
;
mu
[
i
][
0
]
=
xvalue
;
mu
[
i
][
1
]
=
yvalue
;
mu
[
i
][
2
]
=
zvalue
;
double
lensq
=
mu
[
i
][
0
]
*
mu
[
i
][
0
]
+
mu
[
i
][
1
]
*
mu
[
i
][
1
]
+
mu
[
i
][
2
]
*
mu
[
i
][
2
];
double
scale
=
atom
->
dipole
[
atom
->
type
[
i
]]
/
sqrt
(
lensq
);
mu
[
i
][
0
]
*=
scale
;
mu
[
i
][
1
]
*=
scale
;
mu
[
i
][
2
]
*=
scale
;
}
}
else
if
(
keyword
==
QUAT
)
{
double
PI
=
4.0
*
atan
(
1.0
);
double
theta2
=
0.5
*
PI
*
wvalue
/
180.0
;
double
sintheta2
=
sin
(
theta2
);
double
**
quat
=
atom
->
quat
;
quat
[
i
][
0
]
=
cos
(
theta2
);
quat
[
i
][
1
]
=
xvalue
*
sintheta2
;
quat
[
i
][
2
]
=
yvalue
*
sintheta2
;
quat
[
i
][
3
]
=
zvalue
*
sintheta2
;
MathExtra
::
normalize4
(
quat
[
i
]);
}
count
++
;
}
}
/* ----------------------------------------------------------------------
set an owned atom property randomly
set seed based on atom tag
make atom result independent of what proc owns it
------------------------------------------------------------------------- */
void
Set
::
setrandom
(
int
keyword
)
{
int
i
;
if
(
keyword
==
DIPOLE_RANDOM
)
atom
->
check_dipole
();
selection
(
atom
->
nlocal
);
RanPark
*
random
=
new
RanPark
(
lmp
,
1
);
double
**
x
=
atom
->
x
;
int
seed
=
ivalue
;
// set fraction of atom types to newtype
if
(
keyword
==
TYPE_FRACTION
)
{
int
nlocal
=
atom
->
nlocal
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
random
->
reset
(
seed
,
x
[
i
]);
if
(
random
->
uniform
()
>
fraction
)
continue
;
atom
->
type
[
i
]
=
newtype
;
count
++
;
}
// set dipole moments to random orientations in 3d or 2d
// dipole length is determined by dipole type array
}
else
if
(
keyword
==
DIPOLE_RANDOM
)
{
int
*
type
=
atom
->
type
;
double
*
dipole
=
atom
->
dipole
;
double
**
mu
=
atom
->
mu
;
int
nlocal
=
atom
->
nlocal
;
double
msq
,
scale
;
if
(
domain
->
dimension
==
3
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
if
(
dipole
[
type
[
i
]]
>
0.0
)
{
random
->
reset
(
seed
,
x
[
i
]);
mu
[
i
][
0
]
=
random
->
uniform
()
-
0.5
;
mu
[
i
][
1
]
=
random
->
uniform
()
-
0.5
;
mu
[
i
][
2
]
=
random
->
uniform
()
-
0.5
;
msq
=
mu
[
i
][
0
]
*
mu
[
i
][
0
]
+
mu
[
i
][
1
]
*
mu
[
i
][
1
]
+
mu
[
i
][
2
]
*
mu
[
i
][
2
];
scale
=
dipole
[
type
[
i
]]
/
sqrt
(
msq
);
mu
[
i
][
0
]
*=
scale
;
mu
[
i
][
1
]
*=
scale
;
mu
[
i
][
2
]
*=
scale
;
count
++
;
}
}
}
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
if
(
dipole
[
type
[
i
]]
>
0.0
)
{
random
->
reset
(
seed
,
x
[
i
]);
mu
[
i
][
0
]
=
random
->
uniform
()
-
0.5
;
mu
[
i
][
1
]
=
random
->
uniform
()
-
0.5
;
mu
[
i
][
2
]
=
0.0
;
msq
=
mu
[
i
][
0
]
*
mu
[
i
][
0
]
+
mu
[
i
][
1
]
*
mu
[
i
][
1
];
scale
=
dipole
[
type
[
i
]]
/
sqrt
(
msq
);
mu
[
i
][
0
]
*=
scale
;
mu
[
i
][
1
]
*=
scale
;
count
++
;
}
}
}
// set quaternions to random orientations in 3d or 2d
// no need to normalize quats since creations algorithms already do
}
else
if
(
keyword
==
QUAT_RANDOM
)
{
double
**
quat
=
atom
->
quat
;
int
nlocal
=
atom
->
nlocal
;
if
(
domain
->
dimension
==
3
)
{
double
s
,
t1
,
t2
,
theta1
,
theta2
;
double
PI
=
4.0
*
atan
(
1.0
);
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
random
->
reset
(
seed
,
x
[
i
]);
s
=
random
->
uniform
();
t1
=
sqrt
(
1.0
-
s
);
t2
=
sqrt
(
s
);
theta1
=
2.0
*
PI
*
random
->
uniform
();
theta2
=
2.0
*
PI
*
random
->
uniform
();
quat
[
i
][
0
]
=
cos
(
theta2
)
*
t2
;
quat
[
i
][
1
]
=
sin
(
theta1
)
*
t1
;
quat
[
i
][
2
]
=
cos
(
theta1
)
*
t1
;
quat
[
i
][
3
]
=
sin
(
theta2
)
*
t2
;
count
++
;
}
}
else
{
double
theta2
;
double
PI
=
4.0
*
atan
(
1.0
);
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
select
[
i
])
{
random
->
reset
(
seed
,
x
[
i
]);
theta2
=
PI
*
random
->
uniform
();
quat
[
i
][
0
]
=
cos
(
theta2
);
quat
[
i
][
1
]
=
0.0
;
quat
[
i
][
2
]
=
0.0
;
quat
[
i
][
3
]
=
sin
(
theta2
);
count
++
;
}
}
}
delete
random
;
}
/* ---------------------------------------------------------------------- */
void
Set
::
topology
(
int
keyword
)
{
int
m
,
atom1
,
atom2
,
atom3
,
atom4
;
// border swap to acquire ghost atom info
// enforce PBC before in case atoms are outside box
// init entire system since comm->exchange is done
// comm::init needs neighbor::init needs pair::init needs kspace::init, etc
if
(
comm
->
me
==
0
&&
screen
)
fprintf
(
screen
,
" system init for set ...
\n
"
);
lmp
->
init
();
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
domain
->
reset_box
();
comm
->
setup
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
// select both owned and ghost atoms
selection
(
atom
->
nlocal
+
atom
->
nghost
);
// for BOND, each of 2 atoms must be in group
if
(
keyword
==
BOND
)
{
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
for
(
m
=
0
;
m
<
atom
->
num_bond
[
i
];
m
++
)
{
atom1
=
atom
->
map
(
atom
->
bond_atom
[
i
][
m
]);
if
(
atom1
==
-
1
)
error
->
one
(
"Bond atom missing in set command"
);
if
(
select
[
i
]
&&
select
[
atom1
])
{
atom
->
bond_type
[
i
][
m
]
=
ivalue
;
count
++
;
}
}
}
// for ANGLE, each of 3 atoms must be in group
if
(
keyword
==
ANGLE
)
{
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
for
(
m
=
0
;
m
<
atom
->
num_angle
[
i
];
m
++
)
{
atom1
=
atom
->
map
(
atom
->
angle_atom1
[
i
][
m
]);
atom2
=
atom
->
map
(
atom
->
angle_atom2
[
i
][
m
]);
atom3
=
atom
->
map
(
atom
->
angle_atom3
[
i
][
m
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
||
atom3
==
-
1
)
error
->
one
(
"Angle atom missing in set command"
);
if
(
select
[
atom1
]
&&
select
[
atom2
]
&&
select
[
atom3
])
{
atom
->
angle_type
[
i
][
m
]
=
ivalue
;
count
++
;
}
}
}
// for DIHEDRAL, each of 4 atoms must be in group
if
(
keyword
==
DIHEDRAL
)
{
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
for
(
m
=
0
;
m
<
atom
->
num_dihedral
[
i
];
m
++
)
{
atom1
=
atom
->
map
(
atom
->
dihedral_atom1
[
i
][
m
]);
atom2
=
atom
->
map
(
atom
->
dihedral_atom2
[
i
][
m
]);
atom3
=
atom
->
map
(
atom
->
dihedral_atom3
[
i
][
m
]);
atom4
=
atom
->
map
(
atom
->
dihedral_atom4
[
i
][
m
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
||
atom3
==
-
1
||
atom4
==
-
1
)
error
->
one
(
"Dihedral atom missing in set command"
);
if
(
select
[
atom1
]
&&
select
[
atom2
]
&&
select
[
atom3
]
&&
select
[
atom4
])
{
atom
->
dihedral_type
[
i
][
m
]
=
ivalue
;
count
++
;
}
}
}
// for IMPROPER, each of 4 atoms must be in group
if
(
keyword
==
IMPROPER
)
{
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
for
(
m
=
0
;
m
<
atom
->
num_improper
[
i
];
m
++
)
{
atom1
=
atom
->
map
(
atom
->
improper_atom1
[
i
][
m
]);
atom2
=
atom
->
map
(
atom
->
improper_atom2
[
i
][
m
]);
atom3
=
atom
->
map
(
atom
->
improper_atom3
[
i
][
m
]);
atom4
=
atom
->
map
(
atom
->
improper_atom4
[
i
][
m
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
||
atom3
==
-
1
||
atom4
==
-
1
)
error
->
one
(
"Improper atom missing in set command"
);
if
(
select
[
atom1
]
&&
select
[
atom2
]
&&
select
[
atom3
]
&&
select
[
atom4
])
{
atom
->
improper_type
[
i
][
m
]
=
ivalue
;
count
++
;
}
}
}
}
Event Timeline
Log In to Comment