Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86892091
group.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
Wed, Oct 9, 06:05
Size
24 KB
Mime Type
text/x-c
Expires
Fri, Oct 11, 06:05 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
21501187
Attached To
rLAMMPS lammps
group.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 "mpi.h"
#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "group.h"
#include "domain.h"
#include "atom.h"
#include "force.h"
#include "region.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
#define MAX_GROUP 32
enum
{
TYPE
,
MOLECULE
,
ID
};
enum
{
LT
,
LE
,
GT
,
GE
,
EQ
,
NEQ
,
BETWEEN
};
#define BIG 1.0e20
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
/* ----------------------------------------------------------------------
initialize group memory
------------------------------------------------------------------------- */
Group
::
Group
(
LAMMPS
*
lmp
)
:
Pointers
(
lmp
)
{
MPI_Comm_rank
(
world
,
&
me
);
names
=
(
char
**
)
memory
->
smalloc
(
MAX_GROUP
*
sizeof
(
char
*
),
"group:names"
);
bitmask
=
new
int
[
MAX_GROUP
];
inversemask
=
new
int
[
MAX_GROUP
];
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
bitmask
[
i
]
=
1
<<
i
;
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
inversemask
[
i
]
=
bitmask
[
i
]
^
~
0
;
// create "all" group
char
*
str
=
(
char
*
)
"all"
;
int
n
=
strlen
(
str
)
+
1
;
names
[
0
]
=
(
char
*
)
memory
->
smalloc
(
n
*
sizeof
(
char
),
"group:names[]"
);
strcpy
(
names
[
0
],
str
);
ngroup
=
1
;
}
/* ----------------------------------------------------------------------
free all memory
------------------------------------------------------------------------- */
Group
::~
Group
()
{
for
(
int
i
=
0
;
i
<
ngroup
;
i
++
)
memory
->
sfree
(
names
[
i
]);
memory
->
sfree
(
names
);
delete
[]
bitmask
;
delete
[]
inversemask
;
}
/* ----------------------------------------------------------------------
assign atoms to a new or existing group
------------------------------------------------------------------------- */
void
Group
::
assign
(
int
narg
,
char
**
arg
)
{
int
i
;
if
(
domain
->
box_exist
==
0
)
error
->
all
(
"Group command before simulation box is defined"
);
if
(
narg
<
2
)
error
->
all
(
"Illegal group command"
);
// find group in existing list
// igroup = -1 is a new group name, add it
int
igroup
=
find
(
arg
[
0
]);
if
(
igroup
==
-
1
)
{
if
(
ngroup
==
MAX_GROUP
)
error
->
all
(
"Too many groups"
);
igroup
=
ngroup
;
ngroup
++
;
int
n
=
strlen
(
arg
[
0
])
+
1
;
names
[
igroup
]
=
(
char
*
)
memory
->
smalloc
(
n
*
sizeof
(
char
),
"group:names[]"
);
strcpy
(
names
[
igroup
],
arg
[
0
]);
}
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
int
bit
=
bitmask
[
igroup
];
// style = region
// add to group if atom is in region
if
(
strcmp
(
arg
[
1
],
"region"
)
==
0
)
{
if
(
narg
!=
3
)
error
->
all
(
"Illegal group command"
);
int
iregion
=
domain
->
find_region
(
arg
[
2
]);
if
(
iregion
==
-
1
)
error
->
all
(
"Group region ID does not exist"
);
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
mask
[
i
]
|=
bit
;
// style = logical condition
}
else
if
(
narg
>=
3
&&
(
strcmp
(
arg
[
2
],
"<"
)
==
0
||
strcmp
(
arg
[
2
],
">"
)
==
0
||
strcmp
(
arg
[
2
],
"<="
)
==
0
||
strcmp
(
arg
[
2
],
">="
)
==
0
||
strcmp
(
arg
[
2
],
"<>"
)
==
0
))
{
if
(
narg
<
4
||
narg
>
5
)
error
->
all
(
"Illegal group command"
);
int
category
,
condition
,
bound1
,
bound2
;
if
(
strcmp
(
arg
[
1
],
"type"
)
==
0
)
category
=
TYPE
;
else
if
(
strcmp
(
arg
[
1
],
"molecule"
)
==
0
)
category
=
MOLECULE
;
else
if
(
strcmp
(
arg
[
1
],
"id"
)
==
0
)
category
=
ID
;
else
error
->
all
(
"Illegal group command"
);
if
(
strcmp
(
arg
[
2
],
"<"
)
==
0
)
condition
=
LT
;
else
if
(
strcmp
(
arg
[
2
],
"<="
)
==
0
)
condition
=
LE
;
else
if
(
strcmp
(
arg
[
2
],
">"
)
==
0
)
condition
=
GT
;
else
if
(
strcmp
(
arg
[
2
],
">="
)
==
0
)
condition
=
GE
;
else
if
(
strcmp
(
arg
[
2
],
"=="
)
==
0
)
condition
=
EQ
;
else
if
(
strcmp
(
arg
[
2
],
"!="
)
==
0
)
condition
=
NEQ
;
else
if
(
strcmp
(
arg
[
2
],
"<>"
)
==
0
)
condition
=
BETWEEN
;
else
error
->
all
(
"Illegal group command"
);
bound1
=
atoi
(
arg
[
3
]);
bound2
=
-
1
;
if
(
condition
==
BETWEEN
)
{
if
(
narg
!=
5
)
error
->
all
(
"Illegal group command"
);
bound2
=
atoi
(
arg
[
4
]);
}
int
*
attribute
;
if
(
category
==
TYPE
)
attribute
=
atom
->
type
;
else
if
(
category
==
MOLECULE
)
attribute
=
atom
->
molecule
;
else
if
(
category
==
ID
)
attribute
=
atom
->
tag
;
// add to group if meets condition
if
(
condition
==
LT
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
<
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
LE
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
<=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
GT
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
>
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
GE
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
>=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
EQ
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
==
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
NEQ
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
!=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
BETWEEN
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
>=
bound1
&&
attribute
[
i
]
<=
bound2
)
mask
[
i
]
|=
bit
;
}
// style = list of values
}
else
if
(
strcmp
(
arg
[
1
],
"type"
)
==
0
||
strcmp
(
arg
[
1
],
"molecule"
)
==
0
||
strcmp
(
arg
[
1
],
"id"
)
==
0
)
{
if
(
narg
<
3
)
error
->
all
(
"Illegal group command"
);
int
length
=
narg
-
2
;
int
*
list
=
new
int
[
length
];
int
category
;
if
(
strcmp
(
arg
[
1
],
"type"
)
==
0
)
category
=
TYPE
;
else
if
(
strcmp
(
arg
[
1
],
"molecule"
)
==
0
)
category
=
MOLECULE
;
else
if
(
strcmp
(
arg
[
1
],
"id"
)
==
0
)
category
=
ID
;
else
error
->
all
(
"Illegal group command"
);
length
=
narg
-
2
;
for
(
int
iarg
=
2
;
iarg
<
narg
;
iarg
++
)
list
[
iarg
-
2
]
=
atoi
(
arg
[
iarg
]);
int
*
attribute
;
if
(
category
==
TYPE
)
attribute
=
atom
->
type
;
else
if
(
category
==
MOLECULE
)
attribute
=
atom
->
molecule
;
else
if
(
category
==
ID
)
attribute
=
atom
->
tag
;
// add to group if attribute is any in list
for
(
int
ilist
=
0
;
ilist
<
length
;
ilist
++
)
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
==
list
[
ilist
])
mask
[
i
]
|=
bit
;
delete
[]
list
;
// style = subtract
}
else
if
(
strcmp
(
arg
[
1
],
"subtract"
)
==
0
)
{
if
(
narg
<
4
)
error
->
all
(
"Illegal group command"
);
int
length
=
narg
-
2
;
int
*
list
=
new
int
[
length
];
int
jgroup
;
for
(
int
iarg
=
2
;
iarg
<
narg
;
iarg
++
)
{
jgroup
=
find
(
arg
[
iarg
]);
if
(
jgroup
==
-
1
)
error
->
all
(
"Group ID does not exist"
);
list
[
iarg
-
2
]
=
jgroup
;
}
// add to group if in 1st group in list
int
otherbit
=
bitmask
[
list
[
0
]];
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
otherbit
)
mask
[
i
]
|=
bit
;
// remove atoms if they are in any of the other groups
// AND with inverse mask removes the atom from group
int
inverse
=
inversemask
[
igroup
];
for
(
int
ilist
=
1
;
ilist
<
length
;
ilist
++
)
{
otherbit
=
bitmask
[
list
[
ilist
]];
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
otherbit
)
mask
[
i
]
&=
inverse
;
}
delete
[]
list
;
// style = union
}
else
if
(
strcmp
(
arg
[
1
],
"union"
)
==
0
)
{
if
(
narg
<
3
)
error
->
all
(
"Illegal group command"
);
int
length
=
narg
-
2
;
int
*
list
=
new
int
[
length
];
int
jgroup
;
for
(
int
iarg
=
2
;
iarg
<
narg
;
iarg
++
)
{
jgroup
=
find
(
arg
[
iarg
]);
if
(
jgroup
==
-
1
)
error
->
all
(
"Group ID does not exist"
);
list
[
iarg
-
2
]
=
jgroup
;
}
// add to group if in any other group in list
int
otherbit
;
for
(
int
ilist
=
0
;
ilist
<
length
;
ilist
++
)
{
otherbit
=
bitmask
[
list
[
ilist
]];
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
otherbit
)
mask
[
i
]
|=
bit
;
}
delete
[]
list
;
// style = intersect
}
else
if
(
strcmp
(
arg
[
1
],
"intersect"
)
==
0
)
{
if
(
narg
<
4
)
error
->
all
(
"Illegal group command"
);
int
length
=
narg
-
2
;
int
*
list
=
new
int
[
length
];
int
jgroup
;
for
(
int
iarg
=
2
;
iarg
<
narg
;
iarg
++
)
{
jgroup
=
find
(
arg
[
iarg
]);
if
(
jgroup
==
-
1
)
error
->
all
(
"Group ID does not exist"
);
list
[
iarg
-
2
]
=
jgroup
;
}
// add to group if in all groups in list
int
otherbit
,
ok
,
ilist
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
ok
=
1
;
for
(
ilist
=
0
;
ilist
<
length
;
ilist
++
)
{
otherbit
=
bitmask
[
list
[
ilist
]];
if
((
mask
[
i
]
&
otherbit
)
==
0
)
ok
=
0
;
}
if
(
ok
)
mask
[
i
]
|=
bit
;
}
delete
[]
list
;
// not a valid group style
}
else
error
->
all
(
"Illegal group command"
);
// print stats for changed group
int
n
;
n
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
bit
)
n
++
;
double
rlocal
=
n
;
double
all
;
MPI_Allreduce
(
&
rlocal
,
&
all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
me
==
0
)
{
if
(
screen
)
fprintf
(
screen
,
"%.15g atoms in group %s
\n
"
,
all
,
names
[
igroup
]);
if
(
logfile
)
fprintf
(
logfile
,
"%.15g atoms in group %s
\n
"
,
all
,
names
[
igroup
]);
}
}
/* ----------------------------------------------------------------------
add flagged atoms to a new or existing group
------------------------------------------------------------------------- */
void
Group
::
create
(
char
*
name
,
int
*
flag
)
{
int
i
;
// find group in existing list
// igroup = -1 is a new group name, add it
int
igroup
=
find
(
name
);
if
(
igroup
==
-
1
)
{
if
(
ngroup
==
MAX_GROUP
)
error
->
all
(
"Too many groups"
);
igroup
=
ngroup
;
ngroup
++
;
int
n
=
strlen
(
name
)
+
1
;
names
[
igroup
]
=
(
char
*
)
memory
->
smalloc
(
n
*
sizeof
(
char
),
"group:names[]"
);
strcpy
(
names
[
igroup
],
name
);
}
// add atom to group if flag is set
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
int
bit
=
bitmask
[
igroup
];
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
flag
[
i
])
mask
[
i
]
|=
bit
;
}
/* ----------------------------------------------------------------------
return group index if name matches existing group, -1 if no such group
------------------------------------------------------------------------- */
int
Group
::
find
(
const
char
*
name
)
{
for
(
int
igroup
=
0
;
igroup
<
ngroup
;
igroup
++
)
if
(
strcmp
(
name
,
names
[
igroup
])
==
0
)
return
igroup
;
return
-
1
;
}
/* ----------------------------------------------------------------------
write group info to a restart file
only called by proc 0
------------------------------------------------------------------------- */
void
Group
::
write_restart
(
FILE
*
fp
)
{
fwrite
(
&
ngroup
,
sizeof
(
int
),
1
,
fp
);
int
n
;
for
(
int
i
=
0
;
i
<
ngroup
;
i
++
)
{
n
=
strlen
(
names
[
i
])
+
1
;
fwrite
(
&
n
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
names
[
i
],
sizeof
(
char
),
n
,
fp
);
}
}
/* ----------------------------------------------------------------------
read group info from a restart file
proc 0 reads, bcast to all procs
------------------------------------------------------------------------- */
void
Group
::
read_restart
(
FILE
*
fp
)
{
int
i
,
n
;
for
(
i
=
0
;
i
<
ngroup
;
i
++
)
memory
->
sfree
(
names
[
i
]);
if
(
me
==
0
)
fread
(
&
ngroup
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
ngroup
,
1
,
MPI_INT
,
0
,
world
);
for
(
i
=
0
;
i
<
ngroup
;
i
++
)
{
if
(
me
==
0
)
fread
(
&
n
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
names
[
i
]
=
(
char
*
)
memory
->
smalloc
(
n
*
sizeof
(
char
),
"group:names[]"
);
if
(
me
==
0
)
fread
(
names
[
i
],
sizeof
(
char
),
n
,
fp
);
MPI_Bcast
(
names
[
i
],
n
,
MPI_CHAR
,
0
,
world
);
}
}
// ----------------------------------------------------------------------
// computations on a group of atoms
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
count atoms in group
compute in double precision in case system is huge
------------------------------------------------------------------------- */
double
Group
::
count
(
int
igroup
)
{
int
groupbit
=
bitmask
[
igroup
];
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
n
++
;
double
nsingle
=
n
;
double
nall
;
MPI_Allreduce
(
&
nsingle
,
&
nall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
return
nall
;
}
/* ----------------------------------------------------------------------
compute the total mass of group of atoms
use either per-type mass or per-atom rmass
------------------------------------------------------------------------- */
double
Group
::
mass
(
int
igroup
)
{
int
groupbit
=
bitmask
[
igroup
];
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
one
=
0.0
;
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
one
+=
rmass
[
i
];
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
one
+=
mass
[
type
[
i
]];
}
double
all
;
MPI_Allreduce
(
&
one
,
&
all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
return
all
;
}
/* ----------------------------------------------------------------------
compute the total charge of group of atoms
------------------------------------------------------------------------- */
double
Group
::
charge
(
int
igroup
)
{
int
groupbit
=
bitmask
[
igroup
];
int
*
mask
=
atom
->
mask
;
double
*
q
=
atom
->
q
;
int
nlocal
=
atom
->
nlocal
;
double
qone
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
qone
+=
q
[
i
];
double
qall
;
MPI_Allreduce
(
&
qone
,
&
qall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
return
qall
;
}
/* ----------------------------------------------------------------------
compute the coordinate bounds of the group of atoms
periodic images are not considered, so atoms are NOT unwrapped
------------------------------------------------------------------------- */
void
Group
::
bounds
(
int
igroup
,
double
*
minmax
)
{
int
groupbit
=
bitmask
[
igroup
];
double
extent
[
6
];
extent
[
0
]
=
extent
[
2
]
=
extent
[
4
]
=
BIG
;
extent
[
1
]
=
extent
[
3
]
=
extent
[
5
]
=
-
BIG
;
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
extent
[
0
]
=
MIN
(
extent
[
0
],
x
[
i
][
0
]);
extent
[
1
]
=
MAX
(
extent
[
1
],
x
[
i
][
0
]);
extent
[
2
]
=
MIN
(
extent
[
2
],
x
[
i
][
1
]);
extent
[
3
]
=
MAX
(
extent
[
3
],
x
[
i
][
1
]);
extent
[
4
]
=
MIN
(
extent
[
4
],
x
[
i
][
2
]);
extent
[
5
]
=
MAX
(
extent
[
5
],
x
[
i
][
2
]);
}
}
// compute extent across all procs
// flip sign of MIN to do it in one Allreduce MAX
// set box by extent in shrink-wrapped dims
extent
[
0
]
=
-
extent
[
0
];
extent
[
2
]
=
-
extent
[
2
];
extent
[
4
]
=
-
extent
[
4
];
MPI_Allreduce
(
extent
,
minmax
,
6
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
minmax
[
0
]
=
-
minmax
[
0
];
minmax
[
2
]
=
-
minmax
[
2
];
minmax
[
4
]
=
-
minmax
[
4
];
}
/* ----------------------------------------------------------------------
compute the center-of-mass coords of group with total mass = masstotal
return center-of-mass coords in cm[]
must unwrap atoms to compute center-of-mass correctly
------------------------------------------------------------------------- */
void
Group
::
xcm
(
int
igroup
,
double
masstotal
,
double
*
cm
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
int
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
cmone
[
3
];
cmone
[
0
]
=
cmone
[
1
]
=
cmone
[
2
]
=
0.0
;
double
xprd
=
domain
->
xprd
;
double
yprd
=
domain
->
yprd
;
double
zprd
=
domain
->
zprd
;
int
xbox
,
ybox
,
zbox
;
double
massone
;
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
xbox
=
(
image
[
i
]
&
1023
)
-
512
;
ybox
=
(
image
[
i
]
>>
10
&
1023
)
-
512
;
zbox
=
(
image
[
i
]
>>
20
)
-
512
;
massone
=
rmass
[
i
];
cmone
[
0
]
+=
(
x
[
i
][
0
]
+
xbox
*
xprd
)
*
massone
;
cmone
[
1
]
+=
(
x
[
i
][
1
]
+
ybox
*
yprd
)
*
massone
;
cmone
[
2
]
+=
(
x
[
i
][
2
]
+
zbox
*
zprd
)
*
massone
;
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
xbox
=
(
image
[
i
]
&
1023
)
-
512
;
ybox
=
(
image
[
i
]
>>
10
&
1023
)
-
512
;
zbox
=
(
image
[
i
]
>>
20
)
-
512
;
massone
=
mass
[
type
[
i
]];
cmone
[
0
]
+=
(
x
[
i
][
0
]
+
xbox
*
xprd
)
*
massone
;
cmone
[
1
]
+=
(
x
[
i
][
1
]
+
ybox
*
yprd
)
*
massone
;
cmone
[
2
]
+=
(
x
[
i
][
2
]
+
zbox
*
zprd
)
*
massone
;
}
}
MPI_Allreduce
(
cmone
,
cm
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
cm
[
0
]
/=
masstotal
;
cm
[
1
]
/=
masstotal
;
cm
[
2
]
/=
masstotal
;
}
/* ----------------------------------------------------------------------
compute the center-of-mass velocity of group with total mass = masstotal
return center-of-mass velocity in cm[]
------------------------------------------------------------------------- */
void
Group
::
vcm
(
int
igroup
,
double
masstotal
,
double
*
cm
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
v
=
atom
->
v
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
p
[
3
],
massone
;
p
[
0
]
=
p
[
1
]
=
p
[
2
]
=
0.0
;
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
massone
=
rmass
[
i
];
p
[
0
]
+=
v
[
i
][
0
]
*
massone
;
p
[
1
]
+=
v
[
i
][
1
]
*
massone
;
p
[
2
]
+=
v
[
i
][
2
]
*
massone
;
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
massone
=
mass
[
type
[
i
]];
p
[
0
]
+=
v
[
i
][
0
]
*
massone
;
p
[
1
]
+=
v
[
i
][
1
]
*
massone
;
p
[
2
]
+=
v
[
i
][
2
]
*
massone
;
}
}
MPI_Allreduce
(
p
,
cm
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
cm
[
0
]
/=
masstotal
;
cm
[
1
]
/=
masstotal
;
cm
[
2
]
/=
masstotal
;
}
/* ----------------------------------------------------------------------
compute the total force on group of atoms
------------------------------------------------------------------------- */
void
Group
::
fcm
(
int
igroup
,
double
*
cm
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
f
=
atom
->
f
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
double
flocal
[
3
];
flocal
[
0
]
=
flocal
[
1
]
=
flocal
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
flocal
[
0
]
+=
f
[
i
][
0
];
flocal
[
1
]
+=
f
[
i
][
1
];
flocal
[
2
]
+=
f
[
i
][
2
];
}
MPI_Allreduce
(
flocal
,
cm
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
/* ----------------------------------------------------------------------
compute the total kinetic energy of group and return it
------------------------------------------------------------------------- */
double
Group
::
ke
(
int
igroup
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
v
=
atom
->
v
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
one
=
0.0
;
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
one
+=
(
v
[
i
][
0
]
*
v
[
i
][
0
]
+
v
[
i
][
1
]
*
v
[
i
][
1
]
+
v
[
i
][
2
]
*
v
[
i
][
2
])
*
rmass
[
i
];
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
one
+=
(
v
[
i
][
0
]
*
v
[
i
][
0
]
+
v
[
i
][
1
]
*
v
[
i
][
1
]
+
v
[
i
][
2
]
*
v
[
i
][
2
])
*
mass
[
type
[
i
]];
}
double
all
;
MPI_Allreduce
(
&
one
,
&
all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
all
*=
0.5
*
force
->
mvv2e
;
return
all
;
}
/* ----------------------------------------------------------------------
compute the radius-of-gyration of group around center-of-mass cm
must unwrap atoms to compute Rg correctly
------------------------------------------------------------------------- */
double
Group
::
gyration
(
int
igroup
,
double
masstotal
,
double
*
cm
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
int
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
int
xbox
,
ybox
,
zbox
;
double
dx
,
dy
,
dz
,
massone
;
double
xprd
=
domain
->
xprd
;
double
yprd
=
domain
->
yprd
;
double
zprd
=
domain
->
zprd
;
double
rg
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
xbox
=
(
image
[
i
]
&
1023
)
-
512
;
ybox
=
(
image
[
i
]
>>
10
&
1023
)
-
512
;
zbox
=
(
image
[
i
]
>>
20
)
-
512
;
dx
=
(
x
[
i
][
0
]
+
xbox
*
xprd
)
-
cm
[
0
];
dy
=
(
x
[
i
][
1
]
+
ybox
*
yprd
)
-
cm
[
1
];
dz
=
(
x
[
i
][
2
]
+
zbox
*
zprd
)
-
cm
[
2
];
if
(
rmass
)
massone
=
rmass
[
i
];
else
massone
=
mass
[
type
[
i
]];
rg
+=
(
dx
*
dx
+
dy
*
dy
+
dz
*
dz
)
*
massone
;
}
double
rg_all
;
MPI_Allreduce
(
&
rg
,
&
rg_all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
return
sqrt
(
rg_all
/
masstotal
);
}
/* ----------------------------------------------------------------------
compute the angular momentum L (lmom) of group around center-of-mass cm
must unwrap atoms to compute L correctly
------------------------------------------------------------------------- */
void
Group
::
angmom
(
int
igroup
,
double
*
cm
,
double
*
lmom
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
x
=
atom
->
x
;
double
**
v
=
atom
->
v
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
int
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
int
xbox
,
ybox
,
zbox
;
double
dx
,
dy
,
dz
,
massone
;
double
xprd
=
domain
->
xprd
;
double
yprd
=
domain
->
yprd
;
double
zprd
=
domain
->
zprd
;
double
p
[
3
];
p
[
0
]
=
p
[
1
]
=
p
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
xbox
=
(
image
[
i
]
&
1023
)
-
512
;
ybox
=
(
image
[
i
]
>>
10
&
1023
)
-
512
;
zbox
=
(
image
[
i
]
>>
20
)
-
512
;
dx
=
(
x
[
i
][
0
]
+
xbox
*
xprd
)
-
cm
[
0
];
dy
=
(
x
[
i
][
1
]
+
ybox
*
yprd
)
-
cm
[
1
];
dz
=
(
x
[
i
][
2
]
+
zbox
*
zprd
)
-
cm
[
2
];
if
(
rmass
)
massone
=
rmass
[
i
];
else
massone
=
mass
[
type
[
i
]];
p
[
0
]
+=
massone
*
(
dy
*
v
[
i
][
2
]
-
dz
*
v
[
i
][
1
]);
p
[
1
]
+=
massone
*
(
dz
*
v
[
i
][
0
]
-
dx
*
v
[
i
][
2
]);
p
[
2
]
+=
massone
*
(
dx
*
v
[
i
][
1
]
-
dy
*
v
[
i
][
0
]);
}
MPI_Allreduce
(
p
,
lmom
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
/* ----------------------------------------------------------------------
compute moment of inertia tensor around center-of-mass cm
must unwrap atoms to compute itensor correctly
------------------------------------------------------------------------- */
void
Group
::
inertia
(
int
igroup
,
double
*
cm
,
double
itensor
[
3
][
3
])
{
int
i
,
j
;
int
groupbit
=
bitmask
[
igroup
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
int
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
int
xbox
,
ybox
,
zbox
;
double
dx
,
dy
,
dz
,
massone
;
double
xprd
=
domain
->
xprd
;
double
yprd
=
domain
->
yprd
;
double
zprd
=
domain
->
zprd
;
double
ione
[
3
][
3
];
for
(
i
=
0
;
i
<
3
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
ione
[
i
][
j
]
=
0.0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
xbox
=
(
image
[
i
]
&
1023
)
-
512
;
ybox
=
(
image
[
i
]
>>
10
&
1023
)
-
512
;
zbox
=
(
image
[
i
]
>>
20
)
-
512
;
dx
=
(
x
[
i
][
0
]
+
xbox
*
xprd
)
-
cm
[
0
];
dy
=
(
x
[
i
][
1
]
+
ybox
*
yprd
)
-
cm
[
1
];
dz
=
(
x
[
i
][
2
]
+
zbox
*
zprd
)
-
cm
[
2
];
if
(
rmass
)
massone
=
rmass
[
i
];
else
massone
=
mass
[
type
[
i
]];
ione
[
0
][
0
]
+=
massone
*
(
dy
*
dy
+
dz
*
dz
);
ione
[
1
][
1
]
+=
massone
*
(
dx
*
dx
+
dz
*
dz
);
ione
[
2
][
2
]
+=
massone
*
(
dx
*
dx
+
dy
*
dy
);
ione
[
0
][
1
]
-=
massone
*
dx
*
dy
;
ione
[
1
][
2
]
-=
massone
*
dy
*
dz
;
ione
[
0
][
2
]
-=
massone
*
dx
*
dz
;
}
ione
[
1
][
0
]
=
ione
[
0
][
1
];
ione
[
2
][
1
]
=
ione
[
1
][
2
];
ione
[
2
][
0
]
=
ione
[
0
][
2
];
MPI_Allreduce
(
&
ione
[
0
][
0
],
&
itensor
[
0
][
0
],
9
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
/* ----------------------------------------------------------------------
compute angular velocity omega from L = Iw, inverting I to solve for w
really not a group operation, but L and I were computed for a group
------------------------------------------------------------------------- */
void
Group
::
omega
(
double
*
angmom
,
double
inertia
[
3
][
3
],
double
*
w
)
{
double
inverse
[
3
][
3
];
inverse
[
0
][
0
]
=
inertia
[
1
][
1
]
*
inertia
[
2
][
2
]
-
inertia
[
1
][
2
]
*
inertia
[
2
][
1
];
inverse
[
0
][
1
]
=
-
(
inertia
[
0
][
1
]
*
inertia
[
2
][
2
]
-
inertia
[
0
][
2
]
*
inertia
[
2
][
1
]);
inverse
[
0
][
2
]
=
inertia
[
0
][
1
]
*
inertia
[
1
][
2
]
-
inertia
[
0
][
2
]
*
inertia
[
1
][
1
];
inverse
[
1
][
0
]
=
-
(
inertia
[
1
][
0
]
*
inertia
[
2
][
2
]
-
inertia
[
1
][
2
]
*
inertia
[
2
][
0
]);
inverse
[
1
][
1
]
=
inertia
[
0
][
0
]
*
inertia
[
2
][
2
]
-
inertia
[
0
][
2
]
*
inertia
[
2
][
0
];
inverse
[
1
][
2
]
=
-
(
inertia
[
0
][
0
]
*
inertia
[
1
][
2
]
-
inertia
[
0
][
2
]
*
inertia
[
1
][
0
]);
inverse
[
2
][
0
]
=
inertia
[
1
][
0
]
*
inertia
[
2
][
1
]
-
inertia
[
1
][
1
]
*
inertia
[
2
][
0
];
inverse
[
2
][
1
]
=
-
(
inertia
[
0
][
0
]
*
inertia
[
2
][
1
]
-
inertia
[
0
][
1
]
*
inertia
[
2
][
0
]);
inverse
[
2
][
2
]
=
inertia
[
0
][
0
]
*
inertia
[
1
][
1
]
-
inertia
[
0
][
1
]
*
inertia
[
1
][
0
];
double
determinant
=
inertia
[
0
][
0
]
*
inertia
[
1
][
1
]
*
inertia
[
2
][
2
]
+
inertia
[
0
][
1
]
*
inertia
[
1
][
2
]
*
inertia
[
2
][
0
]
+
inertia
[
0
][
2
]
*
inertia
[
1
][
0
]
*
inertia
[
2
][
1
]
-
inertia
[
0
][
0
]
*
inertia
[
1
][
2
]
*
inertia
[
2
][
1
]
-
inertia
[
0
][
1
]
*
inertia
[
1
][
0
]
*
inertia
[
2
][
2
]
-
inertia
[
2
][
0
]
*
inertia
[
1
][
1
]
*
inertia
[
0
][
2
];
int
i
,
j
;
for
(
i
=
0
;
i
<
3
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
inverse
[
i
][
j
]
/=
determinant
;
w
[
0
]
=
inverse
[
0
][
0
]
*
angmom
[
0
]
+
inverse
[
0
][
1
]
*
angmom
[
1
]
+
inverse
[
0
][
2
]
*
angmom
[
2
];
w
[
1
]
=
inverse
[
1
][
0
]
*
angmom
[
0
]
+
inverse
[
1
][
1
]
*
angmom
[
1
]
+
inverse
[
1
][
2
]
*
angmom
[
2
];
w
[
2
]
=
inverse
[
2
][
0
]
*
angmom
[
0
]
+
inverse
[
2
][
1
]
*
angmom
[
1
]
+
inverse
[
2
][
2
]
*
angmom
[
2
];
}
Event Timeline
Log In to Comment