diff --git a/cmake/Modules/ExternalCMakeProject.cmake b/cmake/Modules/ExternalCMakeProject.cmake index 855ce254c9..75c33ab99e 100644 --- a/cmake/Modules/ExternalCMakeProject.cmake +++ b/cmake/Modules/ExternalCMakeProject.cmake @@ -9,8 +9,22 @@ function(ExternalCMakeProject target url hash basedir cmakedir cmakefile) get_filename_component(archive ${url} NAME) file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/_deps/src) - message(STATUS "Downloading ${url}") - file(DOWNLOAD ${url} ${CMAKE_BINARY_DIR}/_deps/${archive} EXPECTED_HASH MD5=${hash} SHOW_PROGRESS) + if(EXISTS ${CMAKE_BINARY_DIR}/_deps/${archive}) + file(MD5 ${CMAKE_BINARY_DIR}/_deps/${archive} DL_MD5) + endif() + if(NOT "${DL_MD5}" STREQUAL "${hash}") + message(STATUS "Downloading ${url}") + file(DOWNLOAD ${url} ${CMAKE_BINARY_DIR}/_deps/${archive} STATUS DL_STATUS SHOW_PROGRESS) + file(MD5 ${CMAKE_BINARY_DIR}/_deps/${archive} DL_MD5) + if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${hash}")) + set(${target}_URL ${url}) + GetFallbackURL(${target}_URL fallback) + message(WARNING "Download from primary URL ${url} failed\nTrying fallback URL ${fallback}") + file(DOWNLOAD ${fallback} ${CMAKE_BINARY_DIR}/_deps/${archive} EXPECTED_HASH MD5=${hash} SHOW_PROGRESS) + endif() + else() + message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/_deps/${archive}") + endif() message(STATUS "Unpacking and configuring ${archive}") execute_process(COMMAND ${CMAKE_COMMAND} -E tar xzf ${CMAKE_BINARY_DIR}/_deps/${archive} WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/_deps/src) diff --git a/cmake/Modules/LAMMPSUtils.cmake b/cmake/Modules/LAMMPSUtils.cmake index 9b42dafc44..bb69104308 100644 --- a/cmake/Modules/LAMMPSUtils.cmake +++ b/cmake/Modules/LAMMPSUtils.cmake @@ -158,3 +158,14 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (EXISTS /etc/os-release)) set(CMAKE_LINUX_DISTRO ${distro}) set(CMAKE_DISTRO_VERSION ${disversion}) endif() + +function(GetFallbackURL input output) + string(REPLACE "_URL" "" _tmp ${input}) + string(TOLOWER ${_tmp} libname) + string(REGEX REPLACE "^https://.*/([^/]+gz)" "${LAMMPS_THIRDPARTY_URL}/${libname}-\\1" newurl "${${input}}") + if ("${newurl}" STREQUAL "${${input}}") + set(${output} "" PARENT_SCOPE) + else() + set(${output} ${newurl} PARENT_SCOPE) + endif() +endfunction(GetFallbackURL) diff --git a/cmake/Modules/Packages/GPU.cmake b/cmake/Modules/Packages/GPU.cmake index 24d9538206..cf321eee9a 100644 --- a/cmake/Modules/Packages/GPU.cmake +++ b/cmake/Modules/Packages/GPU.cmake @@ -428,15 +428,17 @@ elseif(GPU_API STREQUAL "HIP") if(DOWNLOAD_CUB) message(STATUS "CUB download requested") - set(CUB_URL "https://github.com/NVlabs/cub/archive/1.12.0.tar.gz" CACHE STRING "URL for CUB tarball") + # TODO: test update to current version 1.17.2 + set(CUB_URL "https://github.com/nvidia/cub/archive/1.12.0.tar.gz" CACHE STRING "URL for CUB tarball") set(CUB_MD5 "1cf595beacafff104700921bac8519f3" CACHE STRING "MD5 checksum of CUB tarball") mark_as_advanced(CUB_URL) mark_as_advanced(CUB_MD5) + GetFallbackURL(CUB_URL CUB_FALLBACK) include(ExternalProject) ExternalProject_Add(CUB - URL ${CUB_URL} + URL ${CUB_URL} ${CUB_FALLBACK} URL_MD5 ${CUB_MD5} PREFIX "${CMAKE_CURRENT_BINARY_DIR}" CONFIGURE_COMMAND "" diff --git a/cmake/Modules/Packages/KOKKOS.cmake b/cmake/Modules/Packages/KOKKOS.cmake index de64df7268..b4eb227d61 100644 --- a/cmake/Modules/Packages/KOKKOS.cmake +++ b/cmake/Modules/Packages/KOKKOS.cmake @@ -53,8 +53,10 @@ if(DOWNLOAD_KOKKOS) set(KOKKOS_MD5 "f140e02b826223b1045207d9bc10d404" CACHE STRING "MD5 checksum of KOKKOS tarball") mark_as_advanced(KOKKOS_URL) mark_as_advanced(KOKKOS_MD5) + GetFallbackURL(KOKKOS_URL KOKKOS_FALLBACK) + ExternalProject_Add(kokkos_build - URL ${KOKKOS_URL} + URL ${KOKKOS_URL} ${KOKKOS_FALLBACK} URL_MD5 ${KOKKOS_MD5} CMAKE_ARGS ${KOKKOS_LIB_BUILD_ARGS} BUILD_BYPRODUCTS /lib/libkokkoscore.a /lib/libkokkoscontainers.a diff --git a/cmake/Modules/Packages/LATTE.cmake b/cmake/Modules/Packages/LATTE.cmake index d7793fa257..bedd7a64fc 100644 --- a/cmake/Modules/Packages/LATTE.cmake +++ b/cmake/Modules/Packages/LATTE.cmake @@ -19,6 +19,7 @@ if(DOWNLOAD_LATTE) set(LATTE_MD5 "820e73a457ced178c08c71389a385de7" CACHE STRING "MD5 checksum of LATTE tarball") mark_as_advanced(LATTE_URL) mark_as_advanced(LATTE_MD5) + GetFallbackURL(LATTE_URL LATTE_FALLBACK) # CMake cannot pass BLAS or LAPACK library variable to external project if they are a list list(LENGTH BLAS_LIBRARIES} NUM_BLAS) @@ -30,7 +31,7 @@ if(DOWNLOAD_LATTE) include(ExternalProject) ExternalProject_Add(latte_build - URL ${LATTE_URL} + URL ${LATTE_URL} ${LATTE_FALLBACK} URL_MD5 ${LATTE_MD5} SOURCE_SUBDIR cmake CMAKE_ARGS -DCMAKE_INSTALL_PREFIX= ${CMAKE_REQUEST_PIC} -DCMAKE_INSTALL_LIBDIR=lib diff --git a/cmake/Modules/Packages/MDI.cmake b/cmake/Modules/Packages/MDI.cmake index c1e25347af..67c1b82dab 100644 --- a/cmake/Modules/Packages/MDI.cmake +++ b/cmake/Modules/Packages/MDI.cmake @@ -12,6 +12,7 @@ if(DOWNLOAD_MDI) set(MDI_MD5 "7a222353ae8e03961d5365e6cd48baee" CACHE STRING "MD5 checksum for MDI tarball") mark_as_advanced(MDI_URL) mark_as_advanced(MDI_MD5) + GetFallbackURL(MDI_URL MDI_FALLBACK) enable_language(C) # only ON/OFF are allowed for "mpi" flag when building MDI library @@ -63,7 +64,7 @@ if(DOWNLOAD_MDI) # support cross-compilation and ninja-build include(ExternalProject) ExternalProject_Add(mdi_build - URL ${MDI_URL} + URL ${MDI_URL} ${MDI_FALLBACK} URL_MD5 ${MDI_MD5} PREFIX ${CMAKE_CURRENT_BINARY_DIR}/mdi_build_ext CMAKE_ARGS diff --git a/cmake/Modules/Packages/ML-HDNNP.cmake b/cmake/Modules/Packages/ML-HDNNP.cmake index f71d08f3ab..99017f222d 100644 --- a/cmake/Modules/Packages/ML-HDNNP.cmake +++ b/cmake/Modules/Packages/ML-HDNNP.cmake @@ -6,10 +6,11 @@ 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) # adjust settings from detected compiler to compiler platform in n2p2 library # set compiler specific flag to turn on C++11 syntax (required on macOS and CentOS 7) @@ -72,7 +73,7 @@ if(DOWNLOAD_N2P2) # download compile n2p2 library. much patch MPI calls in LAMMPS interface to accommodate MPI-2 (e.g. for cross-compiling) include(ExternalProject) ExternalProject_Add(n2p2_build - URL ${N2P2_URL} + URL ${N2P2_URL} ${N2P2_FALLBACK} URL_MD5 ${N2P2_MD5} UPDATE_COMMAND "" CONFIGURE_COMMAND "" diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index 50581b0676..4517aa77d4 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -3,9 +3,23 @@ set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2 set(PACELIB_MD5 "4f0b3b5b14456fe9a73b447de3765caa" CACHE STRING "MD5 checksum of PACE evaluator library tarball") mark_as_advanced(PACELIB_URL) mark_as_advanced(PACELIB_MD5) +GetFallbackURL(PACELIB_URL PACELIB_FALLBACK) # download library sources to build folder -file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5}) #SHOW_PROGRESS +if(EXISTS ${CMAKE_BINARY_DIR}/libpace.tar.gz) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) +endif() +if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") + message(STATUS "Downloading ${PACELIB_URL}") + file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) + message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") + file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) + else() + message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/libpace.tar.gz") + endif() +endif() # uncompress downloaded sources execute_process( diff --git a/cmake/Modules/Packages/ML-QUIP.cmake b/cmake/Modules/Packages/ML-QUIP.cmake index 10e0f38b78..a90b77190f 100644 --- a/cmake/Modules/Packages/ML-QUIP.cmake +++ b/cmake/Modules/Packages/ML-QUIP.cmake @@ -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") diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index e0e4ff0dde..f231d148bd 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -59,10 +59,11 @@ if(DOWNLOAD_PLUMED) mark_as_advanced(PLUMED_URL) mark_as_advanced(PLUMED_MD5) + GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) include(ExternalProject) ExternalProject_Add(plumed_build - URL ${PLUMED_URL} + URL ${PLUMED_URL} ${PLUMED_FALLBACK} URL_MD5 ${PLUMED_MD5} BUILD_IN_SOURCE 1 CONFIGURE_COMMAND /configure --prefix= diff --git a/cmake/Modules/Packages/SCAFACOS.cmake b/cmake/Modules/Packages/SCAFACOS.cmake index de611a1edb..9a5580163f 100644 --- a/cmake/Modules/Packages/SCAFACOS.cmake +++ b/cmake/Modules/Packages/SCAFACOS.cmake @@ -18,6 +18,8 @@ if(DOWNLOAD_SCAFACOS) set(SCAFACOS_MD5 "bd46d74e3296bd8a444d731bb10c1738" CACHE STRING "MD5 checksum of SCAFACOS tarball") mark_as_advanced(SCAFACOS_URL) mark_as_advanced(SCAFACOS_MD5) + GetFallbackURL(SCAFACOS_URL SCAFACOS_FALLBACK) + # version 1.0.1 needs a patch to compile and linke cleanly with GCC 10 and later. file(DOWNLOAD ${LAMMPS_THIRDPARTY_URL}/scafacos-1.0.1-fix.diff ${CMAKE_CURRENT_BINARY_DIR}/scafacos-1.0.1.fix.diff @@ -30,7 +32,7 @@ if(DOWNLOAD_SCAFACOS) include(ExternalProject) ExternalProject_Add(scafacos_build - URL ${SCAFACOS_URL} + URL ${SCAFACOS_URL} ${SCAFACOS_FALLBACK} URL_MD5 ${SCAFACOS_MD5} PATCH_COMMAND patch -p1 < ${CMAKE_CURRENT_BINARY_DIR}/scafacos-1.0.1.fix.diff CONFIGURE_COMMAND /configure --prefix= --disable-doc diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index 75e4333abd..c4c33cc40c 100644 --- a/cmake/Modules/Tools.cmake +++ b/cmake/Modules/Tools.cmake @@ -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/) diff --git a/doc/src/pair_gauss.rst b/doc/src/pair_gauss.rst index 4eb45a416b..091bd52d45 100644 --- a/doc/src/pair_gauss.rst +++ b/doc/src/pair_gauss.rst @@ -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 ` -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 ` shift option for the energy of the Gauss-potential portion of the pair diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 2d12ce6845..dc2f1879fd 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -363,6 +363,7 @@ Broglie brownian brownw Broyden +Bruskin Brusselle Bryantsev Btarget @@ -2241,6 +2242,7 @@ msd msi MSI msm +msmeam msmflag msse msst diff --git a/examples/PACKAGES/pimd/para-h2/in.scp b/examples/PACKAGES/pimd/para-h2/in.scp index eac3f066b8..eb8debad15 100644 --- a/examples/PACKAGES/pimd/para-h2/in.scp +++ b/examples/PACKAGES/pimd/para-h2/in.scp @@ -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 diff --git a/examples/PACKAGES/pimd/para-h2/run.sh b/examples/PACKAGES/pimd/para-h2/run.sh index ace20f84ea..5304db5917 100644 --- a/examples/PACKAGES/pimd/para-h2/run.sh +++ b/examples/PACKAGES/pimd/para-h2/run.sh @@ -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 diff --git a/examples/PACKAGES/pimd/prot-hairpin/in.scp b/examples/PACKAGES/pimd/prot-hairpin/in.scp index 220d29dc51..00b5150a8a 100644 --- a/examples/PACKAGES/pimd/prot-hairpin/in.scp +++ b/examples/PACKAGES/pimd/prot-hairpin/in.scp @@ -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 diff --git a/examples/PACKAGES/pimd/prot-hairpin/run.sh b/examples/PACKAGES/pimd/prot-hairpin/run.sh index 1423960325..756cc14461 100644 --- a/examples/PACKAGES/pimd/prot-hairpin/run.sh +++ b/examples/PACKAGES/pimd/prot-hairpin/run.sh @@ -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 diff --git a/lib/gpu/lal_atom.cpp b/lib/gpu/lal_atom.cpp index 72cb59a912..3d1a1cc963 100644 --- a/lib/gpu/lal_atom.cpp +++ b/lib/gpu/lal_atom.cpp @@ -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, diff --git a/lib/hdnnp/.gitignore b/lib/hdnnp/.gitignore new file mode 100644 index 0000000000..2a6742ed7a --- /dev/null +++ b/lib/hdnnp/.gitignore @@ -0,0 +1,3 @@ +/includelink +/liblink +/n2p2-* diff --git a/lib/hdnnp/Install.py b/lib/hdnnp/Install.py index 3b232f400b..c5110a0197 100644 --- a/lib/hdnnp/Install.py +++ b/lib/hdnnp/Install.py @@ -10,14 +10,14 @@ import sys, os, platform, subprocess, shutil from argparse import ArgumentParser sys.path.append('..') -from install_helpers import get_cpus, fullpath, geturl, checkmd5sum +from install_helpers import get_cpus, fullpath, geturl, checkmd5sum, getfallback parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") # 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', \ } @@ -76,14 +77,21 @@ if pathflag: if buildflag: url = "https://github.com/CompPhysVienna/n2p2/archive/v%s.tar.gz" % (version) - filename = "n2p2-%s.tar.gz" %version - print("Downloading n2p2 ...") - geturl(url, filename) + filename = "n2p2-%s.tar.gz" % version + fallback = getfallback('n2p2', url) + print("Downloading n2p2 from", url) + try: + geturl(url, filename) + except: + geturl(fallback, filename) # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: if not checkmd5sum(checksums[version], filename): - sys.exit("Checksum for n2p2 library does not match") + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, filename) + if not checkmd5sum(checksums[version], filename): + sys.exit("Checksum for n2p2 library does not match for fallback, too.") print("Unpacking n2p2 source tarball ...") if os.path.exists("%s/n2p2-%s" % (homepath, version)): diff --git a/lib/install_helpers.py b/lib/install_helpers.py index e1ebd1cde7..59171ee4e0 100644 --- a/lib/install_helpers.py +++ b/lib/install_helpers.py @@ -1,4 +1,7 @@ -import hashlib,os,subprocess +import hashlib +import os +import re +import subprocess # try to auto-detect the maximum number of available CPUs def get_cpus(): @@ -32,31 +35,31 @@ def which(program): return None -def geturl(url,fname): +def geturl(url, fname): success = False if which('curl') != None: - cmd = 'curl -L -o "%s" %s' % (fname,url) + cmd = 'curl -L -o "%s" %s' % (fname, url) try: - subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) + subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True) success = True except subprocess.CalledProcessError as e: print("Calling curl failed with: %s" % e.output.decode('UTF-8')) if not success and which('wget') != None: - cmd = 'wget -O "%s" %s' % (fname,url) + cmd = 'wget -O "%s" %s' % (fname, url) try: - subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) + subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True) success = True except subprocess.CalledProcessError as e: print("Calling wget failed with: %s" % e.output.decode('UTF-8')) if not success: - error("Failed to download source code with 'curl' or 'wget'") + raise Exception("Failed to download source code with 'curl' or 'wget' from " + url) return -def checkmd5sum(md5sum,fname): - with open(fname,'rb') as fh: +def checkmd5sum(md5sum, fname): + with open(fname, 'rb') as fh: m = hashlib.md5() while True: data = fh.read(81920) @@ -66,3 +69,6 @@ def checkmd5sum(md5sum,fname): fh.close() return m.hexdigest() == md5sum +def getfallback(lib, url): + archive = re.sub(r'^https://.*/([^/]+gz)', r'-\1', url) + return 'https://download.lammps.org/thirdparty/' + lib + archive diff --git a/lib/kokkos/core/src/Kokkos_Macros.hpp b/lib/kokkos/core/src/Kokkos_Macros.hpp index 9dbd2de0c8..1ba11481e7 100644 --- a/lib/kokkos/core/src/Kokkos_Macros.hpp +++ b/lib/kokkos/core/src/Kokkos_Macros.hpp @@ -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 diff --git a/lib/latte/Install.py b/lib/latte/Install.py index 2e8f9bee8d..4e3b27859c 100644 --- a/lib/latte/Install.py +++ b/lib/latte/Install.py @@ -10,7 +10,7 @@ import sys, os, subprocess, shutil, tarfile from argparse import ArgumentParser sys.path.append('..') -from install_helpers import fullpath, geturl, checkmd5sum +from install_helpers import fullpath, geturl, checkmd5sum, getfallback parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") @@ -86,17 +86,25 @@ if buildflag: url = "https://github.com/lanl/LATTE/archive/v%s.tar.gz" % version lattepath = fullpath(homepath) lattedir = os.path.join(lattepath, homedir) + fallback = getfallback('latte', url) + filename = 'LATTE.tar.gz' # download and unpack LATTE tarball if buildflag: print("Downloading LATTE ...") - geturl(url, "LATTE.tar.gz") + try: + geturl(url, filename) + except: + geturl(fallback, filename) # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: - if not checkmd5sum(checksums[version], 'LATTE.tar.gz'): - sys.exit("Checksum for LATTE library does not match") + if not checkmd5sum(checksums[version], filename): + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, filename) + if not checkmd5sum(checksums[version], filename): + sys.exit("Checksum for LATTE library does not match for fallback, too.") print("Unpacking LATTE ...") if os.path.exists(lattedir): diff --git a/lib/mdi/Install.py b/lib/mdi/Install.py index 59d218b1f9..77d13d3e8b 100644 --- a/lib/mdi/Install.py +++ b/lib/mdi/Install.py @@ -9,7 +9,7 @@ import sys,os,subprocess import glob sys.path.append('..') -from install_helpers import checkmd5sum +from install_helpers import fullpath, geturl, checkmd5sum, getfallback # help message @@ -51,49 +51,6 @@ def error(str=None): # expand to full path name # process leading '~' or relative path -def fullpath(path): - return os.path.abspath(os.path.expanduser(path)) - -def which(program): - def is_exe(fpath): - return os.path.isfile(fpath) and os.access(fpath, os.X_OK) - - fpath, fname = os.path.split(program) - if fpath: - if is_exe(program): - return program - else: - for path in os.environ["PATH"].split(os.pathsep): - path = path.strip('"') - exe_file = os.path.join(path, program) - if is_exe(exe_file): - return exe_file - - return None - -def geturl(url,fname): - success = False - - if which('curl') != None: - cmd = 'curl -L -o "%s" %s' % (fname,url) - try: - subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) - success = True - except subprocess.CalledProcessError as e: - print("Calling curl failed with: %s" % e.output.decode('UTF-8')) - - if not success and which('wget') != None: - cmd = 'wget -O "%s" %s' % (fname,url) - try: - subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) - success = True - except subprocess.CalledProcessError as e: - print("Calling wget failed with: %s" % e.output.decode('UTF-8')) - - if not success: - error("Failed to download source code with 'curl' or 'wget'") - return - # parse args args = sys.argv[1:] @@ -123,17 +80,24 @@ lib = os.path.basename(cwd) # download and unpack MDI_Library tarball -homepath = "." +homepath = fullpath('.') homedir = "%s/MDI_Library" % homepath print("Downloading MDI_Library ...") -mditar = "%s/v%s.tar.gz" % (homepath,version) -geturl(url, mditar) +mditar = "%s/v%s.tar.gz" % (homepath, version) +fallback = getfallback('mdi', url) +try: + geturl(url, mditar) +except: + geturl(fallback, mditar) # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: if not checkmd5sum(checksums[version], mditar): - sys.exit("Checksum for MDI library does not match") + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, mditar) + if not checkmd5sum(checksums[version], mditar): + sys.exit("Checksum for MDI library does not match") print("Unpacking MDI_Library tarball ...") if os.path.exists("%s/v%s" % (homepath,version)): @@ -199,7 +163,6 @@ makefile_lammps = open(str(dir_path) + "/Makefile.lammps", "a") makefile_lammps.write(str(rpath_option) + "\n") makefile_lammps.close() - shared_files = glob.glob( os.path.join( homepath, "liblink", "lib%s.a" % lib) ) if len(shared_files) > 0: print("Build was successful") diff --git a/lib/mscg/Install.py b/lib/mscg/Install.py index 9d18b7c083..df17005dca 100644 --- a/lib/mscg/Install.py +++ b/lib/mscg/Install.py @@ -10,7 +10,7 @@ import sys, os, subprocess, shutil, tarfile from argparse import ArgumentParser sys.path.append('..') -from install_helpers import fullpath, geturl +from install_helpers import fullpath, geturl, checkmd5sum, getfallback parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") @@ -84,7 +84,19 @@ if pathflag: if buildflag: print("Downloading MS-CG ...") tarname = os.path.join(homepath, tarname) - geturl(url, tarname) + fallback = getfallback('mscg', url) + try: + geturl(url, tarname) + except: + geturl(fallback, tarname) + + # verify downloaded archive integrity via md5 checksum, if known. + if mscgver in checksums: + if not checkmd5sum(checksums[mscgver], tarname): + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, tarname) + if not checkmd5sum(checksums[mscgver], tarname): + sys.exit("Checksum for LATTE library does not match for fallback, too.") print("Unpacking MS-CG tarfile ...") diff --git a/lib/pace/Install.py b/lib/pace/Install.py index c72271aa7a..4f3cf299ac 100644 --- a/lib/pace/Install.py +++ b/lib/pace/Install.py @@ -13,7 +13,7 @@ import sys from argparse import ArgumentParser sys.path.append('..') -from install_helpers import fullpath, geturl, checkmd5sum +from install_helpers import fullpath, geturl, checkmd5sum, getfallback # settings @@ -77,14 +77,21 @@ if buildflag: print("Downloading pace tarball ...") archive_filename = "%s.%s" % (version, archive_extension) download_filename = "%s/%s" % (thisdir, archive_filename) + fallback = getfallback('pacelib', url) print("Downloading from ", url, " to ", download_filename, end=" ") - geturl(url, download_filename) + try: + geturl(url, download_filename) + except: + geturl(fallback, download_filename) print(" done") # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: - if not checkmd5sum(checksums[version], archive_filename): - sys.exit("Checksum for pace library does not match") + if not checkmd5sum(checksums[version], download_filename): + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, download_filename) + if not checkmd5sum(checksums[version], download_filename): + sys.exit("Checksum for pace library does not match for fallback, too.") print("Unpacking pace tarball ...") src_folder = thisdir + "/src" diff --git a/lib/plumed/Install.py b/lib/plumed/Install.py index bb3a158bd4..bb6f27fddb 100644 --- a/lib/plumed/Install.py +++ b/lib/plumed/Install.py @@ -10,7 +10,7 @@ import sys, os, platform, subprocess, shutil from argparse import ArgumentParser sys.path.append('..') -from install_helpers import get_cpus, fullpath, geturl, checkmd5sum +from install_helpers import get_cpus, fullpath, geturl, checkmd5sum, getfallback parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") @@ -86,6 +86,7 @@ buildflag = args.build pathflag = args.path is not None plumedpath = args.path mode = args.mode +version = args.version homepath = fullpath('.') homedir = "%s/plumed2" % (homepath) @@ -101,14 +102,21 @@ if pathflag: if buildflag: url = "https://github.com/plumed/plumed2/releases/download/v%s/plumed-src-%s.tgz" % (version, version) - filename = "plumed-src-%s.tar.gz" %version + filename = "plumed-src-%s.tar.gz" % version + fallback = getfallback('plumed', url) print("Downloading plumed ...") - geturl(url, filename) + try: + geturl(url, filename) + except: + geturl(fallback, filename) # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: if not checkmd5sum(checksums[version], filename): - sys.exit("Checksum for plumed2 library does not match") + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, filename) + if not checkmd5sum(checksums[version], filename): + sys.exit("Checksum for plumed2 library does not match for fallback, too.") print("Unpacking plumed2 source tarball ...") if os.path.exists("%s/plumed-%s" % (homepath, version)): diff --git a/lib/scafacos/Install.py b/lib/scafacos/Install.py index 8044215d47..28bcf43086 100644 --- a/lib/scafacos/Install.py +++ b/lib/scafacos/Install.py @@ -10,15 +10,13 @@ import sys, os, subprocess, shutil, tarfile from argparse import ArgumentParser sys.path.append('..') -from install_helpers import fullpath, geturl, get_cpus, checkmd5sum +from install_helpers import fullpath, geturl, get_cpus, checkmd5sum, getfallback -parser = ArgumentParser(prog='Install.py', - description="LAMMPS library build wrapper script") +parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") # settings version = "1.0.1" -url = "https://github.com/scafacos/scafacos/releases/download/v%s/scafacos-%s.tar.gz" % (version, version) # known checksums for different ScaFaCoS versions. used to validate the download. checksums = { \ @@ -59,6 +57,7 @@ if not args.build and not args.path: buildflag = args.build pathflag = args.path is not None version = args.version +url = "https://github.com/scafacos/scafacos/releases/download/v%s/scafacos-%s.tar.gz" % (version, version) homepath = fullpath(".") scafacospath = os.path.join(homepath, "scafacos-%s" % version) @@ -76,12 +75,20 @@ if pathflag: if buildflag: print("Downloading ScaFaCoS ...") - geturl(url, "%s/scafacos-%s.tar.gz" % (homepath, version)) + filename = "%s/scafacos-%s.tar.gz" % (homepath, version) + fallback = getfallback('scafacos', url) + try: + geturl(url, filename) + except: + geturl(fallback, filename) # verify downloaded archive integrity via md5 checksum, if known. if version in checksums: - if not checkmd5sum(checksums[version], '%s/scafacos-%s.tar.gz' % (homepath, version)): - sys.exit("Checksum for ScaFaCoS library does not match") + if not checkmd5sum(checksums[version], filename): + print("Checksum did not match. Trying fallback URL", fallback) + geturl(fallback, filename) + if not checkmd5sum(checksums[version], filename): + sys.exit("Checksum for ScaFaCoS library does not match for fallback, too.") print("Unpacking ScaFaCoS tarball ...") if os.path.exists(scafacospath): diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp index ca87964e76..b9ae1d65ff 100644 --- a/src/DIELECTRIC/pppm_dielectric.cpp +++ b/src/DIELECTRIC/pppm_dielectric.cpp @@ -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; diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index e8a954616b..2c4de6ada1 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -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; diff --git a/src/ELECTRODE/fix_electrode_conp.cpp b/src/ELECTRODE/fix_electrode_conp.cpp index dbb508a033..beecc18f68 100644 --- a/src/ELECTRODE/fix_electrode_conp.cpp +++ b/src/ELECTRODE/fix_electrode_conp.cpp @@ -917,7 +917,7 @@ std::vector FixElectrodeConp::gather_ngroup(std::vector x_local) std::vector FixElectrodeConp::constraint_correction(std::vector x) { - return constraint_projection(x); + return constraint_projection(std::move(x)); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/atom_vec_hybrid_kokkos.cpp b/src/KOKKOS/atom_vec_hybrid_kokkos.cpp index 4dfbf22fae..fce0b3b337 100644 --- a/src/KOKKOS/atom_vec_hybrid_kokkos.cpp +++ b/src/KOKKOS/atom_vec_hybrid_kokkos.cpp @@ -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); diff --git a/src/KOKKOS/atom_vec_hybrid_kokkos.h b/src/KOKKOS/atom_vec_hybrid_kokkos.h index 2a725577ba..f3aad18937 100644 --- a/src/KOKKOS/atom_vec_hybrid_kokkos.h +++ b/src/KOKKOS/atom_vec_hybrid_kokkos.h @@ -32,7 +32,6 @@ namespace LAMMPS_NS { class AtomVecHybridKokkos : public AtomVecKokkos, public AtomVecHybrid { public: AtomVecHybridKokkos(class LAMMPS *); - ~AtomVecHybridKokkos() override; void grow(int) override; diff --git a/src/KOKKOS/kokkos.cpp b/src/KOKKOS/kokkos.cpp index 9bbfb4157f..33890e4040 100644 --- a/src/KOKKOS/kokkos.cpp +++ b/src/KOKKOS/kokkos.cpp @@ -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) diff --git a/src/KOKKOS/kokkos.h b/src/KOKKOS/kokkos.h index 02189e6c4a..4fa294e497 100644 --- a/src/KOKKOS/kokkos.h +++ b/src/KOKKOS/kokkos.h @@ -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); diff --git a/src/KOKKOS/meam_force_kokkos.h b/src/KOKKOS/meam_force_kokkos.h index 5c4244e99b..d086230fc7 100644 --- a/src/KOKKOS/meam_force_kokkos.h +++ b/src/KOKKOS/meam_force_kokkos.h @@ -446,7 +446,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos::operator()(TagMEAMForcev2D[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; diff --git a/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp b/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp index 6518cccaa8..de6404b691 100644 --- a/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp +++ b/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp @@ -46,13 +46,6 @@ MLIAPDescriptorSO3Kokkos::MLIAPDescriptorSO3Kokkos(LAMMPS *lmp, char /* ---------------------------------------------------------------------- */ -template -MLIAPDescriptorSO3Kokkos::~MLIAPDescriptorSO3Kokkos() -{ -} - -/* ---------------------------------------------------------------------- */ - template void MLIAPDescriptorSO3Kokkos::compute_descriptors(class MLIAPData *data_) { diff --git a/src/KOKKOS/mliap_descriptor_so3_kokkos.h b/src/KOKKOS/mliap_descriptor_so3_kokkos.h index aea4540632..2d15924fe8 100644 --- a/src/KOKKOS/mliap_descriptor_so3_kokkos.h +++ b/src/KOKKOS/mliap_descriptor_so3_kokkos.h @@ -29,7 +29,6 @@ class MLIAPDescriptorSO3Kokkos : public MLIAPDescriptorKokkos { public: MLIAPDescriptorSO3Kokkos(LAMMPS *, char *); - ~MLIAPDescriptorSO3Kokkos() override; void compute_descriptors(class MLIAPData *) override; void compute_forces(class MLIAPData *) override; diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 063eeee7e8..2181c0de97 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -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 diff --git a/src/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index 7dd50a2792..a70f3240db 100644 --- a/src/MC/fix_sgcmc.cpp +++ b/src/MC/fix_sgcmc.cpp @@ -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 localSpeciesCounts(atom->ntypes+1, 0); for (int i = 0; i < atom->nlocal; i++, ++type) { diff --git a/src/MC/fix_widom.cpp b/src/MC/fix_widom.cpp index 6f0afdc6ba..e043eaabd6 100644 --- a/src/MC/fix_widom.cpp +++ b/src/MC/fix_widom.cpp @@ -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 diff --git a/src/MEAM/meam_dens_final.cpp b/src/MEAM/meam_dens_final.cpp index ab0ac8c53f..c5dcc6ba00 100644 --- a/src/MEAM/meam_dens_final.cpp +++ b/src/MEAM/meam_dens_final.cpp @@ -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; diff --git a/src/MEAM/meam_dens_init.cpp b/src/MEAM/meam_dens_init.cpp index 00ad276ad7..894bea1117 100644 --- a/src/MEAM/meam_dens_init.cpp +++ b/src/MEAM/meam_dens_init.cpp @@ -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 } } } - diff --git a/src/MEAM/meam_force.cpp b/src/MEAM/meam_force.cpp index 4bc7380898..23230e0fbc 100644 --- a/src/MEAM/meam_force.cpp +++ b/src/MEAM/meam_force.cpp @@ -15,16 +15,16 @@ #include "math_special.h" -#include #include +#include 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]; diff --git a/src/MEAM/meam_funcs.cpp b/src/MEAM/meam_funcs.cpp index f3a6e99772..13caa1e984 100644 --- a/src/MEAM/meam_funcs.cpp +++ b/src/MEAM/meam_funcs.cpp @@ -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); diff --git a/src/MEAM/meam_impl.cpp b/src/MEAM/meam_impl.cpp index 5290647b18..473b491b01 100644 --- a/src/MEAM/meam_impl.cpp +++ b/src/MEAM/meam_impl.cpp @@ -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); } } diff --git a/src/MEAM/meam_setup_done.cpp b/src/MEAM/meam_setup_done.cpp index de1188349c..4adfd68f19 100644 --- a/src/MEAM/meam_setup_done.cpp +++ b/src/MEAM/meam_setup_done.cpp @@ -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) { - 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 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 ineltypes; 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; } } diff --git a/src/MEAM/meam_setup_global.cpp b/src/MEAM/meam_setup_global.cpp index 5d35242e7c..299fc4da61 100644 --- a/src/MEAM/meam_setup_global.cpp +++ b/src/MEAM/meam_setup_global.cpp @@ -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 -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 -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 } diff --git a/src/MEAM/meam_setup_param.cpp b/src/MEAM/meam_setup_param.cpp index 32634004c2..45c8d0ab78 100644 --- a/src/MEAM/meam_setup_param.cpp +++ b/src/MEAM/meam_setup_param.cpp @@ -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: diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index c4a4cfa1d7..2f095754af 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -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 reader; - if (comm->me == 0) { - reader = std::make_shared(lmp, userfile, "MEAM"); - } + if (comm->me == 0) { reader = std::make_shared(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; } diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index ce04be2cc8..29441cd4b3 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -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); diff --git a/src/REPLICA/fix_pimd.cpp b/src/REPLICA/fix_pimd.cpp index 154f3deecd..3968cd0d4b 100644 --- a/src/REPLICA/fix_pimd.cpp +++ b/src/REPLICA/fix_pimd.cpp @@ -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 diff --git a/src/utils.cpp b/src/utils.cpp index 31a14a352a..bc4eee7385 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -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+$"); } /* ---------------------------------------------------------------------- diff --git a/src/variable.cpp b/src/variable.cpp index 78ce8d8758..4285acfb05 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -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 diff --git a/unittest/force-styles/tests/mol-pair-lepton.yaml b/unittest/force-styles/tests/mol-pair-lepton.yaml index 81d8286588..03117a9aa5 100644 --- a/unittest/force-styles/tests/mol-pair-lepton.yaml +++ b/unittest/force-styles/tests/mol-pair-lepton.yaml @@ -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 diff --git a/unittest/force-styles/tests/mol-pair-lepton_coul.yaml b/unittest/force-styles/tests/mol-pair-lepton_coul.yaml index 06dba1ebb1..2092072eb3 100644 --- a/unittest/force-styles/tests/mol-pair-lepton_coul.yaml +++ b/unittest/force-styles/tests/mol-pair-lepton_coul.yaml @@ -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 diff --git a/unittest/force-styles/tests/mol-pair-lepton_coul_long.yaml b/unittest/force-styles/tests/mol-pair-lepton_coul_long.yaml index e9e6e12ec4..8f909c9492 100644 --- a/unittest/force-styles/tests/mol-pair-lepton_coul_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lepton_coul_long.yaml @@ -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 diff --git a/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml b/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml index 883021c26b..8ffb36eac8 100644 --- a/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml +++ b/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml @@ -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 diff --git a/unittest/testing/core.h b/unittest/testing/core.h index b64affc0c7..272d1d21c5 100644 --- a/unittest/testing/core.h +++ b/unittest/testing/core.h @@ -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() {