jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into collected-small-fixes

This commit is contained in:
Axel Kohlmeyer 2020-01-06 15:14:48 -05:00
parent a63792a07a 00dfde2952
commit fbc0b8a881
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
33 changed files with 684 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
doc/src/Errors_messages.rst
View File

@ -518,6 +518,16 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Bond/react: A deleted atom cannot remain bonded to an atom that is not deleted*
Self-explanatory.
*Bond/react: First neighbors of chiral atoms must be of mutually different types*
Self-explanatory.
*Bond/react: Chiral atoms must have exactly four first neighbors*
Self-explanatory.
*Bond/react: Molecule template 'Coords' section required for chiralIDs keyword*
The coordinates of atoms in the pre-reacted template are used to determine
chirality.
*Bond/react special bond generation overflow*
The number of special bonds per-atom created by a reaction exceeds the
system setting. See the read\_data or create\_box command for how to

28
doc/src/fix_bond_react.rst
View File

@ -20,9 +20,9 @@ Syntax
* the common keyword/values may be appended directly after 'bond/react'
* this applies to all reaction specifications (below)
* common\_keyword = *stabilization*
.. parsed-literal::
*stabilization* values = *no* or *yes* *group-ID* *xmax*
*no* = no reaction site stabilization
*yes* = perform reaction site stabilization
@ -40,9 +40,9 @@ Syntax
* map\_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates
* zero or more individual keyword/value pairs may be appended to each react argument
* individual\_keyword = *prob* or *max\_rxn* or *stabilize\_steps* or *update\_edges*
.. parsed-literal::
*prob* values = fraction seed
fraction = initiate reaction with this probability if otherwise eligible
seed = random number seed (positive integer)
@ -253,7 +253,7 @@ A discussion of correctly handling this is also provided on the
The map file is a text document with the following format:
A map file has a header and a body. The header of map file the
contains one mandatory keyword and four optional keywords. The
contains one mandatory keyword and five optional keywords. The
mandatory keyword is 'equivalences':
@ -269,10 +269,11 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and
N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template
N *deleteIDs* = # of atoms N that are specified for deletion
N *chiralIDs* = # of specified chiral centers N
N *customIDs* = # of atoms N that are specified for a custom update
N *constraints* = # of specified reaction constraints N
The body of the map file contains two mandatory sections and four
The body of the map file contains two mandatory sections and five
optional sections. The first mandatory section begins with the keyword
'BondingIDs' and lists the atom IDs of the bonding atom pair in the
pre-reacted molecule template. The second mandatory section begins
@ -284,12 +285,14 @@ molecule template. The first optional section begins with the keyword
'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
molecule template. The second optional section begins with the keyword
'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to
delete. The third optional section begins with the keyword 'Custom
delete. The third optional section begins with the keyword 'ChiralIDs'
lists the atom IDs of chiral atoms whose handedness should be
enforced. The fourth optional section begins with the keyword 'Custom
Edges' and allows for forcing the update of a specific atom's atomic
charge. The first column is the ID of an atom near the edge of the
pre-reacted molecule template, and the value of the second column is
either 'none' or 'charges.' Further details are provided in the
discussion of the 'update\_edges' keyword. The fourth optional section
discussion of the 'update\_edges' keyword. The fifth optional section
begins with the keyword 'Constraints' and lists additional criteria
that must be satisfied in order for the reaction to occur. Currently,
there are three types of constraints available, as discussed below.
@ -332,6 +335,15 @@ A sample map file is given below:
----------
The handedness of atoms that are chiral centers can be enforced by
listing their IDs in the ChiralIDs section. A chiral atom must be
bonded to four atoms with mutually different atom types. This feature
uses the coordinates and types of the involved atoms in the
pre-reaction template to determine handedness. Three atoms bonded to
the chiral center are arbitrarily chosen, to define an oriented plane,
and the relative position of the fourth bonded atom determines the
chiral center's handedness.
Any number of additional constraints may be specified in the
Constraints section of the map file. The constraint of type 'distance'
has syntax as follows:

13
doc/src/pair_meamc.rst
View File

@ -233,15 +233,20 @@ where
Cmin(I,J,K) = Cmin screening parameter when I-J pair is screened
by K (I<=J); default = 2.0
lattce(I,J) = lattice structure of I-J reference structure:
dia = diamond (interlaced fcc for alloy)
fcc = face centered cubic
bcc = body centered cubic
dim = dimer
b1 = rock salt (NaCl structure)
hcp = hexagonal close-packed
dim = dimer
dia = diamond (interlaced fcc for alloy)
dia3= diamond structure with primary 1NN and secondary 3NN interation
b1 = rock salt (NaCl structure)
c11 = MoSi2 structure
l12 = Cu3Au structure (lower case L, followed by 12)
b2 = CsCl structure (interpenetrating simple cubic)
ch4 = methane-like structure, only for binary system
lin = linear structure (180 degree angle)
zig = zigzag structure with a uniform angle
tri = H2O-like structure that has an angle
nn2(I,J) = turn on second-nearest neighbor MEAM formulation for
I-J pair (see for example :ref:`(Lee) <Lee>`).
0 = second-nearest neighbor formulation off
@ -254,6 +259,8 @@ where
zbl(I,J) = blend the MEAM I-J pair potential with the ZBL potential for small
atom separations :ref:`(ZBL) <ZBL>`
default = 1
theta(I,J) = angle between three atoms in line, zigzag, and trimer reference structures in degrees
default = 180
gsmooth_factor = factor determining the length of the G-function smoothing
region; only significant for ibar=0 or ibar=4.
99.0 = short smoothing region, sharp step

2
doc/txt/pair_meamc.txt
View File

@ -139,6 +139,8 @@ potential parameters:
lat = lattice structure of reference configuration
z = number of nearest neighbors in the reference structure
This field is only read for compatibility, the correct
value is inferred from the lattice structure
ielement = atomic number
atwt = atomic weight
alat = lattice constant of reference structure

2
lib/gpu/Makefile.lammps.standard
View File

@ -7,5 +7,5 @@ endif
gpu_SYSINC =
gpu_SYSLIB = -lcudart -lcuda
gpu_SYSPATH = -L$(CUDA_HOME)/lib64
gpu_SYSPATH = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs

2
lib/gpu/Makefile.linux
View File

@ -54,7 +54,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_SINGLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC) -Xcompiler -fPIC
CUDR_CPP = mpicxx -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1 -fPIC

2
lib/gpu/Makefile.linux.double
View File

@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_DOUBLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math
CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK

2
lib/gpu/Makefile.linux.mixed
View File

@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_SINGLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math
CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK

2
lib/gpu/Makefile.linux.single
View File

@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_SINGLE_SINGLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math
CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK

2
lib/gpu/Makefile.linux_multi
View File

@ -40,7 +40,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_SINGLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC) -Xcompiler "-fPIC -std=c++98"
CUDR_CPP = mpicxx -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1 -fPIC

2
lib/gpu/Makefile.mac
View File

@ -14,7 +14,7 @@ NVCC = nvcc -m64
CUDA_ARCH = -arch=sm_11
CUDA_PRECISION = -D_SINGLE_SINGLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib
CUDA_LIB = -L$(CUDA_HOME)/lib -L$(CUDA_HOME)/lib/stubs
CUDA_OPTS = -DUNIX -DUCL_NO_EXIT -O3 --use_fast_math
CUDR_CPP = mpic++ -m64

2
lib/gpu/Makefile.serial
View File

@ -34,7 +34,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_SINGLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L../../src/STUBS -lmpi_stubs
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs -L../../src/STUBS -lmpi_stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC)
CUDR_CPP = g++ -DMPI_GERYON -DUCL_NO_EXIT -fPIC -I../../src/STUBS

2
lib/gpu/Makefile.shannon
View File

@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
CUDA_PRECISION = -D_DOUBLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math
CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK

2
lib/gpu/Makefile.xk7
View File

@ -14,7 +14,7 @@ CUDA_ARCH = -arch=sm_35
CUDA_PRECISION = -D_SINGLE_DOUBLE
CUDA_INCLUDE = -I$(CUDA_HOME)/include
CUDA_LIB = -L$(CUDA_HOME)/lib64
CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
CUDA_OPTS = -DUNIX -O3 --use_fast_math
CUDR_CPP = CC -DCUDA_PROXY -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK

5
lib/gpu/README
View File

@ -3,7 +3,7 @@
--------------------------------
W. Michael Brown (ORNL)
Trung Dac Nguyen (ORNL)
Trung Dac Nguyen (ORNL/Northwestern)
Peng Wang (NVIDIA)
Axel Kohlmeyer (Temple)
Steve Plimpton (SNL)
@ -129,7 +129,8 @@ that ships with the CUDA toolkit, but also with the CUDA driver library
on the head node of a GPU cluster, this library may not be installed,
so you may need to copy it over from one of the compute nodes (best into
this directory). Recent CUDA toolkits starting from CUDA 9 provide a dummy
libcuda.so library, that can be used for linking (but not for running).
libcuda.so library (typically under $(CUDA_HOME)/lib64/stubs), that can be used for
linking.
The gpu library supports 3 precision modes as determined by
the CUDA_PRECISION variable:

4
src/.gitignore vendored
View File

@ -280,6 +280,8 @@
/bond_oxdna_fene.h
/bond_oxdna2_fene.cpp
/bond_oxdna2_fene.h
/bond_oxrna2_fene.cpp
/bond_oxrna2_fene.h
/bond_quartic.cpp
/bond_quartic.h
/bond_table.cpp
@ -982,6 +984,8 @@
/pair_oxdna_*.h
/pair_oxdna2_*.cpp
/pair_oxdna2_*.h
/pair_oxrna2_*.cpp
/pair_oxrna2_*.h
/mf_oxdna.h
/pair_peri_eps.cpp
/pair_peri_eps.h

4
src/USER-INTEL/npair_intel.cpp
View File

@ -360,7 +360,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
if (THREE) ttag[u] = tag[j];
}
if (FULL == 0 || TRI == 1) {
if (FULL == 0 && TRI != 1) {
icount = 0;
istart = ncount;
IP_PRE_edge_align(istart, sizeof(int));
@ -392,7 +392,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
// ---------------------- Loop over i bin
int n = 0;
if (FULL == 0 || TRI == 1) {
if (FULL == 0 && TRI != 1) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep

69
src/USER-MEAMC/meam.h
View File

@ -2,13 +2,15 @@
#define LMP_MEAM_H
#include <cmath>
#include <cstring>
#include "math_const.h"
#define maxelt 5
namespace LAMMPS_NS {
class Memory;
typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;
typedef enum { FCC, BCC, HCP, DIM, DIA, DIA3, B1, C11, L12, B2, CH4, LIN, ZIG, TRI } lattice_t;
class MEAM
{
@ -26,9 +28,7 @@ private:
// Ec_meam = cohesive energy
// re_meam = nearest-neighbor distance
// Omega_meam = atomic volume
// B_meam = bulk modulus
// Z_meam = number of first neighbors for reference structure
// ielt_meam = atomic number of element
// A_meam = adjustable parameter
// alpha_meam = sqrt(9*Omega*B/Ec)
@ -64,8 +64,11 @@ private:
// nr,dr = pair function discretization parameters
// nrar,rdrar = spline coeff array parameters
// theta = angle between three atoms in line, zigzag, and trimer reference structures
// stheta_meam = sin(theta/2) in radian used in line, zigzag, and trimer reference structures
// ctheta_meam = cos(theta/2) in radian used in line, zigzag, and trimer reference structures
double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
double Omega_meam[maxelt], Z_meam[maxelt];
double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];
double delta_meam[maxelt][maxelt];
double beta0_meam[maxelt], beta1_meam[maxelt];
@ -108,6 +111,10 @@ public:
int maxneigh;
double *scrfcn, *dscrfcn, *fcpair;
//angle for trimer, zigzag, line reference structures
double stheta_meam[maxelt][maxelt];
double ctheta_meam[maxelt][maxelt];
protected:
// meam_funcs.cpp
@ -191,9 +198,12 @@ protected:
double embedding(const double A, const double Ec, const double rhobar, double& dF) const;
static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);
static void get_shpfcn(const lattice_t latt, double (&s)[3]);
static int get_Zij(const lattice_t latt);
static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);
static void get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3]);
static int get_Zij2(const lattice_t latt, const double cmin, const double cmax,
const double sthe, double &a, double &S);
static int get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S);
protected:
void meam_checkindex(int, int, int, int*, int*);
void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
@ -204,6 +214,7 @@ protected:
void alloyparams();
void compute_pair_meam();
double phi_meam(double, int, int);
const double phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat);
void compute_reference_density();
void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double,
double, double, double, int, int, lattice_t);
@ -212,7 +223,41 @@ protected:
void interpolate_meam(int);
public:
void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
//-----------------------------------------------------------------------------
// convert lattice spec to lattice_t
// only use single-element lattices if single=true
// return false on failure
// return true and set lat on success
static bool str_to_lat(const char* str, bool single, lattice_t& lat)
{
if (strcmp(str,"fcc") == 0) lat = FCC;
else if (strcmp(str,"bcc") == 0) lat = BCC;
else if (strcmp(str,"hcp") == 0) lat = HCP;
else if (strcmp(str,"dim") == 0) lat = DIM;
else if (strcmp(str,"dia") == 0) lat = DIA;
else if (strcmp(str,"dia3") == 0) lat = DIA3;
else if (strcmp(str,"lin") == 0) lat = LIN;
else if (strcmp(str,"zig") == 0) lat = ZIG;
else if (strcmp(str,"tri") == 0) lat = TRI;
else {
if (single)
return false;
if (strcmp(str,"b1") == 0) lat = B1;
else if (strcmp(str,"c11") == 0) lat = C11;
else if (strcmp(str,"l12") == 0) lat = L12;
else if (strcmp(str,"b2") == 0) lat = B2;
else if (strcmp(str,"ch4") == 0) lat = CH4;
else if (strcmp(str,"lin") == 0) lat =LIN;
else if (strcmp(str,"zig") == 0) lat = ZIG;
else if (strcmp(str,"tri") == 0) lat = TRI;
else return false;
}
return true;
}
static int get_Zij(const lattice_t latt);
void meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* atwt, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar);
@ -222,9 +267,9 @@ public:
void meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset);
void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, int& errorflag);
double* eatom, int ntype, int* type, int* fmap, double** scale, int& errorflag);
void meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* eatom, int ntype, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom);
};
@ -234,6 +279,10 @@ static inline bool iszero(const double f) {
return fabs(f) < 1e-20;
}
static inline bool isone(const double f) {
return fabs(f-1.0) < 1e-20;
}
// Helper functions
static inline double fdiv_zero(const double n, const double d) {

11
src/USER-MEAMC/meam_dens_final.cpp
View File

@ -4,18 +4,20 @@ using namespace LAMMPS_NS;
void
MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
double* eatom, int /*ntype*/, int* type, int* fmap, int& errorflag)
double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, int& errorflag)
{
int i, elti;
int m;
double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;
double denom, rho_bkgd, Fl;
double scaleii;
// Complete the calculation of density
for (i = 0; i < nlocal; i++) {
elti = fmap[type[i]];
if (elti >= 0) {
scaleii = scale[type[i]][type[i]];
rho1[i] = 0.0;
rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];
rho3[i] = 0.0;
@ -52,12 +54,14 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
}
Z = this->Z_meam[elti];
Z = get_Zij(this->lattce_meam[elti][elti]);
G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
if (errorflag != 0)
return;
get_shpfcn(this->lattce_meam[elti][elti], shp);
get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp);
if (this->ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
@ -113,6 +117,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);
if (eflag_either != 0) {
Fl *= scaleii;
if (eflag_global != 0) {
*eng_vdwl = *eng_vdwl + Fl;
}

2
src/USER-MEAMC/meam_dens_init.cpp
View File

@ -197,7 +197,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
// Now compute derivatives
dscrfcn[jn] = 0.0;
sfcij = sij * fcij;
if (!iszero(sfcij) && !iszero(sfcij - 1.0)) {
if (!iszero(sfcij) && !isone(sfcij)) {
for (kn = 0; kn < numneigh_full; kn++) {
k = firstneigh_full[kn];
if (k == j) continue;

35
src/USER-MEAMC/meam_force.cpp
View File

@ -8,7 +8,7 @@ using namespace LAMMPS_NS;
void
MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,
double* eatom, int /*ntype*/, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom)
{
int j, jn, k, kn, kk, m, n, p, q;
@ -17,7 +17,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
double v[6], fi[3], fj[3];
double third, sixth;
double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem;
double r, recip, phi, phip;
double recip, phi, phip;
double sij;
double a1, a1i, a1j, a2, a2i, a2j;
double a3i, a3j;
@ -42,6 +42,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;
double dsij1, dsij2, force1, force2;
double t1i, t2i, t3i, t1j, t2j, t3j;
double scaleij;
third = 1.0 / 3.0;
sixth = 1.0 / 6.0;
@ -59,6 +60,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
for (jn = 0; jn < numneigh; jn++) {
j = firstneigh[jn];
eltj = fmap[type[j]];
scaleij = scale[type[i]][type[j]];
if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {
@ -69,7 +71,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
if (rij2 < this->cutforcesq) {
rij = sqrt(rij2);
r = rij;
recip = 1.0 / rij;
// Compute phi and phip
ind = this->eltind[elti][eltj];
@ -78,17 +80,16 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
kk = std::min(kk, this->nrar - 2);
pp = pp - kk;
pp = std::min(pp, 1.0);
phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
this->phirar[ind][kk];
phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + this->phirar[ind][kk];
phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
recip = 1.0 / r;
if (eflag_either != 0) {
double phi_sc = phi * scaleij;
if (eflag_global != 0)
*eng_vdwl = *eng_vdwl + phi * sij;
*eng_vdwl = *eng_vdwl + phi_sc * sij;
if (eflag_atom != 0) {
eatom[i] = eatom[i] + 0.5 * phi * sij;
eatom[j] = eatom[j] + 0.5 * phi * sij;
eatom[i] = eatom[i] + 0.5 * phi_sc * sij;
eatom[j] = eatom[j] + 0.5 * phi_sc * sij;
}
}
@ -287,8 +288,9 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
}
// Compute derivatives of total density wrt rij, sij and rij(3)
get_shpfcn(this->lattce_meam[elti][elti], shpi);
get_shpfcn(this->lattce_meam[eltj][eltj], shpj);
get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi);
get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj);
drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
dt3dr1 * rho3[i] + t3i * drho3dr1) -
@ -379,6 +381,13 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
for (m = 0; m < 3; m++) {
dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
}
if (!isone(scaleij)) {
dUdrij *= scaleij;
dUdsij *= scaleij;
dUdrijm[0] *= scaleij;
dUdrijm[1] *= scaleij;
dUdrijm[2] *= scaleij;
}
// Add the part of the force due to dUdrij and dUdsij
@ -410,7 +419,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
// Now compute forces on other atoms k due to change in sij
if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop
if (iszero(sij) || isone(sij)) continue; //: cont jn loop
double dxik(0), dyik(0), dzik(0);
double dxjk(0), dyjk(0), dzjk(0);
@ -429,7 +438,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
dsij1 = 0.0;
dsij2 = 0.0;
if (!iszero(sij) && !iszero(sij - 1.0)) {
if (!iszero(sij) && !isone(sij)) {
const double rbound = rij2 * this->ebound_meam[elti][eltj];
delc = Cmax - Cmin;
dxjk = x[k][0] - x[j][0];

96
src/USER-MEAMC/meam_funcs.cpp
View File

@ -198,7 +198,7 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec
// Shape factors for various configurations
//
void
MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3])
{
switch (latt) {
case FCC:
@ -214,7 +214,9 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
s[1] = 0.0;
s[2] = 1.0 / 3.0;
break;
case CH4: // CH4 actually needs shape factor for diamond for C, dimer for H
case DIA:
case DIA3:
s[0] = 0.0;
s[1] = 0.0;
s[2] = 32.0 / 9.0;
@ -222,8 +224,20 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
case DIM:
s[0] = 1.0;
s[1] = 2.0 / 3.0;
// s(3) = 1.d0
s[2] = 0.40;
// s(3) = 1.d0 // this should be 0.4 unless (1-legendre) is multiplied in the density calc.
s[2] = 0.40; // this is (1-legendre) where legendre = 0.6 in dynamo is accounted.
break;
case LIN: //linear, theta being 180
s[0] = 0.0;
s[1] = 8.0 / 3.0; // 4*(co**4 + si**4 - 1.0/3.0) in zig become 4*(1-1/3)
s[2] = 0.0;
break;
case ZIG: //zig-zag
case TRI: //trimer e.g. H2O
s[0] = 4.0*pow(cthe,2);
s[1] = 4.0*(pow(cthe,4) + pow(sthe,4) - 1.0/3.0);
s[2] = 4.0*(pow(cthe,2) * (3*pow(sthe,4) + pow(cthe,4)));
s[2] = s[2] - 0.6*s[0]; //legend in dyn, 0.6 is default value.
break;
default:
s[0] = 0.0;
@ -232,7 +246,7 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
}
//-----------------------------------------------------------------------------
// Number of neighbors for the reference structure
// Number of first neighbors for reference structure
//
int
MEAM::get_Zij(const lattice_t latt)
@ -244,18 +258,25 @@ MEAM::get_Zij(const lattice_t latt)
return 8;
case HCP:
return 12;
case B1:
return 6;
case DIA:
case DIA3:
return 4;
case DIM:
return 1;
case B1:
return 6;
case C11:
return 10;
case L12:
return 12;
case B2:
return 8;
case CH4: // DYNAMO currenly implemented this way while it needs two Z values, 4 and 1
return 4;
case LIN:
case ZIG:
case TRI:
return 2;
// call error('Lattice not defined in get_Zij.')
}
return 0;
@ -263,11 +284,13 @@ MEAM::get_Zij(const lattice_t latt)
//-----------------------------------------------------------------------------
// Number of second neighbors for the reference structure
// a = distance ratio R1/R2
// S = second neighbor screening function
// a = distance ratio R1/R2 (a2nn in dynamo)
// numscr = number of atoms that screen the 2NN bond
// S = second neighbor screening function (xfac, a part of b2nn in dynamo)
//
int
MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, double& a, double& S)
MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax,
const double stheta, double& a, double& S)
{
double C, x, sijk;
@ -299,7 +322,7 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
numscr = 2;
break;
case DIA:
case DIA: // 2NN
Zij2 = 12;
a = sqrt(8.0 / 3.0);
numscr = 1;
@ -308,12 +331,35 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
}
break;
case DIA3: // 3NN
Zij2 = 12;
a = sqrt(11.0 / 3.0);
numscr = 4;
if (cmin < 0.500001) {
// call error('can not do 2NN MEAM for dia')
}
break;
case CH4: //does not have 2nn structure so it returns 0
case LIN: //line
case DIM:
// this really shouldn't be allowed; make sure screening is zero
a = 1.0;
S = 0.0;
return 0;
case TRI: //TRI
Zij2 = 1;
a = 2.0*stheta;
numscr=2;
break;
case ZIG: //zig-zag
Zij2 = 2;
a = 2.0*stheta;
numscr=1;
break;
case L12:
Zij2 = 6;
a = sqrt(2.0);
@ -335,10 +381,38 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
}
// Compute screening for each first neighbor
C = 4.0 / (a * a) - 1.0;
if (latt==DIA3){ // special case for 3NN diamond structure
C = 1.0;
} else {
C = 4.0 / (a * a) - 1.0;
}
x = (C - cmin) / (cmax - cmin);
sijk = fcut(x);
// There are numscr first neighbors screening the second neighbors
S = MathSpecial::powint(sijk, numscr);
return Zij2;
}
int
MEAM::get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S)
{
double x, sijk, C;
int numscr = 0, Zij2 = 0;
switch (latt) {
case ZIG: //zig-zag for b11s and b22s term
case TRI: //trimer for b11s
Zij2 = 2;
numscr=1;
break;
default:
// unknown lattic flag in get Zij
// call error('Lattice not defined in get_Zij.')
break;
}
C = 1.0;
x = (C - cmin) / (cmax - cmin);
sijk = fcut(x);
S = MathSpecial::powint(sijk, numscr);
return Zij2;
}

2
src/USER-MEAMC/meam_impl.cpp
View File

@ -38,7 +38,7 @@ MEAM::MEAM(Memory* mem)
neltypes = 0;
for (int i = 0; i < maxelt; i++) {
Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
A_meam[i] = rho0_meam[i] = beta0_meam[i] =
beta1_meam[i]= beta2_meam[i] = beta3_meam[i] =
t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] =
rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;

365
src/USER-MEAMC/meam_setup_done.cpp
View File

@ -103,6 +103,9 @@ MEAM::alloyparams(void)
this->alpha_meam[i][j] = this->alpha_meam[j][i];
this->lattce_meam[i][j] = this->lattce_meam[j][i];
this->nn2_meam[i][j] = this->nn2_meam[j][i];
// theta for lin,tri,zig references
this->stheta_meam[i][j] = this->stheta_meam[j][i];
this->ctheta_meam[i][j] = this->ctheta_meam[j][i];
// If i<j and term is unset, use default values (e.g. mean of i-i and
// j-j)
} else if (j > i) {
@ -161,10 +164,10 @@ void
MEAM::compute_pair_meam(void)
{
double r /*ununsed:, temp*/;
double r, b2nn, phi_val;
int j, a, b, nv2;
double astar, frac, phizbl;
int n, nmax, Z1, Z2;
int n, Z1, Z2;
double arat, rarat, scrn, scrn2;
double phiaa, phibb /*unused:,phitmp*/;
double C, s111, s112, s221, S11, S22;
@ -215,7 +218,7 @@ MEAM::compute_pair_meam(void)
if (this->nn2_meam[a][b] == 1) {
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b], arat, scrn);
this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);
// The B1, B2, and L12 cases with NN2 have a trick to them; we need to
// compute the contributions from second nearest neighbors, like a-a
@ -229,35 +232,26 @@ MEAM::compute_pair_meam(void)
phiaa = phi_meam(rarat, a, a);
Z1 = get_Zij(this->lattce_meam[a][a]);
Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);
}
}
this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
phiaa+= phi_meam_series(scrn, Z1, Z2, a, a, rarat, arat);
// phi_bb
phibb = phi_meam(rarat, b, b);
Z1 = get_Zij(this->lattce_meam[b][b]);
Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
this->Cmax_meam[b][b][b], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);
}
}
this->Cmax_meam[b][b][b], this->stheta_meam[b][b], arat, scrn);
phibb+= phi_meam_series(scrn, Z1, Z2, b, b, rarat, arat);
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
this->lattce_meam[a][b] == DIA) {
// Add contributions to the B1 or B2 potential
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b], arat, scrn);
this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a],
this->Cmax_meam[b][b][a], arat, scrn2);
this->Cmax_meam[b][b][a], this->stheta_meam[a][b], arat, scrn2);
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
} else if (this->lattce_meam[a][b] == L12) {
@ -279,11 +273,7 @@ MEAM::compute_pair_meam(void)
}
} else {
nmax = 10;
for (n = 1; n <= nmax; n++) {
this->phir[nv2][j] =
this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);
}
this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);
}
}
@ -328,11 +318,12 @@ MEAM::phi_meam(double r, int a, int b)
double rho02, rho12, rho22, rho32;
double scalfac, phiaa, phibb;
double Eu;
double arat, scrn /*unused:,scrn2*/;
double arat, scrn, scrn2;
int Z12, errorflag;
int n, nmax, Z1nn, Z2nn;
int n, Z1nn, Z2nn;
lattice_t latta /*unused:,lattb*/;
double rho_bkgd1, rho_bkgd2;
double b11s, b22s;
double phi_m = 0.0;
@ -341,6 +332,8 @@ MEAM::phi_meam(double r, int a, int b)
// get number of neighbors in the reference structure
// Nref[i][j] = # of i's neighbors of type j
Z1 = get_Zij(this->lattce_meam[a][a]);
Z2 = get_Zij(this->lattce_meam[b][b]);
Z12 = get_Zij(this->lattce_meam[a][b]);
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32);
@ -401,19 +394,17 @@ MEAM::phi_meam(double r, int a, int b)
// If using mixing rule for t, apply to reference structure; else
// use precomputed values
if (this->mix_ref_t == 1) {
Z1 = this->Z_meam[a];
Z2 = this->Z_meam[b];
if (this->ibar_meam[a] <= 0)
G1 = 1.0;
else {
get_shpfcn(this->lattce_meam[a][a], s1);
get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1);
Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);
G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
}
if (this->ibar_meam[b] <= 0)
G2 = 1.0;
else {
get_shpfcn(this->lattce_meam[b][b], s2);
get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b], s2);
Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);
G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
}
@ -439,8 +430,8 @@ MEAM::phi_meam(double r, int a, int b)
rho_bkgd2 = rho0_2;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a];
rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b];
rho_bkgd1 = this->rho0_meam[a] * Z1;
rho_bkgd2 = this->rho0_meam[b] * Z2;
} else {
rho_bkgd1 = this->rho_ref_meam[a];
rho_bkgd2 = this->rho_ref_meam[b];
@ -475,20 +466,44 @@ MEAM::phi_meam(double r, int a, int b)
// account for second neighbor a-a potential here...
Z1nn = get_Zij(this->lattce_meam[a][a]);
Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a);
}
}
this->Cmax_meam[a][a][a], this->stheta_meam[a][b], arat, scrn);
phiaa += phi_meam_series(scrn, Z1nn, Z2nn, a, a, r, arat);
phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
} else if (this->lattce_meam[a][b] == CH4) {
phi_m = (5 * Eu - F1 - 4*F2)/4;
} else if (this->lattce_meam[a][b] == ZIG){
if (a==b){
phi_m = (2 * Eu - F1 - F2) / Z12;
} else{
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
b11s = -Z2/Z1*scrn;
Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], scrn2);
b22s = -Z2/Z1*scrn2;
phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
phibb = phi_meam(2.0*this->stheta_meam[a][b]*r, b, b);
phi_m = (2.0*Eu - F1 - F2 + phiaa*b11s + phibb*b22s) / Z12;
}
} else if (this->lattce_meam[a][b] == TRI) {
if (a==b){
phi_m = (3.0*Eu - 2.0*F1 - F2) / Z12;
} else {
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
b11s = -Z2/Z1*scrn;
phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12;
}
} else {
//
// potential is computed from Rose function and embedding energy
phi_m = (2 * Eu - F1 - F2) / Z12;
//
}
// if r = 0, just return 0
@ -499,6 +514,31 @@ MEAM::phi_meam(double r, int a, int b)
return phi_m;
}
//----------------------------------------------------------------------c
// Compute 2NN series terms for phi
// To avoid nan values of phir due to rapid decrease of b2nn^n or/and
// argument of phi_meam, i.e. r*arat^n, in some cases (3NN dia with low Cmin value)
//
const double
MEAM::phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat)
{
double phi_sum = 0.0;
double b2nn, phi_val;
if (scrn > 0.0) {
b2nn = -Z2*scrn/Z1;
for (int n = 1; n <= 10; n++) {
phi_val = MathSpecial::powint(b2nn,n) * phi_meam(r * MathSpecial::powint(arat, n), a, b);
if (iszero(phi_val)) {
// once either term becomes zero at some point, all folliwng will also be zero
// necessary to avoid numerical error (nan or infty) due to exponential decay in phi_meam
break;
}
phi_sum += phi_val;
}
}
return phi_sum;
}
//----------------------------------------------------------------------c
// Compute background density for reference structure of each element
void
@ -510,11 +550,11 @@ MEAM::compute_reference_density(void)
// loop over element types
for (a = 0; a < this->neltypes; a++) {
Z = (int)this->Z_meam[a];
Z = get_Zij(this->lattce_meam[a][a]);
if (this->ibar_meam[a] <= 0)
Gbar = 1.0;
else {
get_shpfcn(this->lattce_meam[a][a], shp);
get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], shp);
gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);
Gbar = G_gam(gam, this->ibar_meam[a], errorflag);
}
@ -529,7 +569,7 @@ MEAM::compute_reference_density(void)
// screening)
if (this->nn2_meam[a][a] == 1) {
Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
rho0 = rho0 + Z2 * rho0_2nn * scrn;
}
@ -555,30 +595,43 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, dou
*t12av = t12;
*t22av = t22;
*t32av = t32;
} else if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || latt == B1 || latt == DIM || latt == B2) {
// all neighbors are of the opposite type
*t11av = t12;
*t21av = t22;
*t31av = t32;
*t12av = t11;
*t22av = t21;
*t32av = t31;
} else {
a1 = r / this->re_meam[a][a] - 1.0;
a2 = r / this->re_meam[b][b] - 1.0;
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
if (latt == L12) {
rho01 = 8 * rhoa01 + 4 * rhoa02;
*t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
} else switch (latt) {
case FCC:
case BCC:
case DIA:
case DIA3:
case HCP:
case B1:
case DIM:
case B2:
case CH4:
case LIN:
case ZIG:
case TRI:
// all neighbors are of the opposite type
*t11av = t12;
*t21av = t22;
*t31av = t32;
*t12av = t11;
*t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
*t22av = t21;
*t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
*t32av = t31;
} else {
// call error('Lattice not defined in get_tavref.')
}
break;
default:
a1 = r / this->re_meam[a][a] - 1.0;
a2 = r / this->re_meam[b][b] - 1.0;
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
if (latt == L12) {
rho01 = 8 * rhoa01 + 4 * rhoa02;
*t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
*t12av = t11;
*t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
*t22av = t21;
*t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
*t32av = t31;
} else {
// call error('Lattice not defined in get_tavref.')
}
}
}
@ -600,7 +653,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
double a1, a2;
double s[3];
lattice_t lat;
int Zij2nn;
int Zij,Zij2nn;
double rhoa01nn, rhoa02nn;
double rhoa01, rhoa11, rhoa21, rhoa31;
double rhoa02, rhoa12, rhoa22, rhoa32;
@ -621,6 +674,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
lat = this->lattce_meam[a][b];
Zij = get_Zij(lat);
*rho11 = 0.0;
*rho21 = 0.0;
*rho31 = 0.0;
@ -628,64 +683,133 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
*rho22 = 0.0;
*rho32 = 0.0;
if (lat == FCC) {
*rho01 = 12.0 * rhoa02;
*rho02 = 12.0 * rhoa01;
} else if (lat == BCC) {
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
} else if (lat == B1) {
*rho01 = 6.0 * rhoa02;
*rho02 = 6.0 * rhoa01;
} else if (lat == DIA) {
*rho01 = 4.0 * rhoa02;
*rho02 = 4.0 * rhoa01;
*rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
*rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
} else if (lat == HCP) {
*rho01 = 12 * rhoa02;
*rho02 = 12 * rhoa01;
*rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
*rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
} else if (lat == DIM) {
get_shpfcn(DIM, s);
*rho01 = rhoa02;
*rho02 = rhoa01;
*rho11 = s[0] * rhoa12 * rhoa12;
*rho12 = s[0] * rhoa11 * rhoa11;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho31 = s[2] * rhoa32 * rhoa32;
*rho32 = s[2] * rhoa31 * rhoa31;
} else if (lat == C11) {
*rho01 = rhoa01;
*rho02 = rhoa02;
*rho11 = rhoa11;
*rho12 = rhoa12;
*rho21 = rhoa21;
*rho22 = rhoa22;
*rho31 = rhoa31;
*rho32 = rhoa32;
} else if (lat == L12) {
*rho01 = 8 * rhoa01 + 4 * rhoa02;
*rho02 = 12 * rhoa01;
if (this->ialloy == 1) {
*rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]);
denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]);
if (denom > 0.)
*rho21 = *rho21 / denom * *rho01;
} else
*rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
} else if (lat == B2) {
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
} else {
switch (lat) {
case FCC:
*rho01 = 12.0 * rhoa02;
*rho02 = 12.0 * rhoa01;
break;
case BCC:
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
break;
case B1:
*rho01 = 6.0 * rhoa02;
*rho02 = 6.0 * rhoa01;
break;
case DIA:
case DIA3:
*rho01 = 4.0 * rhoa02;
*rho02 = 4.0 * rhoa01;
*rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
*rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
break;
case HCP:
*rho01 = 12 * rhoa02;
*rho02 = 12 * rhoa01;
*rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
*rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
break;
case DIM:
get_shpfcn(DIM, 0, 0, s);
*rho01 = rhoa02;
*rho02 = rhoa01;
*rho11 = s[0] * rhoa12 * rhoa12;
*rho12 = s[0] * rhoa11 * rhoa11;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho31 = s[2] * rhoa32 * rhoa32;
*rho32 = s[2] * rhoa31 * rhoa31;
break;
case C11:
*rho01 = rhoa01;
*rho02 = rhoa02;
*rho11 = rhoa11;
*rho12 = rhoa12;
*rho21 = rhoa21;
*rho22 = rhoa22;
*rho31 = rhoa31;
*rho32 = rhoa32;
break;
case L12:
*rho01 = 8 * rhoa01 + 4 * rhoa02;
*rho02 = 12 * rhoa01;
if (this->ialloy == 1) {
*rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]);
denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]);
if (denom > 0.)
*rho21 = *rho21 / denom * *rho01;
} else
*rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
break;
case B2:
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
break;
case CH4:
*rho01 = 4.0 * rhoa02; //in assumption that 'a' represent carbon
*rho02 = rhoa01; //in assumption that 'b' represent hydrogen
get_shpfcn(DIM, 0, 0, s); //H
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
get_shpfcn(CH4, 0, 0, s); //C
*rho11 = s[0] * rhoa12 * rhoa12;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho31 = s[2] * rhoa32 * rhoa32;
break;
case LIN:
*rho01 = rhoa02*Zij;
*rho02 = rhoa01*Zij;
get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
*rho11 = s[0] * rhoa12 * rhoa12;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho31 = s[2] * rhoa32 * rhoa32;
break;
case ZIG:
*rho01 = rhoa02*Zij;
*rho02 = rhoa01*Zij;
get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
*rho11 = s[0] * rhoa12 * rhoa12;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho31 = s[2] * rhoa32 * rhoa32;
break;
case TRI:
*rho01 = rhoa02;
*rho02 = rhoa01*Zij;
get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
s[0] = 1.0;
s[1] = 2.0/3.0;
s[2] = 1.0 - 0.6*s[0];
*rho11 = s[0] * rhoa12 * rhoa12;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho31 = s[2] * rhoa32 * rhoa32;
break;
// default:
// call error('Lattice not defined in get_densref.')
}
if (this->nn2_meam[a][b] == 1) {
Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], arat, scrn);
Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b],
this->stheta_meam[a][b], arat, scrn);
a1 = arat * r / this->re_meam[a][a] - 1.0;
a2 = arat * r / this->re_meam[b][b] - 1.0;
@ -712,7 +836,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
*rho01 = *rho01 + Zij2nn * scrn * rhoa01nn;
// Assume Zij2nn and arat don't depend on order, but scrn might
Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], arat, scrn);
Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a],
this->stheta_meam[a][b], arat, scrn);
*rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
}
}

37
src/USER-MEAMC/meam_setup_global.cpp
View File

@ -18,7 +18,7 @@ static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
}
void
MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* /*atwt*/, double* alpha,
MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*/, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar)
@ -32,7 +32,6 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub
for (i = 0; i < nelt; i++) {
this->lattce_meam[i][i] = lat[i];
this->Z_meam[i] = z[i];
this->ielt_meam[i] = ielement[i];
this->alpha_meam[i][i] = alpha[i];
this->beta0_meam[i] = b0[i];
@ -49,17 +48,26 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub
this->rho0_meam[i] = rozero[i];
this->ibar_meam[i] = ibar[i];
if (this->lattce_meam[i][i] == FCC)
this->re_meam[i][i] = tmplat[i] / sqrt(2.0);
else if (this->lattce_meam[i][i] == BCC)
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
else if (this->lattce_meam[i][i] == HCP)
this->re_meam[i][i] = tmplat[i];
else if (this->lattce_meam[i][i] == DIM)
this->re_meam[i][i] = tmplat[i];
else if (this->lattce_meam[i][i] == DIA)
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
else {
switch(this->lattce_meam[i][i]) {
case FCC:
this->re_meam[i][i] = tmplat[i] / sqrt(2.0);
break;
case BCC:
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
break;
case HCP:
case DIM:
case CH4:
case LIN:
case ZIG:
case TRI:
this->re_meam[i][i] = tmplat[i];
break;
case DIA:
case DIA3:
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
break;
//default:
// error
}
}
@ -82,4 +90,7 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub
this->emb_lin_neg = 0;
this->bkgd_dyn = 0;
this->erose_form = 0;
// for trimer, zigzag, line refernece structure, sungkwang
setall2d(this->stheta_meam, 1.0); // stheta = sin(theta/2*pi/180) where theta is 180, so 1.0
setall2d(this->ctheta_meam, 0.0); // stheta = cos(theta/2*pi/180) where theta is 180, so 0
}

17
src/USER-MEAMC/meam_setup_param.cpp
View File

@ -1,6 +1,7 @@
#include "meam.h"
#include <algorithm>
using namespace LAMMPS_NS;
using namespace MathConst;
//
// do a sanity check on index parameters
@ -46,6 +47,8 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr
// 18 = zbl_meam
// 19 = emb_lin_neg
// 20 = bkgd_dyn
// 21 = theta
void
MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag)
@ -203,6 +206,20 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3
this->bkgd_dyn = (int)value;
break;
// 21 = theta
// see alloyparams(void) in meam_setup_done.cpp
case 21:
// const double PI = 3.141592653589793238463;
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
i1 = std::min(index[0], index[1]);
i2 = std::max(index[0], index[1]);
// we don't use theta, instead stheta and ctheta
this->stheta_meam[i1][i2] = sin(value/2*MY_PI/180.0);
this->ctheta_meam[i1][i2] = cos(value/2*MY_PI/180.0);
break;
default:
*errorflag = 1;
}

65
src/USER-MEAMC/pair_meamc.cpp
View File

@ -33,13 +33,13 @@ using namespace LAMMPS_NS;
#define MAXLINE 1024
static const int nkeywords = 21;
static const int nkeywords = 22;
static const char *keywords[] = {
"Ec","alpha","rho0","delta","lattce",
"attrac","repuls","nn2","Cmin","Cmax","rc","delr",
"augt1","gsmooth_factor","re","ialloy",
"mixture_ref_t","erose_form","zbl",
"emb_lin_neg","bkgd_dyn"};
"emb_lin_neg","bkgd_dyn", "theta"};
/* ---------------------------------------------------------------------- */
@ -56,6 +56,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
elements = NULL;
mass = NULL;
meam_inst = new MEAM(memory);
scale = NULL;
// set comm size needed by this Pair
@ -79,6 +80,7 @@ PairMEAMC::~PairMEAMC()
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(scale);
delete [] map;
}
}
@ -143,7 +145,7 @@ void PairMEAMC::compute(int eflag, int vflag)
comm->reverse_comm_pair(this);
meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom,
&eng_vdwl,eatom,ntype,type,map,errorflag);
&eng_vdwl,eatom,ntype,type,map,scale,errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
@ -164,7 +166,7 @@ void PairMEAMC::compute(int eflag, int vflag)
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom,
vflag_atom,&eng_vdwl,eatom,ntype,type,map,x,
vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x,
numneigh_half[i],firstneigh_half[i],
numneigh_full[i],firstneigh_full[i],
offset,f,vptr);
@ -183,6 +185,7 @@ void PairMEAMC::allocate()
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(scale,n+1,n+1,"pair:scale");
map = new int[n+1];
}
@ -267,13 +270,16 @@ void PairMEAMC::coeff(int narg, char **arg)
// set mass for i,i in atom class
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
count++;
}
scale[i][j] = 1.0;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
@ -312,8 +318,10 @@ void PairMEAMC::init_list(int id, NeighList *ptr)
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMEAMC::init_one(int /*i*/, int /*j*/)
double PairMEAMC::init_one(int i, int j)
{
if (setflag[i][j] == 0) scale[i][j] = 1.0;
scale[j][i] = scale[i][j];
return cutmax;
}
@ -430,13 +438,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
found[i] = true;
// map lat string to an integer
if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
else if (strcmp(words[1],"dia") == 0) lat[i] = DIA;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
if (!MEAM::str_to_lat(words[1], true, lat[i]))
error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
// store parameters
@ -458,8 +462,12 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
rozero[i] = atof(words[17]);
ibar[i] = atoi(words[18]);
if (!iszero(t0[i]-1.0))
error->all(FLERR,"Unsupported parameter in MEAM potential file");
if (!isone(t0[i]))
error->all(FLERR,"Unsupported parameter in MEAM potential file: t0!=1");
// z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
if (z[i] != MEAM::get_Zij(lat[i]))
error->all(FLERR,"Mismatched parameter in MEAM potential file: z!=lat");
nset++;
}
@ -471,7 +479,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
// pass element parameters to MEAM package
meam_inst->meam_setup_global(nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
meam_inst->meam_setup_global(nelements,lat,ielement,atwt,alpha,b0,b1,b2,b3,
alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
// set element masses
@ -523,6 +531,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
int which;
double value;
lattice_t latt;
int nindex,index[3];
int maxparams = 6;
char **params = new char*[maxparams];
@ -571,16 +580,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
// map lattce_meam value to an integer
if (which == 4) {
if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
else if (strcmp(params[nparams-1],"dia") == 0) value = DIA;
else if (strcmp(params[nparams-1],"b1") == 0) value = B1;
else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
else if (strcmp(params[nparams-1],"b2") == 0) value = B2;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
if (!MEAM::str_to_lat(params[nparams-1], false, latt))
error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
value = latt;
}
else value = atof(params[nparams-1]);
@ -783,3 +785,12 @@ void PairMEAMC::neigh_strip(int inum, int *ilist,
for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
}
}
/* ---------------------------------------------------------------------- */
void *PairMEAMC::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"scale") == 0) return (void *) scale;
return NULL;
}

2
src/USER-MEAMC/pair_meamc.h
View File

@ -35,6 +35,7 @@ class PairMEAMC : public Pair {
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
virtual void *extract(const char *, int &);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
@ -50,6 +51,7 @@ class PairMEAMC : public Pair {
double *mass; // mass of each element
int *map; // mapping from atom types (1-indexed) to elements (1-indexed)
double **scale; // scaling factor for adapt
void allocate();
void read_files(char *, char *);

93
src/USER-MISC/fix_bond_react.cpp
View File

@ -36,6 +36,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
#include "group.h"
#include "citeme.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
@ -297,6 +298,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges");
memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms");
memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms");
for (int j = 0; j < nreacts; j++)
for (int i = 0; i < max_natoms; i++) {
@ -304,6 +306,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (update_edges_flag[j] == 1) custom_edges[i][j] = 1;
else custom_edges[i][j] = 0;
delete_atoms[i][j] = 0;
for (int k = 0; k < 6; k++) {
chiral_atoms[i][k][j] = 0;
}
}
// read all map files afterward
@ -430,6 +435,7 @@ FixBondReact::~FixBondReact()
memory->destroy(landlocked_atoms);
memory->destroy(custom_edges);
memory->destroy(delete_atoms);
memory->destroy(chiral_atoms);
memory->destroy(nevery);
memory->destroy(cutsq);
@ -1675,6 +1681,7 @@ int FixBondReact::check_constraints()
delx1 = x[atom1][0] - x[atom2][0];
dely1 = x[atom1][1] - x[atom2][1];
delz1 = x[atom1][2] - x[atom2][2];
domain->minimum_image(delx1,dely1,delz1); // ghost location fix
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
@ -1682,6 +1689,7 @@ int FixBondReact::check_constraints()
delx2 = x[atom3][0] - x[atom2][0];
dely2 = x[atom3][1] - x[atom2][1];
delz2 = x[atom3][2] - x[atom2][2];
domain->minimum_image(delx2,dely2,delz2); // ghost location fix
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
@ -1699,6 +1707,30 @@ int FixBondReact::check_constraints()
}
}
}
// let's also check chirality within 'check_constraint'
for (int i = 0; i < onemol->natoms; i++) {
if (chiral_atoms[i][0][rxnID] == 1) {
double my4coords[12];
// already ensured, by transitive property, that chiral simulation atom has four neighs
for (int j = 0; j < 4; j++) {
atom1 = atom->map(glove[i][1]);
// loop over known types involved in chiral center
for (int jj = 0; jj < 4; jj++) {
if (atom->type[atom->map(xspecial[atom1][j])] == chiral_atoms[i][jj+2][rxnID]) {
atom2 = atom->map(xspecial[atom1][j]);
atom2 = domain->closest_image(atom1,atom2);
for (int k = 0; k < 3; k++) {
my4coords[3*jj+k] = x[atom2][k];
}
break;
}
}
}
if (get_chirality(my4coords) != chiral_atoms[i][1][rxnID]) return 0;
}
}
return 1;
}
@ -1739,6 +1771,33 @@ double FixBondReact::get_temperature()
return t;
}
/* ----------------------------------------------------------------------
return handedness (1 or -1) of a chiral center, given ordered set of coordinates
------------------------------------------------------------------------- */
int FixBondReact::get_chirality(double four_coords[12])
{
// define oriented plane with first three coordinates
double vec1[3],vec2[3],vec3[3],vec4[3],mean3[3],dot;
for (int i = 0; i < 3; i++) {
vec1[i] = four_coords[i]-four_coords[i+3];
vec2[i] = four_coords[i+3]-four_coords[i+6];
}
MathExtra::cross3(vec1,vec2,vec3);
for (int i = 0; i < 3; i++) {
mean3[i] = (four_coords[i] + four_coords[i+3] +
four_coords[i+6])/3;
vec4[i] = four_coords[i+9] - mean3[i];
}
dot = MathExtra::dot3(vec3,vec4);
dot = dot/fabs(dot);
return (int) dot;
}
/* ----------------------------------------------------------------------
Get xspecials for current molecule templates
------------------------------------------------------------------------- */
@ -2867,6 +2926,7 @@ void FixBondReact::read(int myrxn)
}
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral);
else if (strstr(line,"constraints")) {
sscanf(line,"%d",&nconstr);
memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");
@ -2898,6 +2958,8 @@ void FixBondReact::read(int myrxn)
CustomEdges(line, myrxn);
} else if (strcmp(keyword,"DeleteIDs") == 0) {
DeleteAtoms(line, myrxn);
} else if (strcmp(keyword,"ChiralIDs") == 0) {
ChiralCenters(line, myrxn);
} else if (strcmp(keyword,"Constraints") == 0) {
Constraints(line, myrxn);
} else error->one(FLERR,"Bond/react: Unknown section in map file");
@ -2975,6 +3037,37 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn)
}
}
void FixBondReact::ChiralCenters(char *line, int myrxn)
{
int tmp;
for (int i = 0; i < nchiral; i++) {
readline(line);
sscanf(line,"%d",&tmp);
chiral_atoms[tmp-1][0][myrxn] = 1;
if (onemol->xflag == 0)
error->one(FLERR,"Bond/react: Molecule template 'Coords' section required for chiralIDs keyword");
if ((int) onemol_nxspecial[tmp-1][0] != 4)
error->one(FLERR,"Bond/react: Chiral atoms must have exactly four first neighbors");
for (int j = 0; j < 4; j++) {
for (int k = j+1; k < 4; k++) {
if (onemol->type[onemol_xspecial[tmp-1][j]-1] ==
onemol->type[onemol_xspecial[tmp-1][k]-1])
error->one(FLERR,"Bond/react: First neighbors of chiral atoms must be of mutually different types");
}
}
// record order of atom types, and coords
double my4coords[12];
for (int j = 0; j < 4; j++) {
chiral_atoms[tmp-1][j+2][myrxn] = onemol->type[onemol_xspecial[tmp-1][j]-1];
for (int k = 0; k < 3; k++) {
my4coords[3*j+k] = onemol->x[onemol_xspecial[tmp-1][j]-1][k];
}
}
// get orientation
chiral_atoms[tmp-1][1][myrxn] = get_chirality(my4coords);
}
}
void FixBondReact::Constraints(char *line, int myrxn)
{
double tmp[MAXCONARGS];

17
src/USER-MISC/fix_bond_react.h
View File

@ -110,7 +110,7 @@ class FixBondReact : public Fix {
int *ibonding,*jbonding;
int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge,nequivalent,ncustom,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file
int nedge,nequivalent,ncustom,ndelete,nchiral,nconstr; // # edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -126,6 +126,7 @@ class FixBondReact : public Fix {
int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
int **custom_edges; // atoms in molecule templates with incorrect valences
int **delete_atoms; // atoms in pre-reacted templates to delete
int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types
int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@ -149,6 +150,7 @@ class FixBondReact : public Fix {
void Equivalences(char *,int);
void CustomEdges(char *,int);
void DeleteAtoms(char *,int);
void ChiralCenters(char *,int);
void Constraints(char *,int);
void make_a_guess ();
@ -159,6 +161,7 @@ class FixBondReact : public Fix {
void ring_check();
int check_constraints();
double get_temperature();
int get_chirality(double[12]); // get handedness given an ordered set of coordinates
void open(char *);
void readline(char *);
@ -249,6 +252,18 @@ E: Bond/react: A deleted atom cannot remain bonded to an atom that is not delete
Self-explanatory.
E: Bond/react: First neighbors of chiral atoms must be of mutually different types
Self-explanatory.
E: Bond/react: Chiral atoms must have exactly four first neighbors
Self-explanatory.
E: Bond/react: Molecule template 'Coords' section required for chiralIDs keyword
The coordinates of atoms in the pre-reacted template are used to determine chirality.
E: Bond/react special bond generation overflow
The number of special bonds per-atom created by a reaction exceeds the

7
src/USER-PHONON/dynamical_matrix.cpp
View File

@ -19,6 +19,7 @@
#include "improper.h"
#include "kspace.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "timer.h"
@ -385,6 +386,7 @@ void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
void DynamicalMatrix::update_force()
{
force_clear();
int n_post_force = modify->n_post_force;
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
@ -405,6 +407,11 @@ void DynamicalMatrix::update_force()
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications,
if (n_post_force) modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
++ update->nsteps;
}

2
src/math_special.cpp
View File

@ -538,6 +538,8 @@ double MathSpecial::exp2_x86(double x)
double MathSpecial::fm_exp(double x)
{
#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
if (x < -1022.0/FM_DOUBLE_LOG2OFE) return 0;
if (x > 1023.0/FM_DOUBLE_LOG2OFE) return INFINITY;
return exp2_x86(FM_DOUBLE_LOG2OFE * x);
#else
return ::exp(x);

6
src/min.h
View File

@ -39,9 +39,9 @@ class Min : protected Pointers {
virtual bigint memory_usage() {return 0;}
void modify_params(int, char **);
virtual int modify_param(int, char **) {return 0;}
double fnorm_sqr();
double fnorm_inf();
double fnorm_max();
virtual double fnorm_sqr();
virtual double fnorm_inf();
virtual double fnorm_max();
enum{TWO,MAX,INF};