Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90938804
fix_atom_swap.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, Nov 6, 05:24
Size
23 KB
Mime Type
text/x-c
Expires
Fri, Nov 8, 05:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
22163592
Attached To
rLAMMPS lammps
fix_atom_swap.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Paul Crozier (SNL)
Alexander Stukowski
------------------------------------------------------------------------- */
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include "fix_atom_swap.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_hybrid.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "compute.h"
#include "group.h"
#include "domain.h"
#include "region.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "thermo.h"
#include "output.h"
#include "neighbor.h"
#include <iostream>
using
namespace
std
;
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
using
namespace
MathConst
;
/* ---------------------------------------------------------------------- */
FixAtomSwap
::
FixAtomSwap
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
),
idregion
(
NULL
),
type_list
(
NULL
),
mu
(
NULL
),
qtype
(
NULL
),
sqrt_mass_ratio
(
NULL
),
random_equal
(
NULL
),
random_unequal
(
NULL
)
{
if
(
narg
<
10
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
dynamic_group_allow
=
1
;
vector_flag
=
1
;
size_vector
=
2
;
global_freq
=
1
;
extvector
=
0
;
restart_global
=
1
;
time_depend
=
1
;
// required args
nevery
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
ncycles
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
seed
=
force
->
inumeric
(
FLERR
,
arg
[
5
]);
double
temperature
=
force
->
numeric
(
FLERR
,
arg
[
6
]);
beta
=
1.0
/
(
force
->
boltz
*
temperature
);
if
(
nevery
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
if
(
ncycles
<
0
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
if
(
seed
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
memory
->
create
(
type_list
,
atom
->
ntypes
,
"atom/swap:type_list"
);
memory
->
create
(
mu
,
atom
->
ntypes
+
1
,
"atom/swap:mu"
);
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
mu
[
i
]
=
0.0
;
// read options from end of input line
options
(
narg
-
7
,
&
arg
[
7
]);
// random number generator, same for all procs
random_equal
=
new
RanPark
(
lmp
,
seed
);
// random number generator, not the same for all procs
random_unequal
=
new
RanPark
(
lmp
,
seed
);
// set up reneighboring
force_reneighbor
=
1
;
next_reneighbor
=
update
->
ntimestep
+
1
;
// zero out counters
nswap_attempts
=
0.0
;
nswap_successes
=
0.0
;
atom_swap_nmax
=
0
;
local_swap_atom_list
=
NULL
;
local_swap_iatom_list
=
NULL
;
local_swap_jatom_list
=
NULL
;
// set comm size needed by this Fix
if
(
atom
->
q_flag
)
comm_forward
=
2
;
else
comm_forward
=
1
;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void
FixAtomSwap
::
options
(
int
narg
,
char
**
arg
)
{
if
(
narg
<
0
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
regionflag
=
0
;
conserve_ke_flag
=
1
;
semi_grand_flag
=
0
;
nswaptypes
=
0
;
nmutypes
=
0
;
iregion
=
-
1
;
int
iarg
=
0
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"region"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
iregion
=
domain
->
find_region
(
arg
[
iarg
+
1
]);
if
(
iregion
==
-
1
)
error
->
all
(
FLERR
,
"Region ID for fix atom/swap does not exist"
);
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
idregion
=
new
char
[
n
];
strcpy
(
idregion
,
arg
[
iarg
+
1
]);
regionflag
=
1
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"ke"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"no"
)
==
0
)
conserve_ke_flag
=
0
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"yes"
)
==
0
)
conserve_ke_flag
=
1
;
else
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"semi-grand"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"no"
)
==
0
)
semi_grand_flag
=
0
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"yes"
)
==
0
)
semi_grand_flag
=
1
;
else
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"types"
)
==
0
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
iarg
++
;
while
(
iarg
<
narg
)
{
if
(
isalpha
(
arg
[
iarg
][
0
]))
break
;
if
(
nswaptypes
>=
atom
->
ntypes
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
type_list
[
nswaptypes
]
=
force
->
numeric
(
FLERR
,
arg
[
iarg
]);
nswaptypes
++
;
iarg
++
;
}
}
else
if
(
strcmp
(
arg
[
iarg
],
"mu"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
iarg
++
;
while
(
iarg
<
narg
)
{
if
(
isalpha
(
arg
[
iarg
][
0
]))
break
;
nmutypes
++
;
if
(
nmutypes
>
atom
->
ntypes
)
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
mu
[
nmutypes
]
=
force
->
numeric
(
FLERR
,
arg
[
iarg
]);
iarg
++
;
}
}
else
error
->
all
(
FLERR
,
"Illegal fix atom/swap command"
);
}
}
/* ---------------------------------------------------------------------- */
FixAtomSwap
::~
FixAtomSwap
()
{
memory
->
destroy
(
type_list
);
memory
->
destroy
(
mu
);
memory
->
destroy
(
qtype
);
memory
->
destroy
(
sqrt_mass_ratio
);
if
(
regionflag
)
delete
[]
idregion
;
delete
random_equal
;
delete
random_unequal
;
}
/* ---------------------------------------------------------------------- */
int
FixAtomSwap
::
setmask
()
{
int
mask
=
0
;
mask
|=
PRE_EXCHANGE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixAtomSwap
::
init
()
{
char
*
id_pe
=
(
char
*
)
"thermo_pe"
;
int
ipe
=
modify
->
find_compute
(
id_pe
);
c_pe
=
modify
->
compute
[
ipe
];
int
*
type
=
atom
->
type
;
if
(
nswaptypes
<
2
)
error
->
all
(
FLERR
,
"Must specify at least 2 types in fix atom/swap command"
);
if
(
semi_grand_flag
)
{
if
(
nswaptypes
!=
nmutypes
)
error
->
all
(
FLERR
,
"Need nswaptypes mu values in fix atom/swap command"
);
}
else
{
if
(
nswaptypes
!=
2
)
error
->
all
(
FLERR
,
"Only 2 types allowed when not using semi-grand in fix atom/swap command"
);
if
(
nmutypes
!=
0
)
error
->
all
(
FLERR
,
"Mu not allowed when not using semi-grand in fix atom/swap command"
);
}
for
(
int
iswaptype
=
0
;
iswaptype
<
nswaptypes
;
iswaptype
++
)
if
(
type_list
[
iswaptype
]
<=
0
||
type_list
[
iswaptype
]
>
atom
->
ntypes
)
error
->
all
(
FLERR
,
"Invalid atom type in fix atom/swap command"
);
// this is only required for non-semi-grand
// in which case, nswaptypes = 2
if
(
atom
->
q_flag
&&
!
semi_grand_flag
)
{
double
qmax
,
qmin
;
int
firstall
,
first
;
memory
->
create
(
qtype
,
nswaptypes
,
"atom/swap:qtype"
);
for
(
int
iswaptype
=
0
;
iswaptype
<
nswaptypes
;
iswaptype
++
)
{
first
=
1
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
{
if
(
atom
->
mask
[
i
]
&
groupbit
)
{
if
(
type
[
i
]
==
type_list
[
iswaptype
])
{
if
(
first
)
{
qtype
[
iswaptype
]
=
atom
->
q
[
i
];
first
=
0
;
}
else
if
(
qtype
[
iswaptype
]
!=
atom
->
q
[
i
])
error
->
one
(
FLERR
,
"All atoms of a swapped type must have the same charge."
);
}
}
}
MPI_Allreduce
(
&
first
,
&
firstall
,
1
,
MPI_INT
,
MPI_MIN
,
world
);
if
(
firstall
)
error
->
all
(
FLERR
,
"At least one atom of each swapped type must be present to define charges."
);
if
(
first
)
qtype
[
iswaptype
]
=
-
DBL_MAX
;
MPI_Allreduce
(
&
qtype
[
iswaptype
],
&
qmax
,
1
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
if
(
first
)
qtype
[
iswaptype
]
=
DBL_MAX
;
MPI_Allreduce
(
&
qtype
[
iswaptype
],
&
qmin
,
1
,
MPI_DOUBLE
,
MPI_MIN
,
world
);
if
(
qmax
!=
qmin
)
error
->
all
(
FLERR
,
"All atoms of a swapped type must have same charge."
);
}
}
memory
->
create
(
sqrt_mass_ratio
,
atom
->
ntypes
+
1
,
atom
->
ntypes
+
1
,
"atom/swap:sqrt_mass_ratio"
);
for
(
int
itype
=
1
;
itype
<=
atom
->
ntypes
;
itype
++
)
for
(
int
jtype
=
1
;
jtype
<=
atom
->
ntypes
;
jtype
++
)
sqrt_mass_ratio
[
itype
][
jtype
]
=
sqrt
(
atom
->
mass
[
itype
]
/
atom
->
mass
[
jtype
]);
// check to see if itype and jtype cutoffs are the same
// if not, reneighboring will be needed between swaps
double
**
cutsq
=
force
->
pair
->
cutsq
;
unequal_cutoffs
=
false
;
for
(
int
iswaptype
=
0
;
iswaptype
<
nswaptypes
;
iswaptype
++
)
for
(
int
jswaptype
=
0
;
jswaptype
<
nswaptypes
;
jswaptype
++
)
for
(
int
ktype
=
1
;
ktype
<=
atom
->
ntypes
;
ktype
++
)
if
(
cutsq
[
type_list
[
iswaptype
]][
ktype
]
!=
cutsq
[
type_list
[
jswaptype
]][
ktype
])
unequal_cutoffs
=
true
;
// check that no swappable atoms are in atom->firstgroup
// swapping such an atom might not leave firstgroup atoms first
if
(
atom
->
firstgroup
>=
0
)
{
int
*
mask
=
atom
->
mask
;
int
firstgroupbit
=
group
->
bitmask
[
atom
->
firstgroup
];
int
flag
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
if
((
mask
[
i
]
==
groupbit
)
&&
(
mask
[
i
]
&&
firstgroupbit
))
flag
=
1
;
int
flagall
;
MPI_Allreduce
(
&
flag
,
&
flagall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flagall
)
error
->
all
(
FLERR
,
"Cannot do atom/swap on atoms in atom_modify first group"
);
}
}
/* ----------------------------------------------------------------------
attempt Monte Carlo swaps
------------------------------------------------------------------------- */
void
FixAtomSwap
::
pre_exchange
()
{
// just return if should not be called on this timestep
if
(
next_reneighbor
!=
update
->
ntimestep
)
return
;
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
energy_stored
=
energy_full
();
int
nsuccess
=
0
;
if
(
semi_grand_flag
)
{
update_semi_grand_atoms_list
();
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
nsuccess
+=
attempt_semi_grand
();
}
else
{
update_swap_atoms_list
();
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
nsuccess
+=
attempt_swap
();
}
nswap_attempts
+=
ncycles
;
nswap_successes
+=
nsuccess
;
energy_full
();
next_reneighbor
=
update
->
ntimestep
+
nevery
;
}
/* ----------------------------------------------------------------------
Note: atom charges are assumed equal and so are not updated
------------------------------------------------------------------------- */
int
FixAtomSwap
::
attempt_semi_grand
()
{
if
(
nswap
==
0
)
return
0
;
double
energy_before
=
energy_stored
;
int
itype
,
jtype
,
jswaptype
;
int
i
=
pick_semi_grand_atom
();
if
(
i
>=
0
)
{
jswaptype
=
static_cast
<
int
>
(
nswaptypes
*
random_unequal
->
uniform
());
jtype
=
type_list
[
jswaptype
];
itype
=
atom
->
type
[
i
];
while
(
itype
==
jtype
)
{
jswaptype
=
static_cast
<
int
>
(
nswaptypes
*
random_unequal
->
uniform
());
jtype
=
type_list
[
jswaptype
];
}
atom
->
type
[
i
]
=
jtype
;
}
if
(
unequal_cutoffs
)
{
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
}
else
{
comm
->
forward_comm_fix
(
this
);
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
double
energy_after
=
energy_full
();
int
success
=
0
;
if
(
i
>=
0
)
if
(
random_unequal
->
uniform
()
<
exp
(
-
beta
*
(
energy_after
-
energy_before
+
mu
[
jtype
]
-
mu
[
itype
])))
success
=
1
;
int
success_all
=
0
;
MPI_Allreduce
(
&
success
,
&
success_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
if
(
success_all
)
{
update_semi_grand_atoms_list
();
energy_stored
=
energy_after
;
if
(
conserve_ke_flag
)
{
if
(
i
>=
0
)
{
atom
->
v
[
i
][
0
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
atom
->
v
[
i
][
1
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
atom
->
v
[
i
][
2
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
}
}
return
1
;
}
else
{
if
(
i
>=
0
)
{
atom
->
type
[
i
]
=
itype
;
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
energy_stored
=
energy_before
;
if
(
unequal_cutoffs
)
{
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
}
else
{
comm
->
forward_comm_fix
(
this
);
}
}
return
0
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int
FixAtomSwap
::
attempt_swap
()
{
if
((
niswap
==
0
)
||
(
njswap
==
0
))
return
0
;
double
energy_before
=
energy_stored
;
int
i
=
pick_i_swap_atom
();
int
j
=
pick_j_swap_atom
();
int
itype
=
type_list
[
0
];
int
jtype
=
type_list
[
1
];
if
(
i
>=
0
)
{
atom
->
type
[
i
]
=
jtype
;
if
(
atom
->
q_flag
)
atom
->
q
[
i
]
=
qtype
[
1
];
}
if
(
j
>=
0
)
{
atom
->
type
[
j
]
=
itype
;
if
(
atom
->
q_flag
)
atom
->
q
[
j
]
=
qtype
[
0
];
}
if
(
unequal_cutoffs
)
{
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
}
else
{
comm
->
forward_comm_fix
(
this
);
}
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before
-
energy_after
)))
{
update_swap_atoms_list
();
energy_stored
=
energy_after
;
if
(
conserve_ke_flag
)
{
if
(
i
>=
0
)
{
atom
->
v
[
i
][
0
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
atom
->
v
[
i
][
1
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
atom
->
v
[
i
][
2
]
*=
sqrt_mass_ratio
[
itype
][
jtype
];
}
if
(
j
>=
0
)
{
atom
->
v
[
j
][
0
]
*=
sqrt_mass_ratio
[
jtype
][
itype
];
atom
->
v
[
j
][
1
]
*=
sqrt_mass_ratio
[
jtype
][
itype
];
atom
->
v
[
j
][
2
]
*=
sqrt_mass_ratio
[
jtype
][
itype
];
}
}
return
1
;
}
else
{
if
(
i
>=
0
)
{
atom
->
type
[
i
]
=
type_list
[
0
];
if
(
atom
->
q_flag
)
atom
->
q
[
i
]
=
qtype
[
0
];
}
if
(
j
>=
0
)
{
atom
->
type
[
j
]
=
type_list
[
1
];
if
(
atom
->
q_flag
)
atom
->
q
[
j
]
=
qtype
[
1
];
}
energy_stored
=
energy_before
;
if
(
unequal_cutoffs
)
{
if
(
domain
->
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
comm
->
borders
();
if
(
domain
->
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
}
else
{
comm
->
forward_comm_fix
(
this
);
}
}
return
0
;
}
/* ----------------------------------------------------------------------
compute system potential energy
------------------------------------------------------------------------- */
double
FixAtomSwap
::
energy_full
()
{
int
eflag
=
1
;
int
vflag
=
0
;
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
if
(
modify
->
n_pre_force
)
modify
->
pre_force
(
vflag
);
if
(
force
->
pair
)
force
->
pair
->
compute
(
eflag
,
vflag
);
if
(
atom
->
molecular
)
{
if
(
force
->
bond
)
force
->
bond
->
compute
(
eflag
,
vflag
);
if
(
force
->
angle
)
force
->
angle
->
compute
(
eflag
,
vflag
);
if
(
force
->
dihedral
)
force
->
dihedral
->
compute
(
eflag
,
vflag
);
if
(
force
->
improper
)
force
->
improper
->
compute
(
eflag
,
vflag
);
}
if
(
force
->
kspace
)
force
->
kspace
->
compute
(
eflag
,
vflag
);
if
(
modify
->
n_post_force
)
modify
->
post_force
(
vflag
);
if
(
modify
->
n_end_of_step
)
modify
->
end_of_step
();
update
->
eflag_global
=
update
->
ntimestep
;
double
total_energy
=
c_pe
->
compute_scalar
();
return
total_energy
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int
FixAtomSwap
::
pick_semi_grand_atom
()
{
int
i
=
-
1
;
int
iwhichglobal
=
static_cast
<
int
>
(
nswap
*
random_equal
->
uniform
());
if
((
iwhichglobal
>=
nswap_before
)
&&
(
iwhichglobal
<
nswap_before
+
nswap_local
))
{
int
iwhichlocal
=
iwhichglobal
-
nswap_before
;
i
=
local_swap_atom_list
[
iwhichlocal
];
}
return
i
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int
FixAtomSwap
::
pick_i_swap_atom
()
{
int
i
=
-
1
;
int
iwhichglobal
=
static_cast
<
int
>
(
niswap
*
random_equal
->
uniform
());
if
((
iwhichglobal
>=
niswap_before
)
&&
(
iwhichglobal
<
niswap_before
+
niswap_local
))
{
int
iwhichlocal
=
iwhichglobal
-
niswap_before
;
i
=
local_swap_iatom_list
[
iwhichlocal
];
}
return
i
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int
FixAtomSwap
::
pick_j_swap_atom
()
{
int
j
=
-
1
;
int
jwhichglobal
=
static_cast
<
int
>
(
njswap
*
random_equal
->
uniform
());
if
((
jwhichglobal
>=
njswap_before
)
&&
(
jwhichglobal
<
njswap_before
+
njswap_local
))
{
int
jwhichlocal
=
jwhichglobal
-
njswap_before
;
j
=
local_swap_jatom_list
[
jwhichlocal
];
}
return
j
;
}
/* ----------------------------------------------------------------------
update the list of gas atoms
------------------------------------------------------------------------- */
void
FixAtomSwap
::
update_semi_grand_atoms_list
()
{
int
nlocal
=
atom
->
nlocal
;
double
**
x
=
atom
->
x
;
if
(
atom
->
nmax
>
atom_swap_nmax
)
{
memory
->
sfree
(
local_swap_atom_list
);
atom_swap_nmax
=
atom
->
nmax
;
local_swap_atom_list
=
(
int
*
)
memory
->
smalloc
(
atom_swap_nmax
*
sizeof
(
int
),
"MCSWAP:local_swap_atom_list"
);
}
nswap_local
=
0
;
if
(
regionflag
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
])
==
1
)
{
if
(
atom
->
mask
[
i
]
&
groupbit
)
{
int
itype
=
atom
->
type
[
i
];
int
iswaptype
;
for
(
iswaptype
=
0
;
iswaptype
<
nswaptypes
;
iswaptype
++
)
if
(
itype
==
type_list
[
iswaptype
])
break
;
if
(
iswaptype
==
nswaptypes
)
continue
;
local_swap_atom_list
[
nswap_local
]
=
i
;
nswap_local
++
;
}
}
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
mask
[
i
]
&
groupbit
)
{
int
itype
=
atom
->
type
[
i
];
int
iswaptype
;
for
(
iswaptype
=
0
;
iswaptype
<
nswaptypes
;
iswaptype
++
)
if
(
itype
==
type_list
[
iswaptype
])
break
;
if
(
iswaptype
==
nswaptypes
)
continue
;
local_swap_atom_list
[
nswap_local
]
=
i
;
nswap_local
++
;
}
}
}
MPI_Allreduce
(
&
nswap_local
,
&
nswap
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Scan
(
&
nswap_local
,
&
nswap_before
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
nswap_before
-=
nswap_local
;
}
/* ----------------------------------------------------------------------
update the list of gas atoms
------------------------------------------------------------------------- */
void
FixAtomSwap
::
update_swap_atoms_list
()
{
int
nlocal
=
atom
->
nlocal
;
int
*
type
=
atom
->
type
;
double
**
x
=
atom
->
x
;
if
(
atom
->
nmax
>
atom_swap_nmax
)
{
memory
->
sfree
(
local_swap_iatom_list
);
memory
->
sfree
(
local_swap_jatom_list
);
atom_swap_nmax
=
atom
->
nmax
;
local_swap_iatom_list
=
(
int
*
)
memory
->
smalloc
(
atom_swap_nmax
*
sizeof
(
int
),
"MCSWAP:local_swap_iatom_list"
);
local_swap_jatom_list
=
(
int
*
)
memory
->
smalloc
(
atom_swap_nmax
*
sizeof
(
int
),
"MCSWAP:local_swap_jatom_list"
);
}
niswap_local
=
0
;
njswap_local
=
0
;
if
(
regionflag
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
])
==
1
)
{
if
(
atom
->
mask
[
i
]
&
groupbit
)
{
if
(
type
[
i
]
==
type_list
[
0
])
{
local_swap_iatom_list
[
niswap_local
]
=
i
;
niswap_local
++
;
}
else
if
(
type
[
i
]
==
type_list
[
1
])
{
local_swap_jatom_list
[
njswap_local
]
=
i
;
njswap_local
++
;
}
}
}
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
mask
[
i
]
&
groupbit
)
{
if
(
type
[
i
]
==
type_list
[
0
])
{
local_swap_iatom_list
[
niswap_local
]
=
i
;
niswap_local
++
;
}
else
if
(
type
[
i
]
==
type_list
[
1
])
{
local_swap_jatom_list
[
njswap_local
]
=
i
;
njswap_local
++
;
}
}
}
}
MPI_Allreduce
(
&
niswap_local
,
&
niswap
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Scan
(
&
niswap_local
,
&
niswap_before
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
niswap_before
-=
niswap_local
;
MPI_Allreduce
(
&
njswap_local
,
&
njswap
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Scan
(
&
njswap_local
,
&
njswap_before
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
njswap_before
-=
njswap_local
;
}
/* ---------------------------------------------------------------------- */
int
FixAtomSwap
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
int
i
,
j
,
m
;
int
*
type
=
atom
->
type
;
double
*
q
=
atom
->
q
;
m
=
0
;
if
(
atom
->
q_flag
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
type
[
j
];
buf
[
m
++
]
=
q
[
j
];
}
}
else
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
type
[
j
];
}
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
FixAtomSwap
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
int
*
type
=
atom
->
type
;
double
*
q
=
atom
->
q
;
m
=
0
;
last
=
first
+
n
;
if
(
atom
->
q_flag
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
type
[
i
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
q
[
i
]
=
buf
[
m
++
];
}
}
else
{
for
(
i
=
first
;
i
<
last
;
i
++
)
type
[
i
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
}
}
/* ----------------------------------------------------------------------
return acceptance ratio
------------------------------------------------------------------------- */
double
FixAtomSwap
::
compute_vector
(
int
n
)
{
if
(
n
==
0
)
return
nswap_attempts
;
if
(
n
==
1
)
return
nswap_successes
;
return
0.0
;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double
FixAtomSwap
::
memory_usage
()
{
double
bytes
=
atom_swap_nmax
*
sizeof
(
int
);
return
bytes
;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void
FixAtomSwap
::
write_restart
(
FILE
*
fp
)
{
int
n
=
0
;
double
list
[
4
];
list
[
n
++
]
=
random_equal
->
state
();
list
[
n
++
]
=
random_unequal
->
state
();
list
[
n
++
]
=
next_reneighbor
;
if
(
comm
->
me
==
0
)
{
int
size
=
n
*
sizeof
(
double
);
fwrite
(
&
size
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
list
,
sizeof
(
double
),
n
,
fp
);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void
FixAtomSwap
::
restart
(
char
*
buf
)
{
int
n
=
0
;
double
*
list
=
(
double
*
)
buf
;
seed
=
static_cast
<
int
>
(
list
[
n
++
]);
random_equal
->
reset
(
seed
);
seed
=
static_cast
<
int
>
(
list
[
n
++
]);
random_unequal
->
reset
(
seed
);
next_reneighbor
=
static_cast
<
int
>
(
list
[
n
++
]);
}
Event Timeline
Log In to Comment