ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3620 from akohlmey/collected-small-fixes

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2023-02-02 14:37:39 -05:00
committed by GitHub
parent b8b5e385b6 c9b23b8a03
commit 11871ceb2e
43 changed files with 933 additions and 982 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
cmake/Modules/Packages/ML-HDNNP.cmake
View File

@ -6,8 +6,8 @@ else()
endif()
option(DOWNLOAD_N2P2 "Download n2p2 library instead of using an already installed one)" ${DOWNLOAD_N2P2_DEFAULT})
if(DOWNLOAD_N2P2)
set(N2P2_URL "https://github.com/CompPhysVienna/n2p2/archive/v2.1.4.tar.gz" CACHE STRING "URL for n2p2 tarball")
set(N2P2_MD5 "9595b066636cd6b90b0fef93398297a5" CACHE STRING "MD5 checksum of N2P2 tarball")
set(N2P2_URL "https://github.com/CompPhysVienna/n2p2/archive/v2.2.0.tar.gz" CACHE STRING "URL for n2p2 tarball")
set(N2P2_MD5 "a2d9ab7f676b3a74a324fc1eda0a911d" CACHE STRING "MD5 checksum of N2P2 tarball")
mark_as_advanced(N2P2_URL)
mark_as_advanced(N2P2_MD5)
GetFallbackURL(N2P2_URL N2P2_FALLBACK)

1
cmake/Modules/Packages/ML-QUIP.cmake
View File

@ -16,6 +16,7 @@ if(DOWNLOAD_QUIP)
set(temp "${temp}DEFINES += -DGETARG_F2003 -DFORTRAN_UNDERSCORE\n")
set(temp "${temp}F95FLAGS += -fpp -free -fPIC\n")
set(temp "${temp}F77FLAGS += -fpp -fixed -fPIC\n")
set(temp "${temp}F95_PRE_FILENAME_FLAG = -Tf\n")
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
set(temp "${temp}FPP=${CMAKE_Fortran_COMPILER} -E -x f95-cpp-input\nOPTIM=${CMAKE_Fortran_FLAGS_${BTYPE}}\n")
set(temp "${temp}DEFINES += -DGETARG_F2003 -DGETENV_F2003 -DGFORTRAN -DFORTRAN_UNDERSCORE\n")

6
cmake/Modules/Tools.cmake
View File

@ -50,12 +50,16 @@ if(BUILD_LAMMPS_SHELL)
add_executable(lammps-shell ${LAMMPS_TOOLS_DIR}/lammps-shell/lammps-shell.cpp ${ICON_RC_FILE})
target_include_directories(lammps-shell PRIVATE ${LAMMPS_TOOLS_DIR}/lammps-shell)
target_link_libraries(lammps-shell PRIVATE lammps PkgConfig::READLINE)
# workaround for broken readline pkg-config file on FreeBSD
if(CMAKE_SYSTEM_NAME STREQUAL "FreeBSD")
target_include_directories(lammps-shell PRIVATE /usr/local/include)
endif()
target_link_libraries(lammps-shell PRIVATE lammps PkgConfig::READLINE)
if(CMAKE_SYSTEM_NAME STREQUAL "LinuxMUSL")
pkg_check_modules(TERMCAP IMPORTED_TARGET REQUIRED termcap)
target_link_libraries(lammps-shell PRIVATE lammps PkgConfig::TERMCAP)
endif()
install(TARGETS lammps-shell EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR})
install(DIRECTORY ${LAMMPS_TOOLS_DIR}/lammps-shell/icons DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/)
install(FILES ${LAMMPS_TOOLS_DIR}/lammps-shell/lammps-shell.desktop DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/applications/)

5
doc/src/pair_gauss.rst
View File

@ -127,9 +127,8 @@ before mixing, and converted back after mixing
This way, if either particle is repulsive (if Ai<0 or Aj<0),
then the default interaction between both particles will be repulsive.
The *gauss* style does not support the :doc:`pair_modify <pair_modify>`
shift option. There is no effect due to the Gaussian well beyond the
cutoff; hence reasonable cutoffs need to be specified.
For the *gauss* style there is no effect due to the Gaussian well
beyond the cutoff; hence reasonable cutoffs need to be specified.
The *gauss/cut* style supports the :doc:`pair_modify <pair_modify>` shift
option for the energy of the Gauss-potential portion of the pair

2
doc/utils/sphinx-config/false_positives.txt
View File

@ -363,6 +363,7 @@ Broglie
brownian
brownw
Broyden
Bruskin
Brusselle
Bryantsev
Btarget
@ -2241,6 +2242,7 @@ msd
msi
MSI
msm
msmeam
msmflag
msse
msst

17
examples/PACKAGES/pimd/para-h2/in.scp
View File

@ -1,5 +1,5 @@
variable ibead uloop 32 pad
variable out_freq string 1
variable out_freq string 100
variable job_name string H2
units real
@ -13,17 +13,18 @@ read_data H2.data
pair_coeff 1 * pair.table PAIR_H2
timestep 0.001
thermo_style custom step temp pe etotal pzz
thermo ${out_freq}
timestep 0.001
velocity all create 1.0 1985 rot yes dist gaussian
fix 1 all pimd method nmpimd fmass 1.0 temp 25.0 nhc 4
dump dcd all dcd ${out_freq} dcd/${job_name}_${ibead}.dcd
thermo_style custom step temp pe etotal pzz f_1[1] f_1[2] f_1[3]
thermo_modify colname f_1[1] espring colname f_1[2] T_ring colname f_1[3] virial
thermo ${out_freq}
restart ${out_freq} restart/${job_name}_${ibead}.restart1 restart/${job_name}_${ibead}.restart2
#dump dcd all dcd ${out_freq} ${job_name}_${ibead}.dcd
run 100000
#restart ${out_freq} ${job_name}_${ibead}.restart1 ${job_name}_${ibead}.restart2
run 10000

2
examples/PACKAGES/pimd/para-h2/run.sh
View File

@ -1 +1 @@
mpirun -np 64 lmp_mpi -partition 8x8 -in in.scp -log logfile/log.lammps -screen screen/screen
mpirun -np 64 lmp_mpi -partition 8x8 -in in.scp -log log.lammps -screen none

11
examples/PACKAGES/pimd/prot-hairpin/in.scp
View File

@ -14,18 +14,19 @@ pair_modify mix arithmetic
kspace_style pppm 1e-4
read_data system.data
#read_restart restart/system_${ibead}.rest1
#read_restart system_${ibead}.rest1
special_bonds charmm
fix 1 all pimd method nmpimd fmass 1.0 temp 300.0 nhc 4 sp 2.0
thermo 10
thermo_style custom step temp pe etotal
thermo_style custom step temp pe etotal f_1[1] f_1[2] f_1[3]
thermo_modify colname f_1[1] espring colname f_1[2] T_ring colname f_1[3] virial
timestep 0.08
restart 100 restart/system_${ibead}.rest1 restart/system_${ibead}.rest2
# restart 100 system_${ibead}.rest1 system_${ibead}.rest2
group prot id <= 256
dump 1 prot dcd 100 dcd/prot_${ibead}.dcd
# dump 1 prot dcd 100 prot_${ibead}.dcd
run 2000000
run 200

2
examples/PACKAGES/pimd/prot-hairpin/run.sh
View File

@ -1 +1 @@
mpirun -np 8 lmp_mpi -partition 8x1 -in in.scp -log logfile/log.lammps -screen screen/screen
mpirun -np 8 lmp_mpi -partition 8x1 -in in.scp -log log.lammps -screen none

4
lib/gpu/lal_atom.cpp
View File

@ -202,14 +202,14 @@ bool AtomT::add_fields(const bool charge, const bool rot,
if (extra_fields > 0 && _extra_fields==0) {
_extra_fields=extra_fields;
_other=true;
if (_host_view==false) {
if (!_host_view) {
success=success && (extra.alloc(_max_atoms*_extra_fields,*dev,UCL_WRITE_ONLY,
UCL_READ_ONLY)==UCL_SUCCESS);
gpu_bytes+=extra.device.row_bytes();
}
}
if (bonds && _bonds==false) {
if (bonds && !_bonds) {
_bonds=true;
if (_bonds && _gpu_nbor>0) {
success=success && (dev_tag.alloc(_max_atoms,*dev,

3
lib/hdnnp/Install.py
View File

@ -17,7 +17,7 @@ parser = ArgumentParser(prog='Install.py',
# settings
version = "2.1.4"
version = "2.2.0"
# help message
@ -38,6 +38,7 @@ make lib-hdnnp args="-p $HOME/n2p2" # use existing n2p2 installation in $HOME/n2
# known checksums for different n2p2 versions. used to validate the download.
checksums = { \
'2.2.0' : 'a2d9ab7f676b3a74a324fc1eda0a911d', \
'2.1.4' : '9595b066636cd6b90b0fef93398297a5', \
}

2
lib/kokkos/core/src/Kokkos_Macros.hpp
View File

@ -653,7 +653,9 @@ static constexpr bool kokkos_omp_on_host() { return false; }
#if (defined(KOKKOS_COMPILER_GNU) || defined(KOKKOS_COMPILER_CLANG) || \
defined(KOKKOS_COMPILER_INTEL) || defined(KOKKOS_COMPILER_PGI)) && \
!defined(_WIN32)
#if (defined(__linux__) && defined(__GLIBC_MINOR__))
#define KOKKOS_IMPL_ENABLE_STACKTRACE
#endif
#define KOKKOS_IMPL_ENABLE_CXXABI
#endif

2
src/DIELECTRIC/pppm_dielectric.cpp
View File

@ -356,7 +356,7 @@ void PPPMDielectric::make_rho()
// (mx,my,mz) = global coords of moving stencil pt
double *q = atom->q_scaled;
if (use_qscaled == false) q = atom->q;
if (!use_qscaled) q = atom->q;
double **x = atom->x;
int nlocal = atom->nlocal;

2
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

@ -636,7 +636,7 @@ void PPPMDispDielectric::make_rho_c()
// (mx,my,mz) = global coords of moving stencil pt
double *q = atom->q_scaled;
if (use_qscaled == false) q = atom->q;
if (!use_qscaled) q = atom->q;
double **x = atom->x;
int nlocal = atom->nlocal;

2
src/ELECTRODE/fix_electrode_conp.cpp
View File

@ -917,7 +917,7 @@ std::vector<double> FixElectrodeConp::gather_ngroup(std::vector<double> x_local)
std::vector<double> FixElectrodeConp::constraint_correction(std::vector<double> x)
{
return constraint_projection(x);
return constraint_projection(std::move(x));
}
/* ----------------------------------------------------------------------

7
src/KOKKOS/atom_vec_hybrid_kokkos.cpp
View File

@ -37,13 +37,6 @@ AtomVecKokkos(lmp), AtomVecHybrid(lmp)
/* ---------------------------------------------------------------------- */
AtomVecHybridKokkos::~AtomVecHybridKokkos()
{
}
/* ---------------------------------------------------------------------- */
void AtomVecHybridKokkos::grow(int n)
{
for (int k = 0; k < nstyles; k++) styles[k]->grow(n);

1
src/KOKKOS/atom_vec_hybrid_kokkos.h
View File

@ -32,7 +32,6 @@ namespace LAMMPS_NS {
class AtomVecHybridKokkos : public AtomVecKokkos, public AtomVecHybrid {
public:
AtomVecHybridKokkos(class LAMMPS *);
~AtomVecHybridKokkos() override;
void grow(int) override;

2
src/KOKKOS/kokkos.cpp
View File

@ -342,7 +342,7 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
/* ---------------------------------------------------------------------- */
void KokkosLMP::initialize(Kokkos::InitializationSettings args, Error *error)
void KokkosLMP::initialize(const Kokkos::InitializationSettings& args, Error *error)
{
if (!Kokkos::is_initialized()) {
if (is_finalized)

2
src/KOKKOS/kokkos.h
View File

@ -56,7 +56,7 @@ class KokkosLMP : protected Pointers {
KokkosLMP(class LAMMPS *, int, char **);
static void initialize(Kokkos::InitializationSettings, Error *);
static void initialize(const Kokkos::InitializationSettings&, Error *);
static void finalize();
void accelerator(int, char **);
int neigh_count(int);

2
src/KOKKOS/meam_force_kokkos.h
View File

@ -446,7 +446,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
nv2 = 0;
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
arg = delij[n] * delij[p] * this->v2D[nv2];
arg = delij[n] * delij[p] * v2D[nv2];
drho3mdrm1[m] = drho3mdrm1[m] + d_arho3m(i, vind3D[m][n][p]) * arg;
drho3mdrm2[m] = drho3mdrm2[m] + d_arho3m(j, vind3D[m][n][p]) * arg;
nv2 = nv2 + 1;

7
src/KOKKOS/mliap_descriptor_so3_kokkos.cpp
View File

@ -46,13 +46,6 @@ MLIAPDescriptorSO3Kokkos<DeviceType>::MLIAPDescriptorSO3Kokkos(LAMMPS *lmp, char
/* ---------------------------------------------------------------------- */
template <class DeviceType>
MLIAPDescriptorSO3Kokkos<DeviceType>::~MLIAPDescriptorSO3Kokkos()
{
}
/* ---------------------------------------------------------------------- */
template <class DeviceType>
void MLIAPDescriptorSO3Kokkos<DeviceType>::compute_descriptors(class MLIAPData *data_)
{

1
src/KOKKOS/mliap_descriptor_so3_kokkos.h
View File

@ -29,7 +29,6 @@ class MLIAPDescriptorSO3Kokkos :
public MLIAPDescriptorKokkos<DeviceType> {
public:
MLIAPDescriptorSO3Kokkos(LAMMPS *, char *);
~MLIAPDescriptorSO3Kokkos() override;
void compute_descriptors(class MLIAPData *) override;
void compute_forces(class MLIAPData *) override;

15
src/MC/fix_gcmc.cpp
View File

@ -138,10 +138,17 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
region_zlo = region->extent_zlo;
region_zhi = region->extent_zhi;
if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] ||
region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] ||
region_zlo < domain->boxlo[2] || region_zhi > domain->boxhi[2])
error->all(FLERR,"Fix gcmc region extends outside simulation box");
if (triclinic) {
if ((region_xlo < domain->boxlo_bound[0]) || (region_xhi > domain->boxhi_bound[0]) ||
(region_ylo < domain->boxlo_bound[1]) || (region_yhi > domain->boxhi_bound[1]) ||
(region_zlo < domain->boxlo_bound[2]) || (region_zhi > domain->boxhi_bound[2]))
error->all(FLERR,"Fix gcmc region {} extends outside simulation box", region->id);
} else {
if ((region_xlo < domain->boxlo[0]) || (region_xhi > domain->boxhi[0]) ||
(region_ylo < domain->boxlo[1]) || (region_yhi > domain->boxhi[1]) ||
(region_zlo < domain->boxlo[2]) || (region_zhi > domain->boxhi[2]))
error->all(FLERR,"Fix gcmc region {} extends outside simulation box", region->id);
}
// estimate region volume using MC trials

4
src/MC/fix_sgcmc.cpp
View File

@ -406,7 +406,7 @@ void FixSemiGrandCanonicalMC::doMC()
}
}
if (kappa != 0.0 && serialMode == false) {
if (kappa != 0.0 && !serialMode) {
// What follows is the second rejection test for the variance-constrained
// semi-grandcanonical method.
@ -472,7 +472,7 @@ void FixSemiGrandCanonicalMC::doMC()
// For (parallelized) semi-grandcanonical MC we have to determine the current concentrations now.
// For the serial version and variance-constrained MC it has already been done in the loop.
if (kappa == 0.0 && serialMode == false) {
if (kappa == 0.0 && !serialMode) {
const int *type = atom->type;
std::vector<int> localSpeciesCounts(atom->ntypes+1, 0);
for (int i = 0; i < atom->nlocal; i++, ++type) {

15
src/MC/fix_widom.cpp
View File

@ -111,10 +111,17 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
region_zlo = region->extent_zlo;
region_zhi = region->extent_zhi;
if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] ||
region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] ||
region_zlo < domain->boxlo[2] || region_zhi > domain->boxhi[2])
error->all(FLERR,"Fix widom region {} extends outside simulation box", region->id);
if (triclinic) {
if ((region_xlo < domain->boxlo_bound[0]) || (region_xhi > domain->boxhi_bound[0]) ||
(region_ylo < domain->boxlo_bound[1]) || (region_yhi > domain->boxhi_bound[1]) ||
(region_zlo < domain->boxlo_bound[2]) || (region_zhi > domain->boxhi_bound[2]))
error->all(FLERR,"Fix widom region {} extends outside simulation box", region->id);
} else {
if ((region_xlo < domain->boxlo[0]) || (region_xhi > domain->boxhi[0]) ||
(region_ylo < domain->boxlo[1]) || (region_yhi > domain->boxhi[1]) ||
(region_zlo < domain->boxlo[2]) || (region_zhi > domain->boxhi[2]))
error->all(FLERR,"Fix widom region {} extends outside simulation box", region->id);
}
// estimate region volume using MC trials

100
src/MEAM/meam_dens_final.cpp
View File

@ -15,9 +15,9 @@
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, double** scale, int& errorflag)
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,
double **scale, int &errorflag)
{
int i, elti;
int m;
@ -27,7 +27,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
// Complete the calculation of density
if (this->msmeamflag) {
if (msmeamflag) {
for (i = 0; i < nlocal; i++) {
elti = fmap[type[i]];
if (elti >= 0) {
@ -40,16 +40,16 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
rho1[i] = rho1[i] + arho1[i][m] * arho1[i][m]
- arho1m[i][m] * arho1m[i][m];
rho3[i] = rho3[i] - 3.0 / 5.0 * (arho3b[i][m] * arho3b[i][m]
- arho3mb[i][m] * arho3mb[i][m]);
- arho3mb[i][m] * arho3mb[i][m]);
}
for (m = 0; m < 6; m++) {
rho2[i] = rho2[i] + this->v2D[m] * (arho2[i][m] * arho2[i][m]
- arho2m[i][m] * arho2m[i][m]);
rho2[i] = rho2[i] + v2D[m] * (arho2[i][m] * arho2[i][m]
- arho2m[i][m] * arho2m[i][m]);
}
for (m = 0; m < 10; m++) {
rho3[i] = rho3[i] + this->v3D[m] * (arho3[i][m] * arho3[i][m]
- arho3m[i][m] * arho3m[i][m]);
rho3[i] = rho3[i] + v3D[m] * (arho3[i][m] * arho3[i][m]
- arho3m[i][m] * arho3m[i][m]);
}
// all the t weights are already accounted for with msmeam
@ -59,48 +59,48 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
}
Z = get_Zij(this->lattce_meam[elti][elti]);
Z = get_Zij(lattce_meam[elti][elti]);
G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
G = G_gam(gamma[i], ibar_meam[elti], errorflag);
if (errorflag != 0)
return;
get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp);
get_shpfcn(lattce_meam[elti][elti], stheta_meam[elti][elti], ctheta_meam[elti][elti], shp);
if (this->ibar_meam[elti] <= 0) {
if (ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
if (this->mix_ref_t == 1) {
if (mix_ref_t == 1) {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
} else {
gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /
gam = (t1_meam[elti] * shp[0] + t2_meam[elti] * shp[1] + t3_meam[elti] * shp[2]) /
(Z * Z);
}
Gbar = G_gam(gam, this->ibar_meam[elti], errorflag);
Gbar = G_gam(gam, ibar_meam[elti], errorflag);
}
rho[i] = rho0[i] * G;
if (this->mix_ref_t == 1) {
if (this->ibar_meam[elti] <= 0) {
if (mix_ref_t == 1) {
if (ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
Gbar = dG_gam(gam, this->ibar_meam[elti], dGbar);
Gbar = dG_gam(gam, ibar_meam[elti], dGbar);
}
rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
rho_bkgd = rho0_meam[elti] * Z * Gbar;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd = this->rho0_meam[elti] * Z;
if (bkgd_dyn == 1) {
rho_bkgd = rho0_meam[elti] * Z;
} else {
rho_bkgd = this->rho_ref_meam[elti];
rho_bkgd = rho_ref_meam[elti];
}
}
rhob = rho[i] / rho_bkgd;
denom = 1.0 / rho_bkgd;
G = dG_gam(gamma[i], this->ibar_meam[elti], dG);
G = dG_gam(gamma[i], ibar_meam[elti], dG);
dgamma1[i] = (G - 2 * dG * gamma[i]) * denom;
@ -113,13 +113,13 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
// dgamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (this->mix_ref_t == 1) {
if (mix_ref_t == 1) {
dgamma3[i] = rho0[i] * G * dGbar / (Gbar * Z * Z) * denom;
} else {
dgamma3[i] = 0.0;
}
Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);
Fl = embedding(A_meam[elti], Ec_meam[elti][elti], rhob, frhop[i]);
if (eflag_either != 0) {
Fl *= scaleii;
if (eflag_global != 0) {
@ -144,21 +144,21 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
rho3[i] = rho3[i] - 3.0 / 5.0 * arho3b[i][m] * arho3b[i][m];
}
for (m = 0; m < 6; m++) {
rho2[i] = rho2[i] + this->v2D[m] * arho2[i][m] * arho2[i][m];
rho2[i] = rho2[i] + v2D[m] * arho2[i][m] * arho2[i][m];
}
for (m = 0; m < 10; m++) {
rho3[i] = rho3[i] + this->v3D[m] * arho3[i][m] * arho3[i][m];
rho3[i] = rho3[i] + v3D[m] * arho3[i][m] * arho3[i][m];
}
if (rho0[i] > 0.0) {
if (this->ialloy == 1) {
if (ialloy == 1) {
t_ave[i][0] = fdiv_zero(t_ave[i][0], tsq_ave[i][0]);
t_ave[i][1] = fdiv_zero(t_ave[i][1], tsq_ave[i][1]);
t_ave[i][2] = fdiv_zero(t_ave[i][2], tsq_ave[i][2]);
} else if (this->ialloy == 2) {
t_ave[i][0] = this->t1_meam[elti];
t_ave[i][1] = this->t2_meam[elti];
t_ave[i][2] = this->t3_meam[elti];
} else if (ialloy == 2) {
t_ave[i][0] = t1_meam[elti];
t_ave[i][1] = t2_meam[elti];
t_ave[i][2] = t3_meam[elti];
} else {
t_ave[i][0] = t_ave[i][0] / rho0[i];
t_ave[i][1] = t_ave[i][1] / rho0[i];
@ -172,48 +172,48 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
}
Z = get_Zij(this->lattce_meam[elti][elti]);
Z = get_Zij(lattce_meam[elti][elti]);
G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
G = G_gam(gamma[i], ibar_meam[elti], errorflag);
if (errorflag != 0)
return;
get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp);
get_shpfcn(lattce_meam[elti][elti], stheta_meam[elti][elti], ctheta_meam[elti][elti], shp);
if (this->ibar_meam[elti] <= 0) {
if (ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
if (this->mix_ref_t == 1) {
if (mix_ref_t == 1) {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
} else {
gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /
gam = (t1_meam[elti] * shp[0] + t2_meam[elti] * shp[1] + t3_meam[elti] * shp[2]) /
(Z * Z);
}
Gbar = G_gam(gam, this->ibar_meam[elti], errorflag);
Gbar = G_gam(gam, ibar_meam[elti], errorflag);
}
rho[i] = rho0[i] * G;
if (this->mix_ref_t == 1) {
if (this->ibar_meam[elti] <= 0) {
if (mix_ref_t == 1) {
if (ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
Gbar = dG_gam(gam, this->ibar_meam[elti], dGbar);
Gbar = dG_gam(gam, ibar_meam[elti], dGbar);
}
rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
rho_bkgd = rho0_meam[elti] * Z * Gbar;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd = this->rho0_meam[elti] * Z;
if (bkgd_dyn == 1) {
rho_bkgd = rho0_meam[elti] * Z;
} else {
rho_bkgd = this->rho_ref_meam[elti];
rho_bkgd = rho_ref_meam[elti];
}
}
rhob = rho[i] / rho_bkgd;
denom = 1.0 / rho_bkgd;
G = dG_gam(gamma[i], this->ibar_meam[elti], dG);
G = dG_gam(gamma[i], ibar_meam[elti], dG);
dgamma1[i] = (G - 2 * dG * gamma[i]) * denom;
@ -226,13 +226,13 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
// dgamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (this->mix_ref_t == 1) {
if (mix_ref_t == 1) {
dgamma3[i] = rho0[i] * G * dGbar / (Gbar * Z * Z) * denom;
} else {
dgamma3[i] = 0.0;
}
Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);
Fl = embedding(A_meam[elti], Ec_meam[elti][elti], rhob, frhop[i]);
if (eflag_either != 0) {
Fl *= scaleii;

177
src/MEAM/meam_dens_init.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -20,8 +19,7 @@
using namespace LAMMPS_NS;
void
MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
void MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
{
int i, j;
@ -46,7 +44,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
memory->destroy(t_ave);
memory->destroy(tsq_ave);
// msmeam params
if (this->msmeamflag) {
if (msmeamflag) {
memory->destroy(arho1m);
memory->destroy(arho2m);
memory->destroy(arho3m);
@ -74,7 +72,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
memory->create(t_ave, nmax, 3, "pair:t_ave");
memory->create(tsq_ave, nmax, 3, "pair:tsq_ave");
// msmeam params
if (this->msmeamflag) {
if (msmeamflag) {
memory->create(arho1m, nmax, 3, "pair:arho1m");
memory->create(arho2m, nmax, 6, "pair:arho2m");
memory->create(arho3m, nmax, 10, "pair:arho3m");
@ -99,36 +97,27 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
if (this->msmeamflag) {
if (msmeamflag) {
arho2mb[i] = 0.0;
arho1m[i][0] = arho1m[i][1] = arho1m[i][2] = 0.0;
}
for (j = 0; j < 6; j++) {
arho2[i][j] = 0.0;
if (this->msmeamflag) {
arho2m[i][j] = 0.0;
}
if (msmeamflag) { arho2m[i][j] = 0.0; }
}
for (j = 0; j < 10; j++) {
arho3[i][j] = 0.0;
if (this->msmeamflag) {
arho3m[i][j] = 0.0;
}
if (msmeamflag) { arho3m[i][j] = 0.0; }
}
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
if (this->msmeamflag) {
arho3mb[i][0] = arho3mb[i][1] = arho3mb[i][2] = 0.0;
}
if (msmeamflag) { arho3mb[i][0] = arho3mb[i][1] = arho3mb[i][2] = 0.0; }
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
}
void
MEAM::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::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)
{
// Compute screening function and derivatives
getscreen(i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, numneigh, firstneigh,
@ -140,9 +129,9 @@ MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
int* firstneigh, int numneigh_full, int* firstneigh_full, int /*ntype*/, int* type, int* fmap)
void MEAM::getscreen(int i, double *scrfcn, double *dscrfcn, double *fcpair, double **x,
int numneigh, int *firstneigh, int numneigh_full, int *firstneigh_full,
int /*ntype*/, int *type, int *fmap)
{
int jn, j, kn, k;
int elti, eltj, eltk;
@ -154,7 +143,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
double dCikj;
double rnorm, fc, dfc, drinv;
drinv = 1.0 / this->delr_meam;
drinv = 1.0 / delr_meam;
elti = fmap[type[i]];
if (elti < 0) return;
@ -177,16 +166,16 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
delzij = zjtmp - zitmp;
rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
if (rij2 > this->cutforcesq) {
if (rij2 > cutforcesq) {
dscrfcn[jn] = 0.0;
scrfcn[jn] = 0.0;
fcpair[jn] = 0.0;
continue;
}
const double rbound = this->ebound_meam[elti][eltj] * rij2;
const double rbound = ebound_meam[elti][eltj] * rij2;
rij = sqrt(rij2);
rnorm = (this->cutforce - rij) * drinv;
rnorm = (cutforce - rij) * drinv;
sij = 1.0;
// if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
@ -220,8 +209,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
if (a <= 0.0) continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = this->Cmax_meam[elti][eltj][eltk];
Cmin = this->Cmin_meam[elti][eltj][eltk];
Cmax = Cmax_meam[elti][eltj][eltk];
Cmin = Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax) continue;
// note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case here
@ -271,8 +260,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
if (a <= 0.0) continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = this->Cmax_meam[elti][eltj][eltk];
Cmin = this->Cmin_meam[elti][eltj][eltk];
Cmax = Cmax_meam[elti][eltj][eltk];
Cmin = Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax) {
continue;
// Note that cikj may be slightly negative (within numerical
@ -303,9 +292,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* scrfcn, double* fcpair)
void MEAM::calc_rho1(int i, int /*ntype*/, int *type, int *fmap, double **x, int numneigh,
int *firstneigh, double *scrfcn, double *fcpair)
{
int jn, j, m, n, p, elti, eltj;
int nv2, nv3;
@ -330,58 +318,58 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
delij[1] = x[j][1] - ytmp;
delij[2] = x[j][2] - ztmp;
rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
if (rij2 < this->cutforcesq) {
if (rij2 < cutforcesq) {
eltj = fmap[type[j]];
rij = sqrt(rij2);
ai = rij / this->re_meam[elti][elti] - 1.0;
aj = rij / this->re_meam[eltj][eltj] - 1.0;
ro0i = this->rho0_meam[elti];
ro0j = this->rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj) * sij;
rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj) * sij;
rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj) * sij;
rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj) * sij;
if (this->msmeamflag){
rhoa1mj = ro0j * this->t1m_meam[eltj] * MathSpecial::fm_exp(-this->beta1m_meam[eltj] * aj) * sij;
rhoa2mj = ro0j * this->t2m_meam[eltj] * MathSpecial::fm_exp(-this->beta2m_meam[eltj] * aj) * sij;
rhoa3mj = ro0j * this->t3m_meam[eltj] * MathSpecial::fm_exp(-this->beta3m_meam[eltj] * aj) * sij;
ai = rij / re_meam[elti][elti] - 1.0;
aj = rij / re_meam[eltj][eltj] - 1.0;
ro0i = rho0_meam[elti];
ro0j = rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-beta0_meam[eltj] * aj) * sij;
rhoa1j = ro0j * MathSpecial::fm_exp(-beta1_meam[eltj] * aj) * sij;
rhoa2j = ro0j * MathSpecial::fm_exp(-beta2_meam[eltj] * aj) * sij;
rhoa3j = ro0j * MathSpecial::fm_exp(-beta3_meam[eltj] * aj) * sij;
if (msmeamflag) {
rhoa1mj = ro0j * t1m_meam[eltj] * MathSpecial::fm_exp(-beta1m_meam[eltj] * aj) * sij;
rhoa2mj = ro0j * t2m_meam[eltj] * MathSpecial::fm_exp(-beta2m_meam[eltj] * aj) * sij;
rhoa3mj = ro0j * t3m_meam[eltj] * MathSpecial::fm_exp(-beta3m_meam[eltj] * aj) * sij;
}
rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai) * sij;
rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai) * sij;
rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai) * sij;
rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai) * sij;
if (this->msmeamflag){
rhoa1mi = ro0i * this->t1m_meam[elti] * MathSpecial::fm_exp(-this->beta1m_meam[elti] * ai) * sij;
rhoa2mi = ro0i * this->t2m_meam[elti] * MathSpecial::fm_exp(-this->beta2m_meam[elti] * ai) * sij;
rhoa3mi = ro0i * this->t3m_meam[elti] * MathSpecial::fm_exp(-this->beta3m_meam[elti] * ai) * sij;
rhoa0i = ro0i * MathSpecial::fm_exp(-beta0_meam[elti] * ai) * sij;
rhoa1i = ro0i * MathSpecial::fm_exp(-beta1_meam[elti] * ai) * sij;
rhoa2i = ro0i * MathSpecial::fm_exp(-beta2_meam[elti] * ai) * sij;
rhoa3i = ro0i * MathSpecial::fm_exp(-beta3_meam[elti] * ai) * sij;
if (msmeamflag) {
rhoa1mi = ro0i * t1m_meam[elti] * MathSpecial::fm_exp(-beta1m_meam[elti] * ai) * sij;
rhoa2mi = ro0i * t2m_meam[elti] * MathSpecial::fm_exp(-beta2m_meam[elti] * ai) * sij;
rhoa3mi = ro0i * t3m_meam[elti] * MathSpecial::fm_exp(-beta3m_meam[elti] * ai) * sij;
}
if (this->ialloy == 1) {
rhoa1j = rhoa1j * this->t1_meam[eltj];
rhoa2j = rhoa2j * this->t2_meam[eltj];
rhoa3j = rhoa3j * this->t3_meam[eltj];
rhoa1i = rhoa1i * this->t1_meam[elti];
rhoa2i = rhoa2i * this->t2_meam[elti];
rhoa3i = rhoa3i * this->t3_meam[elti];
if (ialloy == 1) {
rhoa1j = rhoa1j * t1_meam[eltj];
rhoa2j = rhoa2j * t2_meam[eltj];
rhoa3j = rhoa3j * t3_meam[eltj];
rhoa1i = rhoa1i * t1_meam[elti];
rhoa2i = rhoa2i * t2_meam[elti];
rhoa3i = rhoa3i * t3_meam[elti];
}
rho0[i] = rho0[i] + rhoa0j;
rho0[j] = rho0[j] + rhoa0i;
// For ialloy = 2, use single-element value (not average)
// For ialloy = 2, use single-element value (not average)
if (this->ialloy != 2) {
t_ave[i][0] = t_ave[i][0] + this->t1_meam[eltj] * rhoa0j;
t_ave[i][1] = t_ave[i][1] + this->t2_meam[eltj] * rhoa0j;
t_ave[i][2] = t_ave[i][2] + this->t3_meam[eltj] * rhoa0j;
t_ave[j][0] = t_ave[j][0] + this->t1_meam[elti] * rhoa0i;
t_ave[j][1] = t_ave[j][1] + this->t2_meam[elti] * rhoa0i;
t_ave[j][2] = t_ave[j][2] + this->t3_meam[elti] * rhoa0i;
if (ialloy != 2) {
t_ave[i][0] = t_ave[i][0] + t1_meam[eltj] * rhoa0j;
t_ave[i][1] = t_ave[i][1] + t2_meam[eltj] * rhoa0j;
t_ave[i][2] = t_ave[i][2] + t3_meam[eltj] * rhoa0j;
t_ave[j][0] = t_ave[j][0] + t1_meam[elti] * rhoa0i;
t_ave[j][1] = t_ave[j][1] + t2_meam[elti] * rhoa0i;
t_ave[j][2] = t_ave[j][2] + t3_meam[elti] * rhoa0i;
}
if (this->ialloy == 1) {
tsq_ave[i][0] = tsq_ave[i][0] + this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j;
tsq_ave[i][1] = tsq_ave[i][1] + this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j;
tsq_ave[i][2] = tsq_ave[i][2] + this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j;
tsq_ave[j][0] = tsq_ave[j][0] + this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i;
tsq_ave[j][1] = tsq_ave[j][1] + this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i;
tsq_ave[j][2] = tsq_ave[j][2] + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i;
if (ialloy == 1) {
tsq_ave[i][0] = tsq_ave[i][0] + t1_meam[eltj] * t1_meam[eltj] * rhoa0j;
tsq_ave[i][1] = tsq_ave[i][1] + t2_meam[eltj] * t2_meam[eltj] * rhoa0j;
tsq_ave[i][2] = tsq_ave[i][2] + t3_meam[eltj] * t3_meam[eltj] * rhoa0j;
tsq_ave[j][0] = tsq_ave[j][0] + t1_meam[elti] * t1_meam[elti] * rhoa0i;
tsq_ave[j][1] = tsq_ave[j][1] + t2_meam[elti] * t2_meam[elti] * rhoa0i;
tsq_ave[j][2] = tsq_ave[j][2] + t3_meam[elti] * t3_meam[elti] * rhoa0i;
}
arho2b[i] = arho2b[i] + rhoa2j;
arho2b[j] = arho2b[j] + rhoa2i;
@ -394,41 +382,41 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
A3i = rhoa3i / (rij2 * rij);
nv2 = 0;
nv3 = 0;
if (this->msmeamflag) {
if (msmeamflag) {
arho2mb[i] = arho2mb[i] + rhoa2mj;
arho2mb[j] = arho2mb[j] + rhoa2mi;
A1mj = rhoa1mj/rij;
A2mj = rhoa2mj/rij2;
A3mj = rhoa3mj/(rij2*rij);
A1mi = rhoa1mi/rij;
A2mi = rhoa2mi/rij2;
A3mi = rhoa3mi/(rij2*rij);
A1mj = rhoa1mj / rij;
A2mj = rhoa2mj / rij2;
A3mj = rhoa3mj / (rij2 * rij);
A1mi = rhoa1mi / rij;
A2mi = rhoa2mi / rij2;
A3mi = rhoa3mi / (rij2 * rij);
}
for (m = 0; m < 3; m++) {
arho1[i][m] = arho1[i][m] + A1j * delij[m];
arho1[j][m] = arho1[j][m] - A1i * delij[m];
arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij;
arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij;
if (this->msmeamflag) {
arho1m[i][m] = arho1m[i][m] + A1mj*delij[m];
arho1m[j][m] = arho1m[j][m] - A1mi*delij[m];
arho3mb[i][m] = arho3mb[i][m] + rhoa3mj*delij[m] / rij;
arho3mb[j][m] = arho3mb[j][m] - rhoa3mi*delij[m] / rij;
if (msmeamflag) {
arho1m[i][m] = arho1m[i][m] + A1mj * delij[m];
arho1m[j][m] = arho1m[j][m] - A1mi * delij[m];
arho3mb[i][m] = arho3mb[i][m] + rhoa3mj * delij[m] / rij;
arho3mb[j][m] = arho3mb[j][m] - rhoa3mi * delij[m] / rij;
}
for (n = m; n < 3; n++) {
arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n];
arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n];
if (this->msmeamflag) {
arho2m[i][nv2] = arho2m[i][nv2] + A2mj*delij[m] * delij[n];
arho2m[j][nv2] = arho2m[j][nv2] + A2mi*delij[m] * delij[n];
if (msmeamflag) {
arho2m[i][nv2] = arho2m[i][nv2] + A2mj * delij[m] * delij[n];
arho2m[j][nv2] = arho2m[j][nv2] + A2mi * delij[m] * delij[n];
}
nv2 = nv2 + 1;
for (p = n; p < 3; p++) {
arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p];
arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p];
if (this->msmeamflag) {
arho3m[i][nv3] = arho3m[i][nv3] + A3mj*delij[m]*delij[n]*delij[p];
arho3m[j][nv3] = arho3m[j][nv3] - A3mi*delij[m]*delij[n]*delij[p];
if (msmeamflag) {
arho3m[i][nv3] = arho3m[i][nv3] + A3mj * delij[m] * delij[n] * delij[p];
arho3m[j][nv3] = arho3m[j][nv3] - A3mi * delij[m] * delij[n] * delij[p];
}
nv3 = nv3 + 1;
}
@ -438,4 +426,3 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
}
}
}

166
src/MEAM/meam_force.cpp
View File

@ -15,16 +15,16 @@
#include "math_special.h"
#include <cmath>
#include <algorithm>
#include <cmath>
using namespace LAMMPS_NS;
void
MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom,
double* eng_vdwl, 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, double *virial)
void MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom,
double *eng_vdwl, 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,
double *virial)
{
int j, jn, k, kn, kk, m, n, p, q;
int nv2, nv3, elti, eltj, eltk, ind;
@ -98,18 +98,18 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
delij[1] = x[j][1] - yitmp;
delij[2] = x[j][2] - zitmp;
rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
if (rij2 < this->cutforcesq) {
if (rij2 < cutforcesq) {
rij = sqrt(rij2);
recip = 1.0 / rij;
// Compute phi and phip
ind = this->eltind[elti][eltj];
pp = rij * this->rdrar;
ind = eltind[elti][eltj];
pp = rij * rdrar;
kk = (int)pp;
kk = std::min(kk, this->nrar - 2);
kk = std::min(kk, 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];
phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
phi = ((phirar3[ind][kk] * pp + phirar2[ind][kk]) * pp + phirar1[ind][kk]) * pp + phirar[ind][kk];
phip = (phirar6[ind][kk] * pp + phirar5[ind][kk]) * pp + phirar4[ind][kk];
if (eflag_either != 0) {
double phi_sc = phi * scaleij;
@ -126,47 +126,47 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
// Compute pair densities and derivatives
invrei = 1.0 / this->re_meam[elti][elti];
invrei = 1.0 / re_meam[elti][elti];
ai = rij * invrei - 1.0;
ro0i = this->rho0_meam[elti];
rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai);
drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i;
rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai);
drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i;
rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai);
drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i;
rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai);
drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i;
ro0i = rho0_meam[elti];
rhoa0i = ro0i * MathSpecial::fm_exp(-beta0_meam[elti] * ai);
drhoa0i = -beta0_meam[elti] * invrei * rhoa0i;
rhoa1i = ro0i * MathSpecial::fm_exp(-beta1_meam[elti] * ai);
drhoa1i = -beta1_meam[elti] * invrei * rhoa1i;
rhoa2i = ro0i * MathSpecial::fm_exp(-beta2_meam[elti] * ai);
drhoa2i = -beta2_meam[elti] * invrei * rhoa2i;
rhoa3i = ro0i * MathSpecial::fm_exp(-beta3_meam[elti] * ai);
drhoa3i = -beta3_meam[elti] * invrei * rhoa3i;
if (this->msmeamflag) {
rhoa1mi = ro0i * MathSpecial::fm_exp(-this->beta1m_meam[elti] * ai) * t1m_meam[elti];
drhoa1mi = -this->beta1m_meam[elti] * invrei * rhoa1mi;
rhoa2mi = ro0i * MathSpecial::fm_exp(-this->beta2m_meam[elti] * ai) * t2m_meam[elti];
drhoa2mi = -this->beta2m_meam[elti] * invrei * rhoa2mi;
rhoa3mi = ro0i * MathSpecial::fm_exp(-this->beta3m_meam[elti] * ai) * t3m_meam[elti];
drhoa3mi = -this->beta3m_meam[elti] * invrei * rhoa3mi;
if (msmeamflag) {
rhoa1mi = ro0i * MathSpecial::fm_exp(-beta1m_meam[elti] * ai) * t1m_meam[elti];
drhoa1mi = -beta1m_meam[elti] * invrei * rhoa1mi;
rhoa2mi = ro0i * MathSpecial::fm_exp(-beta2m_meam[elti] * ai) * t2m_meam[elti];
drhoa2mi = -beta2m_meam[elti] * invrei * rhoa2mi;
rhoa3mi = ro0i * MathSpecial::fm_exp(-beta3m_meam[elti] * ai) * t3m_meam[elti];
drhoa3mi = -beta3m_meam[elti] * invrei * rhoa3mi;
}
if (elti != eltj) {
invrej = 1.0 / this->re_meam[eltj][eltj];
invrej = 1.0 / re_meam[eltj][eltj];
aj = rij * invrej - 1.0;
ro0j = this->rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj);
drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j;
rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj);
drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j;
rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj);
drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j;
rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj);
drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j;
ro0j = rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-beta0_meam[eltj] * aj);
drhoa0j = -beta0_meam[eltj] * invrej * rhoa0j;
rhoa1j = ro0j * MathSpecial::fm_exp(-beta1_meam[eltj] * aj);
drhoa1j = -beta1_meam[eltj] * invrej * rhoa1j;
rhoa2j = ro0j * MathSpecial::fm_exp(-beta2_meam[eltj] * aj);
drhoa2j = -beta2_meam[eltj] * invrej * rhoa2j;
rhoa3j = ro0j * MathSpecial::fm_exp(-beta3_meam[eltj] * aj);
drhoa3j = -beta3_meam[eltj] * invrej * rhoa3j;
if (this->msmeamflag) {
rhoa1mj = ro0j * t1m_meam[eltj] * MathSpecial::fm_exp(-this->beta1m_meam[eltj] * aj);
drhoa1mj = -this->beta1m_meam[eltj] * invrej * rhoa1mj;
rhoa2mj = ro0j * t2m_meam[eltj] * MathSpecial::fm_exp(-this->beta2m_meam[eltj] * aj);
drhoa2mj = -this->beta2m_meam[eltj] * invrej * rhoa2mj;
rhoa3mj = ro0j * t3m_meam[eltj] * MathSpecial::fm_exp(-this->beta3m_meam[eltj] * aj);
drhoa3mj = -this->beta3m_meam[eltj] * invrej * rhoa3mj;
if (msmeamflag) {
rhoa1mj = ro0j * t1m_meam[eltj] * MathSpecial::fm_exp(-beta1m_meam[eltj] * aj);
drhoa1mj = -beta1m_meam[eltj] * invrej * rhoa1mj;
rhoa2mj = ro0j * t2m_meam[eltj] * MathSpecial::fm_exp(-beta2m_meam[eltj] * aj);
drhoa2mj = -beta2m_meam[eltj] * invrej * rhoa2mj;
rhoa3mj = ro0j * t3m_meam[eltj] * MathSpecial::fm_exp(-beta3m_meam[eltj] * aj);
drhoa3mj = -beta3m_meam[eltj] * invrej * rhoa3mj;
}
} else {
@ -179,7 +179,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
rhoa3j = rhoa3i;
drhoa3j = drhoa3i;
if (this->msmeamflag) {
if (msmeamflag) {
rhoa1mj = rhoa1mi;
drhoa1mj = drhoa1mi;
rhoa2mj = rhoa2mi;
@ -189,17 +189,17 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
}
}
const double t1mi = this->t1_meam[elti];
const double t2mi = this->t2_meam[elti];
const double t3mi = this->t3_meam[elti];
const double t1mj = this->t1_meam[eltj];
const double t2mj = this->t2_meam[eltj];
const double t3mj = this->t3_meam[eltj];
const double t1mi = t1_meam[elti];
const double t2mi = t2_meam[elti];
const double t3mi = t3_meam[elti];
const double t1mj = t1_meam[eltj];
const double t2mj = t2_meam[eltj];
const double t3mj = t3_meam[eltj];
// ialloy mod not needed in MS-MEAM, but similarity here is that we multply rhos by t.
// We did this above with rhoa1mj, rhoa2mj, etc.
if (this->ialloy == 1 || this->msmeamflag) {
if (ialloy == 1 || msmeamflag) {
rhoa1j *= t1mj;
rhoa2j *= t2mj;
rhoa3j *= t3mj;
@ -227,12 +227,12 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
for (q = p; q < 3; q++) {
arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
arg = delij[n] * delij[p] * delij[q] * v3D[nv3];
arg1i3 = arg1i3 + arho3[i][nv3] * arg;
arg1j3 = arg1j3 - arho3[j][nv3] * arg;
nv3 = nv3 + 1;
}
arg = delij[n] * delij[p] * this->v2D[nv2];
arg = delij[n] * delij[p] * v2D[nv2];
arg1i2 = arg1i2 + arho2[i][nv2] * arg;
arg1j2 = arg1j2 + arho2[j][nv2] * arg;
nv2 = nv2 + 1;
@ -255,16 +255,16 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
arg1j3m = 0.0;
arg3i3m = 0.0;
arg3j3m = 0.0;
if (this->msmeamflag) {
if (msmeamflag) {
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
for (q = p; q < 3; q++) {
arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
arg = delij[n] * delij[p] * delij[q] * v3D[nv3];
arg1i3m = arg1i3m - arho3m[i][nv3] * arg;
arg1j3m = arg1j3m + arho3m[j][nv3] * arg;
nv3 = nv3 + 1;
}
arg = delij[n] * delij[p] * this->v2D[nv2];
arg = delij[n] * delij[p] * v2D[nv2];
arg1i2m = arg1i2m + arho2m[i][nv2] * arg;
arg1j2m = arg1j2m + arho2m[j][nv2] * arg;
nv2 = nv2 + 1;
@ -299,8 +299,8 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
drho2drm1[m] = 0.0;
drho2drm2[m] = 0.0;
for (n = 0; n < 3; n++) {
drho2drm1[m] = drho2drm1[m] + arho2[i][this->vind2D[m][n]] * delij[n];
drho2drm2[m] = drho2drm2[m] - arho2[j][this->vind2D[m][n]] * delij[n];
drho2drm1[m] = drho2drm1[m] + arho2[i][vind2D[m][n]] * delij[n];
drho2drm2[m] = drho2drm2[m] - arho2[j][vind2D[m][n]] * delij[n];
}
drho2drm1[m] = a2 * rhoa2j * drho2drm1[m];
drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m];
@ -320,9 +320,9 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
nv2 = 0;
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
arg = delij[n] * delij[p] * this->v2D[nv2];
drho3drm1[m] = drho3drm1[m] + arho3[i][this->vind3D[m][n][p]] * arg;
drho3drm2[m] = drho3drm2[m] + arho3[j][this->vind3D[m][n][p]] * arg;
arg = delij[n] * delij[p] * v2D[nv2];
drho3drm1[m] = drho3drm1[m] + arho3[i][vind3D[m][n][p]] * arg;
drho3drm2[m] = drho3drm2[m] + arho3[j][vind3D[m][n][p]] * arg;
nv2 = nv2 + 1;
}
}
@ -330,7 +330,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
}
if (this->msmeamflag) {
if (msmeamflag) {
// rho1m terms
a1 = 2 * sij / rij;
drho1mdr1 = a1 * (drhoa1mj - rhoa1mj / rij) * arg1i1m;
@ -352,8 +352,8 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
drho2mdrm1[m] = 0.0;
drho2mdrm2[m] = 0.0;
for (n = 0; n < 3; n++) {
drho2mdrm1[m] += arho2m[i][this->vind2D[m][n]] * delij[n];
drho2mdrm2[m] -= arho2m[j][this->vind2D[m][n]] * delij[n];
drho2mdrm1[m] += arho2m[i][vind2D[m][n]] * delij[n];
drho2mdrm2[m] -= arho2m[j][vind2D[m][n]] * delij[n];
}
drho2mdrm1[m] = a2 * rhoa2mj * drho2mdrm1[m];
drho2mdrm2[m] = -a2 * rhoa2mi * drho2mdrm2[m];
@ -376,9 +376,9 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
nv2 = 0;
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
arg = delij[n] * delij[p] * this->v2D[nv2];
drho3mdrm1[m] += arho3m[i][this->vind3D[m][n][p]] * arg;
drho3mdrm2[m] += arho3m[j][this->vind3D[m][n][p]] * arg;
arg = delij[n] * delij[p] * v2D[nv2];
drho3mdrm1[m] += arho3m[i][vind3D[m][n][p]] * arg;
drho3mdrm2[m] += arho3m[j][vind3D[m][n][p]] * arg;
nv2 = nv2 + 1;
}
}
@ -399,7 +399,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
// compute derivatives of weighting functions t wrt rij
// weighting functions t set to unity for MS-MEAM
if (this->msmeamflag) {
if (msmeamflag) {
t1i = 1.0;
t2i = 1.0;
@ -423,7 +423,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
t2j = t_ave[j][1];
t3j = t_ave[j][2];
if (this->ialloy == 1) {
if (ialloy == 1) {
a1i = fdiv_zero(drhoa0j * sij, tsq_ave[i][0]);
a1j = fdiv_zero(drhoa0i * sij, tsq_ave[j][0]);
@ -439,7 +439,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
dt3dr1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
dt3dr2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
} else if (this->ialloy == 2) {
} else if (ialloy == 2) {
dt1dr1 = 0.0;
dt1dr2 = 0.0;
@ -468,10 +468,10 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
}
// Compute derivatives of total density wrt rij, sij and rij(3)
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);
get_shpfcn(lattce_meam[elti][elti], stheta_meam[elti][elti], ctheta_meam[elti][elti], shpi);
get_shpfcn(lattce_meam[eltj][eltj], stheta_meam[elti][elti], ctheta_meam[elti][elti], shpj);
if (this->msmeamflag) {
if (msmeamflag) {
drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * (drho1dr1 - drho1mdr1) +
dt2dr1 * rho2[i] + t2i * (drho2dr1 - drho2mdr1) +
@ -526,7 +526,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
if (this->msmeamflag) {
if (msmeamflag) {
drho1mds1 = a1 * rhoa1mj * arg1i1m;
drho1mds2 = a1 * rhoa1mi * arg1j1m;
drho2mds1 = a2 * rhoa2mj * arg1i2m - 2.0 / 3.0 * arho2mb[i] * rhoa2mj;
@ -544,7 +544,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
drho3mds2 = 0.0;
}
if (this->ialloy == 1) {
if (ialloy == 1) {
a1i = fdiv_zero(rhoa0j, tsq_ave[i][0]);
a1j = fdiv_zero(rhoa0i, tsq_ave[j][0]);
@ -560,7 +560,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
dt3ds1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
dt3ds2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
} else if (this->ialloy == 2) {
} else if (ialloy == 2) {
dt1ds1 = 0.0;
dt1ds2 = 0.0;
@ -586,7 +586,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
dt3ds2 = aj * (t3mi - t3j);
}
if (this->msmeamflag) {
if (msmeamflag) {
drhods1 = dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * (drho1ds1 - drho1mds1) +
dt2ds1 * rho2[i] + t2i * (drho2ds1 - drho2mds1) +
@ -681,13 +681,13 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
double delc, rik2, rjk2;
sij = scrfcn[jn+fnoffset] * fcpair[jn+fnoffset];
const double Cmax = this->Cmax_meam[elti][eltj][eltk];
const double Cmin = this->Cmin_meam[elti][eltj][eltk];
const double Cmax = Cmax_meam[elti][eltj][eltk];
const double Cmin = Cmin_meam[elti][eltj][eltk];
dsij1 = 0.0;
dsij2 = 0.0;
if (!iszero(sij) && !isone(sij)) {
const double rbound = rij2 * this->ebound_meam[elti][eltj];
const double rbound = rij2 * ebound_meam[elti][eltj];
delc = Cmax - Cmin;
dxjk = x[k][0] - x[j][0];
dyjk = x[k][1] - x[j][1];

242
src/MEAM/meam_funcs.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -24,7 +23,6 @@
using namespace LAMMPS_NS;
//-----------------------------------------------------------------------------
// Compute G(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+gamma)
@ -34,8 +32,7 @@ using namespace LAMMPS_NS;
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
//
double
MEAM::G_gam(const double gamma, const int ibar, int& errorflag) const
double MEAM::G_gam(const double gamma, const int ibar, int &errorflag) const
{
double gsmooth_switchpoint;
@ -76,8 +73,7 @@ MEAM::G_gam(const double gamma, const int ibar, int& errorflag) const
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
//
double
MEAM::dG_gam(const double gamma, const int ibar, double& dG) const
double MEAM::dG_gam(const double gamma, const int ibar, double &dG) const
{
double gsmooth_switchpoint;
double G;
@ -125,12 +121,11 @@ MEAM::dG_gam(const double gamma, const int ibar, double& dG) const
//-----------------------------------------------------------------------------
// Compute ZBL potential
//
double
MEAM::zbl(const double r, const int z1, const int z2)
double MEAM::zbl(const double r, const int z1, const int z2)
{
int i;
const double c[] = { 0.028171, 0.28022, 0.50986, 0.18175 };
const double d[] = { 0.20162, 0.40290, 0.94229, 3.1998 };
const double c[] = {0.028171, 0.28022, 0.50986, 0.18175};
const double d[] = {0.20162, 0.40290, 0.94229, 3.1998};
const double azero = 0.4685;
const double cc = 14.3997;
double a, x;
@ -138,33 +133,29 @@ MEAM::zbl(const double r, const int z1, const int z2)
a = azero / (pow(z1, 0.23) + pow(z2, 0.23));
double result = 0.0;
x = r / a;
for (i = 0; i <= 3; i++) {
result = result + c[i] * MathSpecial::fm_exp(-d[i] * x);
}
if (r > 0.0)
result = result * z1 * z2 / r * cc;
for (i = 0; i <= 3; i++) { result = result + c[i] * MathSpecial::fm_exp(-d[i] * x); }
if (r > 0.0) result = result * z1 * z2 / r * cc;
return result;
}
//-----------------------------------------------------------------------------
// Compute embedding function F(rhobar) and derivative F'(rhobar), eqn I.5
//
double
MEAM::embedding(const double A, const double Ec, const double rhobar, double& dF) const
double MEAM::embedding(const double A, const double Ec, const double rhobar, double &dF) const
{
const double AEc = A * Ec;
if (rhobar > 0.0) {
const double lrb = log(rhobar);
dF = AEc * (1.0 + lrb);
return AEc * rhobar * lrb;
const double lrb = log(rhobar);
dF = AEc * (1.0 + lrb);
return AEc * rhobar * lrb;
} else {
if (this->emb_lin_neg == 0) {
if (emb_lin_neg == 0) {
dF = 0.0;
return 0.0;
} else {
dF = - AEc;
return - AEc * rhobar;
dF = -AEc;
return -AEc * rhobar;
}
}
}
@ -172,9 +163,8 @@ MEAM::embedding(const double A, const double Ec, const double rhobar, double& dF
//-----------------------------------------------------------------------------
// Compute Rose energy function, I.16
//
double
MEAM::erose(const double r, const double re, const double alpha, const double Ec, const double repuls,
const double attrac, const int form)
double MEAM::erose(const double r, const double re, const double alpha, const double Ec,
const double repuls, const double attrac, const int form)
{
double astar, a3;
double result = 0.0;
@ -188,11 +178,13 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec
a3 = repuls;
if (form == 1)
result = -Ec * (1 + astar + (-attrac + repuls / r) * MathSpecial::cube(astar)) * MathSpecial::fm_exp(-astar);
result = -Ec * (1 + astar + (-attrac + repuls / r) * MathSpecial::cube(astar)) *
MathSpecial::fm_exp(-astar);
else if (form == 2)
result = -Ec * (1 + astar + a3 * MathSpecial::cube(astar)) * MathSpecial::fm_exp(-astar);
else
result = -Ec * (1 + astar + a3 * MathSpecial::cube(astar) / (r / re)) * MathSpecial::fm_exp(-astar);
result = -Ec * (1 + astar + a3 * MathSpecial::cube(astar) / (r / re)) *
MathSpecial::fm_exp(-astar);
}
return result;
}
@ -200,8 +192,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, const double sthe, const double cthe, double (&s)[3])
void MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3])
{
switch (latt) {
case FCC:
@ -218,7 +209,7 @@ MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, dou
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 CH4: // CH4 actually needs shape factor for diamond for C, dimer for H
case DIA:
case DIA3:
s[0] = 0.0;
@ -229,19 +220,19 @@ MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, dou
s[0] = 1.0;
s[1] = 2.0 / 3.0;
// 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.
s[2] = 0.40; // this is (1-legendre) where legendre = 0.6 in dynamo is accounted.
break;
case LIN: //linear, theta being 180
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[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.
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;
@ -252,8 +243,7 @@ MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, dou
//-----------------------------------------------------------------------------
// Number of first neighbors for reference structure
//
int
MEAM::get_Zij(const lattice_t latt)
int MEAM::get_Zij(const lattice_t latt)
{
switch (latt) {
case FCC:
@ -276,7 +266,7 @@ MEAM::get_Zij(const lattice_t latt)
return 12;
case B2:
return 8;
case CH4: // DYNAMO currently implemented this way while it needs two Z values, 4 and 1
case CH4: // DYNAMO currently implemented this way while it needs two Z values, 4 and 1
return 4;
case LIN:
case ZIG:
@ -293,9 +283,8 @@ MEAM::get_Zij(const lattice_t latt)
// 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,
const double stheta, double& a, double& S)
int MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, const double stheta,
double &a, double &S)
{
double C, x, sijk;
@ -303,91 +292,91 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax,
switch (latt) {
case FCC:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case FCC:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case BCC:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case BCC:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case HCP:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case HCP:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case B1:
case SC:
Zij2 = 12;
a = sqrt(2.0);
numscr = 2;
break;
case B1:
case SC:
Zij2 = 12;
a = sqrt(2.0);
numscr = 2;
break;
case DIA: // 2NN
Zij2 = 12;
a = sqrt(8.0 / 3.0);
numscr = 1;
if (cmin < 0.500001) {
case DIA: // 2NN
Zij2 = 12;
a = sqrt(8.0 / 3.0);
numscr = 1;
if (cmin < 0.500001) {
// call error('can not do 2NN MEAM for dia')
}
break;
}
break;
case DIA3: // 3NN
Zij2 = 12;
a = sqrt(11.0 / 3.0);
numscr = 4;
if (cmin < 0.500001) {
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;
}
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 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 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 ZIG: //zig-zag
Zij2 = 2;
a = 2.0 * stheta;
numscr = 1;
break;
case L12:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case L12:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case B2:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case C11:
// unsupported lattice flag C11 in get_Zij
break;
default:
// unknown lattic flag in get Zij
// call error('Lattice not defined in get_Zij.')
break;
case B2:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case C11:
// unsupported lattice flag C11 in get_Zij
break;
default:
// unknown lattic flag in get Zij
// call error('Lattice not defined in get_Zij.')
break;
}
// Compute screening for each first neighbor
if (latt==DIA3) { // special case for 3NN diamond structure
if (latt == DIA3) { // special case for 3NN diamond structure
C = 1.0;
} else {
C = 4.0 / (a * a) - 1.0;
@ -399,22 +388,21 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax,
return Zij2;
}
int
MEAM::get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S)
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;
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);

85
src/MEAM/meam_impl.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -24,8 +23,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MEAM::MEAM(Memory* mem)
: memory(mem)
MEAM::MEAM(Memory *mem) : memory(mem)
{
phir = phirar = phirar1 = phirar2 = phirar3 = phirar4 = phirar5 = phirar6 = nullptr;
@ -45,15 +43,14 @@ MEAM::MEAM(Memory* mem)
neltypes = 0;
for (int i = 0; i < maxelt; 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] =
t1m_meam[i] = t2m_meam[i] = t3m_meam[i] =
beta1m_meam[i] = beta2m_meam[i] = beta3m_meam[i] = 0.0;
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] = t1m_meam[i] = t2m_meam[i] = t3m_meam[i] = beta1m_meam[i] =
beta2m_meam[i] = beta3m_meam[i] = 0.0;
for (int j = 0; j < maxelt; j++) {
lattce_meam[i][j] = FCC;
Ec_meam[i][j] = re_meam[i][j] = alpha_meam[i][j] = delta_meam[i][j] = ebound_meam[i][j] = attrac_meam[i][j] = repuls_meam[i][j] = 0.0;
Ec_meam[i][j] = re_meam[i][j] = alpha_meam[i][j] = delta_meam[i][j] = ebound_meam[i][j] =
attrac_meam[i][j] = repuls_meam[i][j] = 0.0;
nn2_meam[i][j] = zbl_meam[i][j] = eltind[i][j] = 0;
}
}
@ -63,44 +60,44 @@ MEAM::~MEAM()
{
if (copymode) return;
memory->destroy(this->phirar6);
memory->destroy(this->phirar5);
memory->destroy(this->phirar4);
memory->destroy(this->phirar3);
memory->destroy(this->phirar2);
memory->destroy(this->phirar1);
memory->destroy(this->phirar);
memory->destroy(this->phir);
memory->destroy(phirar6);
memory->destroy(phirar5);
memory->destroy(phirar4);
memory->destroy(phirar3);
memory->destroy(phirar2);
memory->destroy(phirar1);
memory->destroy(phirar);
memory->destroy(phir);
memory->destroy(this->rho);
memory->destroy(this->rho0);
memory->destroy(this->rho1);
memory->destroy(this->rho2);
memory->destroy(this->rho3);
memory->destroy(this->frhop);
memory->destroy(this->gamma);
memory->destroy(this->dgamma1);
memory->destroy(this->dgamma2);
memory->destroy(this->dgamma3);
memory->destroy(this->arho2b);
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(this->arho1);
memory->destroy(this->arho2);
memory->destroy(this->arho3);
memory->destroy(this->arho3b);
memory->destroy(this->t_ave);
memory->destroy(this->tsq_ave);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
memory->destroy(this->scrfcn);
memory->destroy(this->dscrfcn);
memory->destroy(this->fcpair);
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
// msmeam
if (this->msmeamflag){
memory->destroy(this->arho1m);
memory->destroy(this->arho2m);
memory->destroy(this->arho3m);
memory->destroy(this->arho2mb);
memory->destroy(this->arho3mb);
if (msmeamflag) {
memory->destroy(arho1m);
memory->destroy(arho2m);
memory->destroy(arho3m);
memory->destroy(arho2mb);
memory->destroy(arho3mb);
}
}

552
src/MEAM/meam_setup_done.cpp
View File

@ -26,15 +26,15 @@ void MEAM::meam_setup_done(double* cutmax)
int nv2, nv3, m, n, p;
// Force cutoff
this->cutforce = this->rc_meam;
this->cutforcesq = this->cutforce * this->cutforce;
cutforce = rc_meam;
cutforcesq = cutforce * cutforce;
// Pass cutoff back to calling program
*cutmax = this->cutforce;
*cutmax = cutforce;
// Augment t1 term
for (int i = 0; i < maxelt; i++)
this->t1_meam[i] = this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i];
t1_meam[i] = t1_meam[i] + augt1 * 3.0 / 5.0 * t3_meam[i];
// Compute off-diagonal alloy parameters
alloyparams();
@ -44,44 +44,44 @@ void MEAM::meam_setup_done(double* cutmax)
nv3 = 0;
for (m = 0; m < 3; m++) {
for (n = m; n < 3; n++) {
this->vind2D[m][n] = nv2;
this->vind2D[n][m] = nv2;
vind2D[m][n] = nv2;
vind2D[n][m] = nv2;
nv2 = nv2 + 1;
for (p = n; p < 3; p++) {
this->vind3D[m][n][p] = nv3;
this->vind3D[m][p][n] = nv3;
this->vind3D[n][m][p] = nv3;
this->vind3D[n][p][m] = nv3;
this->vind3D[p][m][n] = nv3;
this->vind3D[p][n][m] = nv3;
vind3D[m][n][p] = nv3;
vind3D[m][p][n] = nv3;
vind3D[n][m][p] = nv3;
vind3D[n][p][m] = nv3;
vind3D[p][m][n] = nv3;
vind3D[p][n][m] = nv3;
nv3 = nv3 + 1;
}
}
}
this->v2D[0] = 1;
this->v2D[1] = 2;
this->v2D[2] = 2;
this->v2D[3] = 1;
this->v2D[4] = 2;
this->v2D[5] = 1;
v2D[0] = 1;
v2D[1] = 2;
v2D[2] = 2;
v2D[3] = 1;
v2D[4] = 2;
v2D[5] = 1;
this->v3D[0] = 1;
this->v3D[1] = 3;
this->v3D[2] = 3;
this->v3D[3] = 3;
this->v3D[4] = 6;
this->v3D[5] = 3;
this->v3D[6] = 1;
this->v3D[7] = 3;
this->v3D[8] = 3;
this->v3D[9] = 1;
v3D[0] = 1;
v3D[1] = 3;
v3D[2] = 3;
v3D[3] = 3;
v3D[4] = 6;
v3D[5] = 3;
v3D[6] = 1;
v3D[7] = 3;
v3D[8] = 3;
v3D[9] = 1;
nv2 = 0;
for (m = 0; m < this->neltypes; m++) {
for (n = m; n < this->neltypes; n++) {
this->eltind[m][n] = nv2;
this->eltind[n][m] = nv2;
for (m = 0; m < neltypes; m++) {
for (n = m; n < neltypes; n++) {
eltind[m][n] = nv2;
eltind[n][m] = nv2;
nv2 = nv2 + 1;
}
}
@ -90,8 +90,8 @@ void MEAM::meam_setup_done(double* cutmax)
compute_reference_density();
// Compute pair potentials and setup arrays for interpolation
this->nr = 1000;
this->dr = 1.1 * this->rc_meam / this->nr;
nr = 1000;
dr = 1.1 * rc_meam / nr;
compute_pair_meam();
}
@ -104,52 +104,52 @@ void MEAM::alloyparams()
double eb;
// Loop over pairs
for (i = 0; i < this->neltypes; i++) {
for (j = 0; j < this->neltypes; j++) {
for (i = 0; i < neltypes; i++) {
for (j = 0; j < neltypes; j++) {
// Treat off-diagonal pairs
// If i>j, set all equal to i<j case (which has already been set,
// here or in the input file)
if (i > j) {
this->re_meam[i][j] = this->re_meam[j][i];
this->Ec_meam[i][j] = this->Ec_meam[j][i];
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];
re_meam[i][j] = re_meam[j][i];
Ec_meam[i][j] = Ec_meam[j][i];
alpha_meam[i][j] = alpha_meam[j][i];
lattce_meam[i][j] = lattce_meam[j][i];
nn2_meam[i][j] = 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];
stheta_meam[i][j] = stheta_meam[j][i];
ctheta_meam[i][j] = 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) {
if (iszero(this->Ec_meam[i][j])) {
if (this->lattce_meam[i][j] == L12)
this->Ec_meam[i][j] =
(3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 - this->delta_meam[i][j];
else if (this->lattce_meam[i][j] == C11) {
if (this->lattce_meam[i][i] == DIA)
this->Ec_meam[i][j] =
(2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
if (iszero(Ec_meam[i][j])) {
if (lattce_meam[i][j] == L12)
Ec_meam[i][j] =
(3 * Ec_meam[i][i] + Ec_meam[j][j]) / 4.0 - delta_meam[i][j];
else if (lattce_meam[i][j] == C11) {
if (lattce_meam[i][i] == DIA)
Ec_meam[i][j] =
(2 * Ec_meam[i][i] + Ec_meam[j][j]) / 3.0 - delta_meam[i][j];
else
this->Ec_meam[i][j] =
(this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
Ec_meam[i][j] =
(Ec_meam[i][i] + 2 * Ec_meam[j][j]) / 3.0 - delta_meam[i][j];
} else
this->Ec_meam[i][j] = (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - this->delta_meam[i][j];
Ec_meam[i][j] = (Ec_meam[i][i] + Ec_meam[j][j]) / 2.0 - delta_meam[i][j];
}
if (iszero(this->alpha_meam[i][j]))
this->alpha_meam[i][j] = (this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0;
if (iszero(this->re_meam[i][j]))
this->re_meam[i][j] = (this->re_meam[i][i] + this->re_meam[j][j]) / 2.0;
if (iszero(alpha_meam[i][j]))
alpha_meam[i][j] = (alpha_meam[i][i] + alpha_meam[j][j]) / 2.0;
if (iszero(re_meam[i][j]))
re_meam[i][j] = (re_meam[i][i] + re_meam[j][j]) / 2.0;
}
}
}
// Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets
// where i>j, set equal to the i<j element. Likewise for Cmax.
for (i = 1; i < this->neltypes; i++) {
for (i = 1; i < neltypes; i++) {
for (j = 0; j < i; j++) {
for (k = 0; k < this->neltypes; k++) {
this->Cmin_meam[i][j][k] = this->Cmin_meam[j][i][k];
this->Cmax_meam[i][j][k] = this->Cmax_meam[j][i][k];
for (k = 0; k < neltypes; k++) {
Cmin_meam[i][j][k] = Cmin_meam[j][i][k];
Cmax_meam[i][j][k] = Cmax_meam[j][i][k];
}
}
}
@ -158,11 +158,11 @@ void MEAM::alloyparams()
// atom k definitely lies outside the screening function ellipse (so
// there is no need to calculate its effects). Here, compute it for all
// triplets [i][j][k] so that ebound[i][j] is the maximized over k
for (i = 0; i < this->neltypes; i++) {
for (j = 0; j < this->neltypes; j++) {
for (k = 0; k < this->neltypes; k++) {
eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0));
this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb);
for (i = 0; i < neltypes; i++) {
for (j = 0; j < neltypes; j++) {
for (k = 0; k < neltypes; k++) {
eb = (Cmax_meam[i][j][k] * Cmax_meam[i][j][k]) / (4.0 * (Cmax_meam[i][j][k] - 1.0));
ebound_meam[i][j] = std::max(ebound_meam[i][j], eb);
}
}
}
@ -183,87 +183,87 @@ void MEAM::compute_pair_meam()
double C, s111, s112, s221, S11, S22;
// check for previously allocated arrays and free them
if (this->phir != nullptr)
memory->destroy(this->phir);
if (this->phirar != nullptr)
memory->destroy(this->phirar);
if (this->phirar1 != nullptr)
memory->destroy(this->phirar1);
if (this->phirar2 != nullptr)
memory->destroy(this->phirar2);
if (this->phirar3 != nullptr)
memory->destroy(this->phirar3);
if (this->phirar4 != nullptr)
memory->destroy(this->phirar4);
if (this->phirar5 != nullptr)
memory->destroy(this->phirar5);
if (this->phirar6 != nullptr)
memory->destroy(this->phirar6);
if (phir != nullptr)
memory->destroy(phir);
if (phirar != nullptr)
memory->destroy(phirar);
if (phirar1 != nullptr)
memory->destroy(phirar1);
if (phirar2 != nullptr)
memory->destroy(phirar2);
if (phirar3 != nullptr)
memory->destroy(phirar3);
if (phirar4 != nullptr)
memory->destroy(phirar4);
if (phirar5 != nullptr)
memory->destroy(phirar5);
if (phirar6 != nullptr)
memory->destroy(phirar6);
// allocate memory for array that defines the potential
memory->create(this->phir, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phir");
memory->create(phir, (neltypes * (neltypes + 1)) / 2, nr, "pair:phir");
// allocate coeff memory
memory->create(this->phirar, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar");
memory->create(this->phirar1, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar1");
memory->create(this->phirar2, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar2");
memory->create(this->phirar3, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar3");
memory->create(this->phirar4, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar4");
memory->create(this->phirar5, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar5");
memory->create(this->phirar6, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar6");
memory->create(phirar, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar");
memory->create(phirar1, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar1");
memory->create(phirar2, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar2");
memory->create(phirar3, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar3");
memory->create(phirar4, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar4");
memory->create(phirar5, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar5");
memory->create(phirar6, (neltypes * (neltypes + 1)) / 2, nr, "pair:phirar6");
// loop over pairs of element types
nv2 = 0;
for (a = 0; a < this->neltypes; a++) {
for (b = a; b < this->neltypes; b++) {
for (a = 0; a < neltypes; a++) {
for (b = a; b < neltypes; b++) {
// loop over r values and compute
for (j = 0; j < this->nr; j++) {
r = j * this->dr;
this->phir[nv2][j] = phi_meam(r, a, b);
for (j = 0; j < nr; j++) {
r = j * dr;
phir[nv2][j] = phi_meam(r, a, b);
// if using second-nearest neighbor, solve recursive problem
// (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
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], this->stheta_meam[a][b], arat, scrn);
if (nn2_meam[a][b] == 1) {
Z1 = get_Zij(lattce_meam[a][b]);
Z2 = get_Zij2(lattce_meam[a][b], Cmin_meam[a][a][b],
Cmax_meam[a][a][b], 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
// pairs, but need to include NN2 contributions to those pairs as
// well.
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
this->lattce_meam[a][b] == L12 || this->lattce_meam[a][b] == DIA) {
if (lattce_meam[a][b] == B1 || lattce_meam[a][b] == B2 ||
lattce_meam[a][b] == L12 || lattce_meam[a][b] == DIA) {
rarat = r * arat;
// phi_aa
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], this->stheta_meam[a][a], arat, scrn);
Z1 = get_Zij(lattce_meam[a][a]);
Z2 = get_Zij2(lattce_meam[a][a], Cmin_meam[a][a][a],
Cmax_meam[a][a][a], 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], this->stheta_meam[b][b], arat, scrn);
Z1 = get_Zij(lattce_meam[b][b]);
Z2 = get_Zij2(lattce_meam[b][b], Cmin_meam[b][b][b],
Cmax_meam[b][b][b], 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) {
if (lattce_meam[a][b] == B1 || lattce_meam[a][b] == B2 ||
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], 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], this->stheta_meam[a][b], arat, scrn2);
Z1 = get_Zij(lattce_meam[a][b]);
Z2 = get_Zij2(lattce_meam[a][b], Cmin_meam[a][a][b],
Cmax_meam[a][a][b], stheta_meam[a][b], arat, scrn);
phir[nv2][j] = phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
Z2 = get_Zij2(lattce_meam[a][b], Cmin_meam[b][b][a],
Cmax_meam[b][b][a], stheta_meam[a][b], arat, scrn2);
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
phir[nv2][j] = phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
} else if (this->lattce_meam[a][b] == L12) {
} else if (lattce_meam[a][b] == L12) {
// The L12 case has one last trick; we have to be careful to
// compute
// the correct screening between 2nd-neighbor pairs. 1-1
@ -278,11 +278,11 @@ void MEAM::compute_pair_meam()
get_sijk(C, b, b, a, &s221);
S11 = s111 * s111 * s112 * s112;
S22 = pow(s221, 4);
this->phir[nv2][j] = this->phir[nv2][j] - 0.75 * S11 * phiaa - 0.25 * S22 * phibb;
phir[nv2][j] = phir[nv2][j] - 0.75 * S11 * phiaa - 0.25 * S22 * phibb;
}
} else {
this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);
phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);
}
}
@ -292,14 +292,14 @@ void MEAM::compute_pair_meam()
// else if -3 < astar < -1
// potential is linear combination with zbl potential
// endif
if (this->zbl_meam[a][b] == 1) {
astar = this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0);
if (zbl_meam[a][b] == 1) {
astar = alpha_meam[a][b] * (r / re_meam[a][b] - 1.0);
if (astar <= -3.0)
this->phir[nv2][j] = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
phir[nv2][j] = zbl(r, ielt_meam[a], ielt_meam[b]);
else if (astar > -3.0 && astar < -1.0) {
frac = fcut(1 - (astar + 1.0) / (-3.0 + 1.0));
phizbl = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
this->phir[nv2][j] = frac * this->phir[nv2][j] + (1 - frac) * phizbl;
phizbl = zbl(r, ielt_meam[a], ielt_meam[b]);
phir[nv2][j] = frac * phir[nv2][j] + (1 - frac) * phizbl;
}
}
}
@ -343,12 +343,12 @@ double 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]);
Z1 = get_Zij(lattce_meam[a][a]);
Z2 = get_Zij(lattce_meam[b][b]);
Z12 = get_Zij(lattce_meam[a][b]);
// this function has extra args for msmeam
if (this->msmeamflag) {
if (msmeamflag) {
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32,
&rho1m1, &rho2m1, &rho3m1,
&rho1m2, &rho2m2, &rho3m2);
@ -362,39 +362,39 @@ double MEAM::phi_meam(double r, int a, int b)
return 0.0;
// calculate average weighting factors for the reference structure
if (this->lattce_meam[a][b] == C11) {
if (this->ialloy == 2) {
t11av = this->t1_meam[a];
t12av = this->t1_meam[b];
t21av = this->t2_meam[a];
t22av = this->t2_meam[b];
t31av = this->t3_meam[a];
t32av = this->t3_meam[b];
if (lattce_meam[a][b] == C11) {
if (ialloy == 2) {
t11av = t1_meam[a];
t12av = t1_meam[b];
t21av = t2_meam[a];
t22av = t2_meam[b];
t31av = t3_meam[a];
t32av = t3_meam[b];
} else {
scalfac = 1.0 / (rho01 + rho02);
t11av = scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02);
t11av = scalfac * (t1_meam[a] * rho01 + t1_meam[b] * rho02);
t12av = t11av;
t21av = scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02);
t21av = scalfac * (t2_meam[a] * rho01 + t2_meam[b] * rho02);
t22av = t21av;
t31av = scalfac * (this->t3_meam[a] * rho01 + this->t3_meam[b] * rho02);
t31av = scalfac * (t3_meam[a] * rho01 + t3_meam[b] * rho02);
t32av = t31av;
}
} else {
// average weighting factors for the reference structure, eqn. I.8
get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, this->t1_meam[a], this->t2_meam[a],
this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b,
this->lattce_meam[a][b]);
get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, t1_meam[a], t2_meam[a],
t3_meam[a], t1_meam[b], t2_meam[b], t3_meam[b], r, a, b,
lattce_meam[a][b]);
// with msmeam call twice with different sets of variables
if (this->msmeamflag) {
get_tavref(&t1m1av, &t2m1av, &t3m1av, &t1m2av, &t2m2av, &t3m2av, this->t1m_meam[a], this->t2m_meam[a],
this->t3m_meam[a], this->t1m_meam[b], this->t2m_meam[b], this->t3m_meam[b], r, a, b,
this->lattce_meam[a][b]);
if (msmeamflag) {
get_tavref(&t1m1av, &t2m1av, &t3m1av, &t1m2av, &t2m2av, &t3m2av, t1m_meam[a], t2m_meam[a],
t3m_meam[a], t1m_meam[b], t2m_meam[b], t3m_meam[b], r, a, b,
lattce_meam[a][b]);
}
}
// for c11b structure, calculate background electron densities
if (this->lattce_meam[a][b] == C11) {
latta = this->lattce_meam[a][a];
if (lattce_meam[a][b] == C11) {
latta = lattce_meam[a][a];
if (latta == DIA) {
rhobar1 = MathSpecial::square((Z12 / 2) * (rho02 + rho01))
+ t11av * MathSpecial::square(rho12 - rho11)
@ -418,26 +418,26 @@ double MEAM::phi_meam(double r, int a, int b)
// composition-dependent scaling, equation I.7
// If using mixing rule for t, apply to reference structure; else
// use precomputed values
if (this->mix_ref_t == 1) {
if (this->ibar_meam[a] <= 0)
if (mix_ref_t == 1) {
if (ibar_meam[a] <= 0)
G1 = 1.0;
else {
get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1);
get_shpfcn(lattce_meam[a][a], stheta_meam[a][a], 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);
G1 = G_gam(Gam1, ibar_meam[a], errorflag);
}
if (this->ibar_meam[b] <= 0)
if (ibar_meam[b] <= 0)
G2 = 1.0;
else {
get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b], s2);
get_shpfcn(lattce_meam[b][b], stheta_meam[b][b], 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);
G2 = G_gam(Gam2, ibar_meam[b], errorflag);
}
rho0_1 = this->rho0_meam[a] * Z1 * G1;
rho0_2 = this->rho0_meam[b] * Z2 * G2;
rho0_1 = rho0_meam[a] * Z1 * G1;
rho0_2 = rho0_meam[b] * Z2 * G2;
}
if (this->msmeamflag) {
if (msmeamflag) {
// no additional use of t's here; all included in definitions of rho's for msmeam
Gam1 = rho11 + rho21 + rho31 - (rho1m1 + rho2m1 + rho3m1);
if (rho01 < 1.0e-14)
@ -464,18 +464,18 @@ double MEAM::phi_meam(double r, int a, int b)
Gam2 = Gam2 / (rho02 * rho02);
}
G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
if (this->mix_ref_t == 1) {
G1 = G_gam(Gam1, ibar_meam[a], errorflag);
G2 = G_gam(Gam2, ibar_meam[b], errorflag);
if (mix_ref_t == 1) {
rho_bkgd1 = rho0_1;
rho_bkgd2 = rho0_2;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd1 = this->rho0_meam[a] * Z1;
rho_bkgd2 = this->rho0_meam[b] * Z2;
if (bkgd_dyn == 1) {
rho_bkgd1 = rho0_meam[a] * Z1;
rho_bkgd2 = rho0_meam[b] * Z2;
} else {
rho_bkgd1 = this->rho_ref_meam[a];
rho_bkgd2 = this->rho_ref_meam[b];
rho_bkgd1 = rho_ref_meam[a];
rho_bkgd2 = rho_ref_meam[b];
}
}
rhobar1 = rho01 / rho_bkgd1 * G1;
@ -484,17 +484,17 @@ double MEAM::phi_meam(double r, int a, int b)
// compute embedding functions, eqn I.5
F1 = embedding(this->A_meam[a], this->Ec_meam[a][a], rhobar1, dF);
F2 = embedding(this->A_meam[b], this->Ec_meam[b][b], rhobar2, dF);
F1 = embedding(A_meam[a], Ec_meam[a][a], rhobar1, dF);
F2 = embedding(A_meam[b], Ec_meam[b][b], rhobar2, dF);
// compute Rose function, I.16
Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],
this->attrac_meam[a][b], this->erose_form);
Eu = erose(r, re_meam[a][b], alpha_meam[a][b], Ec_meam[a][b], repuls_meam[a][b],
attrac_meam[a][b], erose_form);
// calculate the pair energy
if (this->lattce_meam[a][b] == C11) {
latta = this->lattce_meam[a][a];
if (lattce_meam[a][b] == C11) {
latta = lattce_meam[a][a];
if (latta == DIA) {
phiaa = phi_meam(r, a, a);
phi_m = (3 * Eu - F2 - 2 * F1 - 5 * phiaa) / Z12;
@ -502,44 +502,44 @@ double MEAM::phi_meam(double r, int a, int b)
phibb = phi_meam(r, b, b);
phi_m = (3 * Eu - F1 - 2 * F2 - 5 * phibb) / Z12;
}
} else if (this->lattce_meam[a][b] == L12) {
} else if (lattce_meam[a][b] == L12) {
phiaa = phi_meam(r, a, a);
// 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], this->stheta_meam[a][b], arat, scrn);
Z1nn = get_Zij(lattce_meam[a][a]);
Z2nn = get_Zij2(lattce_meam[a][a], Cmin_meam[a][a][a],
Cmax_meam[a][a][a], 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) {
} else if (lattce_meam[a][b] == CH4) {
phi_m = (5 * Eu - F1 - 4*F2)/4;
} else if (this->lattce_meam[a][b] == ZIG) {
} else if (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);
Z1 = get_Zij(lattce_meam[a][b]);
Z2 = get_Zij2_b2nn(lattce_meam[a][b], Cmin_meam[a][a][b], 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);
Z2 = get_Zij2_b2nn(lattce_meam[a][b], Cmin_meam[b][b][a], 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);
phiaa = phi_meam(2.0*stheta_meam[a][b]*r, a, a);
phibb = phi_meam(2.0*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) {
} else if (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);
Z1 = get_Zij(lattce_meam[a][b]);
Z2 = get_Zij2_b2nn(lattce_meam[a][b], Cmin_meam[a][a][b], Cmax_meam[a][a][b], scrn);
b11s = -Z2/Z1*scrn;
phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
phiaa = phi_meam(2.0*stheta_meam[a][b]*r, a, a);
phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12;
}
@ -590,32 +590,32 @@ void MEAM::compute_reference_density()
double rho0, rho0_2nn, arat, scrn;
// loop over element types
for (a = 0; a < this->neltypes; a++) {
Z = get_Zij(this->lattce_meam[a][a]);
if (this->ibar_meam[a] <= 0)
for (a = 0; a < neltypes; a++) {
Z = get_Zij(lattce_meam[a][a]);
if (ibar_meam[a] <= 0)
Gbar = 1.0;
else {
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);
get_shpfcn(lattce_meam[a][a], stheta_meam[a][a], ctheta_meam[a][a], shp);
gam = (t1_meam[a] * shp[0] + t2_meam[a] * shp[1] + t3_meam[a] * shp[2]) / (Z * Z);
Gbar = G_gam(gam, ibar_meam[a], errorflag);
}
// The zeroth order density in the reference structure, with
// equilibrium spacing, is just the number of first neighbors times
// the rho0_meam coefficient...
rho0 = this->rho0_meam[a] * Z;
rho0 = rho0_meam[a] * Z;
// ...unless we have unscreened second neighbors, in which case we
// add on the contribution from those (accounting for partial
// 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], this->stheta_meam[a][a], arat, scrn);
rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
if (nn2_meam[a][a] == 1) {
Z2 = get_Zij2(lattce_meam[a][a], Cmin_meam[a][a][a],
Cmax_meam[a][a][a], stheta_meam[a][a], arat, scrn);
rho0_2nn = rho0_meam[a] * MathSpecial::fm_exp(-beta0_meam[a] * (arat - 1));
rho0 = rho0 + Z2 * rho0_2nn * scrn;
}
this->rho_ref_meam[a] = rho0 * Gbar;
rho_ref_meam[a] = rho0 * Gbar;
}
}
@ -628,7 +628,7 @@ void MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av
double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/;
// For ialloy = 2, no averaging is done
if (this->ialloy == 2) {
if (ialloy == 2) {
*t11av = t11;
*t21av = t21;
*t31av = t31;
@ -658,10 +658,10 @@ void MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av
*t32av = t31;
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);
a1 = r / re_meam[a][a] - 1.0;
a2 = r / re_meam[b][b] - 1.0;
rhoa01 = rho0_meam[a] * MathSpecial::fm_exp(-beta0_meam[a] * a1);
rhoa02 = rho0_meam[b] * MathSpecial::fm_exp(-beta0_meam[b] * a2);
if (latt == L12) {
rho01 = 8 * rhoa01 + 4 * rhoa02;
*t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
@ -680,7 +680,7 @@ void MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av
void MEAM::get_sijk(double C, int i, int j, int k, double* sijk)
{
double x;
x = (C - this->Cmin_meam[i][j][k]) / (this->Cmax_meam[i][j][k] - this->Cmin_meam[i][j][k]);
x = (C - Cmin_meam[i][j][k]) / (Cmax_meam[i][j][k] - Cmin_meam[i][j][k]);
*sijk = fcut(x);
}
@ -703,38 +703,38 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
// msmeam
double rhoa1m1, rhoa2m1, rhoa3m1, rhoa1m2, rhoa2m2, rhoa3m2;
a1 = r / this->re_meam[a][a] - 1.0;
a2 = r / this->re_meam[b][b] - 1.0;
a1 = r / re_meam[a][a] - 1.0;
a2 = r / re_meam[b][b] - 1.0;
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa01 = rho0_meam[a] * MathSpecial::fm_exp(-beta0_meam[a] * a1);
if (this->msmeamflag) {
if (msmeamflag) {
// the rho variables are multiplied by t here since ialloy not needed in msmeam
rhoa11 = this->rho0_meam[a] * this->t1_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1);
rhoa21 = this->rho0_meam[a] * this->t2_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1);
rhoa31 = this->rho0_meam[a] * this->t3_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
rhoa12 = this->rho0_meam[b] * this->t1_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2);
rhoa22 = this->rho0_meam[b] * this->t2_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2);
rhoa32 = this->rho0_meam[b] * this->t3_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2);
rhoa11 = rho0_meam[a] * t1_meam[a] * MathSpecial::fm_exp(-beta1_meam[a] * a1);
rhoa21 = rho0_meam[a] * t2_meam[a] * MathSpecial::fm_exp(-beta2_meam[a] * a1);
rhoa31 = rho0_meam[a] * t3_meam[a] * MathSpecial::fm_exp(-beta3_meam[a] * a1);
rhoa02 = rho0_meam[b] * MathSpecial::fm_exp(-beta0_meam[b] * a2);
rhoa12 = rho0_meam[b] * t1_meam[b] * MathSpecial::fm_exp(-beta1_meam[b] * a2);
rhoa22 = rho0_meam[b] * t2_meam[b] * MathSpecial::fm_exp(-beta2_meam[b] * a2);
rhoa32 = rho0_meam[b] * t3_meam[b] * MathSpecial::fm_exp(-beta3_meam[b] * a2);
// msmeam specific rho vars
rhoa1m1 = this->rho0_meam[a] * this->t1m_meam[a] * MathSpecial::fm_exp(-this->beta1m_meam[a] * a1);
rhoa2m1 = this->rho0_meam[a] * this->t2m_meam[a] * MathSpecial::fm_exp(-this->beta2m_meam[a] * a1);
rhoa3m1 = this->rho0_meam[a] * this->t3m_meam[a] * MathSpecial::fm_exp(-this->beta3m_meam[a] * a1);
rhoa1m2 = this->rho0_meam[b] * this->t1m_meam[b] * MathSpecial::fm_exp(-this->beta1m_meam[b] * a2);
rhoa2m2 = this->rho0_meam[b] * this->t2m_meam[b] * MathSpecial::fm_exp(-this->beta2m_meam[b] * a2);
rhoa3m2 = this->rho0_meam[b] * this->t3m_meam[b] * MathSpecial::fm_exp(-this->beta3m_meam[b] * a2);
rhoa1m1 = rho0_meam[a] * t1m_meam[a] * MathSpecial::fm_exp(-beta1m_meam[a] * a1);
rhoa2m1 = rho0_meam[a] * t2m_meam[a] * MathSpecial::fm_exp(-beta2m_meam[a] * a1);
rhoa3m1 = rho0_meam[a] * t3m_meam[a] * MathSpecial::fm_exp(-beta3m_meam[a] * a1);
rhoa1m2 = rho0_meam[b] * t1m_meam[b] * MathSpecial::fm_exp(-beta1m_meam[b] * a2);
rhoa2m2 = rho0_meam[b] * t2m_meam[b] * MathSpecial::fm_exp(-beta2m_meam[b] * a2);
rhoa3m2 = rho0_meam[b] * t3m_meam[b] * MathSpecial::fm_exp(-beta3m_meam[b] * a2);
} else {
rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1);
rhoa21 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1);
rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
rhoa12 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2);
rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2);
rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2);
rhoa11 = rho0_meam[a] * MathSpecial::fm_exp(-beta1_meam[a] * a1);
rhoa21 = rho0_meam[a] * MathSpecial::fm_exp(-beta2_meam[a] * a1);
rhoa31 = rho0_meam[a] * MathSpecial::fm_exp(-beta3_meam[a] * a1);
rhoa02 = rho0_meam[b] * MathSpecial::fm_exp(-beta0_meam[b] * a2);
rhoa12 = rho0_meam[b] * MathSpecial::fm_exp(-beta1_meam[b] * a2);
rhoa22 = rho0_meam[b] * MathSpecial::fm_exp(-beta2_meam[b] * a2);
rhoa32 = rho0_meam[b] * MathSpecial::fm_exp(-beta3_meam[b] * a2);
}
lat = this->lattce_meam[a][b];
lat = lattce_meam[a][b];
Zij = get_Zij(lat);
@ -744,7 +744,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho12 = 0.0;
*rho22 = 0.0;
*rho32 = 0.0;
if (this->msmeamflag) {
if (msmeamflag) {
*rho1m1 = 0.0;
*rho2m1 = 0.0;
*rho3m1 = 0.0;
@ -774,7 +774,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho02 = 4.0 * rhoa01;
*rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
*rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
if (this->msmeamflag) {
if (msmeamflag) {
*rho3m1 = 32.0 / 9.0 * rhoa3m2 * rhoa3m2;
*rho3m2 = 32.0 / 9.0 * rhoa3m1 * rhoa3m1;
}
@ -784,7 +784,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho02 = 12 * rhoa01;
*rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
*rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
if (this->msmeamflag) {
if (msmeamflag) {
*rho3m1 = 1.0 / 3.0 * rhoa3m2 * rhoa3m2;
*rho3m2 = 1.0 / 3.0 * rhoa3m1 * rhoa3m1;
}
@ -799,7 +799,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho22 = s[1] * rhoa21 * rhoa21;
*rho31 = s[2] * rhoa32 * rhoa32;
*rho32 = s[2] * rhoa31 * rhoa31;
if (this->msmeamflag) {
if (msmeamflag) {
*rho1m1 = s[0] * rhoa1m2 * rhoa1m2;
*rho1m2 = s[0] * rhoa1m1 * rhoa1m1;
*rho2m1 = s[1] * rhoa2m2 * rhoa2m2;
@ -817,7 +817,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho22 = rhoa22;
*rho31 = rhoa31;
*rho32 = rhoa32;
if (this->msmeamflag) {
if (msmeamflag) {
*rho1m1 = rhoa1m1;
*rho1m2 = rhoa1m2;
*rho2m1 = rhoa2m1;
@ -829,14 +829,14 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
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 (ialloy ==1){
*rho21 = 8. / 3. * MathSpecial::square(rhoa21 * t2_meam[a] - rhoa22 * t2_meam[b]);
denom = 8 * rhoa01 * MathSpecial::square(t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(t2_meam[b]);
if (denom > 0.)
*rho21 = *rho21 / denom * *rho01;
} else
*rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
if (this->msmeamflag) {
if (msmeamflag) {
*rho2m1 = 8. / 3. * (rhoa2m1 - rhoa2m2) * (rhoa2m1 - rhoa2m2);
}
break;
@ -862,7 +862,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho01 = rhoa02*Zij;
*rho02 = rhoa01*Zij;
get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
get_shpfcn(LIN, stheta_meam[a][b], ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
@ -874,7 +874,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho01 = rhoa02*Zij;
*rho02 = rhoa01*Zij;
get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
get_shpfcn(ZIG, stheta_meam[a][b], ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
@ -886,7 +886,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho01 = rhoa02;
*rho02 = rhoa01*Zij;
get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
get_shpfcn(TRI, stheta_meam[a][b], ctheta_meam[a][b], s);
*rho12 = s[0] * rhoa11 * rhoa11;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho32 = s[2] * rhoa31 * rhoa31;
@ -904,17 +904,17 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
// call error('Lattice not defined in get_densref.')
}
if (this->nn2_meam[a][b] == 1) {
if (nn2_meam[a][b] == 1) {
Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b],
this->stheta_meam[a][b], arat, scrn);
Zij2nn = get_Zij2(lat, Cmin_meam[a][a][b], Cmax_meam[a][a][b],
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;
a1 = arat * r / re_meam[a][a] - 1.0;
a2 = arat * r / re_meam[b][b] - 1.0;
rhoa01nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa02nn = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
rhoa01nn = rho0_meam[a] * MathSpecial::fm_exp(-beta0_meam[a] * a1);
rhoa02nn = rho0_meam[b] * MathSpecial::fm_exp(-beta0_meam[b] * a2);
if (lat == L12) {
// As usual, L12 thinks it's special; we need to be careful computing
@ -935,8 +935,8 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*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],
this->stheta_meam[a][b], arat, scrn);
Zij2nn = get_Zij2(lat, Cmin_meam[b][b][a], Cmax_meam[b][b][a],
stheta_meam[a][b], arat, scrn);
*rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
}
}
@ -950,38 +950,38 @@ void MEAM::interpolate_meam(int ind)
// map to coefficient space
this->nrar = this->nr;
drar = this->dr;
this->rdrar = 1.0 / drar;
nrar = nr;
drar = dr;
rdrar = 1.0 / drar;
// phir interp
for (j = 0; j < this->nrar; j++) {
this->phirar[ind][j] = this->phir[ind][j];
for (j = 0; j < nrar; j++) {
phirar[ind][j] = phir[ind][j];
}
this->phirar1[ind][0] = this->phirar[ind][1] - this->phirar[ind][0];
this->phirar1[ind][1] = 0.5 * (this->phirar[ind][2] - this->phirar[ind][0]);
this->phirar1[ind][this->nrar - 2] =
0.5 * (this->phirar[ind][this->nrar - 1] - this->phirar[ind][this->nrar - 3]);
this->phirar1[ind][this->nrar - 1] = 0.0;
for (j = 2; j < this->nrar - 2; j++) {
this->phirar1[ind][j] = ((this->phirar[ind][j - 2] - this->phirar[ind][j + 2]) +
8.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j - 1])) /
phirar1[ind][0] = phirar[ind][1] - phirar[ind][0];
phirar1[ind][1] = 0.5 * (phirar[ind][2] - phirar[ind][0]);
phirar1[ind][nrar - 2] =
0.5 * (phirar[ind][nrar - 1] - phirar[ind][nrar - 3]);
phirar1[ind][nrar - 1] = 0.0;
for (j = 2; j < nrar - 2; j++) {
phirar1[ind][j] = ((phirar[ind][j - 2] - phirar[ind][j + 2]) +
8.0 * (phirar[ind][j + 1] - phirar[ind][j - 1])) /
12.;
}
for (j = 0; j < this->nrar - 1; j++) {
this->phirar2[ind][j] = 3.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]) -
2.0 * this->phirar1[ind][j] - this->phirar1[ind][j + 1];
this->phirar3[ind][j] = this->phirar1[ind][j] + this->phirar1[ind][j + 1] -
2.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]);
for (j = 0; j < nrar - 1; j++) {
phirar2[ind][j] = 3.0 * (phirar[ind][j + 1] - phirar[ind][j]) -
2.0 * phirar1[ind][j] - phirar1[ind][j + 1];
phirar3[ind][j] = phirar1[ind][j] + phirar1[ind][j + 1] -
2.0 * (phirar[ind][j + 1] - phirar[ind][j]);
}
this->phirar2[ind][this->nrar - 1] = 0.0;
this->phirar3[ind][this->nrar - 1] = 0.0;
phirar2[ind][nrar - 1] = 0.0;
phirar3[ind][nrar - 1] = 0.0;
for (j = 0; j < this->nrar; j++) {
this->phirar4[ind][j] = this->phirar1[ind][j] / drar;
this->phirar5[ind][j] = 2.0 * this->phirar2[ind][j] / drar;
this->phirar6[ind][j] = 3.0 * this->phirar3[ind][j] / drar;
for (j = 0; j < nrar; j++) {
phirar4[ind][j] = phirar1[ind][j] / drar;
phirar5[ind][j] = 2.0 * phirar2[ind][j] / drar;
phirar6[ind][j] = 3.0 * phirar3[ind][j] / drar;
}
}

126
src/MEAM/meam_setup_global.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -18,68 +17,66 @@
using namespace LAMMPS_NS;
template <typename TYPE, int maxi, int maxj>
static inline void setall2d(TYPE (&arr)[maxi][maxj], const TYPE v) {
static inline void setall2d(TYPE (&arr)[maxi][maxj], const TYPE v)
{
for (int i = 0; i < maxi; i++)
for (int j = 0; j < maxj; j++)
arr[i][j] = v;
for (int j = 0; j < maxj; j++) arr[i][j] = v;
}
template <typename TYPE, int maxi, int maxj, int maxk>
static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v)
{
for (int i = 0; i < maxi; i++)
for (int j = 0; j < maxj; j++)
for (int k = 0; k < maxk; k++)
arr[i][j][k] = v;
for (int k = 0; k < maxk; k++) arr[i][j][k] = v;
}
void
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, double* b1m, double *b2m, double *b3m, double *t1m, double *t2m,
double *t3m)
void 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, double *b1m,
double *b2m, double *b3m, double *t1m, double *t2m, double *t3m)
{
int i;
double tmplat[maxelt];
this->neltypes = nelt;
neltypes = nelt;
for (i = 0; i < nelt; i++) {
this->lattce_meam[i][i] = lat[i];
lattce_meam[i][i] = lat[i];
this->ielt_meam[i] = ielement[i];
this->alpha_meam[i][i] = alpha[i];
this->beta0_meam[i] = b0[i];
this->beta1_meam[i] = b1[i];
this->beta2_meam[i] = b2[i];
this->beta3_meam[i] = b3[i];
if (this->msmeamflag){
this->beta1m_meam[i] = b1m[i];
this->beta2m_meam[i] = b2m[i];
this->beta3m_meam[i] = b3m[i];
ielt_meam[i] = ielement[i];
alpha_meam[i][i] = alpha[i];
beta0_meam[i] = b0[i];
beta1_meam[i] = b1[i];
beta2_meam[i] = b2[i];
beta3_meam[i] = b3[i];
if (msmeamflag) {
beta1m_meam[i] = b1m[i];
beta2m_meam[i] = b2m[i];
beta3m_meam[i] = b3m[i];
}
tmplat[i] = alat[i];
this->Ec_meam[i][i] = esub[i];
this->A_meam[i] = asub[i];
this->t0_meam[i] = t0[i];
this->t1_meam[i] = t1[i];
this->t2_meam[i] = t2[i];
this->t3_meam[i] = t3[i];
if (this->msmeamflag){
this->t1m_meam[i] = t1m[i];
this->t2m_meam[i] = t2m[i];
this->t3m_meam[i] = t3m[i];
Ec_meam[i][i] = esub[i];
A_meam[i] = asub[i];
t0_meam[i] = t0[i];
t1_meam[i] = t1[i];
t2_meam[i] = t2[i];
t3_meam[i] = t3[i];
if (msmeamflag) {
t1m_meam[i] = t1m[i];
t2m_meam[i] = t2m[i];
t3m_meam[i] = t3m[i];
}
this->rho0_meam[i] = rozero[i];
this->ibar_meam[i] = ibar[i];
rho0_meam[i] = rozero[i];
ibar_meam[i] = ibar[i];
switch(this->lattce_meam[i][i]) {
switch (lattce_meam[i][i]) {
case FCC:
this->re_meam[i][i] = tmplat[i] / sqrt(2.0);
re_meam[i][i] = tmplat[i] / sqrt(2.0);
break;
case BCC:
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
break;
case HCP:
case DIM:
@ -88,11 +85,11 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*
case ZIG:
case TRI:
case SC:
this->re_meam[i][i] = tmplat[i];
re_meam[i][i] = tmplat[i];
break;
case DIA:
case DIA3:
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
break;
case B1:
case B2:
@ -100,31 +97,30 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*
case L12:
// do nothing
break;
default:
;
// error
default:;
// error
}
}
// Set some defaults
this->rc_meam = 4.0;
this->delr_meam = 0.1;
setall2d(this->attrac_meam, 0.0);
setall2d(this->repuls_meam, 0.0);
setall3d(this->Cmax_meam, 2.8);
setall3d(this->Cmin_meam, 2.0);
setall2d(this->ebound_meam, (2.8*2.8) / (4.0 * (2.8 - 1.0)));
setall2d(this->delta_meam, 0.0);
setall2d(this->nn2_meam, 0);
setall2d(this->zbl_meam, 1);
this->gsmooth_factor = 99.0;
this->augt1 = 1;
this->ialloy = 0;
this->mix_ref_t = 0;
this->emb_lin_neg = 0;
this->bkgd_dyn = 0;
this->erose_form = 0;
rc_meam = 4.0;
delr_meam = 0.1;
setall2d(attrac_meam, 0.0);
setall2d(repuls_meam, 0.0);
setall3d(Cmax_meam, 2.8);
setall3d(Cmin_meam, 2.0);
setall2d(ebound_meam, (2.8 * 2.8) / (4.0 * (2.8 - 1.0)));
setall2d(delta_meam, 0.0);
setall2d(nn2_meam, 0);
setall2d(zbl_meam, 1);
gsmooth_factor = 99.0;
augt1 = 1;
ialloy = 0;
mix_ref_t = 0;
emb_lin_neg = 0;
bkgd_dyn = 0;
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
setall2d(stheta_meam, 1.0); // stheta = sin(theta/2*pi/180) where theta is 180, so 1.0
setall2d(ctheta_meam, 0.0); // stheta = cos(theta/2*pi/180) where theta is 180, so 0
}

96
src/MEAM/meam_setup_param.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -23,8 +22,8 @@ using namespace LAMMPS_NS;
using MathConst::MY_PI;
// do a sanity check on index parameters
void
MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr)
void MEAM::meam_checkindex(int num, int lim, int nidx, int *idx /*idx(3)*/, int *ierr)
{
//: idx[0..2]
*ierr = 0;
@ -74,8 +73,8 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr
// 2 = not enough indices given
// 3 = an element index is out of range
void
MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag)
void MEAM::meam_setup_param(int which, double value, int nindex, int *index /*index(3)*/,
int *errorflag)
{
//: index[0..2]
int i1, i2;
@ -86,161 +85,148 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3
// 0 = Ec_meam
case 0:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Ec_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
Ec_meam[index[0]][index[1]] = value;
break;
// 1 = alpha_meam
case 1:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->alpha_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
alpha_meam[index[0]][index[1]] = value;
break;
// 2 = rho0_meam
case 2:
meam_checkindex(1, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->rho0_meam[index[0]] = value;
if (*errorflag != 0) return;
rho0_meam[index[0]] = value;
break;
// 3 = delta_meam
case 3:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->delta_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
delta_meam[index[0]][index[1]] = value;
break;
// 4 = lattce_meam
case 4:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
vlat = (lattice_t)value;
if (*errorflag != 0) return;
vlat = (lattice_t) value;
this->lattce_meam[index[0]][index[1]] = vlat;
lattce_meam[index[0]][index[1]] = vlat;
break;
// 5 = attrac_meam
case 5:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->attrac_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
attrac_meam[index[0]][index[1]] = value;
break;
// 6 = repuls_meam
case 6:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->repuls_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
repuls_meam[index[0]][index[1]] = value;
break;
// 7 = nn2_meam
case 7:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
if (*errorflag != 0) return;
i1 = std::min(index[0], index[1]);
i2 = std::max(index[0], index[1]);
this->nn2_meam[i1][i2] = (int)value;
nn2_meam[i1][i2] = (int) value;
break;
// 8 = Cmin_meam
case 8:
meam_checkindex(3, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Cmin_meam[index[0]][index[1]][index[2]] = value;
if (*errorflag != 0) return;
Cmin_meam[index[0]][index[1]][index[2]] = value;
break;
// 9 = Cmax_meam
case 9:
meam_checkindex(3, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Cmax_meam[index[0]][index[1]][index[2]] = value;
if (*errorflag != 0) return;
Cmax_meam[index[0]][index[1]][index[2]] = value;
break;
// 10 = rc_meam
case 10:
this->rc_meam = value;
rc_meam = value;
break;
// 11 = delr_meam
case 11:
this->delr_meam = value;
delr_meam = value;
break;
// 12 = augt1
case 12:
this->augt1 = (int)value;
augt1 = (int) value;
break;
// 13 = gsmooth
case 13:
this->gsmooth_factor = value;
gsmooth_factor = value;
break;
// 14 = re_meam
case 14:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->re_meam[index[0]][index[1]] = value;
if (*errorflag != 0) return;
re_meam[index[0]][index[1]] = value;
break;
// 15 = ialloy
case 15:
this->ialloy = (int)value;
ialloy = (int) value;
break;
// 16 = mixture_ref_t
case 16:
this->mix_ref_t = (int)value;
mix_ref_t = (int) value;
break;
// 17 = erose_form
case 17:
this->erose_form = (int)value;
erose_form = (int) value;
break;
// 18 = zbl_meam
case 18:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
if (*errorflag != 0) return;
i1 = std::min(index[0], index[1]);
i2 = std::max(index[0], index[1]);
this->zbl_meam[i1][i2] = (int)value;
zbl_meam[i1][i2] = (int) value;
break;
// 19 = emb_lin_neg
case 19:
this->emb_lin_neg = (int)value;
emb_lin_neg = (int) value;
break;
// 20 = bkgd_dyn
case 20:
this->bkgd_dyn = (int)value;
bkgd_dyn = (int) value;
break;
// 21 = theta
// see alloyparams(void) in meam_setup_done.cpp
case 21:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
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);
stheta_meam[i1][i2] = sin(value / 2 * MY_PI / 180.0);
ctheta_meam[i1][i2] = cos(value / 2 * MY_PI / 180.0);
break;
default:

211
src/MEAM/pair_meam.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -39,11 +38,10 @@ using namespace LAMMPS_NS;
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", "theta"};
"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",
"theta"};
/* ---------------------------------------------------------------------- */
@ -75,8 +73,7 @@ PairMEAM::~PairMEAM()
{
if (copymode) return;
if (meam_inst)
delete meam_inst;
if (meam_inst) delete meam_inst;
if (allocated) {
memory->destroy(setflag);
@ -89,10 +86,10 @@ PairMEAM::~PairMEAM()
void PairMEAM::compute(int eflag, int vflag)
{
int i,ii,n,inum_half,errorflag;
int *ilist_half,*numneigh_half,**firstneigh_half;
int *numneigh_full,**firstneigh_full;
ev_init(eflag,vflag);
int i, ii, n, inum_half, errorflag;
int *ilist_half, *numneigh_half, **firstneigh_half;
int *numneigh_full, **firstneigh_full;
ev_init(eflag, vflag);
// neighbor list info
@ -107,8 +104,8 @@ void PairMEAM::compute(int eflag, int vflag)
// necessary before doing neigh_f2c and neigh_c2f conversions each step
if (neighbor->ago == 0) {
neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
neigh_strip(inum_half, ilist_half, numneigh_half, firstneigh_half);
neigh_strip(inum_half, ilist_half, numneigh_full, firstneigh_full);
}
// check size of scrfcn based on half neighbor list
@ -133,17 +130,14 @@ void PairMEAM::compute(int eflag, int vflag)
errorflag = 0;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
meam_inst->meam_dens_init(i,ntype,type,map,x,
numneigh_half[i],firstneigh_half[i],
numneigh_full[i],firstneigh_full[i],
offset);
meam_inst->meam_dens_init(i, ntype, type, map, x, numneigh_half[i], firstneigh_half[i],
numneigh_full[i], firstneigh_full[i], offset);
offset += numneigh_half[i];
}
comm->reverse_comm(this);
meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom,
&eng_vdwl,eatom,ntype,type,map,scale,errorflag);
if (errorflag)
error->one(FLERR,"MEAM library error {}",errorflag);
meam_inst->meam_dens_final(nlocal, eflag_either, eflag_global, eflag_atom, &eng_vdwl, eatom,
ntype, type, map, scale, errorflag);
if (errorflag) error->one(FLERR, "MEAM library error {}", errorflag);
comm->forward_comm(this);
@ -156,11 +150,9 @@ void PairMEAM::compute(int eflag, int vflag)
if (vflag_atom) vptr = vatom;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
meam_inst->meam_force(i,eflag_global,eflag_atom,vflag_global,
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,virial);
meam_inst->meam_force(i, eflag_global, eflag_atom, vflag_global, 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, virial);
offset += numneigh_half[i];
}
if (vflag_fdotr) virial_fdotr_compute();
@ -171,13 +163,13 @@ void PairMEAM::compute(int eflag, int vflag)
void PairMEAM::allocate()
{
allocated = 1;
int n = atom->ntypes;
int np1 = atom->ntypes + 1;
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");
memory->create(setflag, np1, np1, "pair:setflag");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(scale, np1, np1, "pair:scale");
map = new int[n+1];
map = new int[np1];
}
/* ----------------------------------------------------------------------
@ -186,13 +178,13 @@ void PairMEAM::allocate()
void PairMEAM::settings(int narg, char ** /*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style {} command", myname);
if (narg != 0) error->all(FLERR, "Illegal pair_style {} command", myname);
// set comm size needed by this Pair
if (msmeamflag) {
comm_forward = 38+23; // plus 23 for msmeam
comm_reverse = 30+23; // plus 23 for msmeam
comm_forward = 38 + 23; // plus 23 for msmeam
comm_reverse = 30 + 23; // plus 23 for msmeam
} else {
comm_forward = 38;
comm_reverse = 30;
@ -205,39 +197,38 @@ void PairMEAM::settings(int narg, char ** /*arg*/)
void PairMEAM::coeff(int narg, char **arg)
{
int m,n;
int m, n;
if (!allocated) allocate();
if (narg < 6) error->all(FLERR,"Incorrect args for pair style {} coefficients", myname);
if (narg < 6) error->all(FLERR, "Incorrect args for pair style {} coefficients", myname);
// check for presence of first meam file
std::string lib_file = utils::get_potential_file_path(arg[2]);
if (lib_file.empty())
error->all(FLERR,"Cannot open MEAM library file {}",lib_file);
if (lib_file.empty()) error->all(FLERR, "Cannot open MEAM library file {}", lib_file);
// find meam parameter file in arguments:
// first word that is a file or "NULL" after the MEAM library file
// we need to extract at least one element, so start from index 4
int paridx=-1;
int paridx = -1;
std::string par_file;
for (int i = 4; i < narg; ++i) {
if (strcmp(arg[i],"NULL") == 0) {
if (strcmp(arg[i], "NULL") == 0) {
par_file = "NULL";
paridx = i;
break;
}
par_file = utils::get_potential_file_path(arg[i]);
if (!par_file.empty()) {
paridx=i;
paridx = i;
break;
}
}
if (paridx < 0) error->all(FLERR,"No MEAM parameter file in pair coefficients");
if (paridx < 0) error->all(FLERR, "No MEAM parameter file in pair coefficients");
if ((narg - paridx - 1) != atom->ntypes)
error->all(FLERR,"Incorrect args for pair style {} coefficients", myname);
error->all(FLERR, "Incorrect args for pair style {} coefficients", myname);
// MEAM element names between 2 filenames
// nlibelements = # of MEAM elements
@ -249,17 +240,21 @@ void PairMEAM::coeff(int narg, char **arg)
}
nlibelements = paridx - 3;
if (nlibelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
if (nlibelements < 1) error->all(FLERR, "Incorrect args for pair coefficients");
if (nlibelements > maxelt)
error->all(FLERR,"Too many elements extracted from MEAM library (current limit: {}). "
"Increase 'maxelt' in meam.h and recompile.", maxelt);
error->all(FLERR,
"Too many elements extracted from MEAM library (current limit: {}). "
"Increase 'maxelt' in meam.h and recompile.",
maxelt);
for (int i = 0; i < nlibelements; i++) {
if (std::any_of(libelements.begin(), libelements.end(),
[&](const std::string &elem) { return elem == arg[i+3]; }))
error->all(FLERR,"Must not extract the same element ({}) from MEAM library twice. ", arg[i+3]);
if (std::any_of(libelements.begin(), libelements.end(), [&](const std::string &elem) {
return elem == arg[i + 3];
}))
error->all(FLERR, "Must not extract the same element ({}) from MEAM library twice. ",
arg[i + 3]);
libelements.emplace_back(arg[i+3]);
libelements.emplace_back(arg[i + 3]);
mass.push_back(0.0);
}
@ -267,28 +262,30 @@ void PairMEAM::coeff(int narg, char **arg)
// pass all parameters to MEAM package
// tell MEAM package that setup is done
read_files(lib_file,par_file);
read_files(lib_file, par_file);
meam_inst->meam_setup_done(&cutmax);
// read args that map atom types to MEAM elements
// map[i] = which element the Ith atom type is, -1 if not mapped
for (int i = 4 + nlibelements; i < narg; i++) {
m = i - (4+nlibelements) + 1;
m = i - (4 + nlibelements) + 1;
int j;
for (j = 0; j < nlibelements; j++)
if (libelements[j] == arg[i]) break;
if (j < nlibelements) map[m] = j;
else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
else error->all(FLERR,"Incorrect args for pair style {} coefficients", myname);
if (j < nlibelements)
map[m] = j;
else if (strcmp(arg[i], "NULL") == 0)
map[m] = -1;
else
error->all(FLERR, "Incorrect args for pair style {} coefficients", myname);
}
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
for (int j = i; j <= n; j++) setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass for i,i in atom class
@ -298,14 +295,14 @@ void PairMEAM::coeff(int narg, char **arg)
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]]);
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 style {} coefficients", myname);
if (count == 0) error->all(FLERR, "Incorrect args for pair style {} coefficients", myname);
}
/* ----------------------------------------------------------------------
@ -314,8 +311,7 @@ void PairMEAM::coeff(int narg, char **arg)
void PairMEAM::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style {} requires newton pair on", myname);
if (force->newton_pair == 0) error->all(FLERR, "Pair style {} requires newton pair on", myname);
// need a full and a half neighbor list
@ -330,8 +326,10 @@ void PairMEAM::init_style()
void PairMEAM::init_list(int id, NeighList *ptr)
{
if (id == 1) listfull = ptr;
else if (id == 2) listhalf = ptr;
if (id == 1)
listfull = ptr;
else if (id == 2)
listhalf = ptr;
}
/* ----------------------------------------------------------------------
@ -347,8 +345,7 @@ double PairMEAM::init_one(int i, int j)
/* ---------------------------------------------------------------------- */
void PairMEAM::read_files(const std::string &globalfile,
const std::string &userfile)
void PairMEAM::read_files(const std::string &globalfile, const std::string &userfile)
{
read_global_meam_file(globalfile);
read_user_meam_file(userfile);
@ -394,7 +391,7 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
if (comm->me == 0) {
PotentialFileReader reader(lmp, globalfile, "MEAM", " library");
char * line;
char *line;
constexpr int params_per_line = 19;
int nset = 0;
@ -425,7 +422,7 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
std::string lattice_type = values.next_string();
if (!MEAM::str_to_lat(lattice_type, true, lat[index]))
error->one(FLERR,"Unrecognized lattice type in MEAM library file: {}", lattice_type);
error->one(FLERR, "Unrecognized lattice type in MEAM library file: {}", lattice_type);
// store parameters
@ -458,11 +455,11 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
ibar[index] = values.next_int();
if (!isone(t0[index]))
error->one(FLERR,"Unsupported parameter in MEAM library file: t0 != 1");
error->one(FLERR, "Unsupported parameter in MEAM library file: t0 != 1");
// z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
if (z[index] != MEAM::get_Zij(lat[index]))
error->one(FLERR,"Mismatched parameter in MEAM library file: z != lat");
error->one(FLERR, "Mismatched parameter in MEAM library file: z != lat");
nset++;
} catch (TokenizerException &e) {
@ -479,7 +476,7 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
msg += " ";
msg += libelements[i];
}
error->one(FLERR,msg);
error->one(FLERR, msg);
}
}
@ -514,16 +511,16 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
if (msmeamflag) {
meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(),
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
t2.data(), t3.data(), rozero.data(), ibar.data(), b1m.data(),
b2m.data(), b3m.data(), t1m.data(), t2m.data(), t3m.data());
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
t2.data(), t3.data(), rozero.data(), ibar.data(), b1m.data(),
b2m.data(), b3m.data(), t1m.data(), t2m.data(), t3m.data());
} else {
meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(),
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
t2.data(), t3.data(), rozero.data(), ibar.data(), nullptr,
nullptr, nullptr, nullptr, nullptr, nullptr);
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
t2.data(), t3.data(), rozero.data(), ibar.data(), nullptr, nullptr,
nullptr, nullptr, nullptr, nullptr);
}
// set element masses
@ -543,14 +540,12 @@ void PairMEAM::read_user_meam_file(const std::string &userfile)
std::shared_ptr<PotentialFileReader> reader;
if (comm->me == 0) {
reader = std::make_shared<PotentialFileReader>(lmp, userfile, "MEAM");
}
if (comm->me == 0) { reader = std::make_shared<PotentialFileReader>(lmp, userfile, "MEAM"); }
// read settings
// pass them one at a time to MEAM package
// match strings to list of corresponding ints
char * line = nullptr;
char *line = nullptr;
char buffer[MAXLINE];
while (true) {
@ -562,14 +557,15 @@ void PairMEAM::read_user_meam_file(const std::string &userfile)
line = reader->next_line();
if (line == nullptr) {
nline = -1;
} else nline = strlen(line) + 1;
} else
nline = strlen(line) + 1;
} else {
line = buffer;
}
MPI_Bcast(&nline,1,MPI_INT,0,world);
if (nline<0) break;
MPI_Bcast(line,nline,MPI_CHAR,0,world);
MPI_Bcast(&nline, 1, MPI_INT, 0, world);
if (nline < 0) break;
MPI_Bcast(line, nline, MPI_CHAR, 0, world);
ValueTokenizer values(line, "=(), '\t\n\r\f");
int nparams = values.count();
@ -578,7 +574,7 @@ void PairMEAM::read_user_meam_file(const std::string &userfile)
for (which = 0; which < nkeywords; which++)
if (keyword == keywords[which]) break;
if (which == nkeywords)
error->all(FLERR,"Keyword {} in MEAM parameter file not recognized", keyword);
error->all(FLERR, "Keyword {} in MEAM parameter file not recognized", keyword);
nindex = nparams - 2;
for (int i = 0; i < nindex; i++) index[i] = values.next_int() - 1;
@ -590,30 +586,27 @@ void PairMEAM::read_user_meam_file(const std::string &userfile)
if (!MEAM::str_to_lat(lattice_type, false, latt))
error->all(FLERR, "Unrecognized lattice type in MEAM parameter file: {}", lattice_type);
value = latt;
}
else value = values.next_double();
} else
value = values.next_double();
// pass single setting to MEAM package
int errorflag = 0;
meam_inst->meam_setup_param(which,value,nindex,index,&errorflag);
meam_inst->meam_setup_param(which, value, nindex, index, &errorflag);
if (errorflag) {
const char *descr[] = { "has an unknown error",
"is out of range (please report a bug)",
"expected more indices",
"has out of range element index"};
const char *descr[] = {"has an unknown error", "is out of range (please report a bug)",
"expected more indices", "has out of range element index"};
if ((errorflag < 0) || (errorflag > 3)) errorflag = 0;
error->all(FLERR,"Error in MEAM parameter file: keyword {} {}", keyword, descr[errorflag]);
error->all(FLERR, "Error in MEAM parameter file: keyword {} {}", keyword, descr[errorflag]);
}
}
}
/* ---------------------------------------------------------------------- */
int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int PairMEAM::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
int i, j, k, m;
m = 0;
for (i = 0; i < n; i++) {
@ -663,7 +656,6 @@ int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
buf[m++] = meam_inst->arho3mb[j][1];
buf[m++] = meam_inst->arho3mb[j][2];
}
}
return m;
@ -673,7 +665,7 @@ int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
void PairMEAM::unpack_forward_comm(int n, int first, double *buf)
{
int i,k,m,last;
int i, k, m, last;
m = 0;
last = first + n;
@ -730,7 +722,7 @@ void PairMEAM::unpack_forward_comm(int n, int first, double *buf)
int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,k,m,last;
int i, k, m, last;
m = 0;
last = first + n;
@ -781,7 +773,7 @@ int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,k,m;
int i, j, k, m;
m = 0;
for (i = 0; i < n; i++) {
@ -824,8 +816,6 @@ void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
meam_inst->arho3mb[j][2] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
@ -835,8 +825,8 @@ void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
double PairMEAM::memory_usage()
{
double bytes = 11 * meam_inst->nmax * sizeof(double);
bytes += (double)(3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double);
bytes += (double)3 * meam_inst->maxneigh * sizeof(double);
bytes += (double) (3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double);
bytes += (double) 3 * meam_inst->maxneigh * sizeof(double);
return bytes;
}
@ -847,10 +837,9 @@ double PairMEAM::memory_usage()
done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
------------------------------------------------------------------------- */
void PairMEAM::neigh_strip(int inum, int *ilist,
int *numneigh, int **firstneigh)
void PairMEAM::neigh_strip(int inum, int *ilist, int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int i, j, ii, jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
@ -866,6 +855,6 @@ void PairMEAM::neigh_strip(int inum, int *ilist,
void *PairMEAM::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"scale") == 0) return (void *) scale;
if (strcmp(str, "scale") == 0) return (void *) scale;
return nullptr;
}

4
src/REAXFF/fix_reaxff_species.cpp
View File

@ -61,8 +61,8 @@ static const char cite_reaxff_species_delete[] =
FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Name(nullptr), MolName(nullptr), NMol(nullptr), nd(nullptr),
MolType(nullptr), molmap(nullptr), mark(nullptr), Mol2Spec(nullptr), clusterID(nullptr),
x0(nullptr), BOCut(nullptr), fp(nullptr), pos(nullptr), fdel(nullptr), ele(nullptr),
eletype(nullptr), filepos(nullptr), filedel(nullptr)
x0(nullptr), BOCut(nullptr), fp(nullptr), pos(nullptr), fdel(nullptr), delete_Tcount(nullptr),
ele(nullptr), eletype(nullptr), filepos(nullptr), filedel(nullptr)
{
if (narg < 7) utils::missing_cmd_args(FLERR, "fix reaxff/species", error);

2
src/REPLICA/fix_pimd.cpp
View File

@ -74,6 +74,7 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nhc_temp = 298.15;
nhc_nchain = 2;
sp = 1.0;
np = universe->nworlds;
for (int i = 3; i < narg - 1; i += 2) {
if (strcmp(arg[i], "method") == 0) {
@ -190,7 +191,6 @@ void FixPIMD::init()
// prepare the constants
np = universe->nworlds;
inverse_np = 1.0 / np;
/* The first solution for the force constant, using SI units

7
src/utils.cpp
View File

@ -1360,11 +1360,8 @@ bool utils::is_double(const std::string &str)
{
if (str.empty()) return false;
if (strmatch(str, "^[+-]?\\d+\\.?\\d*$") || strmatch(str, "^[+-]?\\d+\\.?\\d*[eE][+-]?\\d+$") ||
strmatch(str, "^[+-]?\\d*\\.?\\d+$") || strmatch(str, "^[+-]?\\d*\\.?\\d+[eE][+-]?\\d+$"))
return true;
else
return false;
return strmatch(str, "^[+-]?\\d+\\.?\\d*$") || strmatch(str, "^[+-]?\\d+\\.?\\d*[eE][+-]?\\d+$") ||
strmatch(str, "^[+-]?\\d*\\.?\\d+$") || strmatch(str, "^[+-]?\\d*\\.?\\d+[eE][+-]?\\d+$");
}
/* ----------------------------------------------------------------------

14
src/variable.cpp
View File

@ -1907,12 +1907,14 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
char *var = retrieve(word+2);
if (var == nullptr)
print_var_error(FLERR,"Invalid variable evaluation in variable formula",ivar);
if (tree) {
auto newtree = new Tree();
newtree->type = VALUE;
newtree->value = atof(var);
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = atof(var);
if (utils::is_double(var)) {
if (tree) {
auto newtree = new Tree();
newtree->type = VALUE;
newtree->value = atof(var);
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = atof(var);
} else print_var_error(FLERR,"Non-numeric variable value in variable formula",ivar);
// v_name = per-atom vector from atom-style variable
// evaluate the atom-style variable as newtree

2
unittest/force-styles/tests/mol-pair-lepton.yaml
View File

@ -2,7 +2,7 @@
lammps_version: 22 Dec 2022
date_generated: Thu Dec 22 09:57:30 2022
epsilon: 5e-14
skip_tests:
skip_tests: intel
prerequisites: ! |
atom full
pair lepton

2
unittest/force-styles/tests/mol-pair-lepton_coul.yaml
View File

@ -2,7 +2,7 @@
lammps_version: 22 Dec 2022
date_generated: Fri Jan 6 13:07:07 2023
epsilon: 5e-13
skip_tests:
skip_tests: intel
prerequisites: ! |
atom full
pair lepton/coul

4
unittest/force-styles/tests/mol-pair-lepton_coul_long.yaml
View File

@ -1,8 +1,8 @@
---
lammps_version: 22 Dec 2022
date_generated: Fri Jan 6 11:43:19 2023
epsilon: 5e-14
skip_tests:
epsilon: 1e-13
skip_tests: intel
prerequisites: ! |
atom full
pair lepton/coul

2
unittest/force-styles/tests/mol-pair-lepton_zbl.yaml
View File

@ -2,7 +2,7 @@
lammps_version: 22 Dec 2022
date_generated: Sun Jan 8 01:13:49 2023
epsilon: 5e-10
skip_tests:
skip_tests: intel
prerequisites: ! |
atom full
pair lepton

2
unittest/testing/core.h
View File

@ -54,7 +54,7 @@ extern bool verbose;
class LAMMPSTest : public ::testing::Test {
public:
void command(const std::string &line) { lmp->input->one(line.c_str()); }
void command(const std::string &line) { lmp->input->one(line); }
void BEGIN_HIDE_OUTPUT()
{