Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88516594
fix_gcmc.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Oct 19, 06:21
Size
70 KB
Mime Type
text/x-c
Expires
Mon, Oct 21, 06:21 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21784936
Attached To
rLAMMPS lammps
fix_gcmc.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 author: Paul Crozier, Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_gcmc.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_hybrid.h"
#include "molecule.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_extra.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
;
enum
{
ATOM
,
MOLECULE
};
/* ---------------------------------------------------------------------- */
FixGCMC
::
FixGCMC
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
),
idregion
(
NULL
),
full_flag
(
0
),
ngroups
(
0
),
groupstrings
(
NULL
),
ngrouptypes
(
0
),
grouptypestrings
(
NULL
),
grouptypebits
(
NULL
),
grouptypes
(
NULL
),
local_gas_list
(
NULL
),
atom_coord
(
NULL
),
random_equal
(
NULL
),
random_unequal
(
NULL
),
coords
(
NULL
),
imageflags
(
NULL
),
idshake
(
NULL
)
{
if
(
narg
<
11
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
atom
->
molecular
==
2
)
error
->
all
(
FLERR
,
"Fix gcmc does not (yet) work with atom_style template"
);
dynamic_group_allow
=
1
;
vector_flag
=
1
;
size_vector
=
8
;
global_freq
=
1
;
extvector
=
0
;
restart_global
=
1
;
time_depend
=
1
;
// required args
nevery
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
nexchanges
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
nmcmoves
=
force
->
inumeric
(
FLERR
,
arg
[
5
]);
ngcmc_type
=
force
->
inumeric
(
FLERR
,
arg
[
6
]);
seed
=
force
->
inumeric
(
FLERR
,
arg
[
7
]);
reservoir_temperature
=
force
->
numeric
(
FLERR
,
arg
[
8
]);
chemical_potential
=
force
->
numeric
(
FLERR
,
arg
[
9
]);
displace
=
force
->
numeric
(
FLERR
,
arg
[
10
]);
if
(
nevery
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
nexchanges
<
0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
nmcmoves
<
0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
seed
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
reservoir_temperature
<
0.0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
displace
<
0.0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
// read options from end of input line
options
(
narg
-
11
,
&
arg
[
11
]);
// 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
);
// error checks on region and its extent being inside simulation box
region_xlo
=
region_xhi
=
region_ylo
=
region_yhi
=
region_zlo
=
region_zhi
=
0.0
;
if
(
regionflag
)
{
if
(
domain
->
regions
[
iregion
]
->
bboxflag
==
0
)
error
->
all
(
FLERR
,
"Fix gcmc region does not support a bounding box"
);
if
(
domain
->
regions
[
iregion
]
->
dynamic_check
())
error
->
all
(
FLERR
,
"Fix gcmc region cannot be dynamic"
);
region_xlo
=
domain
->
regions
[
iregion
]
->
extent_xlo
;
region_xhi
=
domain
->
regions
[
iregion
]
->
extent_xhi
;
region_ylo
=
domain
->
regions
[
iregion
]
->
extent_ylo
;
region_yhi
=
domain
->
regions
[
iregion
]
->
extent_yhi
;
region_zlo
=
domain
->
regions
[
iregion
]
->
extent_zlo
;
region_zhi
=
domain
->
regions
[
iregion
]
->
extent_zhi
;
if
(
region_xlo
<
domain
->
boxlo
[
0
]
||
region_xhi
>
domain
->
boxhi
[
0
]
||
region_ylo
<
domain
->
boxlo
[
1
]
||
region_yhi
>
domain
->
boxhi
[
1
]
||
region_zlo
<
domain
->
boxlo
[
2
]
||
region_zhi
>
domain
->
boxhi
[
2
])
error
->
all
(
FLERR
,
"Fix gcmc region extends outside simulation box"
);
// estimate region volume using MC trials
double
coord
[
3
];
int
inside
=
0
;
int
attempts
=
10000000
;
for
(
int
i
=
0
;
i
<
attempts
;
i
++
)
{
coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
if
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
!=
0
)
inside
++
;
}
double
max_region_volume
=
(
region_xhi
-
region_xlo
)
*
(
region_yhi
-
region_ylo
)
*
(
region_zhi
-
region_zlo
);
region_volume
=
max_region_volume
*
static_cast
<
double
>
(
inside
)
/
static_cast
<
double
>
(
attempts
);
}
// error check and further setup for mode = MOLECULE
if
(
mode
==
MOLECULE
)
{
if
(
onemols
[
imol
]
->
xflag
==
0
)
error
->
all
(
FLERR
,
"Fix gcmc molecule must have coordinates"
);
if
(
onemols
[
imol
]
->
typeflag
==
0
)
error
->
all
(
FLERR
,
"Fix gcmc molecule must have atom types"
);
if
(
ngcmc_type
!=
0
)
error
->
all
(
FLERR
,
"Atom type must be zero in fix gcmc mol command"
);
if
(
onemols
[
imol
]
->
qflag
==
1
&&
atom
->
q
==
NULL
)
error
->
all
(
FLERR
,
"Fix gcmc molecule has charges, but atom style does not"
);
if
(
atom
->
molecular
==
2
&&
onemols
!=
atom
->
avec
->
onemols
)
error
->
all
(
FLERR
,
"Fix gcmc molecule template ID must be same "
"as atom_style template ID"
);
onemols
[
imol
]
->
check_attributes
(
0
);
}
if
(
charge_flag
&&
atom
->
q
==
NULL
)
error
->
all
(
FLERR
,
"Fix gcmc atom has charge, but atom style does not"
);
if
(
shakeflag
&&
mode
==
ATOM
)
error
->
all
(
FLERR
,
"Cannot use fix gcmc shake and not molecule"
);
// setup of coords and imageflags array
if
(
mode
==
ATOM
)
natoms_per_molecule
=
1
;
else
natoms_per_molecule
=
onemols
[
imol
]
->
natoms
;
memory
->
create
(
coords
,
natoms_per_molecule
,
3
,
"gcmc:coords"
);
memory
->
create
(
imageflags
,
natoms_per_molecule
,
"gcmc:imageflags"
);
memory
->
create
(
atom_coord
,
natoms_per_molecule
,
3
,
"gcmc:atom_coord"
);
// compute the number of MC cycles that occur nevery timesteps
ncycles
=
nexchanges
+
nmcmoves
;
// set up reneighboring
force_reneighbor
=
1
;
next_reneighbor
=
update
->
ntimestep
+
1
;
// zero out counters
ntranslation_attempts
=
0.0
;
ntranslation_successes
=
0.0
;
nrotation_attempts
=
0.0
;
nrotation_successes
=
0.0
;
ndeletion_attempts
=
0.0
;
ndeletion_successes
=
0.0
;
ninsertion_attempts
=
0.0
;
ninsertion_successes
=
0.0
;
gcmc_nmax
=
0
;
local_gas_list
=
NULL
;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void
FixGCMC
::
options
(
int
narg
,
char
**
arg
)
{
if
(
narg
<
0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
// defaults
mode
=
ATOM
;
max_rotation_angle
=
10
*
MY_PI
/
180
;
regionflag
=
0
;
iregion
=
-
1
;
region_volume
=
0
;
max_region_attempts
=
1000
;
molecule_group
=
0
;
molecule_group_bit
=
0
;
molecule_group_inversebit
=
0
;
exclusion_group
=
0
;
exclusion_group_bit
=
0
;
pressure_flag
=
false
;
pressure
=
0.0
;
fugacity_coeff
=
1.0
;
shakeflag
=
0
;
charge
=
0.0
;
charge_flag
=
false
;
full_flag
=
false
;
idshake
=
NULL
;
ngroups
=
0
;
int
ngroupsmax
=
0
;
groupstrings
=
NULL
;
ngrouptypes
=
0
;
int
ngrouptypesmax
=
0
;
grouptypestrings
=
NULL
;
grouptypes
=
NULL
;
grouptypebits
=
NULL
;
energy_intra
=
0.0
;
tfac_insert
=
1.0
;
int
iarg
=
0
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"mol"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
imol
=
atom
->
find_molecule
(
arg
[
iarg
+
1
]);
if
(
imol
==
-
1
)
error
->
all
(
FLERR
,
"Molecule template ID for fix gcmc does not exist"
);
if
(
atom
->
molecules
[
imol
]
->
nset
>
1
&&
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"Molecule template for "
"fix gcmc has multiple molecules"
);
mode
=
MOLECULE
;
onemols
=
atom
->
molecules
;
nmol
=
onemols
[
imol
]
->
nset
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"region"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
iregion
=
domain
->
find_region
(
arg
[
iarg
+
1
]);
if
(
iregion
==
-
1
)
error
->
all
(
FLERR
,
"Region ID for fix gcmc 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
],
"maxangle"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
max_rotation_angle
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
max_rotation_angle
*=
MY_PI
/
180
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"pressure"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
pressure
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
pressure_flag
=
true
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"fugacity_coeff"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
fugacity_coeff
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"charge"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
charge
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
charge_flag
=
true
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"shake"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
delete
[]
idshake
;
idshake
=
new
char
[
n
];
strcpy
(
idshake
,
arg
[
iarg
+
1
]);
shakeflag
=
1
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"full_energy"
)
==
0
)
{
full_flag
=
true
;
iarg
+=
1
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"group"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
ngroups
>=
ngroupsmax
)
{
ngroupsmax
=
ngroups
+
1
;
groupstrings
=
(
char
**
)
memory
->
srealloc
(
groupstrings
,
ngroupsmax
*
sizeof
(
char
*
),
"fix_gcmc:groupstrings"
);
}
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
groupstrings
[
ngroups
]
=
new
char
[
n
];
strcpy
(
groupstrings
[
ngroups
],
arg
[
iarg
+
1
]);
ngroups
++
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"grouptype"
)
==
0
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
if
(
ngrouptypes
>=
ngrouptypesmax
)
{
ngrouptypesmax
=
ngrouptypes
+
1
;
grouptypes
=
(
int
*
)
memory
->
srealloc
(
grouptypes
,
ngrouptypesmax
*
sizeof
(
int
),
"fix_gcmc:grouptypes"
);
grouptypestrings
=
(
char
**
)
memory
->
srealloc
(
grouptypestrings
,
ngrouptypesmax
*
sizeof
(
char
*
),
"fix_gcmc:grouptypestrings"
);
}
grouptypes
[
ngrouptypes
]
=
atoi
(
arg
[
iarg
+
1
]);
int
n
=
strlen
(
arg
[
iarg
+
2
])
+
1
;
grouptypestrings
[
ngrouptypes
]
=
new
char
[
n
];
strcpy
(
grouptypestrings
[
ngrouptypes
],
arg
[
iarg
+
2
]);
ngrouptypes
++
;
iarg
+=
3
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"intra_energy"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
energy_intra
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"tfac_insert"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
tfac_insert
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
error
->
all
(
FLERR
,
"Illegal fix gcmc command"
);
}
}
/* ---------------------------------------------------------------------- */
FixGCMC
::~
FixGCMC
()
{
if
(
regionflag
)
delete
[]
idregion
;
delete
random_equal
;
delete
random_unequal
;
memory
->
destroy
(
local_gas_list
);
memory
->
destroy
(
atom_coord
);
memory
->
destroy
(
coords
);
memory
->
destroy
(
imageflags
);
delete
[]
idshake
;
if
(
ngroups
>
0
)
{
for
(
int
igroup
=
0
;
igroup
<
ngroups
;
igroup
++
)
delete
[]
groupstrings
[
igroup
];
memory
->
sfree
(
groupstrings
);
}
if
(
ngrouptypes
>
0
)
{
memory
->
destroy
(
grouptypes
);
memory
->
destroy
(
grouptypebits
);
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
delete
[]
grouptypestrings
[
igroup
];
memory
->
sfree
(
grouptypestrings
);
}
if
(
full_flag
&&
group
)
{
int
igroupall
=
group
->
find
(
"all"
);
neighbor
->
exclusion_group_group_delete
(
exclusion_group
,
igroupall
);
}
}
/* ---------------------------------------------------------------------- */
int
FixGCMC
::
setmask
()
{
int
mask
=
0
;
mask
|=
PRE_EXCHANGE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixGCMC
::
init
()
{
triclinic
=
domain
->
triclinic
;
// decide whether to switch to the full_energy option
if
(
!
full_flag
)
{
if
((
force
->
kspace
)
||
(
force
->
pair
==
NULL
)
||
(
force
->
pair
->
single_enable
==
0
)
||
(
force
->
pair_match
(
"hybrid"
,
0
))
||
(
force
->
pair_match
(
"eam"
,
0
))
)
{
full_flag
=
true
;
if
(
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"Fix gcmc using full_energy option"
);
}
}
if
(
full_flag
)
{
char
*
id_pe
=
(
char
*
)
"thermo_pe"
;
int
ipe
=
modify
->
find_compute
(
id_pe
);
c_pe
=
modify
->
compute
[
ipe
];
}
int
*
type
=
atom
->
type
;
if
(
mode
==
ATOM
)
{
if
(
ngcmc_type
<=
0
||
ngcmc_type
>
atom
->
ntypes
)
error
->
all
(
FLERR
,
"Invalid atom type in fix gcmc command"
);
}
// if mode == ATOM, warn if any deletable atom has a mol ID
if
((
mode
==
ATOM
)
&&
atom
->
molecule_flag
)
{
tagint
*
molecule
=
atom
->
molecule
;
int
flag
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
if
(
type
[
i
]
==
ngcmc_type
)
if
(
molecule
[
i
])
flag
=
1
;
int
flagall
;
MPI_Allreduce
(
&
flag
,
&
flagall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flagall
&&
comm
->
me
==
0
)
error
->
all
(
FLERR
,
"Fix gcmc cannot exchange individual atoms belonging to a molecule"
);
}
// if mode == MOLECULE, check for unset mol IDs
if
(
mode
==
MOLECULE
)
{
tagint
*
molecule
=
atom
->
molecule
;
int
*
mask
=
atom
->
mask
;
int
flag
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
if
(
mask
[
i
]
==
groupbit
)
if
(
molecule
[
i
]
==
0
)
flag
=
1
;
int
flagall
;
MPI_Allreduce
(
&
flag
,
&
flagall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flagall
&&
comm
->
me
==
0
)
error
->
all
(
FLERR
,
"All mol IDs should be set for fix gcmc group atoms"
);
}
if
(((
mode
==
MOLECULE
)
&&
(
atom
->
molecule_flag
==
0
))
||
((
mode
==
MOLECULE
)
&&
(
!
atom
->
tag_enable
||
!
atom
->
map_style
)))
error
->
all
(
FLERR
,
"Fix gcmc molecule command requires that "
"atoms have molecule attributes"
);
// if shakeflag defined, check for SHAKE fix
// its molecule template must be same as this one
fixshake
=
NULL
;
if
(
shakeflag
)
{
int
ifix
=
modify
->
find_fix
(
idshake
);
if
(
ifix
<
0
)
error
->
all
(
FLERR
,
"Fix gcmc shake fix does not exist"
);
fixshake
=
modify
->
fix
[
ifix
];
int
tmp
;
if
(
onemols
!=
(
Molecule
**
)
fixshake
->
extract
(
"onemol"
,
tmp
))
error
->
all
(
FLERR
,
"Fix gcmc and fix shake not using "
"same molecule template ID"
);
}
// check for fix rigid
for
(
int
irigid
=
0
;
irigid
<
modify
->
nfix
;
irigid
++
)
{
if
(
strncmp
(
modify
->
fix
[
irigid
]
->
style
,
"rigid"
,
5
)
==
0
)
error
->
all
(
FLERR
,
"Fix gcmc can not currently be used with any rigid fix"
);
}
if
(
domain
->
dimension
==
2
)
error
->
all
(
FLERR
,
"Cannot use fix gcmc in a 2d simulation"
);
// create a new group for interaction exclusions
// used for attempted atom or molecule deletions
// skip if already exists from previous init()
if
(
full_flag
&&
!
exclusion_group_bit
)
{
char
**
group_arg
=
new
char
*
[
4
];
// create unique group name for atoms to be excluded
int
len
=
strlen
(
id
)
+
30
;
group_arg
[
0
]
=
new
char
[
len
];
sprintf
(
group_arg
[
0
],
"FixGCMC:gcmc_exclusion_group:%s"
,
id
);
group_arg
[
1
]
=
(
char
*
)
"subtract"
;
group_arg
[
2
]
=
(
char
*
)
"all"
;
group_arg
[
3
]
=
(
char
*
)
"all"
;
group
->
assign
(
4
,
group_arg
);
exclusion_group
=
group
->
find
(
group_arg
[
0
]);
if
(
exclusion_group
==
-
1
)
error
->
all
(
FLERR
,
"Could not find fix gcmc exclusion group ID"
);
exclusion_group_bit
=
group
->
bitmask
[
exclusion_group
];
// neighbor list exclusion setup
// turn off interactions between group all and the exclusion group
int
narg
=
4
;
char
**
arg
=
new
char
*
[
narg
];;
arg
[
0
]
=
(
char
*
)
"exclude"
;
arg
[
1
]
=
(
char
*
)
"group"
;
arg
[
2
]
=
group_arg
[
0
];
arg
[
3
]
=
(
char
*
)
"all"
;
neighbor
->
modify_params
(
narg
,
arg
);
delete
[]
group_arg
[
0
];
delete
[]
group_arg
;
delete
[]
arg
;
}
// create a new group for temporary use with selected molecules
if
(
mode
==
MOLECULE
)
{
char
**
group_arg
=
new
char
*
[
3
];
// create unique group name for atoms to be rotated
int
len
=
strlen
(
id
)
+
30
;
group_arg
[
0
]
=
new
char
[
len
];
sprintf
(
group_arg
[
0
],
"FixGCMC:rotation_gas_atoms:%s"
,
id
);
group_arg
[
1
]
=
(
char
*
)
"molecule"
;
char
digits
[
12
];
sprintf
(
digits
,
"%d"
,
-
1
);
group_arg
[
2
]
=
digits
;
group
->
assign
(
3
,
group_arg
);
molecule_group
=
group
->
find
(
group_arg
[
0
]);
if
(
molecule_group
==
-
1
)
error
->
all
(
FLERR
,
"Could not find fix gcmc rotation group ID"
);
molecule_group_bit
=
group
->
bitmask
[
molecule_group
];
molecule_group_inversebit
=
molecule_group_bit
^
~
0
;
delete
[]
group_arg
[
0
];
delete
[]
group_arg
;
}
// get all of the needed molecule data if mode == MOLECULE,
// otherwise just get the gas mass
if
(
mode
==
MOLECULE
)
{
onemols
[
imol
]
->
compute_mass
();
onemols
[
imol
]
->
compute_com
();
gas_mass
=
onemols
[
imol
]
->
masstotal
;
for
(
int
i
=
0
;
i
<
onemols
[
imol
]
->
natoms
;
i
++
)
{
onemols
[
imol
]
->
x
[
i
][
0
]
-=
onemols
[
imol
]
->
com
[
0
];
onemols
[
imol
]
->
x
[
i
][
1
]
-=
onemols
[
imol
]
->
com
[
1
];
onemols
[
imol
]
->
x
[
i
][
2
]
-=
onemols
[
imol
]
->
com
[
2
];
}
}
else
gas_mass
=
atom
->
mass
[
ngcmc_type
];
if
(
gas_mass
<=
0.0
)
error
->
all
(
FLERR
,
"Illegal fix gcmc gas mass <= 0"
);
// check that no deletable atoms are in atom->firstgroup
// deleting such an atom would 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 GCMC on atoms in atom_modify first group"
);
}
// compute beta, lambda, sigma, and the zz factor
beta
=
1.0
/
(
force
->
boltz
*
reservoir_temperature
);
double
lambda
=
sqrt
(
force
->
hplanck
*
force
->
hplanck
/
(
2.0
*
MY_PI
*
gas_mass
*
force
->
mvv2e
*
force
->
boltz
*
reservoir_temperature
));
sigma
=
sqrt
(
force
->
boltz
*
reservoir_temperature
*
tfac_insert
/
gas_mass
/
force
->
mvv2e
);
zz
=
exp
(
beta
*
chemical_potential
)
/
(
pow
(
lambda
,
3.0
));
if
(
pressure_flag
)
zz
=
pressure
*
fugacity_coeff
*
beta
/
force
->
nktv2p
;
imagezero
=
((
imageint
)
IMGMAX
<<
IMG2BITS
)
|
((
imageint
)
IMGMAX
<<
IMGBITS
)
|
IMGMAX
;
// construct group bitmask for all new atoms
// aggregated over all group keywords
groupbitall
=
1
|
groupbit
;
for
(
int
igroup
=
0
;
igroup
<
ngroups
;
igroup
++
)
{
int
jgroup
=
group
->
find
(
groupstrings
[
igroup
]);
if
(
jgroup
==
-
1
)
error
->
all
(
FLERR
,
"Could not find specified fix gcmc group ID"
);
groupbitall
|=
group
->
bitmask
[
jgroup
];
}
// construct group type bitmasks
// not aggregated over all group keywords
if
(
ngrouptypes
>
0
)
{
memory
->
create
(
grouptypebits
,
ngrouptypes
,
"fix_gcmc:grouptypebits"
);
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
{
int
jgroup
=
group
->
find
(
grouptypestrings
[
igroup
]);
if
(
jgroup
==
-
1
)
error
->
all
(
FLERR
,
"Could not find specified fix gcmc group ID"
);
grouptypebits
[
igroup
]
=
group
->
bitmask
[
jgroup
];
}
}
}
/* ----------------------------------------------------------------------
attempt Monte Carlo translations, rotations, insertions, and deletions
done before exchange, borders, reneighbor
so that ghost atoms and neighbor lists will be correct
------------------------------------------------------------------------- */
void
FixGCMC
::
pre_exchange
()
{
// just return if should not be called on this timestep
if
(
next_reneighbor
!=
update
->
ntimestep
)
return
;
xlo
=
domain
->
boxlo
[
0
];
xhi
=
domain
->
boxhi
[
0
];
ylo
=
domain
->
boxlo
[
1
];
yhi
=
domain
->
boxhi
[
1
];
zlo
=
domain
->
boxlo
[
2
];
zhi
=
domain
->
boxhi
[
2
];
if
(
triclinic
)
{
sublo
=
domain
->
sublo_lamda
;
subhi
=
domain
->
subhi_lamda
;
}
else
{
sublo
=
domain
->
sublo
;
subhi
=
domain
->
subhi
;
}
if
(
regionflag
)
volume
=
region_volume
;
else
volume
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
if
(
full_flag
)
{
energy_stored
=
energy_full
();
if
(
mode
==
MOLECULE
)
{
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
{
int
random_int_fraction
=
static_cast
<
int
>
(
random_equal
->
uniform
()
*
ncycles
)
+
1
;
if
(
random_int_fraction
<=
nmcmoves
)
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_molecule_translation_full
();
else
attempt_molecule_rotation_full
();
}
else
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_molecule_deletion_full
();
else
attempt_molecule_insertion_full
();
}
}
}
else
{
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
{
int
random_int_fraction
=
static_cast
<
int
>
(
random_equal
->
uniform
()
*
ncycles
)
+
1
;
if
(
random_int_fraction
<=
nmcmoves
)
{
attempt_atomic_translation_full
();
}
else
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_atomic_deletion_full
();
else
attempt_atomic_insertion_full
();
}
}
}
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
}
else
{
if
(
mode
==
MOLECULE
)
{
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
{
int
random_int_fraction
=
static_cast
<
int
>
(
random_equal
->
uniform
()
*
ncycles
)
+
1
;
if
(
random_int_fraction
<=
nmcmoves
)
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_molecule_translation
();
else
attempt_molecule_rotation
();
}
else
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_molecule_deletion
();
else
attempt_molecule_insertion
();
}
}
}
else
{
for
(
int
i
=
0
;
i
<
ncycles
;
i
++
)
{
int
random_int_fraction
=
static_cast
<
int
>
(
random_equal
->
uniform
()
*
ncycles
)
+
1
;
if
(
random_int_fraction
<=
nmcmoves
)
{
attempt_atomic_translation
();
}
else
{
if
(
random_equal
->
uniform
()
<
0.5
)
attempt_atomic_deletion
();
else
attempt_atomic_insertion
();
}
}
}
}
next_reneighbor
=
update
->
ntimestep
+
nevery
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_translation
()
{
ntranslation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
int
i
=
pick_random_gas_atom
();
int
success
=
0
;
if
(
i
>=
0
)
{
double
**
x
=
atom
->
x
;
double
energy_before
=
energy
(
i
,
ngcmc_type
,
-
1
,
x
[
i
]);
double
rsq
=
1.1
;
double
rx
,
ry
,
rz
;
rx
=
ry
=
rz
=
0.0
;
double
coord
[
3
];
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_unequal
->
uniform
()
-
1.0
;
ry
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rz
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
x
[
i
][
0
]
+
displace
*
rx
;
coord
[
1
]
=
x
[
i
][
1
]
+
displace
*
ry
;
coord
[
2
]
=
x
[
i
][
2
]
+
displace
*
rz
;
if
(
regionflag
)
{
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_unequal
->
uniform
()
-
1.0
;
ry
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rz
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
x
[
i
][
0
]
+
displace
*
rx
;
coord
[
1
]
=
x
[
i
][
1
]
+
displace
*
ry
;
coord
[
2
]
=
x
[
i
][
2
]
+
displace
*
rz
;
}
}
if
(
!
domain
->
inside_nonperiodic
(
coord
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
double
energy_after
=
energy
(
i
,
ngcmc_type
,
-
1
,
coord
);
if
(
random_unequal
->
uniform
()
<
exp
(
beta
*
(
energy_before
-
energy_after
)))
{
x
[
i
][
0
]
=
coord
[
0
];
x
[
i
][
1
]
=
coord
[
1
];
x
[
i
][
2
]
=
coord
[
2
];
success
=
1
;
}
}
int
success_all
=
0
;
MPI_Allreduce
(
&
success
,
&
success_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
if
(
success_all
)
{
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ntranslation_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_deletion
()
{
ndeletion_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
int
i
=
pick_random_gas_atom
();
int
success
=
0
;
if
(
i
>=
0
)
{
double
deletion_energy
=
energy
(
i
,
ngcmc_type
,
-
1
,
atom
->
x
[
i
]);
if
(
random_unequal
->
uniform
()
<
ngas
*
exp
(
beta
*
deletion_energy
)
/
(
zz
*
volume
))
{
atom
->
avec
->
copy
(
atom
->
nlocal
-
1
,
i
,
1
);
atom
->
nlocal
--
;
success
=
1
;
}
}
int
success_all
=
0
;
MPI_Allreduce
(
&
success
,
&
success_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
if
(
success_all
)
{
atom
->
natoms
--
;
if
(
atom
->
tag_enable
)
{
if
(
atom
->
map_style
)
atom
->
map_init
();
}
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ndeletion_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_insertion
()
{
double
lamda
[
3
];
ninsertion_attempts
+=
1.0
;
// pick coordinates for insertion point
double
coord
[
3
];
if
(
regionflag
)
{
int
region_attempt
=
0
;
coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
region_attempt
++
;
if
(
region_attempt
>=
max_region_attempts
)
return
;
}
if
(
triclinic
)
domain
->
x2lamda
(
coord
,
lamda
);
}
else
{
if
(
triclinic
==
0
)
{
coord
[
0
]
=
xlo
+
random_equal
->
uniform
()
*
(
xhi
-
xlo
);
coord
[
1
]
=
ylo
+
random_equal
->
uniform
()
*
(
yhi
-
ylo
);
coord
[
2
]
=
zlo
+
random_equal
->
uniform
()
*
(
zhi
-
zlo
);
}
else
{
lamda
[
0
]
=
random_equal
->
uniform
();
lamda
[
1
]
=
random_equal
->
uniform
();
lamda
[
2
]
=
random_equal
->
uniform
();
// wasteful, but necessary
if
(
lamda
[
0
]
==
1.0
)
lamda
[
0
]
=
0.0
;
if
(
lamda
[
1
]
==
1.0
)
lamda
[
1
]
=
0.0
;
if
(
lamda
[
2
]
==
1.0
)
lamda
[
2
]
=
0.0
;
domain
->
lamda2x
(
lamda
,
coord
);
}
}
int
proc_flag
=
0
;
if
(
triclinic
==
0
)
{
domain
->
remap
(
coord
);
if
(
!
domain
->
inside
(
coord
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
if
(
coord
[
0
]
>=
sublo
[
0
]
&&
coord
[
0
]
<
subhi
[
0
]
&&
coord
[
1
]
>=
sublo
[
1
]
&&
coord
[
1
]
<
subhi
[
1
]
&&
coord
[
2
]
>=
sublo
[
2
]
&&
coord
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
else
{
if
(
lamda
[
0
]
>=
sublo
[
0
]
&&
lamda
[
0
]
<
subhi
[
0
]
&&
lamda
[
1
]
>=
sublo
[
1
]
&&
lamda
[
1
]
<
subhi
[
1
]
&&
lamda
[
2
]
>=
sublo
[
2
]
&&
lamda
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
int
success
=
0
;
if
(
proc_flag
)
{
int
ii
=
-
1
;
if
(
charge_flag
)
{
ii
=
atom
->
nlocal
+
atom
->
nghost
;
if
(
ii
>=
atom
->
nmax
)
atom
->
avec
->
grow
(
0
);
atom
->
q
[
ii
]
=
charge
;
}
double
insertion_energy
=
energy
(
ii
,
ngcmc_type
,
-
1
,
coord
);
if
(
random_unequal
->
uniform
()
<
zz
*
volume
*
exp
(
-
beta
*
insertion_energy
)
/
(
ngas
+
1
))
{
atom
->
avec
->
create_atom
(
ngcmc_type
,
coord
);
int
m
=
atom
->
nlocal
-
1
;
// add to groups
// optionally add to type-based groups
atom
->
mask
[
m
]
=
groupbitall
;
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
{
if
(
ngcmc_type
==
grouptypes
[
igroup
])
atom
->
mask
[
m
]
|=
grouptypebits
[
igroup
];
}
atom
->
v
[
m
][
0
]
=
random_unequal
->
gaussian
()
*
sigma
;
atom
->
v
[
m
][
1
]
=
random_unequal
->
gaussian
()
*
sigma
;
atom
->
v
[
m
][
2
]
=
random_unequal
->
gaussian
()
*
sigma
;
modify
->
create_attribute
(
m
);
success
=
1
;
}
}
int
success_all
=
0
;
MPI_Allreduce
(
&
success
,
&
success_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
if
(
success_all
)
{
atom
->
natoms
++
;
if
(
atom
->
tag_enable
)
{
atom
->
tag_extend
();
if
(
atom
->
map_style
)
atom
->
map_init
();
}
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ninsertion_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_translation
()
{
ntranslation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
translation_molecule
=
pick_random_gas_molecule
();
if
(
translation_molecule
==
-
1
)
return
;
double
energy_before_sum
=
molecule_energy
(
translation_molecule
);
double
**
x
=
atom
->
x
;
double
rx
,
ry
,
rz
;
double
com_displace
[
3
],
coord
[
3
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_equal
->
uniform
()
-
1.0
;
ry
=
2
*
random_equal
->
uniform
()
-
1.0
;
rz
=
2
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
com_displace
[
0
]
=
displace
*
rx
;
com_displace
[
1
]
=
displace
*
ry
;
com_displace
[
2
]
=
displace
*
rz
;
int
nlocal
=
atom
->
nlocal
;
if
(
regionflag
)
{
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
mask
[
i
]
|=
molecule_group_bit
;
}
else
{
mask
[
i
]
&=
molecule_group_inversebit
;
}
}
double
com
[
3
];
com
[
0
]
=
com
[
1
]
=
com
[
2
]
=
0.0
;
group
->
xcm
(
molecule_group
,
gas_mass
,
com
);
coord
[
0
]
=
com
[
0
]
+
displace
*
rx
;
coord
[
1
]
=
com
[
1
]
+
displace
*
ry
;
coord
[
2
]
=
com
[
2
]
+
displace
*
rz
;
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_equal
->
uniform
()
-
1.0
;
ry
=
2
*
random_equal
->
uniform
()
-
1.0
;
rz
=
2
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
com
[
0
]
+
displace
*
rx
;
coord
[
1
]
=
com
[
1
]
+
displace
*
ry
;
coord
[
2
]
=
com
[
2
]
+
displace
*
rz
;
}
com_displace
[
0
]
=
displace
*
rx
;
com_displace
[
1
]
=
displace
*
ry
;
com_displace
[
2
]
=
displace
*
rz
;
}
double
energy_after
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
coord
[
0
]
=
x
[
i
][
0
]
+
com_displace
[
0
];
coord
[
1
]
=
x
[
i
][
1
]
+
com_displace
[
1
];
coord
[
2
]
=
x
[
i
][
2
]
+
com_displace
[
2
];
if
(
!
domain
->
inside_nonperiodic
(
coord
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
energy_after
+=
energy
(
i
,
atom
->
type
[
i
],
translation_molecule
,
coord
);
}
}
double
energy_after_sum
=
0.0
;
MPI_Allreduce
(
&
energy_after
,
&
energy_after_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before_sum
-
energy_after_sum
)))
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
x
[
i
][
0
]
+=
com_displace
[
0
];
x
[
i
][
1
]
+=
com_displace
[
1
];
x
[
i
][
2
]
+=
com_displace
[
2
];
}
}
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ntranslation_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_rotation
()
{
nrotation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
rotation_molecule
=
pick_random_gas_molecule
();
if
(
rotation_molecule
==
-
1
)
return
;
double
energy_before_sum
=
molecule_energy
(
rotation_molecule
);
int
nlocal
=
atom
->
nlocal
;
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
rotation_molecule
)
{
mask
[
i
]
|=
molecule_group_bit
;
}
else
{
mask
[
i
]
&=
molecule_group_inversebit
;
}
}
double
com
[
3
];
com
[
0
]
=
com
[
1
]
=
com
[
2
]
=
0.0
;
group
->
xcm
(
molecule_group
,
gas_mass
,
com
);
// generate point in unit cube
// then restrict to unit sphere
double
r
[
3
],
rotmat
[
3
][
3
],
quat
[
4
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
r
[
0
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
1
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
2
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
MathExtra
::
dot3
(
r
,
r
);
}
double
theta
=
random_equal
->
uniform
()
*
max_rotation_angle
;
MathExtra
::
norm3
(
r
);
MathExtra
::
axisangle_to_quat
(
r
,
theta
,
quat
);
MathExtra
::
quat_to_mat
(
quat
,
rotmat
);
double
**
x
=
atom
->
x
;
imageint
*
image
=
atom
->
image
;
double
energy_after
=
0.0
;
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
molecule_group_bit
)
{
double
xtmp
[
3
];
domain
->
unmap
(
x
[
i
],
image
[
i
],
xtmp
);
xtmp
[
0
]
-=
com
[
0
];
xtmp
[
1
]
-=
com
[
1
];
xtmp
[
2
]
-=
com
[
2
];
MathExtra
::
matvec
(
rotmat
,
xtmp
,
atom_coord
[
n
]);
atom_coord
[
n
][
0
]
+=
com
[
0
];
atom_coord
[
n
][
1
]
+=
com
[
1
];
atom_coord
[
n
][
2
]
+=
com
[
2
];
xtmp
[
0
]
=
atom_coord
[
n
][
0
];
xtmp
[
1
]
=
atom_coord
[
n
][
1
];
xtmp
[
2
]
=
atom_coord
[
n
][
2
];
domain
->
remap
(
xtmp
);
if
(
!
domain
->
inside
(
xtmp
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
energy_after
+=
energy
(
i
,
atom
->
type
[
i
],
rotation_molecule
,
xtmp
);
n
++
;
}
}
double
energy_after_sum
=
0.0
;
MPI_Allreduce
(
&
energy_after
,
&
energy_after_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before_sum
-
energy_after_sum
)))
{
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
molecule_group_bit
)
{
image
[
i
]
=
imagezero
;
x
[
i
][
0
]
=
atom_coord
[
n
][
0
];
x
[
i
][
1
]
=
atom_coord
[
n
][
1
];
x
[
i
][
2
]
=
atom_coord
[
n
][
2
];
domain
->
remap
(
x
[
i
],
image
[
i
]);
n
++
;
}
}
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
nrotation_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_deletion
()
{
ndeletion_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
deletion_molecule
=
pick_random_gas_molecule
();
if
(
deletion_molecule
==
-
1
)
return
;
double
deletion_energy_sum
=
molecule_energy
(
deletion_molecule
);
if
(
random_equal
->
uniform
()
<
ngas
*
exp
(
beta
*
deletion_energy_sum
)
/
(
zz
*
volume
*
natoms_per_molecule
))
{
int
i
=
0
;
while
(
i
<
atom
->
nlocal
)
{
if
(
atom
->
molecule
[
i
]
==
deletion_molecule
)
{
atom
->
avec
->
copy
(
atom
->
nlocal
-
1
,
i
,
1
);
atom
->
nlocal
--
;
}
else
i
++
;
}
atom
->
natoms
-=
natoms_per_molecule
;
if
(
atom
->
map_style
)
atom
->
map_init
();
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ndeletion_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_insertion
()
{
double
lamda
[
3
];
ninsertion_attempts
+=
1.0
;
double
com_coord
[
3
];
if
(
regionflag
)
{
int
region_attempt
=
0
;
com_coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
com_coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
com_coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
while
(
domain
->
regions
[
iregion
]
->
match
(
com_coord
[
0
],
com_coord
[
1
],
com_coord
[
2
])
==
0
)
{
com_coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
com_coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
com_coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
region_attempt
++
;
if
(
region_attempt
>=
max_region_attempts
)
return
;
}
if
(
triclinic
)
domain
->
x2lamda
(
com_coord
,
lamda
);
}
else
{
if
(
triclinic
==
0
)
{
com_coord
[
0
]
=
xlo
+
random_equal
->
uniform
()
*
(
xhi
-
xlo
);
com_coord
[
1
]
=
ylo
+
random_equal
->
uniform
()
*
(
yhi
-
ylo
);
com_coord
[
2
]
=
zlo
+
random_equal
->
uniform
()
*
(
zhi
-
zlo
);
}
else
{
lamda
[
0
]
=
random_equal
->
uniform
();
lamda
[
1
]
=
random_equal
->
uniform
();
lamda
[
2
]
=
random_equal
->
uniform
();
// wasteful, but necessary
if
(
lamda
[
0
]
==
1.0
)
lamda
[
0
]
=
0.0
;
if
(
lamda
[
1
]
==
1.0
)
lamda
[
1
]
=
0.0
;
if
(
lamda
[
2
]
==
1.0
)
lamda
[
2
]
=
0.0
;
domain
->
lamda2x
(
lamda
,
com_coord
);
}
}
// generate point in unit cube
// then restrict to unit sphere
double
r
[
3
],
rotmat
[
3
][
3
],
quat
[
4
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
r
[
0
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
1
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
2
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
MathExtra
::
dot3
(
r
,
r
);
}
double
theta
=
random_equal
->
uniform
()
*
MY_2PI
;
MathExtra
::
norm3
(
r
);
MathExtra
::
axisangle_to_quat
(
r
,
theta
,
quat
);
MathExtra
::
quat_to_mat
(
quat
,
rotmat
);
double
insertion_energy
=
0.0
;
bool
procflag
[
natoms_per_molecule
];
for
(
int
i
=
0
;
i
<
natoms_per_molecule
;
i
++
)
{
MathExtra
::
matvec
(
rotmat
,
onemols
[
imol
]
->
x
[
i
],
atom_coord
[
i
]);
atom_coord
[
i
][
0
]
+=
com_coord
[
0
];
atom_coord
[
i
][
1
]
+=
com_coord
[
1
];
atom_coord
[
i
][
2
]
+=
com_coord
[
2
];
// use temporary variable for remapped position
// so unmapped position is preserved in atom_coord
double
xtmp
[
3
];
xtmp
[
0
]
=
atom_coord
[
i
][
0
];
xtmp
[
1
]
=
atom_coord
[
i
][
1
];
xtmp
[
2
]
=
atom_coord
[
i
][
2
];
domain
->
remap
(
xtmp
);
if
(
!
domain
->
inside
(
xtmp
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
procflag
[
i
]
=
false
;
if
(
triclinic
==
0
)
{
if
(
xtmp
[
0
]
>=
sublo
[
0
]
&&
xtmp
[
0
]
<
subhi
[
0
]
&&
xtmp
[
1
]
>=
sublo
[
1
]
&&
xtmp
[
1
]
<
subhi
[
1
]
&&
xtmp
[
2
]
>=
sublo
[
2
]
&&
xtmp
[
2
]
<
subhi
[
2
])
procflag
[
i
]
=
true
;
}
else
{
domain
->
x2lamda
(
xtmp
,
lamda
);
if
(
lamda
[
0
]
>=
sublo
[
0
]
&&
lamda
[
0
]
<
subhi
[
0
]
&&
lamda
[
1
]
>=
sublo
[
1
]
&&
lamda
[
1
]
<
subhi
[
1
]
&&
lamda
[
2
]
>=
sublo
[
2
]
&&
lamda
[
2
]
<
subhi
[
2
])
procflag
[
i
]
=
true
;
}
if
(
procflag
[
i
])
{
int
ii
=
-
1
;
if
(
onemols
[
imol
]
->
qflag
==
1
)
{
ii
=
atom
->
nlocal
+
atom
->
nghost
;
if
(
ii
>=
atom
->
nmax
)
atom
->
avec
->
grow
(
0
);
atom
->
q
[
ii
]
=
onemols
[
imol
]
->
q
[
i
];
}
insertion_energy
+=
energy
(
ii
,
onemols
[
imol
]
->
type
[
i
],
-
1
,
xtmp
);
}
}
double
insertion_energy_sum
=
0.0
;
MPI_Allreduce
(
&
insertion_energy
,
&
insertion_energy_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
if
(
random_equal
->
uniform
()
<
zz
*
volume
*
natoms_per_molecule
*
exp
(
-
beta
*
insertion_energy_sum
)
/
(
ngas
+
natoms_per_molecule
))
{
tagint
maxmol
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
maxmol
=
MAX
(
maxmol
,
atom
->
molecule
[
i
]);
tagint
maxmol_all
;
MPI_Allreduce
(
&
maxmol
,
&
maxmol_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
maxmol_all
++
;
if
(
maxmol_all
>=
MAXTAGINT
)
error
->
all
(
FLERR
,
"Fix gcmc ran out of available molecule IDs"
);
tagint
maxtag
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
maxtag
=
MAX
(
maxtag
,
atom
->
tag
[
i
]);
tagint
maxtag_all
;
MPI_Allreduce
(
&
maxtag
,
&
maxtag_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
int
nlocalprev
=
atom
->
nlocal
;
double
vnew
[
3
];
vnew
[
0
]
=
random_equal
->
gaussian
()
*
sigma
;
vnew
[
1
]
=
random_equal
->
gaussian
()
*
sigma
;
vnew
[
2
]
=
random_equal
->
gaussian
()
*
sigma
;
for
(
int
i
=
0
;
i
<
natoms_per_molecule
;
i
++
)
{
if
(
procflag
[
i
])
{
atom
->
avec
->
create_atom
(
onemols
[
imol
]
->
type
[
i
],
atom_coord
[
i
]);
int
m
=
atom
->
nlocal
-
1
;
// add to groups
// optionally add to type-based groups
atom
->
mask
[
m
]
=
groupbitall
;
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
{
if
(
ngcmc_type
==
grouptypes
[
igroup
])
atom
->
mask
[
m
]
|=
grouptypebits
[
igroup
];
}
atom
->
image
[
m
]
=
imagezero
;
domain
->
remap
(
atom
->
x
[
m
],
atom
->
image
[
m
]);
atom
->
molecule
[
m
]
=
maxmol_all
;
if
(
maxtag_all
+
i
+
1
>=
MAXTAGINT
)
error
->
all
(
FLERR
,
"Fix gcmc ran out of available atom IDs"
);
atom
->
tag
[
m
]
=
maxtag_all
+
i
+
1
;
atom
->
v
[
m
][
0
]
=
vnew
[
0
];
atom
->
v
[
m
][
1
]
=
vnew
[
1
];
atom
->
v
[
m
][
2
]
=
vnew
[
2
];
atom
->
add_molecule_atom
(
onemols
[
imol
],
i
,
m
,
maxtag_all
);
modify
->
create_attribute
(
m
);
}
}
if
(
shakeflag
)
fixshake
->
set_molecule
(
nlocalprev
,
maxtag_all
,
imol
,
com_coord
,
vnew
,
quat
);
atom
->
natoms
+=
natoms_per_molecule
;
if
(
atom
->
natoms
<
0
)
error
->
all
(
FLERR
,
"Too many total atoms"
);
atom
->
nbonds
+=
onemols
[
imol
]
->
nbonds
;
atom
->
nangles
+=
onemols
[
imol
]
->
nangles
;
atom
->
ndihedrals
+=
onemols
[
imol
]
->
ndihedrals
;
atom
->
nimpropers
+=
onemols
[
imol
]
->
nimpropers
;
if
(
atom
->
map_style
)
atom
->
map_init
();
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
update_gas_atoms_list
();
ninsertion_successes
+=
1.0
;
}
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_translation_full
()
{
ntranslation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
double
energy_before
=
energy_stored
;
int
i
=
pick_random_gas_atom
();
double
**
x
=
atom
->
x
;
double
xtmp
[
3
];
xtmp
[
0
]
=
xtmp
[
1
]
=
xtmp
[
2
]
=
0.0
;
tagint
tmptag
=
-
1
;
if
(
i
>=
0
)
{
double
rsq
=
1.1
;
double
rx
,
ry
,
rz
;
rx
=
ry
=
rz
=
0.0
;
double
coord
[
3
];
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_unequal
->
uniform
()
-
1.0
;
ry
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rz
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
x
[
i
][
0
]
+
displace
*
rx
;
coord
[
1
]
=
x
[
i
][
1
]
+
displace
*
ry
;
coord
[
2
]
=
x
[
i
][
2
]
+
displace
*
rz
;
if
(
regionflag
)
{
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_unequal
->
uniform
()
-
1.0
;
ry
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rz
=
2
*
random_unequal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
x
[
i
][
0
]
+
displace
*
rx
;
coord
[
1
]
=
x
[
i
][
1
]
+
displace
*
ry
;
coord
[
2
]
=
x
[
i
][
2
]
+
displace
*
rz
;
}
}
if
(
!
domain
->
inside_nonperiodic
(
coord
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
xtmp
[
0
]
=
x
[
i
][
0
];
xtmp
[
1
]
=
x
[
i
][
1
];
xtmp
[
2
]
=
x
[
i
][
2
];
x
[
i
][
0
]
=
coord
[
0
];
x
[
i
][
1
]
=
coord
[
1
];
x
[
i
][
2
]
=
coord
[
2
];
tmptag
=
atom
->
tag
[
i
];
}
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before
-
energy_after
)))
{
energy_stored
=
energy_after
;
ntranslation_successes
+=
1.0
;
}
else
{
tagint
tmptag_all
;
MPI_Allreduce
(
&
tmptag
,
&
tmptag_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
double
xtmp_all
[
3
];
MPI_Allreduce
(
&
xtmp
,
&
xtmp_all
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
{
if
(
tmptag_all
==
atom
->
tag
[
i
])
{
x
[
i
][
0
]
=
xtmp_all
[
0
];
x
[
i
][
1
]
=
xtmp_all
[
1
];
x
[
i
][
2
]
=
xtmp_all
[
2
];
}
}
energy_stored
=
energy_before
;
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_deletion_full
()
{
double
q_tmp
;
const
int
q_flag
=
atom
->
q_flag
;
ndeletion_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
double
energy_before
=
energy_stored
;
const
int
i
=
pick_random_gas_atom
();
int
tmpmask
;
if
(
i
>=
0
)
{
tmpmask
=
atom
->
mask
[
i
];
atom
->
mask
[
i
]
=
exclusion_group_bit
;
if
(
q_flag
)
{
q_tmp
=
atom
->
q
[
i
];
atom
->
q
[
i
]
=
0.0
;
}
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
ngas
*
exp
(
beta
*
(
energy_before
-
energy_after
))
/
(
zz
*
volume
))
{
if
(
i
>=
0
)
{
atom
->
avec
->
copy
(
atom
->
nlocal
-
1
,
i
,
1
);
atom
->
nlocal
--
;
}
atom
->
natoms
--
;
if
(
atom
->
map_style
)
atom
->
map_init
();
ndeletion_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
if
(
i
>=
0
)
{
atom
->
mask
[
i
]
=
tmpmask
;
if
(
q_flag
)
atom
->
q
[
i
]
=
q_tmp
;
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
energy_stored
=
energy_before
;
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_atomic_insertion_full
()
{
double
lamda
[
3
];
ninsertion_attempts
+=
1.0
;
double
energy_before
=
energy_stored
;
double
coord
[
3
];
if
(
regionflag
)
{
int
region_attempt
=
0
;
coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
region_attempt
++
;
if
(
region_attempt
>=
max_region_attempts
)
return
;
}
if
(
triclinic
)
domain
->
x2lamda
(
coord
,
lamda
);
}
else
{
if
(
triclinic
==
0
)
{
coord
[
0
]
=
xlo
+
random_equal
->
uniform
()
*
(
xhi
-
xlo
);
coord
[
1
]
=
ylo
+
random_equal
->
uniform
()
*
(
yhi
-
ylo
);
coord
[
2
]
=
zlo
+
random_equal
->
uniform
()
*
(
zhi
-
zlo
);
}
else
{
lamda
[
0
]
=
random_equal
->
uniform
();
lamda
[
1
]
=
random_equal
->
uniform
();
lamda
[
2
]
=
random_equal
->
uniform
();
// wasteful, but necessary
if
(
lamda
[
0
]
==
1.0
)
lamda
[
0
]
=
0.0
;
if
(
lamda
[
1
]
==
1.0
)
lamda
[
1
]
=
0.0
;
if
(
lamda
[
2
]
==
1.0
)
lamda
[
2
]
=
0.0
;
domain
->
lamda2x
(
lamda
,
coord
);
}
}
int
proc_flag
=
0
;
if
(
triclinic
==
0
)
{
domain
->
remap
(
coord
);
if
(
!
domain
->
inside
(
coord
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
if
(
coord
[
0
]
>=
sublo
[
0
]
&&
coord
[
0
]
<
subhi
[
0
]
&&
coord
[
1
]
>=
sublo
[
1
]
&&
coord
[
1
]
<
subhi
[
1
]
&&
coord
[
2
]
>=
sublo
[
2
]
&&
coord
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
else
{
if
(
lamda
[
0
]
>=
sublo
[
0
]
&&
lamda
[
0
]
<
subhi
[
0
]
&&
lamda
[
1
]
>=
sublo
[
1
]
&&
lamda
[
1
]
<
subhi
[
1
]
&&
lamda
[
2
]
>=
sublo
[
2
]
&&
lamda
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
if
(
proc_flag
)
{
atom
->
avec
->
create_atom
(
ngcmc_type
,
coord
);
int
m
=
atom
->
nlocal
-
1
;
// add to groups
// optionally add to type-based groups
atom
->
mask
[
m
]
=
groupbitall
;
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
{
if
(
ngcmc_type
==
grouptypes
[
igroup
])
atom
->
mask
[
m
]
|=
grouptypebits
[
igroup
];
}
atom
->
v
[
m
][
0
]
=
random_unequal
->
gaussian
()
*
sigma
;
atom
->
v
[
m
][
1
]
=
random_unequal
->
gaussian
()
*
sigma
;
atom
->
v
[
m
][
2
]
=
random_unequal
->
gaussian
()
*
sigma
;
if
(
charge_flag
)
atom
->
q
[
m
]
=
charge
;
modify
->
create_attribute
(
m
);
}
atom
->
natoms
++
;
if
(
atom
->
tag_enable
)
{
atom
->
tag_extend
();
if
(
atom
->
map_style
)
atom
->
map_init
();
}
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
zz
*
volume
*
exp
(
beta
*
(
energy_before
-
energy_after
))
/
(
ngas
+
1
))
{
ninsertion_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
atom
->
natoms
--
;
if
(
proc_flag
)
atom
->
nlocal
--
;
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
energy_stored
=
energy_before
;
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_translation_full
()
{
ntranslation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
translation_molecule
=
pick_random_gas_molecule
();
if
(
translation_molecule
==
-
1
)
return
;
double
energy_before
=
energy_stored
;
double
**
x
=
atom
->
x
;
double
rx
,
ry
,
rz
;
double
com_displace
[
3
],
coord
[
3
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_equal
->
uniform
()
-
1.0
;
ry
=
2
*
random_equal
->
uniform
()
-
1.0
;
rz
=
2
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
com_displace
[
0
]
=
displace
*
rx
;
com_displace
[
1
]
=
displace
*
ry
;
com_displace
[
2
]
=
displace
*
rz
;
int
nlocal
=
atom
->
nlocal
;
if
(
regionflag
)
{
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
mask
[
i
]
|=
molecule_group_bit
;
}
else
{
mask
[
i
]
&=
molecule_group_inversebit
;
}
}
double
com
[
3
];
com
[
0
]
=
com
[
1
]
=
com
[
2
]
=
0.0
;
group
->
xcm
(
molecule_group
,
gas_mass
,
com
);
coord
[
0
]
=
com
[
0
]
+
displace
*
rx
;
coord
[
1
]
=
com
[
1
]
+
displace
*
ry
;
coord
[
2
]
=
com
[
2
]
+
displace
*
rz
;
while
(
domain
->
regions
[
iregion
]
->
match
(
coord
[
0
],
coord
[
1
],
coord
[
2
])
==
0
)
{
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
rx
=
2
*
random_equal
->
uniform
()
-
1.0
;
ry
=
2
*
random_equal
->
uniform
()
-
1.0
;
rz
=
2
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
rx
*
rx
+
ry
*
ry
+
rz
*
rz
;
}
coord
[
0
]
=
com
[
0
]
+
displace
*
rx
;
coord
[
1
]
=
com
[
1
]
+
displace
*
ry
;
coord
[
2
]
=
com
[
2
]
+
displace
*
rz
;
}
com_displace
[
0
]
=
displace
*
rx
;
com_displace
[
1
]
=
displace
*
ry
;
com_displace
[
2
]
=
displace
*
rz
;
}
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
x
[
i
][
0
]
+=
com_displace
[
0
];
x
[
i
][
1
]
+=
com_displace
[
1
];
x
[
i
][
2
]
+=
com_displace
[
2
];
if
(
!
domain
->
inside_nonperiodic
(
x
[
i
]))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
}
}
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before
-
energy_after
)))
{
ntranslation_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
energy_stored
=
energy_before
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
translation_molecule
)
{
x
[
i
][
0
]
-=
com_displace
[
0
];
x
[
i
][
1
]
-=
com_displace
[
1
];
x
[
i
][
2
]
-=
com_displace
[
2
];
}
}
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_rotation_full
()
{
nrotation_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
rotation_molecule
=
pick_random_gas_molecule
();
if
(
rotation_molecule
==
-
1
)
return
;
double
energy_before
=
energy_stored
;
int
nlocal
=
atom
->
nlocal
;
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
rotation_molecule
)
{
mask
[
i
]
|=
molecule_group_bit
;
}
else
{
mask
[
i
]
&=
molecule_group_inversebit
;
}
}
double
com
[
3
];
com
[
0
]
=
com
[
1
]
=
com
[
2
]
=
0.0
;
group
->
xcm
(
molecule_group
,
gas_mass
,
com
);
// generate point in unit cube
// then restrict to unit sphere
double
r
[
3
],
rotmat
[
3
][
3
],
quat
[
4
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
r
[
0
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
1
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
2
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
MathExtra
::
dot3
(
r
,
r
);
}
double
theta
=
random_equal
->
uniform
()
*
max_rotation_angle
;
MathExtra
::
norm3
(
r
);
MathExtra
::
axisangle_to_quat
(
r
,
theta
,
quat
);
MathExtra
::
quat_to_mat
(
quat
,
rotmat
);
double
**
x
=
atom
->
x
;
imageint
*
image
=
atom
->
image
;
imageint
image_orig
[
natoms_per_molecule
];
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
molecule_group_bit
)
{
atom_coord
[
n
][
0
]
=
x
[
i
][
0
];
atom_coord
[
n
][
1
]
=
x
[
i
][
1
];
atom_coord
[
n
][
2
]
=
x
[
i
][
2
];
image_orig
[
n
]
=
image
[
i
];
double
xtmp
[
3
];
domain
->
unmap
(
x
[
i
],
image
[
i
],
xtmp
);
xtmp
[
0
]
-=
com
[
0
];
xtmp
[
1
]
-=
com
[
1
];
xtmp
[
2
]
-=
com
[
2
];
MathExtra
::
matvec
(
rotmat
,
xtmp
,
x
[
i
]);
x
[
i
][
0
]
+=
com
[
0
];
x
[
i
][
1
]
+=
com
[
1
];
x
[
i
][
2
]
+=
com
[
2
];
image
[
i
]
=
imagezero
;
domain
->
remap
(
x
[
i
],
image
[
i
]);
if
(
!
domain
->
inside
(
x
[
i
]))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
n
++
;
}
}
double
energy_after
=
energy_full
();
if
(
random_equal
->
uniform
()
<
exp
(
beta
*
(
energy_before
-
energy_after
)))
{
nrotation_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
energy_stored
=
energy_before
;
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
molecule_group_bit
)
{
x
[
i
][
0
]
=
atom_coord
[
n
][
0
];
x
[
i
][
1
]
=
atom_coord
[
n
][
1
];
x
[
i
][
2
]
=
atom_coord
[
n
][
2
];
image
[
i
]
=
image_orig
[
n
];
n
++
;
}
}
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_deletion_full
()
{
ndeletion_attempts
+=
1.0
;
if
(
ngas
==
0
)
return
;
tagint
deletion_molecule
=
pick_random_gas_molecule
();
if
(
deletion_molecule
==
-
1
)
return
;
double
energy_before
=
energy_stored
;
int
m
=
0
;
double
q_tmp
[
natoms_per_molecule
];
int
tmpmask
[
atom
->
nlocal
];
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
deletion_molecule
)
{
tmpmask
[
i
]
=
atom
->
mask
[
i
];
atom
->
mask
[
i
]
=
exclusion_group_bit
;
toggle_intramolecular
(
i
);
if
(
atom
->
q_flag
)
{
q_tmp
[
m
]
=
atom
->
q
[
i
];
m
++
;
atom
->
q
[
i
]
=
0.0
;
}
}
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
double
energy_after
=
energy_full
();
// energy_before corrected by energy_intra
double
deltaphi
=
ngas
*
exp
(
beta
*
((
energy_before
-
energy_intra
)
-
energy_after
))
/
(
zz
*
volume
*
natoms_per_molecule
);
if
(
random_equal
->
uniform
()
<
deltaphi
)
{
int
i
=
0
;
while
(
i
<
atom
->
nlocal
)
{
if
(
atom
->
molecule
[
i
]
==
deletion_molecule
)
{
atom
->
avec
->
copy
(
atom
->
nlocal
-
1
,
i
,
1
);
atom
->
nlocal
--
;
}
else
i
++
;
}
atom
->
natoms
-=
natoms_per_molecule
;
if
(
atom
->
map_style
)
atom
->
map_init
();
ndeletion_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
energy_stored
=
energy_before
;
int
m
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
{
if
(
atom
->
molecule
[
i
]
==
deletion_molecule
)
{
atom
->
mask
[
i
]
=
tmpmask
[
i
];
toggle_intramolecular
(
i
);
if
(
atom
->
q_flag
)
{
atom
->
q
[
i
]
=
q_tmp
[
m
];
m
++
;
}
}
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
attempt_molecule_insertion_full
()
{
double
lamda
[
3
];
ninsertion_attempts
+=
1.0
;
double
energy_before
=
energy_stored
;
tagint
maxmol
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
maxmol
=
MAX
(
maxmol
,
atom
->
molecule
[
i
]);
tagint
maxmol_all
;
MPI_Allreduce
(
&
maxmol
,
&
maxmol_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
maxmol_all
++
;
if
(
maxmol_all
>=
MAXTAGINT
)
error
->
all
(
FLERR
,
"Fix gcmc ran out of available molecule IDs"
);
int
insertion_molecule
=
maxmol_all
;
tagint
maxtag
=
0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
maxtag
=
MAX
(
maxtag
,
atom
->
tag
[
i
]);
tagint
maxtag_all
;
MPI_Allreduce
(
&
maxtag
,
&
maxtag_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
int
nlocalprev
=
atom
->
nlocal
;
double
com_coord
[
3
];
if
(
regionflag
)
{
int
region_attempt
=
0
;
com_coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
com_coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
com_coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
while
(
domain
->
regions
[
iregion
]
->
match
(
com_coord
[
0
],
com_coord
[
1
],
com_coord
[
2
])
==
0
)
{
com_coord
[
0
]
=
region_xlo
+
random_equal
->
uniform
()
*
(
region_xhi
-
region_xlo
);
com_coord
[
1
]
=
region_ylo
+
random_equal
->
uniform
()
*
(
region_yhi
-
region_ylo
);
com_coord
[
2
]
=
region_zlo
+
random_equal
->
uniform
()
*
(
region_zhi
-
region_zlo
);
region_attempt
++
;
if
(
region_attempt
>=
max_region_attempts
)
return
;
}
if
(
triclinic
)
domain
->
x2lamda
(
com_coord
,
lamda
);
}
else
{
if
(
triclinic
==
0
)
{
com_coord
[
0
]
=
xlo
+
random_equal
->
uniform
()
*
(
xhi
-
xlo
);
com_coord
[
1
]
=
ylo
+
random_equal
->
uniform
()
*
(
yhi
-
ylo
);
com_coord
[
2
]
=
zlo
+
random_equal
->
uniform
()
*
(
zhi
-
zlo
);
}
else
{
lamda
[
0
]
=
random_equal
->
uniform
();
lamda
[
1
]
=
random_equal
->
uniform
();
lamda
[
2
]
=
random_equal
->
uniform
();
// wasteful, but necessary
if
(
lamda
[
0
]
==
1.0
)
lamda
[
0
]
=
0.0
;
if
(
lamda
[
1
]
==
1.0
)
lamda
[
1
]
=
0.0
;
if
(
lamda
[
2
]
==
1.0
)
lamda
[
2
]
=
0.0
;
domain
->
lamda2x
(
lamda
,
com_coord
);
}
}
// generate point in unit cube
// then restrict to unit sphere
double
r
[
3
],
rotmat
[
3
][
3
],
quat
[
4
];
double
rsq
=
1.1
;
while
(
rsq
>
1.0
)
{
r
[
0
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
1
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
r
[
2
]
=
2.0
*
random_equal
->
uniform
()
-
1.0
;
rsq
=
MathExtra
::
dot3
(
r
,
r
);
}
double
theta
=
random_equal
->
uniform
()
*
MY_2PI
;
MathExtra
::
norm3
(
r
);
MathExtra
::
axisangle_to_quat
(
r
,
theta
,
quat
);
MathExtra
::
quat_to_mat
(
quat
,
rotmat
);
double
vnew
[
3
];
vnew
[
0
]
=
random_equal
->
gaussian
()
*
sigma
;
vnew
[
1
]
=
random_equal
->
gaussian
()
*
sigma
;
vnew
[
2
]
=
random_equal
->
gaussian
()
*
sigma
;
for
(
int
i
=
0
;
i
<
natoms_per_molecule
;
i
++
)
{
double
xtmp
[
3
];
MathExtra
::
matvec
(
rotmat
,
onemols
[
imol
]
->
x
[
i
],
xtmp
);
xtmp
[
0
]
+=
com_coord
[
0
];
xtmp
[
1
]
+=
com_coord
[
1
];
xtmp
[
2
]
+=
com_coord
[
2
];
// need to adjust image flags in remap()
imageint
imagetmp
=
imagezero
;
domain
->
remap
(
xtmp
,
imagetmp
);
if
(
!
domain
->
inside
(
xtmp
))
error
->
one
(
FLERR
,
"Fix gcmc put atom outside box"
);
int
proc_flag
=
0
;
if
(
triclinic
==
0
)
{
if
(
xtmp
[
0
]
>=
sublo
[
0
]
&&
xtmp
[
0
]
<
subhi
[
0
]
&&
xtmp
[
1
]
>=
sublo
[
1
]
&&
xtmp
[
1
]
<
subhi
[
1
]
&&
xtmp
[
2
]
>=
sublo
[
2
]
&&
xtmp
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
else
{
domain
->
x2lamda
(
xtmp
,
lamda
);
if
(
lamda
[
0
]
>=
sublo
[
0
]
&&
lamda
[
0
]
<
subhi
[
0
]
&&
lamda
[
1
]
>=
sublo
[
1
]
&&
lamda
[
1
]
<
subhi
[
1
]
&&
lamda
[
2
]
>=
sublo
[
2
]
&&
lamda
[
2
]
<
subhi
[
2
])
proc_flag
=
1
;
}
if
(
proc_flag
)
{
atom
->
avec
->
create_atom
(
onemols
[
imol
]
->
type
[
i
],
xtmp
);
int
m
=
atom
->
nlocal
-
1
;
// add to groups
// optionally add to type-based groups
atom
->
mask
[
m
]
=
groupbitall
;
for
(
int
igroup
=
0
;
igroup
<
ngrouptypes
;
igroup
++
)
{
if
(
ngcmc_type
==
grouptypes
[
igroup
])
atom
->
mask
[
m
]
|=
grouptypebits
[
igroup
];
}
atom
->
image
[
m
]
=
imagetmp
;
atom
->
molecule
[
m
]
=
insertion_molecule
;
if
(
maxtag_all
+
i
+
1
>=
MAXTAGINT
)
error
->
all
(
FLERR
,
"Fix gcmc ran out of available atom IDs"
);
atom
->
tag
[
m
]
=
maxtag_all
+
i
+
1
;
atom
->
v
[
m
][
0
]
=
vnew
[
0
];
atom
->
v
[
m
][
1
]
=
vnew
[
1
];
atom
->
v
[
m
][
2
]
=
vnew
[
2
];
atom
->
add_molecule_atom
(
onemols
[
imol
],
i
,
m
,
maxtag_all
);
modify
->
create_attribute
(
m
);
}
}
if
(
shakeflag
)
fixshake
->
set_molecule
(
nlocalprev
,
maxtag_all
,
imol
,
com_coord
,
vnew
,
quat
);
atom
->
natoms
+=
natoms_per_molecule
;
if
(
atom
->
natoms
<
0
)
error
->
all
(
FLERR
,
"Too many total atoms"
);
atom
->
nbonds
+=
onemols
[
imol
]
->
nbonds
;
atom
->
nangles
+=
onemols
[
imol
]
->
nangles
;
atom
->
ndihedrals
+=
onemols
[
imol
]
->
ndihedrals
;
atom
->
nimpropers
+=
onemols
[
imol
]
->
nimpropers
;
if
(
atom
->
map_style
)
atom
->
map_init
();
atom
->
nghost
=
0
;
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
double
energy_after
=
energy_full
();
// energy_after corrected by energy_intra
double
deltaphi
=
zz
*
volume
*
natoms_per_molecule
*
exp
(
beta
*
(
energy_before
-
(
energy_after
-
energy_intra
)))
/
(
ngas
+
natoms_per_molecule
);
if
(
random_equal
->
uniform
()
<
deltaphi
)
{
ninsertion_successes
+=
1.0
;
energy_stored
=
energy_after
;
}
else
{
atom
->
nbonds
-=
onemols
[
imol
]
->
nbonds
;
atom
->
nangles
-=
onemols
[
imol
]
->
nangles
;
atom
->
ndihedrals
-=
onemols
[
imol
]
->
ndihedrals
;
atom
->
nimpropers
-=
onemols
[
imol
]
->
nimpropers
;
atom
->
natoms
-=
natoms_per_molecule
;
energy_stored
=
energy_before
;
int
i
=
0
;
while
(
i
<
atom
->
nlocal
)
{
if
(
atom
->
molecule
[
i
]
==
insertion_molecule
)
{
atom
->
avec
->
copy
(
atom
->
nlocal
-
1
,
i
,
1
);
atom
->
nlocal
--
;
}
else
i
++
;
}
if
(
force
->
kspace
)
force
->
kspace
->
qsum_qsq
();
}
update_gas_atoms_list
();
}
/* ----------------------------------------------------------------------
compute particle's interaction energy with the rest of the system
------------------------------------------------------------------------- */
double
FixGCMC
::
energy
(
int
i
,
int
itype
,
tagint
imolecule
,
double
*
coord
)
{
double
delx
,
dely
,
delz
,
rsq
;
double
**
x
=
atom
->
x
;
int
*
type
=
atom
->
type
;
tagint
*
molecule
=
atom
->
molecule
;
int
nall
=
atom
->
nlocal
+
atom
->
nghost
;
pair
=
force
->
pair
;
cutsq
=
force
->
pair
->
cutsq
;
double
fpair
=
0.0
;
double
factor_coul
=
1.0
;
double
factor_lj
=
1.0
;
double
total_energy
=
0.0
;
for
(
int
j
=
0
;
j
<
nall
;
j
++
)
{
if
(
i
==
j
)
continue
;
if
(
mode
==
MOLECULE
)
if
(
imolecule
==
molecule
[
j
])
continue
;
delx
=
coord
[
0
]
-
x
[
j
][
0
];
dely
=
coord
[
1
]
-
x
[
j
][
1
];
delz
=
coord
[
2
]
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
int
jtype
=
type
[
j
];
if
(
rsq
<
cutsq
[
itype
][
jtype
])
total_energy
+=
pair
->
single
(
i
,
j
,
itype
,
jtype
,
rsq
,
factor_coul
,
factor_lj
,
fpair
);
}
return
total_energy
;
}
/* ----------------------------------------------------------------------
compute the energy of the given gas molecule in its current position
sum across all procs that own atoms of the given molecule
------------------------------------------------------------------------- */
double
FixGCMC
::
molecule_energy
(
tagint
gas_molecule_id
)
{
double
mol_energy
=
0.0
;
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
if
(
atom
->
molecule
[
i
]
==
gas_molecule_id
)
{
mol_energy
+=
energy
(
i
,
atom
->
type
[
i
],
gas_molecule_id
,
atom
->
x
[
i
]);
}
double
mol_energy_sum
=
0.0
;
MPI_Allreduce
(
&
mol_energy
,
&
mol_energy_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
return
mol_energy_sum
;
}
/* ----------------------------------------------------------------------
compute system potential energy
------------------------------------------------------------------------- */
double
FixGCMC
::
energy_full
()
{
if
(
triclinic
)
domain
->
x2lamda
(
atom
->
nlocal
);
domain
->
pbc
();
comm
->
exchange
();
atom
->
nghost
=
0
;
comm
->
borders
();
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
+
atom
->
nghost
);
if
(
modify
->
n_pre_neighbor
)
modify
->
pre_neighbor
();
neighbor
->
build
();
int
eflag
=
1
;
int
vflag
=
0
;
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
FixGCMC
::
pick_random_gas_atom
()
{
int
i
=
-
1
;
int
iwhichglobal
=
static_cast
<
int
>
(
ngas
*
random_equal
->
uniform
());
if
((
iwhichglobal
>=
ngas_before
)
&&
(
iwhichglobal
<
ngas_before
+
ngas_local
))
{
int
iwhichlocal
=
iwhichglobal
-
ngas_before
;
i
=
local_gas_list
[
iwhichlocal
];
}
return
i
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
tagint
FixGCMC
::
pick_random_gas_molecule
()
{
int
iwhichglobal
=
static_cast
<
int
>
(
ngas
*
random_equal
->
uniform
());
tagint
gas_molecule_id
=
0
;
if
((
iwhichglobal
>=
ngas_before
)
&&
(
iwhichglobal
<
ngas_before
+
ngas_local
))
{
int
iwhichlocal
=
iwhichglobal
-
ngas_before
;
int
i
=
local_gas_list
[
iwhichlocal
];
gas_molecule_id
=
atom
->
molecule
[
i
];
}
tagint
gas_molecule_id_all
=
0
;
MPI_Allreduce
(
&
gas_molecule_id
,
&
gas_molecule_id_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
return
gas_molecule_id_all
;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void
FixGCMC
::
toggle_intramolecular
(
int
i
)
{
if
(
atom
->
avec
->
bonds_allow
)
for
(
int
m
=
0
;
m
<
atom
->
num_bond
[
i
];
m
++
)
atom
->
bond_type
[
i
][
m
]
=
-
atom
->
bond_type
[
i
][
m
];
if
(
atom
->
avec
->
angles_allow
)
for
(
int
m
=
0
;
m
<
atom
->
num_angle
[
i
];
m
++
)
atom
->
angle_type
[
i
][
m
]
=
-
atom
->
angle_type
[
i
][
m
];
if
(
atom
->
avec
->
dihedrals_allow
)
for
(
int
m
=
0
;
m
<
atom
->
num_dihedral
[
i
];
m
++
)
atom
->
dihedral_type
[
i
][
m
]
=
-
atom
->
dihedral_type
[
i
][
m
];
if
(
atom
->
avec
->
impropers_allow
)
for
(
int
m
=
0
;
m
<
atom
->
num_improper
[
i
];
m
++
)
atom
->
improper_type
[
i
][
m
]
=
-
atom
->
improper_type
[
i
][
m
];
}
/* ----------------------------------------------------------------------
update the list of gas atoms
------------------------------------------------------------------------- */
void
FixGCMC
::
update_gas_atoms_list
()
{
int
nlocal
=
atom
->
nlocal
;
int
*
mask
=
atom
->
mask
;
tagint
*
molecule
=
atom
->
molecule
;
double
**
x
=
atom
->
x
;
if
(
atom
->
nmax
>
gcmc_nmax
)
{
memory
->
sfree
(
local_gas_list
);
gcmc_nmax
=
atom
->
nmax
;
local_gas_list
=
(
int
*
)
memory
->
smalloc
(
gcmc_nmax
*
sizeof
(
int
),
"GCMC:local_gas_list"
);
}
ngas_local
=
0
;
if
(
regionflag
)
{
if
(
mode
==
MOLECULE
)
{
tagint
maxmol
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
maxmol
=
MAX
(
maxmol
,
molecule
[
i
]);
tagint
maxmol_all
;
MPI_Allreduce
(
&
maxmol
,
&
maxmol_all
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
double
comx
[
maxmol_all
];
double
comy
[
maxmol_all
];
double
comz
[
maxmol_all
];
for
(
int
imolecule
=
0
;
imolecule
<
maxmol_all
;
imolecule
++
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
molecule
[
i
]
==
imolecule
)
{
mask
[
i
]
|=
molecule_group_bit
;
}
else
{
mask
[
i
]
&=
molecule_group_inversebit
;
}
}
double
com
[
3
];
com
[
0
]
=
com
[
1
]
=
com
[
2
]
=
0.0
;
group
->
xcm
(
molecule_group
,
gas_mass
,
com
);
comx
[
imolecule
]
=
com
[
0
];
comy
[
imolecule
]
=
com
[
1
];
comz
[
imolecule
]
=
com
[
2
];
}
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
if
(
domain
->
regions
[
iregion
]
->
match
(
comx
[
molecule
[
i
]],
comy
[
molecule
[
i
]],
comz
[
molecule
[
i
]])
==
1
)
{
local_gas_list
[
ngas_local
]
=
i
;
ngas_local
++
;
}
}
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
if
(
domain
->
regions
[
iregion
]
->
match
(
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
])
==
1
)
{
local_gas_list
[
ngas_local
]
=
i
;
ngas_local
++
;
}
}
}
}
}
else
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
local_gas_list
[
ngas_local
]
=
i
;
ngas_local
++
;
}
}
}
MPI_Allreduce
(
&
ngas_local
,
&
ngas
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Scan
(
&
ngas_local
,
&
ngas_before
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
ngas_before
-=
ngas_local
;
}
/* ----------------------------------------------------------------------
return acceptance ratios
------------------------------------------------------------------------- */
double
FixGCMC
::
compute_vector
(
int
n
)
{
if
(
n
==
0
)
return
ntranslation_attempts
;
if
(
n
==
1
)
return
ntranslation_successes
;
if
(
n
==
2
)
return
ninsertion_attempts
;
if
(
n
==
3
)
return
ninsertion_successes
;
if
(
n
==
4
)
return
ndeletion_attempts
;
if
(
n
==
5
)
return
ndeletion_successes
;
if
(
n
==
6
)
return
nrotation_attempts
;
if
(
n
==
7
)
return
nrotation_successes
;
return
0.0
;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double
FixGCMC
::
memory_usage
()
{
double
bytes
=
gcmc_nmax
*
sizeof
(
int
);
return
bytes
;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void
FixGCMC
::
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
FixGCMC
::
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