diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index 71de83afa9..a1795a6776 100644 --- a/doc/src/Errors_messages.rst +++ b/doc/src/Errors_messages.rst @@ -518,6 +518,16 @@ Doc page with :doc:`WARNING messages ` *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 diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index ab63aa232c..b2866eb9c7 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -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: diff --git a/doc/src/pair_meamc.rst b/doc/src/pair_meamc.rst index daa10906c7..ac858ce24b 100644 --- a/doc/src/pair_meamc.rst +++ b/doc/src/pair_meamc.rst @@ -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) `). 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) ` 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 diff --git a/doc/txt/pair_meamc.txt b/doc/txt/pair_meamc.txt index 7c42e9d2f2..2831600e08 100644 --- a/doc/txt/pair_meamc.txt +++ b/doc/txt/pair_meamc.txt @@ -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 diff --git a/lib/gpu/Makefile.lammps.standard b/lib/gpu/Makefile.lammps.standard index 1ceb93b342..9526e8e373 100644 --- a/lib/gpu/Makefile.lammps.standard +++ b/lib/gpu/Makefile.lammps.standard @@ -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 diff --git a/lib/gpu/Makefile.linux b/lib/gpu/Makefile.linux index 7001c6d8b9..a26fbe114c 100644 --- a/lib/gpu/Makefile.linux +++ b/lib/gpu/Makefile.linux @@ -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 diff --git a/lib/gpu/Makefile.linux.double b/lib/gpu/Makefile.linux.double index e65647f160..05f083697d 100644 --- a/lib/gpu/Makefile.linux.double +++ b/lib/gpu/Makefile.linux.double @@ -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 diff --git a/lib/gpu/Makefile.linux.mixed b/lib/gpu/Makefile.linux.mixed index a036b84ee3..ca414f1fc1 100644 --- a/lib/gpu/Makefile.linux.mixed +++ b/lib/gpu/Makefile.linux.mixed @@ -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 diff --git a/lib/gpu/Makefile.linux.single b/lib/gpu/Makefile.linux.single index 808647cea7..1b349faac2 100644 --- a/lib/gpu/Makefile.linux.single +++ b/lib/gpu/Makefile.linux.single @@ -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 diff --git a/lib/gpu/Makefile.linux_multi b/lib/gpu/Makefile.linux_multi index 2f75ca0e2b..ba50170f39 100644 --- a/lib/gpu/Makefile.linux_multi +++ b/lib/gpu/Makefile.linux_multi @@ -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 diff --git a/lib/gpu/Makefile.mac b/lib/gpu/Makefile.mac index f9f8d5179a..b571e8d85a 100644 --- a/lib/gpu/Makefile.mac +++ b/lib/gpu/Makefile.mac @@ -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 diff --git a/lib/gpu/Makefile.serial b/lib/gpu/Makefile.serial index 99153fc471..b0cfb3c86b 100644 --- a/lib/gpu/Makefile.serial +++ b/lib/gpu/Makefile.serial @@ -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 diff --git a/lib/gpu/Makefile.shannon b/lib/gpu/Makefile.shannon index 22c2dc89d7..b14bd01ec2 100644 --- a/lib/gpu/Makefile.shannon +++ b/lib/gpu/Makefile.shannon @@ -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 diff --git a/lib/gpu/Makefile.xk7 b/lib/gpu/Makefile.xk7 index 0b9f029399..f3050575e1 100644 --- a/lib/gpu/Makefile.xk7 +++ b/lib/gpu/Makefile.xk7 @@ -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 diff --git a/lib/gpu/README b/lib/gpu/README index 6cc386270d..2d98749a40 100644 --- a/lib/gpu/README +++ b/lib/gpu/README @@ -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: diff --git a/src/.gitignore b/src/.gitignore index 5848874d94..1b7dc7e6b2 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/USER-INTEL/npair_intel.cpp b/src/USER-INTEL/npair_intel.cpp index a82d3f29e5..4256e03b3c 100644 --- a/src/USER-INTEL/npair_intel.cpp +++ b/src/USER-INTEL/npair_intel.cpp @@ -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 diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 4e5b298ffb..719e1af6f0 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -2,13 +2,15 @@ #define LMP_MEAM_H #include +#include +#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) { diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 867106df88..fe9e74ca7c 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -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; } diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index 155941422c..39e3ff9467 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -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; diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 09aad90111..2b6832e155 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -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]; diff --git a/src/USER-MEAMC/meam_funcs.cpp b/src/USER-MEAMC/meam_funcs.cpp index 312caff686..da984cfeb7 100644 --- a/src/USER-MEAMC/meam_funcs.cpp +++ b/src/USER-MEAMC/meam_funcs.cpp @@ -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; +} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp index 74e8af1cde..0a0b95e14a 100644 --- a/src/USER-MEAMC/meam_impl.cpp +++ b/src/USER-MEAMC/meam_impl.cpp @@ -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; diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index 37bfce5873..99b330a83a 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -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 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; } } diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index b66aafa06c..ac5718914c 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -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 } diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index 585b5b5712..70d8b63598 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,6 +1,7 @@ #include "meam.h" #include 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; } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index a9034a1af3..9733ee2dad 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -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; +} diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index 31dd8ba022..102985f0ca 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -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 *); diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 33aab5b4ad..2995a596a5 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -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]; diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index f59e051ed1..a5bf2bc587 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -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 diff --git a/src/USER-PHONON/dynamical_matrix.cpp b/src/USER-PHONON/dynamical_matrix.cpp index 1495219124..6bb843c16e 100644 --- a/src/USER-PHONON/dynamical_matrix.cpp +++ b/src/USER-PHONON/dynamical_matrix.cpp @@ -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; } diff --git a/src/math_special.cpp b/src/math_special.cpp index bf11a1ad45..d4abc36f25 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -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); diff --git a/src/min.h b/src/min.h index fe7f0be8c6..ba47f95410 100644 --- a/src/min.h +++ b/src/min.h @@ -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};