Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87797655
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
Tue, Oct 15, 00:49
Size
45 KB
Mime Type
text/x-c
Expires
Thu, Oct 17, 00:49 (2 d)
Engine
blob
Format
Raw Data
Handle
21656780
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 "mpi.h"
#include "math.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 "modify.h"
#include "fix.h"
#include "compute.h"
#include "output.h"
#include "input.h"
#include "variable.h"
#include "dump.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
/* ----------------------------------------------------------------------
initialize group memory
------------------------------------------------------------------------- */
Group
::
Group
(
LAMMPS
*
lmp
)
:
Pointers
(
lmp
)
{
MPI_Comm_rank
(
world
,
&
me
);
names
=
new
char
*
[
MAX_GROUP
];
bitmask
=
new
int
[
MAX_GROUP
];
inversemask
=
new
int
[
MAX_GROUP
];
dynamic
=
new
int
[
MAX_GROUP
];
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
names
[
i
]
=
NULL
;
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
;
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
dynamic
[
i
]
=
0
;
// create "all" group
char
*
str
=
(
char
*
)
"all"
;
int
n
=
strlen
(
str
)
+
1
;
names
[
0
]
=
new
char
[
n
];
strcpy
(
names
[
0
],
str
);
ngroup
=
1
;
}
/* ----------------------------------------------------------------------
free all memory
------------------------------------------------------------------------- */
Group
::~
Group
()
{
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
delete
[]
names
[
i
];
delete
[]
names
;
delete
[]
bitmask
;
delete
[]
inversemask
;
delete
[]
dynamic
;
}
/* ----------------------------------------------------------------------
assign atoms to a new or existing group
------------------------------------------------------------------------- */
void
Group
::
assign
(
int
narg
,
char
**
arg
)
{
int
i
;
if
(
domain
->
box_exist
==
0
)
error
->
all
(
FLERR
,
"Group command before simulation box is defined"
);
if
(
narg
<
2
)
error
->
all
(
FLERR
,
"Illegal group command"
);
// delete the group if not being used elsewhere
// clear mask of each atom assigned to this group
if
(
strcmp
(
arg
[
1
],
"delete"
)
==
0
)
{
int
igroup
=
find
(
arg
[
0
]);
if
(
igroup
==
-
1
)
error
->
all
(
FLERR
,
"Could not find group delete group ID"
);
if
(
igroup
==
0
)
error
->
all
(
FLERR
,
"Cannot delete group all"
);
for
(
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
modify
->
fix
[
i
]
->
igroup
==
igroup
)
error
->
all
(
FLERR
,
"Cannot delete group currently used by a fix"
);
for
(
i
=
0
;
i
<
modify
->
ncompute
;
i
++
)
if
(
modify
->
compute
[
i
]
->
igroup
==
igroup
)
error
->
all
(
FLERR
,
"Cannot delete group currently used by a compute"
);
for
(
i
=
0
;
i
<
output
->
ndump
;
i
++
)
if
(
output
->
dump
[
i
]
->
igroup
==
igroup
)
error
->
all
(
FLERR
,
"Cannot delete group currently used by a dump"
);
if
(
atom
->
firstgroupname
&&
strcmp
(
arg
[
0
],
atom
->
firstgroupname
)
==
0
)
error
->
all
(
FLERR
,
"Cannot delete group currently used by atom_modify first"
);
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
int
bits
=
inversemask
[
igroup
];
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
mask
[
i
]
&=
bits
;
if
(
dynamic
[
igroup
])
{
int
n
=
strlen
(
"GROUP_"
)
+
strlen
(
names
[
igroup
])
+
1
;
char
*
fixID
=
new
char
[
n
];
sprintf
(
fixID
,
"GROUP_%s"
,
names
[
igroup
]);
modify
->
delete_fix
(
fixID
);
delete
[]
fixID
;
}
delete
[]
names
[
igroup
];
names
[
igroup
]
=
NULL
;
dynamic
[
igroup
]
=
0
;
ngroup
--
;
return
;
}
// find group in existing list
// add a new group if igroup = -1
int
igroup
=
find
(
arg
[
0
]);
if
(
igroup
==
-
1
)
{
if
(
ngroup
==
MAX_GROUP
)
error
->
all
(
FLERR
,
"Too many groups"
);
igroup
=
find_unused
();
int
n
=
strlen
(
arg
[
0
])
+
1
;
names
[
igroup
]
=
new
char
[
n
];
strcpy
(
names
[
igroup
],
arg
[
0
]);
ngroup
++
;
}
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
(
FLERR
,
"Illegal group command"
);
int
iregion
=
domain
->
find_region
(
arg
[
2
]);
if
(
iregion
==
-
1
)
error
->
all
(
FLERR
,
"Group region ID does not exist"
);
domain
->
regions
[
iregion
]
->
init
();
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
mask
[
i
]
|=
bit
;
// style = type, molecule, id
// add to group if atom matches type/molecule/id or condition
}
else
if
(
strcmp
(
arg
[
1
],
"type"
)
==
0
||
strcmp
(
arg
[
1
],
"molecule"
)
==
0
||
strcmp
(
arg
[
1
],
"id"
)
==
0
)
{
if
(
narg
<
3
)
error
->
all
(
FLERR
,
"Illegal group command"
);
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
;
// args = logical condition
if
(
narg
>
3
&&
(
strcmp
(
arg
[
2
],
"<"
)
==
0
||
strcmp
(
arg
[
2
],
">"
)
==
0
||
strcmp
(
arg
[
2
],
"<="
)
==
0
||
strcmp
(
arg
[
2
],
">="
)
==
0
||
strcmp
(
arg
[
2
],
"<>"
)
==
0
))
{
int
condition
=
-
1
;
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
(
FLERR
,
"Illegal group command"
);
tagint
bound1
,
bound2
;
if
(
sizeof
(
tagint
)
==
sizeof
(
smallint
))
bound1
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
else
bound1
=
force
->
bnumeric
(
FLERR
,
arg
[
3
]);
bound2
=
-
1
;
if
(
condition
==
BETWEEN
)
{
if
(
narg
!=
5
)
error
->
all
(
FLERR
,
"Illegal group command"
);
if
(
sizeof
(
tagint
)
==
sizeof
(
smallint
))
bound2
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
else
bound2
=
force
->
bnumeric
(
FLERR
,
arg
[
4
]);
}
else
if
(
narg
!=
4
)
error
->
all
(
FLERR
,
"Illegal group command"
);
int
*
attribute
=
NULL
;
tagint
*
tattribute
=
NULL
;
if
(
category
==
TYPE
)
attribute
=
atom
->
type
;
else
if
(
category
==
MOLECULE
)
tattribute
=
atom
->
molecule
;
else
if
(
category
==
ID
)
tattribute
=
atom
->
tag
;
// add to group if meets condition
if
(
attribute
)
{
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
;
}
}
else
{
if
(
condition
==
LT
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
<
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
LE
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
<=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
GT
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
>
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
GE
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
>=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
EQ
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
==
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
NEQ
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
!=
bound1
)
mask
[
i
]
|=
bit
;
}
else
if
(
condition
==
BETWEEN
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
>=
bound1
&&
tattribute
[
i
]
<=
bound2
)
mask
[
i
]
|=
bit
;
}
}
// args = list of values
}
else
{
int
*
attribute
=
NULL
;
tagint
*
tattribute
=
NULL
;
if
(
category
==
TYPE
)
attribute
=
atom
->
type
;
else
if
(
category
==
MOLECULE
)
tattribute
=
atom
->
molecule
;
else
if
(
category
==
ID
)
tattribute
=
atom
->
tag
;
char
*
ptr
;
tagint
start
,
stop
,
delta
;
for
(
int
iarg
=
2
;
iarg
<
narg
;
iarg
++
)
{
if
(
strchr
(
arg
[
iarg
],
':'
))
{
start
=
ATOTAGINT
(
strtok
(
arg
[
iarg
],
":"
));
stop
=
ATOTAGINT
(
strtok
(
NULL
,
":"
));
ptr
=
strtok
(
NULL
,
":"
);
if
(
ptr
)
delta
=
ATOTAGINT
(
ptr
);
else
delta
=
1
;
}
else
{
start
=
stop
=
ATOTAGINT
(
arg
[
iarg
]);
delta
=
1
;
}
// add to group if attribute matches value or sequence
if
(
attribute
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
attribute
[
i
]
>=
start
&&
attribute
[
i
]
<=
stop
&&
(
attribute
[
i
]
-
start
)
%
delta
==
0
)
mask
[
i
]
|=
bit
;
}
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
tattribute
[
i
]
>=
start
&&
tattribute
[
i
]
<=
stop
&&
(
tattribute
[
i
]
-
start
)
%
delta
==
0
)
mask
[
i
]
|=
bit
;
}
}
}
// style = variable
// add to group if atom-atyle variable is non-zero
}
else
if
(
strcmp
(
arg
[
1
],
"variable"
)
==
0
)
{
int
ivar
=
input
->
variable
->
find
(
arg
[
2
]);
if
(
ivar
<
0
)
error
->
all
(
FLERR
,
"Variable name for group does not exist"
);
if
(
!
input
->
variable
->
atomstyle
(
ivar
))
error
->
all
(
FLERR
,
"Variable for group is invalid style"
);
double
*
aflag
;
// aflag = evaluation of per-atom variable
memory
->
create
(
aflag
,
nlocal
,
"group:aflag"
);
input
->
variable
->
compute_atom
(
ivar
,
0
,
aflag
,
1
,
0
);
// add to group if per-atom variable evaluated to non-zero
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
aflag
[
i
]
!=
0.0
)
mask
[
i
]
|=
bit
;
memory
->
destroy
(
aflag
);
// style = subtract
}
else
if
(
strcmp
(
arg
[
1
],
"subtract"
)
==
0
)
{
if
(
narg
<
4
)
error
->
all
(
FLERR
,
"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
(
FLERR
,
"Group ID does not exist"
);
if
(
dynamic
[
jgroup
])
error
->
all
(
FLERR
,
"Cannot subtract groups using a dynamic group"
);
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
(
FLERR
,
"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
(
FLERR
,
"Group ID does not exist"
);
if
(
dynamic
[
jgroup
])
error
->
all
(
FLERR
,
"Cannot union groups using a dynamic group"
);
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
(
FLERR
,
"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
(
FLERR
,
"Group ID does not exist"
);
if
(
dynamic
[
jgroup
])
error
->
all
(
FLERR
,
"Cannot intersect groups using a dynamic group"
);
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
;
// style = dynamic
// create a new FixGroup to dynamically determine atoms in group
}
else
if
(
strcmp
(
arg
[
1
],
"dynamic"
)
==
0
)
{
if
(
narg
<
4
)
error
->
all
(
FLERR
,
"Illegal group command"
);
if
(
strcmp
(
arg
[
0
],
arg
[
2
])
==
0
)
error
->
all
(
FLERR
,
"Group dynamic cannot reference itself"
);
if
(
find
(
arg
[
2
])
<
0
)
error
->
all
(
FLERR
,
"Group dynamic parent group does not exist"
);
if
(
igroup
==
0
)
error
->
all
(
FLERR
,
"Group all cannot be made dynamic"
);
// if group is already dynamic, delete existing FixGroup
if
(
dynamic
[
igroup
])
{
int
n
=
strlen
(
"GROUP_"
)
+
strlen
(
names
[
igroup
])
+
1
;
char
*
fixID
=
new
char
[
n
];
sprintf
(
fixID
,
"GROUP_%s"
,
names
[
igroup
]);
modify
->
delete_fix
(
fixID
);
delete
[]
fixID
;
}
dynamic
[
igroup
]
=
1
;
int
n
=
strlen
(
"GROUP_"
)
+
strlen
(
names
[
igroup
])
+
1
;
char
*
fixID
=
new
char
[
n
];
sprintf
(
fixID
,
"GROUP_%s"
,
names
[
igroup
]);
char
**
newarg
=
new
char
*
[
narg
];
newarg
[
0
]
=
fixID
;
newarg
[
1
]
=
arg
[
2
];
newarg
[
2
]
=
(
char
*
)
"GROUP"
;
for
(
int
i
=
3
;
i
<
narg
;
i
++
)
newarg
[
i
]
=
arg
[
i
];
modify
->
add_fix
(
narg
,
newarg
);
delete
[]
newarg
;
delete
[]
fixID
;
// style = static
// remove dynamic FixGroup if necessary
}
else
if
(
strcmp
(
arg
[
1
],
"static"
)
==
0
)
{
if
(
narg
!=
2
)
error
->
all
(
FLERR
,
"Illegal group command"
);
if
(
dynamic
[
igroup
])
{
int
n
=
strlen
(
"GROUP_"
)
+
strlen
(
names
[
igroup
])
+
1
;
char
*
fixID
=
new
char
[
n
];
sprintf
(
fixID
,
"GROUP_%s"
,
names
[
igroup
]);
modify
->
delete_fix
(
fixID
);
delete
[]
fixID
;
}
dynamic
[
igroup
]
=
0
;
// not a valid group style
}
else
error
->
all
(
FLERR
,
"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
(
dynamic
[
igroup
])
{
if
(
screen
)
fprintf
(
screen
,
"dynamic group %s
\n
"
,
names
[
igroup
]);
if
(
logfile
)
fprintf
(
logfile
,
"dynamic group %s
\n
"
,
names
[
igroup
]);
}
else
{
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
// add a new group if igroup = -1
int
igroup
=
find
(
name
);
if
(
igroup
==
-
1
)
{
if
(
ngroup
==
MAX_GROUP
)
error
->
all
(
FLERR
,
"Too many groups"
);
igroup
=
find_unused
();
int
n
=
strlen
(
name
)
+
1
;
names
[
igroup
]
=
new
char
[
n
];
strcpy
(
names
[
igroup
],
name
);
ngroup
++
;
}
// add atoms to group whose flags are 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
<
MAX_GROUP
;
igroup
++
)
if
(
names
[
igroup
]
&&
strcmp
(
name
,
names
[
igroup
])
==
0
)
return
igroup
;
return
-
1
;
}
/* ----------------------------------------------------------------------
return index of first available group
should never be called when group limit has been reached
------------------------------------------------------------------------- */
int
Group
::
find_unused
()
{
for
(
int
igroup
=
0
;
igroup
<
MAX_GROUP
;
igroup
++
)
if
(
names
[
igroup
]
==
NULL
)
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
);
// use count to not change restart format with deleted groups
// remove this on next major release
int
n
;
int
count
=
0
;
for
(
int
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
{
if
(
names
[
i
])
n
=
strlen
(
names
[
i
])
+
1
;
else
n
=
0
;
fwrite
(
&
n
,
sizeof
(
int
),
1
,
fp
);
if
(
n
)
{
fwrite
(
names
[
i
],
sizeof
(
char
),
n
,
fp
);
count
++
;
}
if
(
count
==
ngroup
)
break
;
}
}
/* ----------------------------------------------------------------------
read group info from a restart file
proc 0 reads, bcast to all procs
------------------------------------------------------------------------- */
void
Group
::
read_restart
(
FILE
*
fp
)
{
int
i
,
n
;
// delete existing group names
// atom masks will be overwritten by reading of restart file
for
(
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
delete
[]
names
[
i
];
if
(
me
==
0
)
fread
(
&
ngroup
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
ngroup
,
1
,
MPI_INT
,
0
,
world
);
// use count to not change restart format with deleted groups
// remove this on next major release
int
count
=
0
;
for
(
i
=
0
;
i
<
MAX_GROUP
;
i
++
)
{
if
(
count
==
ngroup
)
{
names
[
i
]
=
NULL
;
continue
;
}
if
(
me
==
0
)
fread
(
&
n
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
if
(
n
)
{
names
[
i
]
=
new
char
[
n
];
if
(
me
==
0
)
fread
(
names
[
i
],
sizeof
(
char
),
n
,
fp
);
MPI_Bcast
(
names
[
i
],
n
,
MPI_CHAR
,
0
,
world
);
count
++
;
}
else
names
[
i
]
=
NULL
;
}
}
// ----------------------------------------------------------------------
// computations on a group of atoms
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
count atoms in group
------------------------------------------------------------------------- */
bigint
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
++
;
bigint
nsingle
=
n
;
bigint
nall
;
MPI_Allreduce
(
&
nsingle
,
&
nall
,
1
,
MPI_LMP_BIGINT
,
MPI_SUM
,
world
);
return
nall
;
}
/* ----------------------------------------------------------------------
count atoms in group and region
------------------------------------------------------------------------- */
bigint
Group
::
count
(
int
igroup
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
n
++
;
bigint
nsingle
=
n
;
bigint
nall
;
MPI_Allreduce
(
&
nsingle
,
&
nall
,
1
,
MPI_LMP_BIGINT
,
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
];
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
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 mass of group of atoms in region
use either per-type mass or per-atom rmass
------------------------------------------------------------------------- */
double
Group
::
mass
(
int
igroup
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
double
one
=
0.0
;
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
one
+=
rmass
[
i
];
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
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
];
double
*
q
=
atom
->
q
;
int
*
mask
=
atom
->
mask
;
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 total charge of group of atoms in region
------------------------------------------------------------------------- */
double
Group
::
charge
(
int
igroup
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
double
*
q
=
atom
->
q
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
double
qone
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
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 coordinate bounds of the group of atoms in region
periodic images are not considered, so atoms are NOT unwrapped
------------------------------------------------------------------------- */
void
Group
::
bounds
(
int
igroup
,
double
*
minmax
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
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 of atoms
masstotal = total mass
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
;
imageint
*
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
massone
;
double
unwrap
[
3
];
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
massone
=
rmass
[
i
];
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
cmone
[
0
]
+=
unwrap
[
0
]
*
massone
;
cmone
[
1
]
+=
unwrap
[
1
]
*
massone
;
cmone
[
2
]
+=
unwrap
[
2
]
*
massone
;
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
massone
=
mass
[
type
[
i
]];
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
cmone
[
0
]
+=
unwrap
[
0
]
*
massone
;
cmone
[
1
]
+=
unwrap
[
1
]
*
massone
;
cmone
[
2
]
+=
unwrap
[
2
]
*
massone
;
}
}
MPI_Allreduce
(
cmone
,
cm
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
masstotal
>
0.0
)
{
cm
[
0
]
/=
masstotal
;
cm
[
1
]
/=
masstotal
;
cm
[
2
]
/=
masstotal
;
}
}
/* ----------------------------------------------------------------------
compute the center-of-mass coords of group of atoms in region
mastotal = total mass
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
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
imageint
*
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
massone
;
double
unwrap
[
3
];
if
(
rmass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
massone
=
rmass
[
i
];
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
cmone
[
0
]
+=
unwrap
[
0
]
*
massone
;
cmone
[
1
]
+=
unwrap
[
1
]
*
massone
;
cmone
[
2
]
+=
unwrap
[
2
]
*
massone
;
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
massone
=
mass
[
type
[
i
]];
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
cmone
[
0
]
+=
unwrap
[
0
]
*
massone
;
cmone
[
1
]
+=
unwrap
[
1
]
*
massone
;
cmone
[
2
]
+=
unwrap
[
2
]
*
massone
;
}
}
MPI_Allreduce
(
cmone
,
cm
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
masstotal
>
0.0
)
{
cm
[
0
]
/=
masstotal
;
cm
[
1
]
/=
masstotal
;
cm
[
2
]
/=
masstotal
;
}
}
/* ----------------------------------------------------------------------
compute the center-of-mass velocity of group of atoms
masstotal = total mass
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
);
if
(
masstotal
>
0.0
)
{
cm
[
0
]
/=
masstotal
;
cm
[
1
]
/=
masstotal
;
cm
[
2
]
/=
masstotal
;
}
}
/* ----------------------------------------------------------------------
compute the center-of-mass velocity of group of atoms in region
masstotal = total mass
return center-of-mass velocity in cm[]
------------------------------------------------------------------------- */
void
Group
::
vcm
(
int
igroup
,
double
masstotal
,
double
*
cm
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
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
);
if
(
masstotal
>
0.0
)
{
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 force on group of atoms in region
------------------------------------------------------------------------- */
void
Group
::
fcm
(
int
igroup
,
double
*
cm
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
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 of atoms 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 total kinetic energy of group of atoms in region and return it
------------------------------------------------------------------------- */
double
Group
::
ke
(
int
igroup
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
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 of atoms
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
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
double
rg
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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
);
if
(
masstotal
>
0.0
)
return
sqrt
(
rg_all
/
masstotal
);
return
0.0
;
}
/* ----------------------------------------------------------------------
compute the radius-of-gyration of group of atoms in region
around center-of-mass cm
must unwrap atoms to compute Rg correctly
------------------------------------------------------------------------- */
double
Group
::
gyration
(
int
igroup
,
double
masstotal
,
double
*
cm
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
double
rg
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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
);
if
(
masstotal
>
0.0
)
return
sqrt
(
rg_all
/
masstotal
);
return
0.0
;
}
/* ----------------------------------------------------------------------
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
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
double
p
[
3
];
p
[
0
]
=
p
[
1
]
=
p
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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 the angular momentum L (lmom) of group of atoms in region
around center-of-mass cm
must unwrap atoms to compute L correctly
------------------------------------------------------------------------- */
void
Group
::
angmom
(
int
igroup
,
double
*
cm
,
double
*
lmom
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
double
**
v
=
atom
->
v
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
double
p
[
3
];
p
[
0
]
=
p
[
1
]
=
p
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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 the torque T (tq) on group
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void
Group
::
torque
(
int
igroup
,
double
*
cm
,
double
*
tq
)
{
int
groupbit
=
bitmask
[
igroup
];
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
int
*
mask
=
atom
->
mask
;
imageint
*
image
=
atom
->
image
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
;
double
unwrap
[
3
];
double
tlocal
[
3
];
tlocal
[
0
]
=
tlocal
[
1
]
=
tlocal
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
cm
[
2
];
tlocal
[
0
]
+=
dy
*
f
[
i
][
2
]
-
dz
*
f
[
i
][
1
];
tlocal
[
1
]
+=
dz
*
f
[
i
][
0
]
-
dx
*
f
[
i
][
2
];
tlocal
[
2
]
+=
dx
*
f
[
i
][
1
]
-
dy
*
f
[
i
][
0
];
}
MPI_Allreduce
(
tlocal
,
tq
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
/* ----------------------------------------------------------------------
compute the torque T (tq) on group of atoms in region
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void
Group
::
torque
(
int
igroup
,
double
*
cm
,
double
*
tq
,
int
iregion
)
{
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
int
*
mask
=
atom
->
mask
;
imageint
*
image
=
atom
->
image
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
;
double
unwrap
[
3
];
double
tlocal
[
3
];
tlocal
[
0
]
=
tlocal
[
1
]
=
tlocal
[
2
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
cm
[
2
];
tlocal
[
0
]
+=
dy
*
f
[
i
][
2
]
-
dz
*
f
[
i
][
1
];
tlocal
[
1
]
+=
dz
*
f
[
i
][
0
]
-
dx
*
f
[
i
][
2
];
tlocal
[
2
]
+=
dx
*
f
[
i
][
1
]
-
dy
*
f
[
i
][
0
];
}
MPI_Allreduce
(
tlocal
,
tq
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
/* ----------------------------------------------------------------------
compute moment of inertia tensor around center-of-mass cm of group
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
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
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
)
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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 moment of inertia tensor around cm of group of atoms in region
must unwrap atoms to compute itensor correctly
------------------------------------------------------------------------- */
void
Group
::
inertia
(
int
igroup
,
double
*
cm
,
double
itensor
[
3
][
3
],
int
iregion
)
{
int
i
,
j
;
int
groupbit
=
bitmask
[
igroup
];
Region
*
region
=
domain
->
regions
[
iregion
];
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
*
type
=
atom
->
type
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
dx
,
dy
,
dz
,
massone
;
double
unwrap
[
3
];
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
&&
region
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
]))
{
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
dx
=
unwrap
[
0
]
-
cm
[
0
];
dy
=
unwrap
[
1
]
-
cm
[
1
];
dz
=
unwrap
[
2
]
-
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
];
if
(
determinant
>
0.0
)
for
(
int
i
=
0
;
i
<
3
;
i
++
)
for
(
int
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