diff --git a/cmake/Modules/Packages/KIM.cmake b/cmake/Modules/Packages/KIM.cmake index d74b82b036..42b7bb884d 100644 --- a/cmake/Modules/Packages/KIM.cmake +++ b/cmake/Modules/Packages/KIM.cmake @@ -34,8 +34,8 @@ if(DOWNLOAD_KIM) enable_language(C) enable_language(Fortran) ExternalProject_Add(kim_build - URL https://s3.openkim.org/kim-api/kim-api-2.1.3.txz - URL_MD5 6ee829a1bbba5f8b9874c88c4c4ebff8 + URL https://s3.openkim.org/kim-api/kim-api-2.2.0.txz + URL_MD5 e7f944e1593cffd7444679a660607f6c BINARY_DIR build CMAKE_ARGS ${CMAKE_REQUEST_PIC} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} @@ -53,11 +53,28 @@ if(DOWNLOAD_KIM) add_library(LAMMPS::KIM UNKNOWN IMPORTED) set_target_properties(LAMMPS::KIM PROPERTIES IMPORTED_LOCATION "${INSTALL_DIR}/lib/libkim-api${CMAKE_SHARED_LIBRARY_SUFFIX}" - INTERFACE_INCLUDE_DIRECTORIES "${INSTALL_DIR}/include/kim-api") - target_link_libraries(lammps PRIVATE LAMMPS::KIM) + INTERFACE_INCLUDE_DIRECTORIES "${INSTALL_DIR}/include/kim-api" + ) add_dependencies(LAMMPS::KIM kim_build) + target_link_libraries(lammps PRIVATE LAMMPS::KIM) + # Set rpath so lammps build directory is relocatable + if("${CMAKE_SYSTEM_NAME}" STREQUAL "Darwin") + set(_rpath_prefix "@loader_path") + else() + set(_rpath_prefix "$ORIGIN") + endif() + set_target_properties(lmp PROPERTIES + BUILD_RPATH "${_rpath_prefix}/kim_build-prefix/lib" + ) else() - find_package(PkgConfig REQUIRED) - pkg_check_modules(KIM-API REQUIRED IMPORTED_TARGET libkim-api>=${KIM-API_MIN_VERSION}) - target_link_libraries(lammps PRIVATE PkgConfig::KIM-API) + if(KIM-API_FOUND AND KIM_API_VERSION VERSION_GREATER_EQUAL 2.2.0) + # For kim-api >= 2.2.0 + find_package(KIM-API ${KIM-API_MIN_VERSION} CONFIG REQUIRED) + target_link_libraries(lammps PRIVATE KIM-API::kim-api) + else() + # For kim-api 2.1.3 (consistent with previous version of this file) + find_package(PkgConfig REQUIRED) + pkg_check_modules(KIM-API REQUIRED IMPORTED_TARGET libkim-api>=KIM-API_MIN_VERSION) + target_link_libraries(lammps PRIVATE PkgConfig::KIM-API) + endif() endif() diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index dd56938773..a826228219 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -35,6 +35,7 @@ OPT. * :doc:`class2 (ko) ` * :doc:`fene (iko) ` * :doc:`fene/expand (o) ` + * :doc:`gaussian ` * :doc:`gromos (o) ` * :doc:`harmonic (iko) ` * :doc:`harmonic/shift (o) ` @@ -84,6 +85,7 @@ OPT. * :doc:`dipole (o) ` * :doc:`fourier (o) ` * :doc:`fourier/simple (o) ` + * :doc:`gaussian ` - multicentered Gaussian-based angle potential * :doc:`harmonic (iko) ` * :doc:`mm3 ` * :doc:`quartic (o) ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 227409e45c..d1bc20c9ca 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -62,6 +62,7 @@ OPT. * :doc:`efield ` * :doc:`ehex ` * :doc:`electron/stopping ` + * :doc:`electron/stopping/fit ` * :doc:`enforce2d (k) ` * :doc:`eos/cv ` * :doc:`eos/table ` diff --git a/doc/src/angle_gaussian.rst b/doc/src/angle_gaussian.rst new file mode 100644 index 0000000000..1caae3b20c --- /dev/null +++ b/doc/src/angle_gaussian.rst @@ -0,0 +1,69 @@ +.. index:: angle_style gaussian + +angle_style gaussian command +================================ + +Syntax +"""""" + +.. code-block:: LAMMPS + + angle_style gaussian + +Examples +"""""""" + +.. code-block:: LAMMPS + + angle_style gaussian + angle_coeff 1 300.0 2 0.0128 0.375 80.0 0.0730 0.148 123.0 + +Description +""""""""""" + +The *gaussian* angle style uses the potential: + +.. math:: + + E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-(\theta-\theta_{i})^2}{w_i^2})\right) \right) + +This analytical form is a suitable potential for obtaining +mesoscale effective force fields which can reproduce target atomistic distributions :ref:`(Milano) ` +The following coefficients must be defined for each angle type via the +:doc:`angle_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands: + +* T temperature at which the potential was derived +* :math:`n` (integer >=1) +* :math:`A_1` (-) +* :math:`w_1` (-) +* :math:`\theta_1` (degrees) +* ... +* :math:`A_n` (-) +* :math:`w_n` (-) +* :math:`\theta_n` (degrees) + + +Restrictions +"""""""""""" + +This angle style can only be used if LAMMPS was built with the +USER-MISC package. See the :doc:`Build package ` doc +page for more info. + +Related commands +"""""""""""""""" + +:doc:`angle_coeff ` + +Default +""""""" + +none + +---------- + +.. _Milano1: + +**(Milano)** G. Milano, S. Goudeau, F. Mueller-Plathe, J. Polym. Sci. B Polym. Phys. 43, 871 (2005). diff --git a/doc/src/angle_style.rst b/doc/src/angle_style.rst index 91232a6b6c..540d15fee2 100644 --- a/doc/src/angle_style.rst +++ b/doc/src/angle_style.rst @@ -87,6 +87,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist. * :doc:`dipole ` - angle that controls orientation of a point dipole * :doc:`fourier ` - angle with multiple cosine terms * :doc:`fourier/simple ` - angle with a single cosine term +* :doc:`gaussian ` - multicentered Gaussian-based angle potential * :doc:`harmonic ` - harmonic angle * :doc:`mm3 ` - anharmonic angle * :doc:`quartic ` - angle with cubic and quartic terms diff --git a/doc/src/bond_gaussian.rst b/doc/src/bond_gaussian.rst new file mode 100644 index 0000000000..dc9bea1cce --- /dev/null +++ b/doc/src/bond_gaussian.rst @@ -0,0 +1,70 @@ +.. index:: bond_style gaussian + +bond_style gaussian command +================================ + +Syntax +"""""" + +.. code-block:: LAMMPS + + bond_style gaussian + +Examples +"""""""" + +.. code-block:: LAMMPS + + bond_style gaussian + bond_coeff 1 300.0 2 0.0128 0.375 3.37 0.0730 0.148 3.63 + +Description +""""""""""" + +The *gaussian* bond style uses the potential: + +.. math:: + + E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-(r-r_{i})^2}{w_i^2})\right) \right) + +This analytical form is a suitable potential for obtaining +mesoscale effective force fields which can reproduce target atomistic distributions :ref:`(Milano) ` + +The following coefficients must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands: + +* T temperature at which the potential was derived +* :math:`n` (integer >=1) +* :math:`A_1` (-) +* :math:`w_1` (-) +* :math:`r_1` (length) +* ... +* :math:`A_n` (-) +* :math:`w_n` (-) +* :math:`r_n` (length) + + +Restrictions +"""""""""""" + +This bond style can only be used if LAMMPS was built with the +USER-MISC package. See the :doc:`Build package ` doc +page for more info. + +Related commands +"""""""""""""""" + +:doc:`bond_coeff ` + +Default +""""""" + +none + +---------- + +.. _Milano0: + +**(Milano)** G. Milano, S. Goudeau, F. Mueller-Plathe, J. Polym. Sci. B Polym. Phys. 43, 871 (2005). diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index 64fdeb95c9..e302ea792d 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -87,6 +87,7 @@ accelerated styles exist. * :doc:`class2 ` - COMPASS (class 2) bond * :doc:`fene ` - FENE (finite-extensible non-linear elastic) bond * :doc:`fene/expand ` - FENE bonds with variable size particles +* :doc:`gaussian ` - multicentered Gaussian-based bond potential * :doc:`gromos ` - GROMOS force field bond * :doc:`harmonic ` - harmonic bond * :doc:`harmonic/shift ` - shifted harmonic bond diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 4d026f462f..37ddbbeb3e 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -205,6 +205,7 @@ accelerated styles exist. * :doc:`efield ` - impose electric field on system * :doc:`ehex ` - enhanced heat exchange algorithm * :doc:`electron/stopping ` - electronic stopping power as a friction force +* :doc:`electron/stopping/fit ` - electronic stopping power as a friction force * :doc:`enforce2d ` - zero out z-dimension velocity and force * :doc:`eos/cv ` - * :doc:`eos/table ` - diff --git a/doc/src/fix_electron_stopping.rst b/doc/src/fix_electron_stopping.rst index 76455e88df..4a61afbc09 100644 --- a/doc/src/fix_electron_stopping.rst +++ b/doc/src/fix_electron_stopping.rst @@ -1,28 +1,41 @@ .. index:: fix electron/stopping +.. index:: fix electron/stopping/fit fix electron/stopping command ============================= +fix electron/stopping/fit command +================================= + Syntax """""" .. parsed-literal:: - fix ID group-ID electron/stopping Ecut file keyword value ... + fix ID group-ID style args * ID, group-ID are documented in :doc:`fix ` command -* electron/stopping = style name of this fix command -* Ecut = minimum kinetic energy for electronic stopping (energy units) -* file = name of the file containing the electronic stopping power table -* zero or more keyword/value pairs may be appended to args -* keyword = *region* or *minneigh* +* style = *electron/stopping* or *electron/stopping/fit* .. parsed-literal:: + *electron/stopping* args = Ecut file keyword value ... + Ecut = minimum kinetic energy for electronic stopping (energy units) + file = name of the file containing the electronic stopping power table + + *electron/stopping/fit* args = Ecut c1 c2 ... + Ecut = minimum kinetic energy for electronic stopping (energy units) + c1 c2 = linear and quadratic coefficients for the fitted quadratic polynomial + +* zero or more keyword/value pairs may be appended to args for style = *electron/stopping* + + .. parsed-literal:: + + keyword = *region* or *minneigh* *region* value = region-ID - region-ID = region, whose atoms will be affected by this fix + region-ID = region whose atoms will be affected by this fix *minneigh* value = minneigh - minneigh = minimum number of neighbors an atom to have stopping applied + minneigh = minimum number of neighbors an atom to have stopping applied Examples """""""" @@ -32,6 +45,8 @@ Examples fix el all electron/stopping 10.0 elstop-table.txt fix el all electron/stopping 10.0 elstop-table.txt minneigh 3 fix el mygroup electron/stopping 1.0 elstop-table.txt region bulk + fix 1 all electron/stopping/fit 4.63 3.3e-3 4.0e-8 + fix 1 all electron/stopping/fit 3.49 1.8e-3 9.0e-8 7.57 4.2e-3 5.0e-8 Description """"""""""" @@ -129,6 +144,19 @@ scientific publications, experimental databases or by using of the impact parameter of the ion; these results can be used to derive the stopping power. +---------- + +Style *electron/stopping/fit* calculates the electronic stopping power +and cumulative energy lost to the electron gas via a quadratic functional +and applies a drag force to the classical equations-of-motion for all +atoms moving above some minimum cutoff velocity (i.e., kinetic energy). +These coefficients can be determined by fitting a quadratic polynomial to +electronic stopping data predicted by, for example, SRIM or TD-DFT. Multiple +'Ecut c1 c2' values can be provided for multi-species simulations in the order +of the atom types. There is an examples/USER/misc/electron_stopping/ directory, +which illustrates uses of this command. Details of this implementation are +further described in :ref:`Stewart2018 ` and :ref:`Lee2020 `. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -181,3 +209,11 @@ The default is no limitation by region, and minneigh = 1. .. _PASS: **(PASS)** PASS webpage: https://www.sdu.dk/en/DPASS + +.. _Stewart2018: + +**(Stewart2018)** J.A. Stewart, et al. (2018) Journal of Applied Physics, 123(16), 165902. + +.. _Lee2020: + +**(Lee2020)** C.W. Lee, et al. (2020) Physical Review B, 102(2), 024107. diff --git a/doc/src/package.rst b/doc/src/package.rst index 238dc50115..725536310a 100644 --- a/doc/src/package.rst +++ b/doc/src/package.rst @@ -448,8 +448,7 @@ does not require atomic operations in the calculation of pair forces. For that reason, *full* is the default setting for GPUs. However, when running on CPUs, a *half* neighbor list is the default because it are often faster, just as it is for non-accelerated pair styles. Similarly, -the *neigh/qeq* keyword determines how neighbor lists are built for :doc:`fix qeq/reax/kk `. If not explicitly set, the value of -*neigh/qeq* will match *neigh*\ . +the *neigh/qeq* keyword determines how neighbor lists are built for :doc:`fix qeq/reax/kk `. If the *neigh/thread* keyword is set to *off*\ , then the KOKKOS package threads only over atoms. However, for small systems, this may not expose diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index d559425b86..eb57526b78 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1128,6 +1128,7 @@ Gonzalez-Melchor googlemail googletest Gordan +Goudeau GPa gpu gpuID @@ -1887,6 +1888,7 @@ mie Mie Mij Mikami +Milano Militzer Minary mincap @@ -1990,6 +1992,7 @@ multi multibody Multibody multicenter +multicentered multicmd multicomponent multicore @@ -2010,6 +2013,7 @@ muVT mux muy muz +Müller mv mV Mvapich @@ -2457,6 +2461,7 @@ polydispersity polyelectrolyte polyhedra polymorphism +Polym popen Popov popstore @@ -2788,6 +2793,7 @@ Schulten Schunk Schuring Schwen +Sci screenshot screenshots Scripps diff --git a/examples/USER/misc/electron_stopping/in.cascade_AlCu b/examples/USER/misc/electron_stopping/in.cascade_AlCu new file mode 100644 index 0000000000..26a8852f24 --- /dev/null +++ b/examples/USER/misc/electron_stopping/in.cascade_AlCu @@ -0,0 +1,38 @@ +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Al lattice with a single incident Cu atom - multiple species simulation +# *** + +units metal +boundary p p p + +lattice fcc 4.0495 + +region box block -10 10 -10 10 -10 10 +create_box 2 box +create_atoms 1 box + +pair_style eam/alloy +pair_coeff * * ../../../../potentials/AlCu.eam.alloy Al Cu + +mass 1 26.982 +mass 2 63.546 + +velocity all create 300 42534 mom yes rot yes + +set atom 1 type 2 +group pka id 1 +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 3.49 1.8e-3 9.0e-8 7.57 4.2e-3 5.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 diff --git a/examples/USER/misc/electron_stopping/in.cascade_SiSi b/examples/USER/misc/electron_stopping/in.cascade_SiSi new file mode 100644 index 0000000000..4596f51904 --- /dev/null +++ b/examples/USER/misc/electron_stopping/in.cascade_SiSi @@ -0,0 +1,36 @@ +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Si lattice with one primary knock-on atom (PKA) - single species simulation +# *** + +units metal +boundary p p p + +lattice diamond 5.431 + +region box block -10 10 -10 10 -10 10 +create_box 1 box +create_atoms 1 box + +pair_style tersoff/zbl +pair_coeff * * ../../../../potentials/SiC.tersoff.zbl Si + +mass 1 28.0855 + +velocity all create 300 42534 mom yes rot yes + +group pka id 1 +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 4.63 3.3e-3 4.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 diff --git a/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.1 b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.1 new file mode 100644 index 0000000000..942a65ea86 --- /dev/null +++ b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.1 @@ -0,0 +1,115 @@ +LAMMPS (18 Sep 2020) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94) + using 1 OpenMP thread(s) per MPI task +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Al lattice with a single incident Cu atom - multiple species simulation +# *** + +units metal +boundary p p p + +lattice fcc 4.0495 +Lattice spacing in x,y,z = 4.0495000 4.0495000 4.0495000 + +region box block -10 10 -10 10 -10 10 +create_box 2 box +Created orthogonal box = (-40.495000 -40.495000 -40.495000) to (40.495000 40.495000 40.495000) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 32000 atoms + create_atoms CPU = 0.004 seconds + +pair_style eam/alloy +pair_coeff * * ../../../../potentials/AlCu.eam.alloy Al Cu +Reading eam/alloy potential file ../../../../potentials/AlCu.eam.alloy with DATE: 2008-10-01 + +mass 1 26.982 +mass 2 63.546 + +velocity all create 300 42534 mom yes rot yes + +set atom 1 type 2 +Setting atom values ... + 1 settings made for type +group pka id 1 +1 atoms in group pka +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 3.49 1.8e-3 9.0e-8 7.57 4.2e-3 5.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 8.6825 + ghost atom cutoff = 8.6825 + binsize = 4.34125, bins = 19 19 19 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair eam/alloy, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 23.27 | 23.27 | 23.27 Mbytes +Step Dt Time Temp PotEng KinEng f_3 + 0 2.4879625e-05 0 53029.167 -106186.96 219339 0 + 5 2.4881895e-05 0.00012440247 53019.542 -106184.2 219299.18 36.968257 + 10 1.0341799e-05 0.00021830163 53006.803 -106159.67 219246.49 64.866504 + 15 5.6753071e-06 0.0002583636 52996.09 -106127.26 219202.18 76.766363 + 20 4.0862476e-06 0.0002830119 52987.566 -106099.31 219166.93 84.086344 + 25 3.3677164e-06 0.00030179992 52980.99 -106077.67 219139.73 89.665096 + 30 3.0218093e-06 0.00031784053 52976.302 -106063.03 219120.34 94.427505 + 35 2.8935922e-06 0.00033262008 52973.489 -106055.77 219108.7 98.815168 + 40 2.9394209e-06 0.00034711037 52972.541 -106056.15 219104.78 103.11678 + 45 3.1822694e-06 0.00036220043 52973.425 -106064.3 219108.44 107.59642 + 50 3.743784e-06 0.00037905999 52976.072 -106080.26 219119.38 112.60152 + 55 5.0685177e-06 0.0003999287 52980.344 -106104.15 219137.05 118.79722 + 60 1.0189784e-05 0.00043198334 52985.861 -106136.52 219159.88 128.31484 + 65 1.8636384e-05 0.00052946777 52985.275 -106162.53 219157.45 157.2625 + 70 1.844772e-05 0.00061001061 52977.927 -106155.89 219127.06 181.17691 + 75 2.4893022e-05 0.00072690857 52972.391 -106168.08 219104.16 215.88136 + 80 7.390618e-06 0.00081149431 52959.379 -106139.89 219050.34 240.98969 + 85 4.1547853e-06 0.00084037647 52948.078 -106101.74 219003.59 249.56079 + 90 2.9763749e-06 0.00085843347 52938.03 -106065.54 218962.03 254.91825 + 95 2.3727508e-06 0.00087197081 52929.043 -106032.37 218924.86 258.93397 + 100 2.0138478e-06 0.00088304936 52921.103 -106002.81 218892.02 262.21977 +Loop time of 9.53154 on 1 procs for 100 steps with 32000 atoms + +Performance: 0.002 ns/day, 13147.213 hours/ns, 10.491 timesteps/s +99.7% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 8.815 | 8.815 | 8.815 | 0.0 | 92.48 +Neigh | 0.38408 | 0.38408 | 0.38408 | 0.0 | 4.03 +Comm | 0.029049 | 0.029049 | 0.029049 | 0.0 | 0.30 +Output | 0.0025912 | 0.0025912 | 0.0025912 | 0.0 | 0.03 +Modify | 0.28624 | 0.28624 | 0.28624 | 0.0 | 3.00 +Other | | 0.01456 | | | 0.15 + +Nlocal: 32000.0 ave 32000 max 32000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 26825.0 ave 26825 max 26825 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 2.81599e+06 ave 2.81599e+06 max 2.81599e+06 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 2815993 +Ave neighs/atom = 87.999781 +Neighbor list builds = 5 +Dangerous builds = 3 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:00:09 diff --git a/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.4 b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.4 new file mode 100644 index 0000000000..aa9d91d43d --- /dev/null +++ b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_AlCu.intel.4 @@ -0,0 +1,115 @@ +LAMMPS (18 Sep 2020) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94) + using 1 OpenMP thread(s) per MPI task +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Al lattice with a single incident Cu atom - multiple species simulation +# *** + +units metal +boundary p p p + +lattice fcc 4.0495 +Lattice spacing in x,y,z = 4.0495000 4.0495000 4.0495000 + +region box block -10 10 -10 10 -10 10 +create_box 2 box +Created orthogonal box = (-40.495000 -40.495000 -40.495000) to (40.495000 40.495000 40.495000) + 1 by 2 by 2 MPI processor grid +create_atoms 1 box +Created 32000 atoms + create_atoms CPU = 0.001 seconds + +pair_style eam/alloy +pair_coeff * * ../../../../potentials/AlCu.eam.alloy Al Cu +Reading eam/alloy potential file ../../../../potentials/AlCu.eam.alloy with DATE: 2008-10-01 + +mass 1 26.982 +mass 2 63.546 + +velocity all create 300 42534 mom yes rot yes + +set atom 1 type 2 +Setting atom values ... + 1 settings made for type +group pka id 1 +1 atoms in group pka +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 3.49 1.8e-3 9.0e-8 7.57 4.2e-3 5.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 8.6825 + ghost atom cutoff = 8.6825 + binsize = 4.34125, bins = 19 19 19 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair eam/alloy, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d/newton + bin: standard +Per MPI rank memory allocation (min/avg/max) = 9.014 | 9.014 | 9.014 Mbytes +Step Dt Time Temp PotEng KinEng f_3 + 0 2.4879625e-05 0 53029.167 -106186.96 219339 0 + 5 2.4881895e-05 0.00012440247 53019.542 -106184.2 219299.18 36.968257 + 10 1.0341742e-05 0.00021830123 53006.803 -106159.67 219246.49 64.866388 + 15 5.6752732e-06 0.00025836298 52996.09 -106127.26 219202.18 76.76618 + 20 4.0862174e-06 0.00028301112 52987.566 -106099.31 219166.93 84.086114 + 25 3.3676848e-06 0.00030179899 52980.99 -106077.67 219139.73 89.664821 + 30 3.021773e-06 0.00031783944 52976.302 -106063.02 219120.33 94.427181 + 35 2.8935472e-06 0.00033261879 52973.489 -106055.77 219108.7 98.814784 + 40 2.9393606e-06 0.00034710883 52972.541 -106056.15 219104.78 103.11632 + 45 3.1821803e-06 0.00036219854 52973.425 -106064.29 219108.43 107.59586 + 50 3.7436309e-06 0.00037905755 52976.071 -106080.26 219119.38 112.60079 + 55 5.0681667e-06 0.0003999252 52980.343 -106104.15 219137.05 118.79618 + 60 1.0187808e-05 0.00043197649 52985.861 -106136.51 219159.87 128.3128 + 65 1.8643099e-05 0.00052944037 52985.278 -106162.53 219157.46 157.25436 + 70 1.8445045e-05 0.00060999223 52977.928 -106155.89 219127.06 181.17146 + 75 2.4893021e-05 0.00072688076 52972.393 -106168.08 219104.17 215.8731 + 80 7.3916674e-06 0.0008114874 52959.382 -106139.9 219050.35 240.98764 + 85 4.1550998e-06 0.00084037284 52948.08 -106101.75 219003.6 249.55971 + 90 2.976545e-06 0.00085843108 52938.032 -106065.55 218962.04 254.91754 + 95 2.3728646e-06 0.00087196913 52929.045 -106032.38 218924.87 258.93348 + 100 2.0139362e-06 0.00088304819 52921.106 -106002.82 218892.03 262.21943 +Loop time of 2.45676 on 4 procs for 100 steps with 32000 atoms + +Performance: 0.007 ns/day, 3388.559 hours/ns, 40.704 timesteps/s +99.8% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 2.2127 | 2.2257 | 2.2399 | 0.8 | 90.59 +Neigh | 0.095856 | 0.097842 | 0.10197 | 0.8 | 3.98 +Comm | 0.03626 | 0.054908 | 0.069787 | 6.3 | 2.23 +Output | 0.00088538 | 0.0011005 | 0.0017236 | 1.1 | 0.04 +Modify | 0.072449 | 0.072553 | 0.072611 | 0.0 | 2.95 +Other | | 0.004684 | | | 0.19 + +Nlocal: 8000.00 ave 8033 max 7977 min +Histogram: 1 0 1 1 0 0 0 0 0 1 +Nghost: 12605.0 ave 12628 max 12572 min +Histogram: 1 0 0 0 0 0 1 1 0 1 +Neighs: 703998.0 ave 706570 max 702282 min +Histogram: 1 1 0 0 1 0 0 0 0 1 + +Total # of neighbors = 2815992 +Ave neighs/atom = 87.999750 +Neighbor list builds = 5 +Dangerous builds = 3 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:00:02 diff --git a/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.1 b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.1 new file mode 100644 index 0000000000..526d88977c --- /dev/null +++ b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.1 @@ -0,0 +1,113 @@ +LAMMPS (18 Sep 2020) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94) + using 1 OpenMP thread(s) per MPI task +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Si lattice with one primary knock-on atom (PKA) - single species simulation +# *** + +units metal +boundary p p p + +lattice diamond 5.431 +Lattice spacing in x,y,z = 5.4310000 5.4310000 5.4310000 + +region box block -10 10 -10 10 -10 10 +create_box 1 box +Created orthogonal box = (-54.310000 -54.310000 -54.310000) to (54.310000 54.310000 54.310000) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 64000 atoms + create_atoms CPU = 0.008 seconds + +pair_style tersoff/zbl +pair_coeff * * ../../../../potentials/SiC.tersoff.zbl Si +Reading tersoff/zbl potential file ../../../../potentials/SiC.tersoff.zbl with DATE: 2009-04-15 + +mass 1 28.0855 + +velocity all create 300 42534 mom yes rot yes + +group pka id 1 +1 atoms in group pka +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 4.63 3.3e-3 4.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 5 + ghost atom cutoff = 5 + binsize = 2.5, bins = 44 44 44 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair tersoff/zbl, perpetual + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 23.91 | 23.91 | 23.91 Mbytes +Step Dt Time Temp PotEng KinEng f_3 + 0 2.4879625e-05 0 21258.724 -296293.96 175863.22 0 + 5 9.2462185e-06 9.0354398e-05 21252.059 -296277.65 175808.08 38.834393 + 10 6.1558479e-06 0.00013003059 21246.736 -296250.63 175764.05 55.881868 + 15 4.9002335e-06 0.00015681379 21242.649 -296228.27 175730.24 67.386915 + 20 5.853687e-06 0.00018286317 21239.571 -296213.99 175704.78 78.574676 + 25 7.0182823e-06 0.00021483214 21237.286 -296208.81 175685.87 92.303028 + 30 8.2570048e-06 0.00025083646 21236.034 -296213.99 175675.52 107.76296 + 35 6.4734302e-06 0.00029194185 21233.473 -296210.33 175654.33 125.41212 + 40 7.3445302e-06 0.00032561085 21231.196 -296205.94 175635.5 139.86641 + 45 6.9480705e-06 0.00036268325 21230.21 -296213.68 175627.33 155.78046 + 50 7.2224188e-06 0.00039655436 21230.512 -296230.74 175629.84 170.32001 + 55 1.0773409e-05 0.00044221823 21230.023 -296246.37 175625.79 189.92217 + 60 5.7527075e-06 0.00048339879 21226.064 -296231.33 175593.04 207.5982 + 65 5.8568503e-06 0.0005110075 21222.544 -296213.97 175563.92 219.44643 + 70 6.7430644e-06 0.00054252027 21220.179 -296207.92 175544.35 232.96808 + 75 7.0523029e-06 0.00057648256 21219.781 -296219.19 175541.06 247.53974 + 80 1.784394e-05 0.00062210154 21221.276 -296251.35 175553.43 267.11364 + 85 2.1885193e-05 0.0007395532 21218.037 -296274.94 175526.64 317.50995 + 90 8.233509e-06 0.00081518257 21211.247 -296251.53 175470.47 349.95382 + 95 5.1490725e-06 0.00084789415 21205.33 -296216.55 175421.52 363.982 + 100 5.7628664e-06 0.0008764946 21200.168 -296186.27 175378.81 376.24357 +Loop time of 20.4868 on 1 procs for 100 steps with 64000 atoms + +Performance: 0.002 ns/day, 9874.909 hours/ns, 4.881 timesteps/s +99.8% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 19.397 | 19.397 | 19.397 | 0.0 | 94.68 +Neigh | 0.45332 | 0.45332 | 0.45332 | 0.0 | 2.21 +Comm | 0.035144 | 0.035144 | 0.035144 | 0.0 | 0.17 +Output | 0.0040397 | 0.0040397 | 0.0040397 | 0.0 | 0.02 +Modify | 0.5739 | 0.5739 | 0.5739 | 0.0 | 2.80 +Other | | 0.0229 | | | 0.11 + +Nlocal: 64000.0 ave 64000 max 64000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 18341.0 ave 18341 max 18341 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0.00000 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 1.79199e+06 ave 1.79199e+06 max 1.79199e+06 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 1791990 +Ave neighs/atom = 27.999844 +Neighbor list builds = 7 +Dangerous builds = 2 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:00:20 diff --git a/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.4 b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.4 new file mode 100644 index 0000000000..ed02f0edb5 --- /dev/null +++ b/examples/USER/misc/electron_stopping/log.18Sep2020.cascade_SiSi.intel.4 @@ -0,0 +1,113 @@ +LAMMPS (18 Sep 2020) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94) + using 1 OpenMP thread(s) per MPI task +# *** +# Example input for including electronic stopping effects using fix electron/stopping/fit +# Si lattice with one primary knock-on atom (PKA) - single species simulation +# *** + +units metal +boundary p p p + +lattice diamond 5.431 +Lattice spacing in x,y,z = 5.4310000 5.4310000 5.4310000 + +region box block -10 10 -10 10 -10 10 +create_box 1 box +Created orthogonal box = (-54.310000 -54.310000 -54.310000) to (54.310000 54.310000 54.310000) + 1 by 2 by 2 MPI processor grid +create_atoms 1 box +Created 64000 atoms + create_atoms CPU = 0.003 seconds + +pair_style tersoff/zbl +pair_coeff * * ../../../../potentials/SiC.tersoff.zbl Si +Reading tersoff/zbl potential file ../../../../potentials/SiC.tersoff.zbl with DATE: 2009-04-15 + +mass 1 28.0855 + +velocity all create 300 42534 mom yes rot yes + +group pka id 1 +1 atoms in group pka +velocity pka set 1120 1620 400 + +fix 1 all nve +fix 2 all dt/reset 1 NULL 0.001 0.05 emax 10.0 +fix 3 all electron/stopping/fit 4.63 3.3e-3 4.0e-8 + +thermo 5 +thermo_style custom step dt time temp pe ke f_3 +thermo_modify lost warn flush yes + +#dump 0 all custom 10 dump.pka_* id type x y z vx vy vz fx fy fz +#dump_modify 0 first yes + +run 100 +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 5 + ghost atom cutoff = 5 + binsize = 2.5, bins = 44 44 44 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair tersoff/zbl, perpetual + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 7.777 | 7.777 | 7.777 Mbytes +Step Dt Time Temp PotEng KinEng f_3 + 0 2.4879625e-05 0 21258.724 -296293.96 175863.22 0 + 5 9.2462211e-06 9.0354422e-05 21252.059 -296277.66 175808.08 38.834403 + 10 6.1568847e-06 0.00013003069 21246.736 -296250.63 175764.05 55.881913 + 15 4.8997821e-06 0.00015681555 21242.649 -296228.27 175730.24 67.387669 + 20 5.8536714e-06 0.00018286203 21239.571 -296213.99 175704.78 78.574187 + 25 7.0160073e-06 0.00021483347 21237.285 -296208.8 175685.86 92.303599 + 30 8.2556386e-06 0.00025082046 21236.035 -296213.99 175675.53 107.75609 + 35 6.4735852e-06 0.00029193023 21233.474 -296210.34 175654.34 125.40713 + 40 7.3441556e-06 0.00032559783 21231.197 -296205.94 175635.5 139.86082 + 45 6.9483099e-06 0.00036267022 21230.21 -296213.68 175627.34 155.77487 + 50 7.2213562e-06 0.0003965413 21230.513 -296230.74 175629.84 170.31441 + 55 1.0776037e-05 0.00044219672 21230.024 -296246.38 175625.8 189.91293 + 60 5.7538246e-06 0.0004833796 21226.067 -296231.34 175593.06 207.58996 + 65 5.856213e-06 0.00051099409 21222.546 -296213.98 175563.94 219.44067 + 70 6.7431217e-06 0.00054250526 21220.18 -296207.92 175544.37 232.96164 + 75 7.0518411e-06 0.00057646788 21219.781 -296219.18 175541.07 247.53344 + 80 1.7829072e-05 0.00062207162 21221.276 -296251.34 175553.43 267.1008 + 85 2.1894958e-05 0.0007395084 21218.04 -296274.95 175526.66 317.49073 + 90 8.2365472e-06 0.00081516502 21211.25 -296251.55 175470.49 349.94629 + 95 5.1493496e-06 0.00084788428 21205.333 -296216.57 175421.54 363.97777 + 100 5.7652664e-06 0.00087648406 21200.171 -296186.3 175378.84 376.23905 +Loop time of 5.23182 on 4 procs for 100 steps with 64000 atoms + +Performance: 0.010 ns/day, 2520.759 hours/ns, 19.114 timesteps/s +99.7% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.8069 | 4.8618 | 4.9229 | 2.0 | 92.93 +Neigh | 0.11442 | 0.11856 | 0.12948 | 1.8 | 2.27 +Comm | 0.040943 | 0.098988 | 0.15807 | 14.7 | 1.89 +Output | 0.0013075 | 0.0014801 | 0.0019936 | 0.8 | 0.03 +Modify | 0.14496 | 0.14502 | 0.1451 | 0.0 | 2.77 +Other | | 0.005981 | | | 0.11 + +Nlocal: 16000.0 ave 16017 max 15987 min +Histogram: 1 0 1 0 0 1 0 0 0 1 +Nghost: 8044.00 ave 8058 max 8026 min +Histogram: 1 0 0 0 1 0 0 0 1 1 +Neighs: 0.00000 ave 0 max 0 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +FullNghs: 447998.0 ave 448471 max 447634 min +Histogram: 1 0 1 0 0 1 0 0 0 1 + +Total # of neighbors = 1791990 +Ave neighs/atom = 27.999844 +Neighbor list builds = 7 +Dangerous builds = 2 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:00:05 diff --git a/lib/kim/Install.py b/lib/kim/Install.py index 33d0071dde..3ebcb8e237 100644 --- a/lib/kim/Install.py +++ b/lib/kim/Install.py @@ -18,12 +18,13 @@ parser = ArgumentParser(prog='Install.py', # settings thisdir = fullpath('.') -version = "2.1.3" +version = "2.2.0" # known checksums for different KIM-API versions. used to validate the download. checksums = { \ '2.1.2' : '6ac52e14ef52967fc7858220b208cba5', \ '2.1.3' : '6ee829a1bbba5f8b9874c88c4c4ebff8', \ + '2.2.0' : 'e7f944e1593cffd7444679a660607f6c', \ } diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp index 605d5cd7ee..5da513d972 100644 --- a/src/KIM/kim_init.cpp +++ b/src/KIM/kim_init.cpp @@ -92,7 +92,12 @@ void KimInit::command(int narg, char **arg) strcpy(user_units,arg[1]); if (narg == 3) { if (strcmp(arg[2],"unit_conversion_mode")==0) unit_conversion_mode = true; - else { error->all(FLERR,"Illegal kim_init command"); } + else { + error->all(FLERR,fmt::format("Illegal kim_init command.\nThe argument " + "followed by unit_style {} is an optional " + "argument and when is used must " + "be unit_conversion_mode", user_units)); + } } else unit_conversion_mode = false; char *model_units; @@ -122,38 +127,41 @@ void get_kim_unit_names( KIM_TimeUnit & timeUnit, Error * error) { - if ((strcmp(system,"real")==0)) { + if (strcmp(system,"real") == 0) { lengthUnit = KIM_LENGTH_UNIT_A; energyUnit = KIM_ENERGY_UNIT_kcal_mol; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_fs; - } else if ((strcmp(system,"metal")==0)) { + } else if (strcmp(system,"metal") == 0) { lengthUnit = KIM_LENGTH_UNIT_A; energyUnit = KIM_ENERGY_UNIT_eV; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_ps; - } else if ((strcmp(system,"si")==0)) { + } else if (strcmp(system,"si") == 0) { lengthUnit = KIM_LENGTH_UNIT_m; energyUnit = KIM_ENERGY_UNIT_J; chargeUnit = KIM_CHARGE_UNIT_C; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_s; - } else if ((strcmp(system,"cgs")==0)) { + } else if (strcmp(system,"cgs") == 0) { lengthUnit = KIM_LENGTH_UNIT_cm; energyUnit = KIM_ENERGY_UNIT_erg; chargeUnit = KIM_CHARGE_UNIT_statC; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_s; - } else if ((strcmp(system,"electron")==0)) { + } else if (strcmp(system,"electron") == 0) { lengthUnit = KIM_LENGTH_UNIT_Bohr; energyUnit = KIM_ENERGY_UNIT_Hartree; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_fs; - } else if ((strcmp(system,"lj")==0)) { - error->all(FLERR,"LAMMPS unit_style lj not supported by KIM models"); + } else if (strcmp(system,"lj") == 0 || + strcmp(system,"micro") ==0 || + strcmp(system,"nano")==0) { + error->all(FLERR,fmt::format("LAMMPS unit_style {} not supported " + "by KIM models", system)); } else { error->all(FLERR,"Unknown unit_style"); } @@ -171,19 +179,19 @@ void KimInit::determine_model_type_and_units(char * model_name, KIM_TemperatureUnit temperatureUnit; KIM_TimeUnit timeUnit; int units_accepted; - KIM_Collections * kim_Coll; + KIM_Collections * collections; KIM_CollectionItemType itemType; - int kim_error = KIM_Collections_Create(&kim_Coll); - if (kim_error) { - error->all(FLERR,"Unable to access KIM Collections to find Model."); - } + int kim_error = KIM_Collections_Create(&collections); + if (kim_error) + error->all(FLERR,"Unable to access KIM Collections to find Model"); - kim_error = KIM_Collections_GetItemType(kim_Coll, model_name, &itemType); - if (kim_error) { - error->all(FLERR,"KIM Model name not found."); - } - KIM_Collections_Destroy(&kim_Coll); + auto logID = fmt::format("{}_Collections", comm->me); + KIM_Collections_SetLogID(collections, logID.c_str()); + + kim_error = KIM_Collections_GetItemType(collections, model_name, &itemType); + if (kim_error) error->all(FLERR,"KIM Model name not found"); + KIM_Collections_Destroy(&collections); if (KIM_CollectionItemType_Equal(itemType, KIM_COLLECTION_ITEM_TYPE_portableModel)) { @@ -199,12 +207,13 @@ void KimInit::determine_model_type_and_units(char * model_name, &units_accepted, &pkim); - if (kim_error) - error->all(FLERR,"Unable to load KIM Simulator Model."); + if (kim_error) error->all(FLERR,"Unable to load KIM Simulator Model"); model_type = MO; if (units_accepted) { + logID = fmt::format("{}_Model", comm->me); + KIM_Model_SetLogID(pkim, logID.c_str()); *model_units = new char[strlen(user_units)+1]; strcpy(*model_units,user_units); return; @@ -226,6 +235,8 @@ void KimInit::determine_model_type_and_units(char * model_name, &units_accepted, &pkim); if (units_accepted) { + logID = fmt::format("{}_Model", comm->me); + KIM_Model_SetLogID(pkim, logID.c_str()); *model_units = new char[strlen(systems[i])+1]; strcpy(*model_units,systems[i]); return; @@ -238,37 +249,40 @@ void KimInit::determine_model_type_and_units(char * model_name, error->all(FLERR,"KIM Model does not support the requested unit system"); } } else if (KIM_CollectionItemType_Equal( - itemType, KIM_COLLECTION_ITEM_TYPE_simulatorModel)) { - KIM_SimulatorModel * kim_SM; - kim_error = KIM_SimulatorModel_Create(model_name, &kim_SM); + itemType, KIM_COLLECTION_ITEM_TYPE_simulatorModel)) { + KIM_SimulatorModel * simulatorModel; + kim_error = KIM_SimulatorModel_Create(model_name, &simulatorModel); if (kim_error) - error->all(FLERR,"Unable to load KIM Simulator Model."); + error->all(FLERR,"Unable to load KIM Simulator Model"); model_type = SM; + logID = fmt::format("{}_SimulatorModel", comm->me); + KIM_SimulatorModel_SetLogID(simulatorModel, logID.c_str()); + int sim_fields; int sim_lines; char const * sim_field; char const * sim_value; - KIM_SimulatorModel_GetNumberOfSimulatorFields(kim_SM, &sim_fields); - KIM_SimulatorModel_CloseTemplateMap(kim_SM); + KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields); + KIM_SimulatorModel_CloseTemplateMap(simulatorModel); for (int i=0; i < sim_fields; ++i) { KIM_SimulatorModel_GetSimulatorFieldMetadata( - kim_SM,i,&sim_lines,&sim_field); + simulatorModel, i, &sim_lines, &sim_field); if (0 == strcmp(sim_field,"units")) { - KIM_SimulatorModel_GetSimulatorFieldLine(kim_SM,i,0,&sim_value); + KIM_SimulatorModel_GetSimulatorFieldLine( + simulatorModel, i, 0, &sim_value); int len=strlen(sim_value)+1; - *model_units = new char[len]; strcpy(*model_units,sim_value); + *model_units = new char[len]; + strcpy(*model_units,sim_value); break; } } - KIM_SimulatorModel_Destroy(&kim_SM); + KIM_SimulatorModel_Destroy(&simulatorModel); if ((! unit_conversion_mode) && (strcmp(*model_units, user_units)!=0)) { - std::string mesg("Incompatible units for KIM Simulator Model, " - "required units = "); - mesg += *model_units; - error->all(FLERR,mesg); + error->all(FLERR,fmt::format("Incompatible units for KIM Simulator Model" + ", required units = {}", *model_units)); } } } @@ -295,13 +309,16 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM KIM_SimulatorModel * simulatorModel; if (model_type == SM) { int kim_error = - KIM_SimulatorModel_Create(model_name,&simulatorModel); + KIM_SimulatorModel_Create(model_name, &simulatorModel); if (kim_error) - error->all(FLERR,"Unable to load KIM Simulator Model."); + error->all(FLERR,"Unable to load KIM Simulator Model"); + + auto logID = fmt::format("{}_SimulatorModel", comm->me); + KIM_SimulatorModel_SetLogID(simulatorModel, logID.c_str()); char const *sim_name, *sim_version; KIM_SimulatorModel_GetSimulatorNameAndVersion( - simulatorModel,&sim_name, &sim_version); + simulatorModel, &sim_name, &sim_version); if (0 != strcmp(sim_name,"LAMMPS")) error->all(FLERR,"Incompatible KIM Simulator Model"); @@ -389,7 +406,7 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM mesg += fmt::format(" {:<8} | {:<{}} | {:<10} | {}\n",i+1,str_name, max_len,data_type,extent); } - } else mesg += "No mutable parameters. \n"; + } else mesg += "No mutable parameters.\n"; KIM_Model_Destroy(&pkim); input->write_echo(mesg); @@ -447,7 +464,7 @@ void KimInit::do_variables(const std::string &from, const std::string &to) conversion_factor); if (ier != 0) error->all(FLERR,fmt::format("Unable to obtain conversion factor: " - "unit = {}; from = {}; to = {}.", + "unit = {}; from = {}; to = {}", units[i], from, to)); variable->internal_set(v_unit,conversion_factor); @@ -461,23 +478,28 @@ void KimInit::do_variables(const std::string &from, const std::string &to) void KimInit::write_log_cite(char *model_name) { - KIM_Collections * coll; - int err = KIM_Collections_Create(&coll); + KIM_Collections * collections; + int err = KIM_Collections_Create(&collections); if (err) return; + auto logID = fmt::format("{}_Collections", comm->me); + KIM_Collections_SetLogID(collections, logID.c_str()); + int extent; if (model_type == MO) { err = KIM_Collections_CacheListOfItemMetadataFiles( - coll,KIM_COLLECTION_ITEM_TYPE_portableModel,model_name,&extent); + collections, KIM_COLLECTION_ITEM_TYPE_portableModel, + model_name,&extent); } else if (model_type == SM) { err = KIM_Collections_CacheListOfItemMetadataFiles( - coll,KIM_COLLECTION_ITEM_TYPE_simulatorModel,model_name,&extent); + collections, KIM_COLLECTION_ITEM_TYPE_simulatorModel, + model_name, &extent); } else { - error->all(FLERR,"Unknown model type."); + error->all(FLERR,"Unknown model type"); } if (err) { - KIM_Collections_Destroy(&coll); + KIM_Collections_Destroy(&collections); return; } @@ -486,7 +508,8 @@ void KimInit::write_log_cite(char *model_name) int availableAsString; char const * fileString; err = KIM_Collections_GetItemMetadataFile( - coll,i,&fileName,nullptr,nullptr,&availableAsString,&fileString); + collections, i, &fileName, nullptr, nullptr, + &availableAsString, &fileString); if (err) continue; if (0 == strncmp("kimcite",fileName,7)) { @@ -494,5 +517,5 @@ void KimInit::write_log_cite(char *model_name) } } - KIM_Collections_Destroy(&coll); + KIM_Collections_Destroy(&collections); } diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp index 73e0f5b275..713f91b56c 100644 --- a/src/KIM/kim_interactions.cpp +++ b/src/KIM/kim_interactions.cpp @@ -16,6 +16,7 @@ Ryan S. Elliott (UMN) Ellad B. Tadmor (UMN) Ronald Miller (Carleton U) + Yaser Afshar (UMN) ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- @@ -87,6 +88,7 @@ void KimInteractions::command(int narg, char **arg) if (!domain->box_exist) error->all(FLERR,"Must use 'kim_interactions' command after " "simulation box is defined"); + do_setup(narg,arg); } @@ -98,7 +100,11 @@ void KimInteractions::do_setup(int narg, char **arg) if ((narg == 1) && (0 == strcmp("fixed_types",arg[0]))) { fixed_types = true; } else if (narg != atom->ntypes) { - error->all(FLERR,"Illegal kim_interactions command"); + error->all(FLERR,fmt::format("Illegal kim_interactions command.\nThe " + "LAMMPS simulation has {} atom type(s), but " + "{} chemical species passed to the " + "kim_interactions command", + atom->ntypes, narg)); } else { fixed_types = false; } @@ -121,22 +127,20 @@ void KimInteractions::do_setup(int narg, char **arg) input->write_echo("#=== BEGIN kim_interactions ==================================\n"); if (simulatorModel) { - if (!fixed_types) { - std::string delimiter(""); - std::string atom_type_sym_list; - std::string atom_type_num_list; + std::string atom_type_sym_list = + fmt::format("{}", fmt::join(arg, arg + narg, " ")); + + std::string atom_type_num_list = + fmt::format("{}", species_to_atomic_no(arg[0])); - for (int i = 0; i < narg; i++) { - atom_type_sym_list += delimiter + arg[i]; - atom_type_num_list += delimiter + std::to_string(species_to_atomic_no(arg[i])); - delimiter = " "; - } + for (int i = 1; i < narg; ++i) + atom_type_num_list += fmt::format(" {}", species_to_atomic_no(arg[i])); KIM_SimulatorModel_AddTemplateMap( - simulatorModel,"atom-type-sym-list",atom_type_sym_list.c_str()); + simulatorModel, "atom-type-sym-list", atom_type_sym_list.c_str()); KIM_SimulatorModel_AddTemplateMap( - simulatorModel,"atom-type-num-list",atom_type_num_list.c_str()); + simulatorModel, "atom-type-num-list", atom_type_num_list.c_str()); KIM_SimulatorModel_CloseTemplateMap(simulatorModel); // validate species selection @@ -150,14 +154,14 @@ void KimInteractions::do_setup(int narg, char **arg) for (auto atom_type_sym : utils::split_words(atom_type_sym_list)) { species_is_supported = false; if (atom_type_sym == "NULL") continue; - for (int i=0; i < sim_num_species; ++i) { - KIM_SimulatorModel_GetSupportedSpecies(simulatorModel,i,&sim_species); + for (int i = 0; i < sim_num_species; ++i) { + KIM_SimulatorModel_GetSupportedSpecies( + simulatorModel, i, &sim_species); if (atom_type_sym == sim_species) species_is_supported = true; } if (!species_is_supported) { - std::string msg = "Species '"; - msg += atom_type_sym + "' is not supported by this KIM Simulator Model"; - error->all(FLERR,msg); + error->all(FLERR,fmt::format("Species '{}' is not supported by this " + "KIM Simulator Model", atom_type_sym)); } } } else { @@ -169,29 +173,39 @@ void KimInteractions::do_setup(int narg, char **arg) int sim_fields, sim_lines; const char *sim_field, *sim_value; KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields); - for (int i=0; i < sim_fields; ++i) { + for (int i = 0; i < sim_fields; ++i) { KIM_SimulatorModel_GetSimulatorFieldMetadata( - simulatorModel,i,&sim_lines,&sim_field); + simulatorModel, i, &sim_lines, &sim_field); - if (0 == strcmp(sim_field,"units")) { - KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel,i,0,&sim_value); - if (0 != strcmp(sim_value,update->unit_style)) + if (strcmp(sim_field,"units") == 0) { + KIM_SimulatorModel_GetSimulatorFieldLine( + simulatorModel, i, 0, &sim_value); + if (strcmp(sim_value,update->unit_style) != 0) error->all(FLERR,"Incompatible units for KIM Simulator Model"); } } - int sim_model_idx=-1; - for (int i=0; i < sim_fields; ++i) { + bool no_model_definition = true; + for (int i = 0; i < sim_fields; ++i) { KIM_SimulatorModel_GetSimulatorFieldMetadata( - simulatorModel,i,&sim_lines,&sim_field); - if (0 == strcmp(sim_field,"model-defn")) { - if (domain->periodicity[0]&&domain->periodicity[1]&&domain->periodicity[2]) input->one("variable kim_periodic equal 1"); - else if (domain->periodicity[0]&&domain->periodicity[1]&&!domain->periodicity[2]) input->one("variable kim_periodic equal 2"); + simulatorModel, i, &sim_lines, &sim_field); + if (strcmp(sim_field,"model-defn") == 0) { + if (domain->periodicity[0]&& + domain->periodicity[1]&& + domain->periodicity[2]) + input->one("variable kim_periodic equal 1"); + else if (domain->periodicity[0]&& + domain->periodicity[1]&& + !domain->periodicity[2]) + input->one("variable kim_periodic equal 2"); else input->one("variable kim_periodic equal 0"); - sim_model_idx = i; - for (int j=0; j < sim_lines; ++j) { + + // KIM Simulator Model has a Model definition + no_model_definition = false; + + for (int j = 0; j < sim_lines; ++j) { KIM_SimulatorModel_GetSimulatorFieldLine( - simulatorModel,sim_model_idx,j,&sim_value); + simulatorModel, i, j, &sim_value); if (utils::strmatch(sim_value,"^KIM_SET_TYPE_PARAMETERS")) { // Notes regarding the KIM_SET_TYPE_PARAMETERS command // * This is an INTERNAL command. @@ -212,7 +226,7 @@ void KimInteractions::do_setup(int narg, char **arg) } } - if (sim_model_idx < 0) + if (no_model_definition) error->all(FLERR,"KIM Simulator Model has no Model definition"); KIM_SimulatorModel_OpenAndInitializeTemplateMap(simulatorModel); @@ -227,14 +241,9 @@ void KimInteractions::do_setup(int narg, char **arg) // NOTE: all references to arg must appear before calls to input->one() // as that will reset the argument vector. - std::string cmd1("pair_style kim "); - cmd1 += model_name; - - std::string cmd2("pair_coeff * * "); - for (int i=0; i < narg; ++i) { - cmd2 += arg[i]; - cmd2 += " "; - } + auto cmd1 = fmt::format("pair_style kim {}", model_name); + auto cmd2 = + fmt::format("pair_coeff * * {}", fmt::join(arg, arg + narg, " ")); input->one(cmd1); input->one(cmd2); @@ -248,22 +257,25 @@ void KimInteractions::do_setup(int narg, char **arg) void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) const { - int nocomment; auto words = utils::split_words(input_line); - std::string key = words[1]; + const std::string key = words[1]; + if (key != "pair" && key != "charge") + error->one(FLERR,fmt::format("Unrecognized KEY {} for " + "KIM_SET_TYPE_PARAMETERS command", key)); + std::string filename = words[2]; std::vector species(words.begin()+3,words.end()); if ((int)species.size() != atom->ntypes) error->one(FLERR,"Incorrect args for KIM_SET_TYPE_PARAMETERS command"); - FILE *fp; - fp = fopen(filename.c_str(),"r"); - if (fp == nullptr) { - error->one(FLERR,"Parameter file not found"); + FILE *fp = nullptr; + if (comm->me == 0) { + fp = fopen(filename.c_str(),"r"); + if (fp == nullptr) error->one(FLERR,"Parameter file not found"); } - char line[MAXLINE],*ptr; + char line[MAXLINE], *ptr; int n, eof = 0; while (1) { @@ -279,27 +291,24 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); - nocomment = line[0] != '#'; + if (ptr = strchr(line,'#')) *ptr = '\0'; + if (strspn(line," \t\n\r") == strlen(line)) continue; - if(nocomment) { - words = utils::split_words(line); - if (key == "pair") { - for (int ia = 0; ia < atom->ntypes; ++ia) { - for (int ib = ia; ib < atom->ntypes; ++ib) - if (((species[ia] == words[0]) && (species[ib] == words[1])) - || ((species[ib] == words[0]) && (species[ia] == words[1]))) - input->one(fmt::format("pair_coeff {} {} {}",ia+1,ib+1,fmt::join(words.begin()+2,words.end()," "))); - } - } else if (key == "charge") { - for (int ia = 0; ia < atom->ntypes; ++ia) - if (species[ia] == words[0]) - input->one(fmt::format("set type {} charge {}",ia+1,words[1])); - } else { - error->one(FLERR,"Unrecognized KEY for KIM_SET_TYPE_PARAMETERS command"); + words = utils::split_words(line); + if (key == "pair") { + for (int ia = 0; ia < atom->ntypes; ++ia) { + for (int ib = ia; ib < atom->ntypes; ++ib) + if (((species[ia] == words[0]) && (species[ib] == words[1])) + || ((species[ib] == words[0]) && (species[ia] == words[1]))) + input->one(fmt::format("pair_coeff {} {} {}",ia+1,ib+1, + fmt::join(words.begin()+2,words.end()," "))); } + } else { + for (int ia = 0; ia < atom->ntypes; ++ia) + if (species[ia] == words[0]) + input->one(fmt::format("set type {} charge {}",ia+1,words[1])); } } - fclose(fp); } /* ---------------------------------------------------------------------- */ diff --git a/src/KIM/kim_param.cpp b/src/KIM/kim_param.cpp index 900e13e09c..0927a467bf 100644 --- a/src/KIM/kim_param.cpp +++ b/src/KIM/kim_param.cpp @@ -57,6 +57,7 @@ #include "kim_param.h" +#include "comm.h" #include "error.h" #include "fix_store_kim.h" #include "force.h" @@ -65,8 +66,9 @@ #include "pair_kim.h" #include "variable.h" +#include #include -#include +#include extern "C" { @@ -118,8 +120,11 @@ void get_kim_unit_names( chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_fs; - } else if ((strcmp(system, "lj") == 0)) { - error->all(FLERR, "LAMMPS unit_style lj not supported by KIM models"); + } else if (strcmp(system,"lj") == 0 || + strcmp(system,"micro") ==0 || + strcmp(system,"nano")==0) { + error->all(FLERR,fmt::format("LAMMPS unit_style {} not supported " + "by KIM models", system)); } else error->all(FLERR, "Unknown unit_style"); } @@ -240,6 +245,9 @@ void KimParam::command(int narg, char **arg) &pkim); if (kim_error) error->all(FLERR, "Unable to create KIM Portable Model"); + + auto logID = fmt::format("{}_Model", comm->me); + KIM_Model_SetLogID(pkim, logID.c_str()); } } @@ -253,7 +261,7 @@ void KimParam::command(int narg, char **arg) // Parameter name char *paramname = nullptr; // Variable name - char *varname = nullptr; + std::string varname; // Loop over all the arguments for (int i = 1; i < narg;) { @@ -277,15 +285,13 @@ void KimParam::command(int narg, char **arg) if (kim_error) error->all(FLERR, "KIM GetParameterMetadata returned error"); - if (strcmp(paramname, str_name) == 0) - break; + if (strcmp(paramname, str_name) == 0) break; } if (param_index >= numberOfParameters) { - std::string msg("Wrong argument in kim_param get command.\n"); - msg += "This Model does not have the requested '"; - msg += paramname; - msg += "' parameter"; + auto msg = fmt::format("Wrong argument in kim_param get command.\n" + "This Model does not have the requested '{}' " + "parameter", paramname); error->all(FLERR, msg); } @@ -299,60 +305,53 @@ void KimParam::command(int narg, char **arg) // Check to see if the indices range contains // only integer numbers and/or range : if (argtostr.find_first_not_of("0123456789:") != std::string::npos) { - std::string msg("Illegal index_range.\n"); - msg += "Expected integer parameter(s) instead of '"; - msg += argtostr; - msg += "' in index_range"; + auto msg = fmt::format("Illegal index_range.\nExpected integer " + "parameter(s) instead of '{}' in " + "index_range", argtostr); error->all(FLERR, msg); } std::string::size_type npos = argtostr.find(':'); if (npos != std::string::npos) { argtostr[npos] = ' '; - std::stringstream str(argtostr); - str >> nlbound >> nubound; + auto words = utils::split_words(argtostr); + nlbound = atoi(words[0].c_str()); + nubound = atoi(words[1].c_str()); + if (nubound < 1 || nubound > extent || nlbound < 1 || nlbound > nubound) { - std::string msg("Illegal index_range '"); - msg += std::to_string(nlbound) + "-"; - msg += std::to_string(nubound) + "' for '"; - msg += paramname; - msg += "' parameter with the extent of '"; - msg += std::to_string(extent); - msg += "'"; + auto msg = fmt::format("Illegal index_range '{}-{}' for '{}' " + "parameter with the extent of '{}'", + nlbound, nubound, paramname, extent); error->all(FLERR, msg); } } else { - std::stringstream str(argtostr); - str >> nlbound; + nlbound = atoi(argtostr.c_str()); + if (nlbound < 1 || nlbound > extent) { - std::string msg("Illegal index '"); - msg += std::to_string(nlbound) + "' for '"; - msg += paramname; - msg += "' parameter with the extent of '"; - msg += std::to_string(extent); - msg += "'"; + auto msg = fmt::format("Illegal index '{}' for '{}' parameter " + "with the extent of '{}'", nlbound, + paramname, extent); error->all(FLERR, msg); } + nubound = nlbound; } } else { - std::string msg("Wrong number of arguments in "); - msg += "'kim_param get' command.\n"; - msg += "Index range after parameter name is mandatory"; + std::string msg("Wrong number of arguments in 'kim_param get' "); + msg += "command.\nIndex range after parameter name is mandatory"; error->all(FLERR, msg); } int const nvars = nubound - nlbound + 1; - char **varsname = nullptr; + std::vector varsname; if (i < narg) { // Get the variable/variable_base name varname = arg[i++]; } else { - std::string msg("Wrong number of arguments in "); - msg += "'kim_param get' command.\n"; - msg += "The LAMMPS variable name is mandatory"; + std::string msg("Wrong number of arguments in 'kim_param get' "); + msg += "command.\nThe LAMMPS variable name is mandatory"; error->all(FLERR, msg); } @@ -362,56 +361,47 @@ void KimParam::command(int narg, char **arg) if (nvars > 1) { if (i < narg) { if (strcmp(arg[i], "split") == 0) { - varsname = new char *[nvars]; + varsname.resize(nvars); for (int j = 0, k = nlbound; j < nvars; ++j, ++k) { - std::stringstream str; - str << varname << "_" << k; - varsname[j] = const_cast(str.str().c_str()); + varsname[j] = fmt::format("{}_{}", varname, k); } } else if (strcmp(arg[i], "list") == 0) { list_requested = true; - varsname = new char *[1]; + varsname.resize(1); varsname[0] = varname; // Default explicit (optional) formatarg } else if (i - 1 + nvars < narg) { - varsname = new char *[nvars]; + varsname.resize(nvars); --i; - for (int j = 0; j < nvars; ++j, ++i) - varsname[j] = arg[i]; + for (int j = 0; j < nvars; ++j, ++i) varsname[j] = arg[i]; if (i < narg) { - if (strcmp(arg[i], "explicit") == 0) - ++i; + if (strcmp(arg[i], "explicit") == 0) ++i; } } else { - std::string msg("Wrong number of arguments in "); - msg += "'kim_param get' command.\nThe LAMMPS '"; - msg += std::to_string(nvars); - msg += "' variable names or '"; - msg += varname; - msg += " split' is mandatory"; + auto msg = + fmt::format("Wrong number of arguments in 'kim_param get' " + "command.\nThe LAMMPS '{}' variable names or " + "'{} split' is mandatory", nvars, varname); error->all(FLERR, msg); } } else { - std::string msg("Wrong number of arguments in "); - msg += "'kim_param get' command.\nThe LAMMPS '"; - msg += std::to_string(nvars); - msg += "' variable names or '"; - msg += varname; - msg += " split/list' is mandatory"; + auto msg = + fmt::format("Wrong number of arguments in 'kim_param get' " + "command.\nThe LAMMPS '{}' variable names or " + "'{} split/list' is mandatory", nvars, varname); error->all(FLERR, msg); } } else { - varsname = new char *[1]; + varsname.resize(1); if (i < narg) { if (strcmp(arg[i], "split") == 0) { - std::stringstream str; - str << varname << "_" << nlbound; - varsname[0] = const_cast(str.str().c_str()); + varsname[0] = fmt::format("{}_{}", varname, nlbound); ++i; } else { if ((strcmp(arg[i], "list") == 0) || - (strcmp(arg[i], "explicit") == 0)) + (strcmp(arg[i], "explicit") == 0)) ++i; + varsname[0] = varname; } } else { @@ -419,98 +409,90 @@ void KimParam::command(int narg, char **arg) } } - char **varcmd = new char *[3]; - varcmd[1] = const_cast("string"); - if (KIM_DataType_Equal(kim_DataType, KIM_DATA_TYPE_Double)) { if (list_requested) { - std::stringstream str; + std::string str; double V; { kim_error = KIM_Model_GetParameterDouble(pkim, param_index, nlbound - 1, &V); if (kim_error) error->all(FLERR, "KIM GetParameterDouble returned error"); - str << V; + + str = fmt::format("{}", V); } for (int j = 1; j < nvars; ++j) { kim_error = KIM_Model_GetParameterDouble(pkim, param_index, nlbound - 1 + j, &V); if (kim_error) error->all(FLERR, "KIM GetParameterDouble returned error"); - str << " " << V; + + str += fmt::format(" {}", V); } - varcmd[0] = varsname[0]; - varcmd[2] = const_cast(str.str().c_str()); - input->variable->set(3, varcmd); - echo_var_assign(varcmd[0], varcmd[2]); + + auto setcmd = fmt::format("{} string {}", varsname[0], str); + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); + } else { + double V; for (int j = 0; j < nvars; ++j) { - varcmd[0] = varsname[j]; - double V; kim_error = KIM_Model_GetParameterDouble(pkim, param_index, nlbound - 1 + j, &V); if (kim_error) error->all(FLERR, "KIM GetParameterDouble returned error"); - std::stringstream str; - str << V; - varcmd[2] = const_cast(str.str().c_str()); - input->variable->set(3, varcmd); - echo_var_assign(varcmd[0], varcmd[2]); + + auto setcmd = fmt::format("{} string {}", varsname[j], V); + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); } } } else if (KIM_DataType_Equal(kim_DataType, KIM_DATA_TYPE_Integer)) { if (list_requested) { - std::stringstream str; + std::string str; int V; { kim_error = KIM_Model_GetParameterInteger(pkim, param_index, nlbound - 1, &V); if (kim_error) error->all(FLERR, "KIM GetParameterInteger returned error"); - str << V; + + str = fmt::format("{}", V); } for (int j = 1; j < nvars; ++j) { kim_error = KIM_Model_GetParameterInteger(pkim, param_index, nlbound - 1 + j, &V); if (kim_error) error->all(FLERR, "KIM GetParameterInteger returned error"); - str << " " << V; + + str += fmt::format(" {}", V); } - varcmd[0] = varsname[0]; - varcmd[2] = const_cast(str.str().c_str()); - input->variable->set(3, varcmd); - echo_var_assign(varcmd[0], varcmd[2]); + + auto setcmd = fmt::format("{} string {}", varsname[0], str); + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); + } else { + int V; for (int j = 0; j < nvars; ++j) { - varcmd[0] = varsname[j]; - int V; kim_error = KIM_Model_GetParameterInteger(pkim, param_index, nlbound - 1 + j, &V); if (kim_error) error->all(FLERR, "KIM GetParameterInteger returned error"); - std::stringstream str; - str << V; - varcmd[2] = const_cast(str.str().c_str()); - input->variable->set(3, varcmd); - echo_var_assign(varcmd[0], varcmd[2]); + + auto setcmd = fmt::format("{} string {}", varsname[j], V); + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); } } } else error->all(FLERR, "Wrong parameter type"); - - delete[] varcmd; - delete[] varsname; } // End of loop over all the arguments // Set the parameters } else { - std::string set_cmd("pair_coeff * * "); - set_cmd += atom_type_list; - for (int i = 1; i < narg; ++i) { - set_cmd += " "; - set_cmd += arg[i]; - } - input->one(set_cmd); + auto setcmd = fmt::format("pair_coeff * * {} {}", atom_type_list, + fmt::join(arg + 1, arg + narg, " ")); + input->one(setcmd); } } else error->all(FLERR, "This model has No mutable parameters"); @@ -521,12 +503,3 @@ void KimParam::command(int narg, char **arg) input->write_echo(fmt::format("#=== END kim-param {} =====================" "==================\n",kim_param_get_set)); } - -/* ---------------------------------------------------------------------- */ - -void KimParam::echo_var_assign(const std::string &name, - const std::string &value) const -{ - input->write_echo(fmt::format("variable {} string {}\n", - name, value)); -} diff --git a/src/KIM/kim_param.h b/src/KIM/kim_param.h index e4fc5ce59a..3e20207cca 100644 --- a/src/KIM/kim_param.h +++ b/src/KIM/kim_param.h @@ -77,9 +77,6 @@ public: ~KimParam(); void command(int, char **); - -private: - void echo_var_assign(const std::string &name, const std::string &value) const; }; } // namespace LAMMPS_NS diff --git a/src/KIM/kim_property.cpp b/src/KIM/kim_property.cpp index 88fa16e60a..17d8778c7a 100644 --- a/src/KIM/kim_property.cpp +++ b/src/KIM/kim_property.cpp @@ -100,11 +100,7 @@ void kimProperty::command(int narg, char **arg) error->all(FLERR, msg); } - if (comm->me == 0) { - std::string msg; - msg = "#=== kim-property ===========================================\n"; - input->write_echo(msg); - } + input->write_echo("#=== kim-property ===========================================\n"); // Get the kim_str ptr to the data associated with a kim_property_str // variable diff --git a/src/KIM/kim_query.cpp b/src/KIM/kim_query.cpp index 0b0d5e97ad..e7b308625b 100644 --- a/src/KIM/kim_query.cpp +++ b/src/KIM/kim_query.cpp @@ -63,11 +63,15 @@ #include "info.h" #include "input.h" #include "modify.h" +#include "utils.h" #include "variable.h" #include "version.h" +#include "tokenizer.h" +#include "fmt/format.h" +#include #include -#include +#include #if defined(LMP_KIM_CURL) #include @@ -92,26 +96,28 @@ static size_t write_callback(void *, size_t, size_t, void *); void KimQuery::command(int narg, char **arg) { - char *varname, *function, *value; - if (narg < 2) error->all(FLERR,"Illegal kim_query command"); // check if we had a kim_init command by finding fix STORE/KIM // retrieve model name. - char * model_name; + char *model_name; - int ifix = modify->find_fix("KIM_MODEL_STORE"); + const int ifix = modify->find_fix("KIM_MODEL_STORE"); if (ifix >= 0) { FixStoreKIM *fix_store = (FixStoreKIM *) modify->fix[ifix]; model_name = (char *)fix_store->getptr("model_name"); } else error->all(FLERR,"Must use 'kim_init' before 'kim_query'"); - varname = arg[0]; + char *varname = arg[0]; + bool split = false; - if (0 == strcmp("split",arg[1])) { - if (narg == 2) error->all(FLERR,"Illegal kim_query command"); - if (0 == strcmp("list",arg[2])) - error->all(FLERR,"Illegal kim_query command"); + if (strcmp("split",arg[1]) == 0) { + if (narg == 2) error->all(FLERR,"Illegal kim_query command.\nThe keyword " + "'split' must be followed by the name of " + "the query function"); + if (strcmp("list",arg[2]) == 0) + error->all(FLERR,"Illegal kim_query command.\nThe 'list' keyword " + "can not be used after 'split'"); split = true; arg++; narg--; @@ -119,78 +125,69 @@ void KimQuery::command(int narg, char **arg) // The “list” is the default setting // the result is returned as a space-separated list of values in variable - if (0 == strcmp("list",arg[1])) { - if (narg == 2) error->all(FLERR,"Illegal kim_query command"); - if (split) error->all(FLERR,"Illegal kim_query command"); + if (strcmp("list",arg[1]) == 0) { + if (narg == 2) error->all(FLERR,"Illegal kim_query command.\nThe 'list' " + "keyword must be followed by ('split' " + "and) the name of the query function"); arg++; narg--; } - function = arg[1]; - for (int i = 2; i < narg; ++i) { - if (0 == strncmp("model=",arg[i], 6)) { - error->all(FLERR,"Illegal 'model' key in kim_query command"); - } - } + char *function = arg[1]; + for (int i = 2; i < narg; ++i) { + if (strncmp("model=",arg[i],6) == 0) + error->all(FLERR,"Illegal 'model' key in kim_query command"); + + if (!strchr(arg[i], '=') || !strchr(arg[i], '[') || !strchr(arg[i], ']')) + error->all(FLERR,fmt::format("Illegal query format.\nInput argument of " + "`{}` to kim_query is wrong. The query " + "format is the keyword=[value], where value " + "is always an array of one or more " + "comma-separated items", arg[i])); + } #if defined(LMP_KIM_CURL) - value = do_query(function, model_name, narg-2, arg+2, comm->me, world); + char *value = do_query(function, model_name, narg-2, arg+2, comm->me, world); // check for valid result // on error the content of "value" is a '\0' byte // as the first element, and then the error message // that was returned by the web server - char errmsg[1024]; - if (0 == strlen(value)) { - sprintf(errmsg,"OpenKIM query failed: %s",value+1); - error->all(FLERR,errmsg); - } else if (0 == strcmp(value,"EMPTY")) { - sprintf(errmsg,"OpenKIM query returned no results"); - error->all(FLERR,errmsg); + if (strlen(value) == 0) { + error->all(FLERR,fmt::format("OpenKIM query failed: {}", value+1)); + } else if (strcmp(value,"EMPTY") == 0) { + error->all(FLERR,fmt::format("OpenKIM query returned no results")); } input->write_echo("#=== BEGIN kim-query =========================================\n"); - char **varcmd = new char*[3]; - varcmd[1] = (char *) "string"; - - std::stringstream ss(value); - std::string token; - + ValueTokenizer values(value, ","); if (split) { int counter = 1; - while(std::getline(ss, token, ',')) { - token.erase(0,token.find_first_not_of(" \n\r\t")); // ltrim - token.erase(token.find_last_not_of(" \n\r\t") + 1); // rtrim - std::stringstream splitname; - splitname << varname << "_" << counter++; - varcmd[0] = const_cast(splitname.str().c_str()); - varcmd[2] = const_cast(token.c_str()); - input->variable->set(3,varcmd); - echo_var_assign(splitname.str(), varcmd[2]); + while (values.has_next()) { + auto svalue = values.next_string(); + auto setcmd = fmt::format("{}_{} string {}", varname, counter++, svalue); + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); } } else { - varcmd[0] = varname; - std::string value_string; - while(std::getline(ss, token, ',')) { - token.erase(0,token.find_first_not_of(" \n\r\t")); // ltrim - token.erase(token.find_last_not_of(" \n\r\t") + 1); // rtrim - if (value_string.size() && token.size()) - value_string += " "; - value_string += token; + auto svalue = values.next_string(); + std::string setcmd = fmt::format("{} string \"{}", varname, svalue); + while (values.has_next()) { + svalue = values.next_string(); + setcmd += fmt::format(" {}", svalue); } - varcmd[2] = const_cast(value_string.c_str());; - input->variable->set(3,varcmd); - echo_var_assign(varname, value_string); + setcmd += "\""; + input->variable->set(setcmd); + input->write_echo(fmt::format("variable {}\n", setcmd)); } input->write_echo("#=== END kim-query ===========================================\n\n"); - delete[] varcmd; delete[] value; #else error->all(FLERR,"Cannot use 'kim_query' command when KIM package " - "is compiled without support for libcurl"); + "is compiled without support for libcurl"); #endif } @@ -201,17 +198,18 @@ void KimQuery::command(int narg, char **arg) size_t write_callback(void *data, size_t size, size_t nmemb, void *userp) { struct WriteBuf *buf = (struct WriteBuf *)userp; - size_t buffer_size = size*nmemb; // copy chunks into the buffer for as long as there is space left if (buf->sizeleft) { - size_t copy_this_much = buf->sizeleft; - if (copy_this_much > buffer_size) - copy_this_much = buffer_size; + const size_t buffer_size = size * nmemb; + const size_t copy_this_much = + buf->sizeleft > buffer_size ? buffer_size : buf->sizeleft; + memcpy(buf->dataptr, data, copy_this_much); buf->dataptr += copy_this_much; buf->sizeleft -= copy_this_much; + return copy_this_much; } return 0; // done @@ -220,16 +218,12 @@ size_t write_callback(void *data, size_t size, size_t nmemb, void *userp) char *do_query(char *qfunction, char * model_name, int narg, char **arg, int rank, MPI_Comm comm) { - char value[512], *retval; + char value[512]; // run the web query from rank 0 only if (rank == 0) { - CURL *handle; - CURLcode res; - // set up and clear receive buffer - struct WriteBuf buf; buf.dataptr = value; buf.sizeleft = 511; @@ -237,19 +231,43 @@ char *do_query(char *qfunction, char * model_name, int narg, char **arg, // create curl web query instance curl_global_init(CURL_GLOBAL_DEFAULT); - handle = curl_easy_init(); + CURL *handle = curl_easy_init(); if (handle) { - std::string url("https://query.openkim.org/api/"); - url += qfunction; - - std::string query(arg[0]); - query += "&model=[\""; - query += model_name; - query += "\"]"; - for (int i=1; i < narg; ++i) { - query += '&'; - query += arg[i]; + auto url = fmt::format("https://query.openkim.org/api/{}", qfunction); + auto query = fmt::format("model=[\"{}\"]", model_name); + for (int i = 0; i < narg; ++i) { + ValueTokenizer values(arg[i], "=[]"); + std::string key = values.next_string(); + std::string val = values.next_string(); + std::string::size_type n = val.find(","); + if (n == std::string::npos) { + if (utils::is_integer(val) || + utils::is_double(val) || + (val.front() == '"' && + val.back() == '"')) { + query += fmt::format("&{}", arg[i]); + } else { + query += fmt::format("&{}=[\"{}\"]", key, val); + } + } else { + query += fmt::format("&{}=[", key); + while (n != std::string::npos){ + std::string sval = val.substr(0, n); + if (utils::is_integer(sval) || + utils::is_double(sval) || + (val.front() == '"' && + val.back() == '"')) { + query += fmt::format("{},", sval); + } else { + query += fmt::format("\"{}\",", sval); + } + val = val.substr(n + 1); + n = val.find(","); + } + if (val.size()) query += fmt::format("\"{}\"]", val); + else query[query.size() - 1]=']'; + } } #if LMP_DEBUG_CURL @@ -274,23 +292,21 @@ char *do_query(char *qfunction, char * model_name, int narg, char **arg, } } - std::string user_agent = std::string("kim_query--LAMMPS/") - + LAMMPS_VERSION - + " (" + Info::get_os_info() + ")"; + std::string user_agent = fmt::format("kim_query--LAMMPS/{} ({})", + LAMMPS_VERSION, Info::get_os_info()); + curl_easy_setopt(handle, CURLOPT_USERAGENT, user_agent.c_str()); - curl_easy_setopt(handle, CURLOPT_URL, url.c_str()); curl_easy_setopt(handle, CURLOPT_FOLLOWLOCATION, 1L); curl_easy_setopt(handle, CURLOPT_POSTFIELDS, query.c_str()); - - curl_easy_setopt(handle, CURLOPT_WRITEFUNCTION,write_callback); + curl_easy_setopt(handle, CURLOPT_WRITEFUNCTION, write_callback); curl_easy_setopt(handle, CURLOPT_WRITEDATA,&buf); // perform OpenKIM query and check for errors - res = curl_easy_perform(handle); + CURLcode res = curl_easy_perform(handle); if (res != CURLE_OK) { // on error we return an "empty" string but add error message after it - value[0]= '\0'; + value[0] = '\0'; strcpy(value+1,curl_easy_strerror(res)); } curl_easy_cleanup(handle); @@ -302,30 +318,31 @@ char *do_query(char *qfunction, char * model_name, int narg, char **arg, // we must make a proper copy of the query, as the stack allocation // for "value" will go out of scope. a valid query has a '[' as // the first character. skip over it (and the last character ']', too) - // an error messages starts with a '\0' character. copy that and + // an error message starts with a '\0' character. copy that and // the remaining string, as that is the error message. + char *retval; + // a valid query has a '[' as the first character if (value[0] == '[') { - int len = strlen(value)-1; + int len = strlen(value) - 1; if (value[len] == ']') { - retval = new char[len]; value[len] = '\0'; - if (0 == strcmp(value+1, "")) { - strcpy(retval,"EMPTY"); - } else - strcpy(retval,value+1); + retval = new char[len]; + if (strcmp(value+1, "") == 0) strcpy(retval,"EMPTY"); + else strcpy(retval,value+1); } else { retval = new char[len+2]; retval[0] = '\0'; strcpy(retval+1,value); } + // an error message starts with a '\0' character } else if (value[0] == '\0') { int len = strlen(value+1)+2; retval = new char[len]; retval[0] = '\0'; strcpy(retval+1,value+1); + // unknown response type. we should not get here. } else { - // unknown response type. we should not get here. // we return an "empty" string but add error message after it int len = strlen(value)+2; retval = new char[len]; @@ -335,11 +352,3 @@ char *do_query(char *qfunction, char * model_name, int narg, char **arg, return retval; } #endif - -/* ---------------------------------------------------------------------- */ - -void KimQuery::echo_var_assign(const std::string &name, - const std::string &value) const -{ - input->write_echo(fmt::format("variable {} string {}\n",name,value)); -} diff --git a/src/KIM/kim_query.h b/src/KIM/kim_query.h index e1c49ed71d..f2523f5a98 100644 --- a/src/KIM/kim_query.h +++ b/src/KIM/kim_query.h @@ -13,7 +13,8 @@ /* ---------------------------------------------------------------------- Contributing authors: Axel Kohlmeyer (Temple U), - Ryan S. Elliott (UMN) + Ryan S. Elliott (UMN), + Yaser Afshar (UMN) ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- @@ -71,9 +72,6 @@ class KimQuery : protected Pointers { public: KimQuery(class LAMMPS *lmp) : Pointers(lmp) {}; void command(int, char **); - private: - void echo_var_assign(const std::string &name, const std::string &value) - const; }; } diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp index 2a5fafefef..7ff067dcac 100644 --- a/src/KIM/pair_kim.cpp +++ b/src/KIM/pair_kim.cpp @@ -67,8 +67,9 @@ #include "neighbor.h" #include "update.h" +#include #include -#include +#include using namespace LAMMPS_NS; @@ -204,8 +205,8 @@ void PairKIM::compute(int eflag, int vflag) KIM_COMPUTE_ARGUMENT_NAME_particleContributing, kim_particleContributing); if (kimerror) - error->all(FLERR, - "Unable to set KIM particle species codes and/or contributing"); + error->all(FLERR,"Unable to set KIM particle species " + "codes and/or contributing"); } // kim_particleSpecies = KIM atom species for each LAMMPS atom @@ -339,7 +340,9 @@ void PairKIM::coeff(int narg, char **arg) // insure I,J args are * * if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) - error->all(FLERR,"Incorrect args for pair coefficients"); + error->all(FLERR,"Incorrect args for pair coefficients.\nThe first two " + "arguments of pair_coeff command must be * * to span " + "all LAMMPS atom types"); int ilo,ihi,jlo,jhi; utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); @@ -410,9 +413,8 @@ void PairKIM::coeff(int narg, char **arg) if (supported) { kim_particle_codes[i] = code; } else { - std::string msg("create_kim_particle_codes: symbol not found: "); - msg += lmps_unique_elements[i]; - error->all(FLERR, msg); + error->all(FLERR,fmt::format("GetSpeciesSupportAndCode: symbol not " + "found: {}",lmps_unique_elements[i])); } } // Set the new values for PM parameters @@ -421,11 +423,9 @@ void PairKIM::coeff(int narg, char **arg) int numberOfParameters(0); KIM_Model_GetNumberOfParameters(pkim, &numberOfParameters); - if (!numberOfParameters) { - std::string msg("Incorrect args for pair coefficients \n"); - msg += "This model has No mutable parameters"; - error->all(FLERR, msg); - } + if (!numberOfParameters) + error->all(FLERR,"Incorrect args for pair coefficients\n" + "This model has No mutable parameters"); int kimerror; @@ -456,11 +456,10 @@ void PairKIM::coeff(int narg, char **arg) } if (param_index >= numberOfParameters) { - std::string msg("Wrong argument for pair coefficients.\n"); - msg += "This Model does not have the requested '"; - msg += paramname; - msg += "' parameter"; - error->all(FLERR, msg); + auto msg = fmt::format("Wrong argument for pair coefficients.\n" + "This Model does not have the requested " + "'{}' parameter", paramname); + error->all(FLERR,msg); } // Get the index_range for the requested parameter @@ -472,48 +471,41 @@ void PairKIM::coeff(int narg, char **arg) // Check to see if the indices range contains only integer numbers & : if (argtostr.find_first_not_of("0123456789:") != std::string::npos) { - std::string msg("Illegal index_range.\n"); - msg += "Expected integer parameter(s) instead of '"; - msg += argtostr; - msg += "' in index_range"; - error->all(FLERR, msg); + auto msg = fmt::format("Illegal index_range.\nExpected integer " + "parameter(s) instead of '{}' in " + "index_range", argtostr); + error->all(FLERR,msg); } std::string::size_type npos = argtostr.find(':'); if (npos != std::string::npos) { argtostr[npos] = ' '; - std::stringstream str(argtostr); - str >> nlbound >> nubound; + auto words = utils::split_words(argtostr); + nlbound = atoi(words[0].c_str()); + nubound = atoi(words[1].c_str()); + if (nubound < 1 || nubound > extent || nlbound < 1 || nlbound > nubound) { - std::string msg("Illegal index_range '"); - msg += std::to_string(nlbound) + "-" ; - msg += std::to_string(nubound) + "' for '"; - msg += paramname; - msg += "' parameter with the extent of '"; - msg += std::to_string(extent); - msg += "'"; - error->all(FLERR, msg); + auto msg = fmt::format("Illegal index_range '{}-{}' for '{}' " + "parameter with the extent of '{}'", + nlbound, nubound, paramname, extent); + error->all(FLERR,msg); } } else { - std::stringstream str(argtostr); - str >> nlbound; + nlbound = atoi(argtostr.c_str()); + if (nlbound < 1 || nlbound > extent) { - std::string msg("Illegal index '"); - msg += std::to_string(nlbound) + "' for '"; - msg += paramname; - msg += "' parameter with the extent of '"; - msg += std::to_string(extent); - msg += "'"; - error->all(FLERR, msg); + auto msg = fmt::format("Illegal index '{}' for '{}' parameter " + "with the extent of '{}'", nlbound, + paramname, extent); + error->all(FLERR,msg); } + nubound = nlbound; } } else { - std::string msg = - "Wrong number of arguments for pair coefficients.\n"; - msg += "Index range after parameter name is mandatory"; - error->all(FLERR, msg); + error->all(FLERR,"Wrong number of arguments for pair coefficients.\n" + "Index range after parameter name is mandatory"); } // Parameter values @@ -524,7 +516,7 @@ void PairKIM::coeff(int narg, char **arg) kimerror = KIM_Model_SetParameterDouble(pkim, param_index, nlbound - 1 + j, V); if (kimerror) - error->all(FLERR, "KIM SetParameterDouble returned error"); + error->all(FLERR,"KIM SetParameterDouble returned error"); } } else if (KIM_DataType_Equal(kim_DataType, KIM_DATA_TYPE_Integer)) { for (int j = 0; j < nubound - nlbound + 1; ++j) { @@ -532,25 +524,22 @@ void PairKIM::coeff(int narg, char **arg) kimerror = KIM_Model_SetParameterInteger(pkim, param_index, nlbound - 1 + j, V); if (kimerror) - error->all(FLERR, "KIM SetParameterInteger returned error"); + error->all(FLERR,"KIM SetParameterInteger returned error"); } } else - error->all(FLERR, "Wrong parameter type to update"); + error->all(FLERR,"Wrong parameter type to update"); } else { - std::string msg = - "Wrong number of variable values for pair coefficients.\n"; - msg += "'"; - msg += std::to_string(nubound - nlbound + 1); - msg += "' values are requested for '"; - msg += paramname; - msg += "' parameter."; - error->all(FLERR, msg); + auto msg = fmt::format("Wrong number of variable values for pair " + "coefficients.\n'{}' values are requested " + "for '{}' parameter", nubound - nlbound + 1, + paramname); + error->all(FLERR,msg); } } kimerror = KIM_Model_ClearThenRefresh(pkim); if (kimerror) - error->all(FLERR, "KIM KIM_Model_ClearThenRefresh returned error"); + error->all(FLERR,"KIM KIM_Model_ClearThenRefresh returned error"); } } @@ -842,11 +831,13 @@ void PairKIM::kim_init() else if (!requestedUnitsAccepted) error->all(FLERR,"KIM Model did not accept the requested unit system"); + auto logID = fmt::format("{}_Model", comm->me); + KIM_Model_SetLogID(pkim, logID.c_str()); + // check that the model does not require unknown capabilities kimerror = check_for_routine_compatibility(); if (kimerror) - error->all(FLERR, - "KIM Model requires unknown Routines. Unable to proceed"); + error->all(FLERR,"KIM Model requires unknown Routines. Unable to proceed"); kimerror = KIM_Model_ComputeArgumentsCreate(pkim, &pargs); if (kimerror) { @@ -854,6 +845,9 @@ void PairKIM::kim_init() error->all(FLERR,"KIM ComputeArgumentsCreate failed"); } else kim_init_ok = true; + logID = fmt::format("{}_ComputeArguments", comm->me); + KIM_ComputeArguments_SetLogID(pargs, logID.c_str()); + // determine KIM Model capabilities (used in this function below) set_kim_model_has_flags(); @@ -985,43 +979,46 @@ void PairKIM::set_lmps_flags() error->all(FLERR,"pair_kim does not support hybrid"); // determine unit system and set lmps_units flag - if ((strcmp(update->unit_style,"real")==0)) { + if (strcmp(update->unit_style,"real") == 0) { lmps_units = REAL; lengthUnit = KIM_LENGTH_UNIT_A; energyUnit = KIM_ENERGY_UNIT_kcal_mol; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_fs; - } else if ((strcmp(update->unit_style,"metal")==0)) { + } else if (strcmp(update->unit_style,"metal") == 0) { lmps_units = METAL; lengthUnit = KIM_LENGTH_UNIT_A; energyUnit = KIM_ENERGY_UNIT_eV; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_ps; - } else if ((strcmp(update->unit_style,"si")==0)) { + } else if (strcmp(update->unit_style,"si") == 0) { lmps_units = SI; lengthUnit = KIM_LENGTH_UNIT_m; energyUnit = KIM_ENERGY_UNIT_J; chargeUnit = KIM_CHARGE_UNIT_C; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_s; - } else if ((strcmp(update->unit_style,"cgs")==0)) { + } else if (strcmp(update->unit_style,"cgs") == 0) { lmps_units = CGS; lengthUnit = KIM_LENGTH_UNIT_cm; energyUnit = KIM_ENERGY_UNIT_erg; chargeUnit = KIM_CHARGE_UNIT_statC; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_s; - } else if ((strcmp(update->unit_style,"electron")==0)) { + } else if (strcmp(update->unit_style,"electron") == 0) { lmps_units = ELECTRON; lengthUnit = KIM_LENGTH_UNIT_Bohr; energyUnit = KIM_ENERGY_UNIT_Hartree; chargeUnit = KIM_CHARGE_UNIT_e; temperatureUnit = KIM_TEMPERATURE_UNIT_K; timeUnit = KIM_TIME_UNIT_fs; - } else if ((strcmp(update->unit_style,"lj")==0)) { - error->all(FLERR,"LAMMPS unit_style lj not supported by KIM models"); + } else if (strcmp(update->unit_style,"lj") == 0 || + strcmp(update->unit_style,"micro") == 0 || + strcmp(update->unit_style,"nano") == 0) { + error->all(FLERR,fmt::format("LAMMPS unit_style {} not supported " + "by KIM models", update->unit_style)); } else { error->all(FLERR,"Unknown unit_style"); } @@ -1102,29 +1099,33 @@ void PairKIM::set_kim_model_has_flags() KIM_SUPPORT_STATUS_required)) { std::string msg("KIM Model requires unsupported compute argument: "); msg += KIM_ComputeArgumentName_ToString(computeArgumentName); - error->all(FLERR, msg); + error->all(FLERR,msg); } } - if (KIM_SupportStatus_Equal(kim_model_support_for_energy, - KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide `partialEnergy'; " - "Potential energy will be zero"); + if (comm->me == 0) { + if (KIM_SupportStatus_Equal(kim_model_support_for_energy, + KIM_SUPPORT_STATUS_notSupported)) + error->warning(FLERR,"KIM Model does not provide 'partialEnergy'; " + "Potential energy will be zero"); - if (KIM_SupportStatus_Equal(kim_model_support_for_forces, - KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide `partialForce'; " - "Forces will be zero"); + if (KIM_SupportStatus_Equal(kim_model_support_for_forces, + KIM_SUPPORT_STATUS_notSupported)) + error->warning(FLERR,"KIM Model does not provide 'partialForce'; " + "Forces will be zero"); - if (KIM_SupportStatus_Equal(kim_model_support_for_particleEnergy, - KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide `partialParticleEnergy'; " - "energy per atom will be zero"); + if (KIM_SupportStatus_Equal(kim_model_support_for_particleEnergy, + KIM_SUPPORT_STATUS_notSupported)) + error->warning(FLERR,"KIM Model does not provide " + "'partialParticleEnergy'; " + "energy per atom will be zero"); - if (KIM_SupportStatus_Equal(kim_model_support_for_particleVirial, - KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide `partialParticleVirial'; " - "virial per atom will be zero"); + if (KIM_SupportStatus_Equal(kim_model_support_for_particleVirial, + KIM_SUPPORT_STATUS_notSupported)) + error->warning(FLERR,"KIM Model does not provide " + "'partialParticleVirial'; " + "virial per atom will be zero"); + } int numberOfComputeCallbackNames; KIM_COMPUTE_CALLBACK_NAME_GetNumberOfComputeCallbackNames( diff --git a/src/KOKKOS/kokkos.cpp b/src/KOKKOS/kokkos.cpp index fe72981d50..64455fef4f 100644 --- a/src/KOKKOS/kokkos.cpp +++ b/src/KOKKOS/kokkos.cpp @@ -199,7 +199,6 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) #endif neigh_thread = 0; neigh_thread_set = 0; - neighflag_qeq_set = 0; if (ngpus > 0) { neighflag = FULL; neighflag_qeq = FULL; @@ -314,7 +313,6 @@ void KokkosLMP::accelerator(int narg, char **arg) neighflag = HALF; } else error->all(FLERR,"Illegal package kokkos command"); - if (!neighflag_qeq_set) neighflag_qeq = neighflag; iarg += 2; } else if (strcmp(arg[iarg],"neigh/qeq") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command"); @@ -326,7 +324,6 @@ void KokkosLMP::accelerator(int narg, char **arg) neighflag_qeq = HALF; } else error->all(FLERR,"Illegal package kokkos command"); - neighflag_qeq_set = 1; iarg += 2; } else if (strcmp(arg[iarg],"binsize") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command"); diff --git a/src/KOKKOS/kokkos_type.h b/src/KOKKOS/kokkos_type.h index 5930a9e207..580b22d35f 100644 --- a/src/KOKKOS/kokkos_type.h +++ b/src/KOKKOS/kokkos_type.h @@ -405,7 +405,7 @@ struct s_EV_FLOAT_REAX { E_FLOAT evdwl; E_FLOAT ecoul; E_FLOAT v[6]; - E_FLOAT ereax[10]; + E_FLOAT ereax[9]; KOKKOS_INLINE_FUNCTION s_EV_FLOAT_REAX() { evdwl = 0; diff --git a/src/USER-MISC/README b/src/USER-MISC/README index f98268762f..538fdb7952 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -22,10 +22,12 @@ angle_style cosine/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 angle_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 angle_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 angle_style fourier/simple, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 +angle_style gaussian, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 Nov 2020 angle_style dipole, Mario Orsi, orsimario at gmail.com, 10 Jan 12 angle_style quartic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 +bond_style gaussian, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 Nov 2020 bond_style special, David Nicholson, davidanich at gmail.com, 31 Jan 2020 compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007 compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013 @@ -51,6 +53,7 @@ dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18 fix accelerate/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015 +fix electron/stopping/fit, James Stewart (SNL), jstewa .at. sandia.gov, 23 Sep 2020 fix electron/stopping, Konstantin Avchaciov, k.avchachov at gmail.com, 26 Feb 2019 fix ffl, David Wilkins (EPFL Lausanne), david.wilkins @ epfl.ch, 28 Sep 2018 fix filter/corotate, Lukas Fath (KIT), lukas.fath at kit.edu, 15 Mar 2017 diff --git a/src/USER-MISC/angle_gaussian.cpp b/src/USER-MISC/angle_gaussian.cpp new file mode 100644 index 0000000000..640be38433 --- /dev/null +++ b/src/USER-MISC/angle_gaussian.cpp @@ -0,0 +1,345 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "angle_gaussian.h" + +#include +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMAL 0.001 +#define SMALL 1.0e-8 + +/* ---------------------------------------------------------------------- */ + +AngleGaussian::AngleGaussian(LAMMPS *lmp) : Angle(lmp) +{ +} + +/* ---------------------------------------------------------------------- */ + +AngleGaussian::~AngleGaussian() +{ + if (allocated && !copymode) { + memory->destroy(setflag); + memory->destroy(nterms); + memory->destroy(angle_temperature); + for (int i = 1; i <= atom->nangletypes; i++) { + if (alpha[i]) delete [] alpha[i]; + if (width[i]) delete [] width[i]; + if (theta0[i]) delete [] theta0[i]; + } + delete [] alpha; + delete [] width; + delete [] theta0; + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleGaussian::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double dtheta,tk; + double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22; + double prefactor, exponent, g_i, sum_g_i, sum_numerator; + + eangle = 0.0; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + s = sqrt(1.0 - c*c); + if (s < SMAL) s = SMAL; + s = 1.0/s; + + // force & energy + double theta = acos(c); + + sum_g_i = 0.0; + sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + dtheta = theta - theta0[type][i]; + prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2))); + exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]); + g_i = prefactor*exp(exponent); + sum_g_i += g_i; + sum_numerator += g_i*dtheta/(width[type][i]*width[type][i]); + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + if (eflag) eangle = -(force->boltz*angle_temperature[type])*log(sum_g_i); + + // I should check about the sign of this expression + a = -4.0*(force->boltz*angle_temperature[type])*(sum_numerator/sum_g_i)*s; + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + + // apply force to each of 3 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleGaussian::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(nterms,n+1,"angle:nterms"); + memory->create(angle_temperature,n+1,"angle:angle_temperature"); + + alpha = new double *[n+1]; + width = new double *[n+1]; + theta0 = new double *[n+1]; + for (int i = 1; i <= n; i++) { + alpha[i] = 0; + width[i] = 0; + theta0[i] = 0; + } + + memory->create(setflag,n+1,"angle:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void AngleGaussian::coeff(int narg, char **arg) +{ + if (narg < 6) error->all(FLERR,"Incorrect args for angle coefficients"); + + int ilo,ihi; + utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error); + + double angle_temperature_one = utils::numeric(FLERR,arg[1],false,lmp); + int n = utils::inumeric(FLERR,arg[2],false,lmp); + if (narg != 3*n + 3) + error->all(FLERR,"Incorrect args for angle coefficients"); + + if (!allocated) allocate(); + + // convert theta0 from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + angle_temperature[i] = angle_temperature_one; + nterms[i] = n; + alpha[i] = new double [n]; + width[i] = new double [n]; + theta0[i] = new double [n]; + for (int j = 0; j < n; j++ ) { + alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp); + width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp); + theta0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp)* MY_PI / 180.0; + setflag[i] = 1; + } + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleGaussian::equilibrium_angle(int i) +{ + return theta0[i][0]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleGaussian::write_restart(FILE *fp) +{ + fwrite(&angle_temperature[1],sizeof(double),atom->nangletypes,fp); + fwrite(&nterms[1],sizeof(int),atom->nangletypes,fp); + for(int i = 1; i <= atom->nangletypes; i++) { + fwrite(alpha[i],sizeof(double),nterms[i],fp); + fwrite(width[i],sizeof(double),nterms[i],fp); + fwrite(theta0[i],sizeof(double),nterms[i],fp); + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleGaussian::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR,&angle_temperature[1],sizeof(double),atom->nangletypes,fp,nullptr,error); + utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nangletypes,fp,nullptr,error); + } + MPI_Bcast(&angle_temperature[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&nterms[1],atom->nangletypes,MPI_INT,0,world); + + // allocate + for(int i = 1; i <= atom->nangletypes; i++) { + alpha[i] = new double [nterms[i]]; + width[i] = new double [nterms[i]]; + theta0[i] = new double [nterms[i]]; + } + + if (comm->me == 0) { + for(int i = 1; i <= atom->nangletypes; i++) { + utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error); + utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error); + utils::sfread(FLERR,theta0[i],sizeof(double),nterms[i],fp,nullptr,error); + } + } + + for(int i = 1; i <= atom->nangletypes; i++) { + MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world); + MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world); + MPI_Bcast(theta0[i],nterms[i],MPI_DOUBLE,0,world); + } + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void AngleGaussian::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nangletypes; i++) { + fprintf(fp,"%d %g %d",i,angle_temperature[i],nterms[i]); + for (int j = 0; j < nterms[i]; j++) { + fprintf(fp," %g %g %g",alpha[i][j],width[i][j],(theta0[i][j]/MY_PI)*180.0); + } + fprintf(fp, "\n"); + } + +} + +/* ---------------------------------------------------------------------- */ + +double AngleGaussian::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c) ; + + double sum_g_i = 0.0; + double sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dtheta = theta - theta0[type][i]; + double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2))); + double exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]); + double g_i = prefactor*exp(exponent); + sum_g_i += g_i; + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + return -(force->boltz*angle_temperature[type])*log(sum_g_i); +} diff --git a/src/USER-MISC/angle_gaussian.h b/src/USER-MISC/angle_gaussian.h new file mode 100644 index 0000000000..cd97d62670 --- /dev/null +++ b/src/USER-MISC/angle_gaussian.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(gaussian,AngleGaussian) + +#else + +#ifndef LMP_ANGLE_GAUSSIAN_H +#define LMP_ANGLE_GAUSSIAN_H + +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleGaussian : public Angle { + public: + AngleGaussian(class LAMMPS *); + virtual ~AngleGaussian(); + virtual void compute(int, int); + virtual void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + virtual void read_restart(FILE *); + void write_data(FILE *); + double single(int, int, int, int); + + protected: + int *nterms; + double *angle_temperature; + double **alpha,**width,**theta0; + + virtual void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for angle coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/src/USER-MISC/bond_gaussian.cpp b/src/USER-MISC/bond_gaussian.cpp new file mode 100644 index 0000000000..b1e9d3d66b --- /dev/null +++ b/src/USER-MISC/bond_gaussian.cpp @@ -0,0 +1,288 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "bond_gaussian.h" + +#include +#include +#include "atom.h" +#include "neighbor.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 1.0e-10 + +/* ---------------------------------------------------------------------- */ + +BondGaussian::BondGaussian(LAMMPS *lmp) : Bond(lmp) +{ + reinitflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +BondGaussian::~BondGaussian() +{ + if (allocated && !copymode) { + memory->destroy(setflag); + memory->destroy(nterms); + memory->destroy(bond_temperature); + for (int i = 1; i <= atom->nbondtypes; i++) { + if (alpha[i]) delete [] alpha[i]; + if (width[i]) delete [] width[i]; + if (r0[i]) delete [] r0[i]; + } + delete [] alpha; + delete [] width; + delete [] r0; + } +} + +/* ---------------------------------------------------------------------- */ + +void BondGaussian::compute(int eflag, int vflag) +{ + int i1,i2,n,type; + double delx,dely,delz,ebond,fbond; + double rsq,r,dr; + double prefactor, exponent, g_i, sum_g_i, sum_numerator; + + ebond = 0.0; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + rsq = delx*delx + dely*dely + delz*delz; + r = sqrt(rsq); + + sum_g_i = 0.0; + sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + dr = r - r0[type][i]; + prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2))); + exponent = -2*dr*dr/(width[type][i]*width[type][i]); + g_i = prefactor*exp(exponent); + sum_g_i += g_i; + sum_numerator += g_i*dr/(width[type][i]*width[type][i]); + } + + // force & energy + if (sum_g_i < SMALL) sum_g_i = SMALL; + + if (r > 0.0) fbond = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r; + else fbond = 0.0; + + if (eflag) ebond = -(force->boltz*bond_temperature[type])*log(sum_g_i); + + // apply force to each of 2 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx*fbond; + f[i1][1] += dely*fbond; + f[i1][2] += delz*fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx*fbond; + f[i2][1] -= dely*fbond; + f[i2][2] -= delz*fbond; + } + + if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondGaussian::allocate() +{ + allocated = 1; + int n = atom->nbondtypes; + + memory->create(nterms,n+1,"bond:nterms"); + memory->create(bond_temperature,n+1,"bond:bond_temperature"); + + alpha = new double *[n+1]; + width = new double *[n+1]; + r0 = new double *[n+1]; + for (int i = 1; i <= n; i++) { + alpha[i] = 0; + width[i] = 0; + r0[i] = 0; + } + + memory->create(setflag,n+1,"bond:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondGaussian::coeff(int narg, char **arg) +{ + if (narg < 6) error->all(FLERR,"Incorrect args for bond coefficients"); + + int ilo,ihi; + utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error); + + double bond_temp_one = utils::numeric(FLERR,arg[1],false,lmp); + int n = utils::inumeric(FLERR,arg[2],false,lmp); + if (narg != 3*n + 3) + error->all(FLERR,"Incorrect args for bond coefficients"); + + if (!allocated) allocate(); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + bond_temperature[i] = bond_temp_one; + nterms[i] = n; + alpha[i] = new double [n]; + width[i] = new double [n]; + r0[i] = new double [n]; + for (int j = 0; j < n; j++ ) { + alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp); + width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp); + r0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp); + setflag[i] = 1; + } + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + return an equilbrium bond length +------------------------------------------------------------------------- */ + +double BondGaussian::equilibrium_distance(int i) +{ + return r0[i][0]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondGaussian::write_restart(FILE *fp) +{ + fwrite(&bond_temperature[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&nterms[1],sizeof(int),atom->nbondtypes,fp); + for(int i = 1; i <= atom->nbondtypes; i++) { + fwrite(alpha[i],sizeof(double),nterms[i],fp); + fwrite(width[i],sizeof(double),nterms[i],fp); + fwrite(r0[i],sizeof(double),nterms[i],fp); + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondGaussian::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR,&bond_temperature[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nbondtypes,fp,nullptr,error); + } + MPI_Bcast(&bond_temperature[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&nterms[1],atom->nbondtypes,MPI_INT,0,world); + + // allocate + for(int i = 1; i <= atom->nbondtypes; i++) { + alpha[i] = new double [nterms[i]]; + width[i] = new double [nterms[i]]; + r0[i] = new double [nterms[i]]; + } + + if (comm->me == 0) { + for(int i = 1; i <= atom->nbondtypes; i++) { + utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error); + utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error); + utils::sfread(FLERR,r0[i],sizeof(double),nterms[i],fp,nullptr,error); + } + } + + for(int i = 1; i <= atom->nbondtypes; i++) { + MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world); + MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world); + MPI_Bcast(r0[i],nterms[i],MPI_DOUBLE,0,world); + } + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondGaussian::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) { + fprintf(fp,"%d %g %d",i,bond_temperature[i],nterms[i]); + for (int j = 0; j < nterms[i]; j++) { + fprintf(fp," %g %g %g",alpha[i][j],width[i][j],r0[i][j]); + } + fprintf(fp, "\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/, + double &fforce) +{ + double r = sqrt(rsq); + fforce = 0; + + double sum_g_i = 0.0; + double sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dr = r - r0[type][i]; + double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2))); + double exponent = -2*dr*dr/(width[type][i]*width[type][i]); + double g_i = prefactor*exp(exponent); + sum_g_i += g_i; + sum_numerator += g_i*dr/(width[type][i]*width[type][i]); + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + if (r > 0.0) fforce = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r; + + return -(force->boltz*bond_temperature[type])*log(sum_g_i); +} + diff --git a/src/USER-MISC/bond_gaussian.h b/src/USER-MISC/bond_gaussian.h new file mode 100644 index 0000000000..6fb587ff2b --- /dev/null +++ b/src/USER-MISC/bond_gaussian.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef BOND_CLASS + +BondStyle(gaussian,BondGaussian) + +#else + +#ifndef LMP_BOND_GAUSSIAN_H +#define LMP_BOND_GAUSSIAN_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondGaussian : public Bond { + public: + BondGaussian(class LAMMPS *); + virtual ~BondGaussian(); + virtual void compute(int, int); + virtual void coeff(int, char **); + double equilibrium_distance(int); + void write_restart(FILE *); + virtual void read_restart(FILE *); + void write_data(FILE *); + double single(int, double, int, int, double &); + + protected: + int *nterms; + double *bond_temperature; + double **alpha,**width,**r0; + + virtual void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for bond coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/src/USER-MISC/fix_electron_stopping_fit.cpp b/src/USER-MISC/fix_electron_stopping_fit.cpp new file mode 100644 index 0000000000..c58bffe83b --- /dev/null +++ b/src/USER-MISC/fix_electron_stopping_fit.cpp @@ -0,0 +1,208 @@ +/* --------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------ */ + +/* --------------------------------------------------------------------- + Contributing authors: Stephen M. Foiles (SNL) + James A. Stewart (SNL) + ------------------------------------------------------------------ */ + +#include "fix_electron_stopping_fit.h" + +#include "atom.h" +#include "citeme.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_special.h" +#include "modify.h" +#include "region.h" +#include "respa.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +// --------------------------------------------------------------------- + +static const char cite_fix_electron_stopping_fit_c[] = + "fix electron/stopping/fit command:\n\n" + "@Article{Stewart2018,\n" + " author = { J.A. Stewart and G. Brookman and P. Price and M. Franco and W. Ji and K. Hattar and R. Dingreville },\n" + " title = { Characterizing single isolated radiation-damage events from molecular dynamics via virtual diffraction methods },\n" + " journal = { Journal of Applied Physics },\n" + " year = { 2018 },\n" + " volume = { 123 },\n" + " number = { 16 },\n" + " pages = { 165902 }\n" + "}\n\n" + "@Article{Lee2020,\n" + " author = { C.W. Lee and J.A. Stewart and S.M. Foiles and R. Dingreville and A. Schleife },\n" + " title = { Multiscale simulations of electron and ion dynamics in self-irradiated silicon },\n" + " journal = { Physical Review B },\n" + " year = { 2020 },\n" + " volume = { 102 },\n" + " number = { 2 },\n" + " pages = { 024107 }\n" + "}\n\n"; + +// --------------------------------------------------------------------- + +FixElectronStoppingFit::FixElectronStoppingFit(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp,narg,arg), energy_coh_in(nullptr), drag_fac_in_1(nullptr), + drag_fac_in_2(nullptr), drag_fac_1(nullptr), drag_fac_2(nullptr), + v_min_sq(nullptr), v_max_sq(nullptr) +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_electron_stopping_fit_c); + + if (narg < 3 + 3*atom->ntypes) { + error->all(FLERR,"Incorrect number of fix electron/stopping/fit arguments"); + } + + scalar_flag = 1; + global_freq = 1; + + energy_coh_in = new double[atom->ntypes+1]; + + drag_fac_in_1 = new double[atom->ntypes+1]; + drag_fac_in_2 = new double[atom->ntypes+1]; + + for (int i = 1; i <= atom->ntypes; i++) { + energy_coh_in[i] = utils::numeric(FLERR,arg[3*i],false,lmp); + drag_fac_in_1[i] = utils::numeric(FLERR,arg[3*i+1],false,lmp); + drag_fac_in_2[i] = utils::numeric(FLERR,arg[3*i+2],false,lmp); + }; + + v_min_sq = new double[atom->ntypes+1]; + v_max_sq = new double[atom->ntypes+1]; + + drag_fac_1 = new double[atom->ntypes+1]; + drag_fac_2 = new double[atom->ntypes+1]; + + for (int i = 1; i <= atom->ntypes; i++) { + double mvv; + mvv = 2.0*energy_coh_in[i]/force->mvv2e; + v_min_sq[i] = 1.0*mvv/atom->mass[i]; + v_max_sq[i] = 2.0*mvv/atom->mass[i]; + drag_fac_1[i] = drag_fac_in_1[i]; + drag_fac_2[i] = drag_fac_in_2[i]; + }; +}; + +// --------------------------------------------------------------------- + +FixElectronStoppingFit::~FixElectronStoppingFit() +{ + delete [] energy_coh_in; + delete [] drag_fac_in_1; + delete [] drag_fac_in_2; + delete [] drag_fac_1; + delete [] drag_fac_2; + delete [] v_min_sq; + delete [] v_max_sq; +}; + +// --------------------------------------------------------------------- + +int FixElectronStoppingFit::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +}; + +// --------------------------------------------------------------------- + +void FixElectronStoppingFit::init() +{ + electronic_loss_this_node = 0.; + electronic_loss = 0.; + f_dot_v_prior = 0.; + f_dot_v_current = 0.; + last_step = update->ntimestep; +}; + +// --------------------------------------------------------------------- + +void FixElectronStoppingFit::setup(int vflag) +{ + if (strcmp(update->integrate_style,"verlet") == 0) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + }; +}; + +// --------------------------------------------------------------------- + +void FixElectronStoppingFit::post_force(int vflag) +{ + double **v = atom->v; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + + f_dot_v_current = 0.0; + for (int i = 0; i < nlocal; i++) { + double vv = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; + if (vv > v_min_sq[type[i]]) { + double gamma_x; + double gamma_y; + double gamma_z; + double v_mag = sqrt(vv); + if (vv < v_max_sq[type[i]]) { + double frac = (vv - v_min_sq[type[i]])/(v_max_sq[type[i]] - v_min_sq[type[i]]); + gamma_x = frac*(drag_fac_2[type[i]]*v[i][0] + drag_fac_1[type[i]]); + gamma_y = frac*(drag_fac_2[type[i]]*v[i][1] + drag_fac_1[type[i]]); + gamma_z = frac*(drag_fac_2[type[i]]*v[i][2] + drag_fac_1[type[i]]); + } else { + gamma_x = drag_fac_2[type[i]]*v[i][0] + drag_fac_1[type[i]]; + gamma_y = drag_fac_2[type[i]]*v[i][1] + drag_fac_1[type[i]]; + gamma_z = drag_fac_2[type[i]]*v[i][2] + drag_fac_1[type[i]]; + }; + f[i][0] -= gamma_x*v[i][0]; + f[i][1] -= gamma_y*v[i][1]; + f[i][2] -= gamma_z*v[i][2]; + f_dot_v_current += v_mag*sqrt( MathSpecial::square(gamma_x*v[i][0]) + + MathSpecial::square(gamma_y*v[i][1]) + + MathSpecial::square(gamma_z*v[i][2]) ); + }; + }; + this_step = update->ntimestep; + electronic_loss_this_node += (this_step - last_step)*update->dt*0.5*(f_dot_v_prior + f_dot_v_current); + last_step = this_step; + f_dot_v_prior = f_dot_v_current; +}; + +// --------------------------------------------------------------------- + +void FixElectronStoppingFit::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +}; + +// --------------------------------------------------------------------- + +double FixElectronStoppingFit::compute_scalar() +{ + MPI_Allreduce(&electronic_loss_this_node,&electronic_loss,1,MPI_DOUBLE,MPI_SUM,world); + return electronic_loss; +}; + +// --------------------------------------------------------------------- diff --git a/src/USER-MISC/fix_electron_stopping_fit.h b/src/USER-MISC/fix_electron_stopping_fit.h new file mode 100644 index 0000000000..3237e93f99 --- /dev/null +++ b/src/USER-MISC/fix_electron_stopping_fit.h @@ -0,0 +1,54 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Stephen M. Foiles (SNL) + James A. Stewart (SNL) +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(electron/stopping/fit,FixElectronStoppingFit) + +#else + +#ifndef LMP_FIX_ELECTRON_STOPPING_FIT_H +#define LMP_FIX_ELECTRON_STOPPING_FIT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixElectronStoppingFit : public Fix { + public: + FixElectronStoppingFit(class LAMMPS *, int, char **); + ~FixElectronStoppingFit(); + int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double compute_scalar(); + + private: + double *energy_coh_in,*v_min_sq,*v_max_sq,*drag_fac_in_1,*drag_fac_in_2,*drag_fac_1,*drag_fac_2; + double electronic_loss,electronic_loss_this_node; + double f_dot_v_prior,f_dot_v_current; + int last_step,this_step; + int nlevels_respa; +}; + +} + +#endif +#endif diff --git a/unittest/commands/test_kim_commands.cpp b/unittest/commands/test_kim_commands.cpp index dd7ac09b19..42e8a441bd 100644 --- a/unittest/commands/test_kim_commands.cpp +++ b/unittest/commands/test_kim_commands.cpp @@ -17,6 +17,7 @@ #include "lammps.h" #include "lmppython.h" #include "modify.h" +#include "output.h" #include "utils.h" #include "variable.h" #include "gmock/gmock.h" @@ -82,17 +83,24 @@ TEST_F(KimCommandsTest, kim_init) { if (!LAMMPS::is_installed_pkg("KIM")) GTEST_SKIP(); - TEST_FAILURE(".*ERROR: Illegal kim_init command.*", lmp->input->one("kim_init");); + TEST_FAILURE(".*ERROR: Illegal kim_init command.*", + lmp->input->one("kim_init");); TEST_FAILURE(".*ERROR: Illegal kim_init command.*", lmp->input->one("kim_init LennardJones_Ar real si");); TEST_FAILURE(".*ERROR: LAMMPS unit_style lj not supported by KIM models.*", lmp->input->one("kim_init LennardJones_Ar lj");); + TEST_FAILURE(".*ERROR: LAMMPS unit_style micro not supported by KIM models.*", + lmp->input->one("kim_init LennardJones_Ar micro");); + TEST_FAILURE(".*ERROR: LAMMPS unit_style nano not supported by KIM models.*", + lmp->input->one("kim_init LennardJones_Ar nano");); TEST_FAILURE(".*ERROR: Unknown unit_style.*", lmp->input->one("kim_init LennardJones_Ar new_style");); TEST_FAILURE(".*ERROR: KIM Model name not found.*", lmp->input->one("kim_init Unknown_Model real");); TEST_FAILURE(".*ERROR: Incompatible units for KIM Simulator Model, required units = metal.*", lmp->input->one("kim_init Sim_LAMMPS_LJcut_AkersonElliott_Alchemy_PbAu real");); + // TEST_FAILURE(".*ERROR: KIM Model does not support the requested unit system.*", + // lmp->input->one("kim_init ex_model_Ar_P_Morse real");); if (!verbose) ::testing::internal::CaptureStdout(); lmp->input->one("kim_init LennardJones_Ar real"); @@ -125,6 +133,17 @@ TEST_F(KimCommandsTest, kim_interactions) lmp->input->one("create_atoms 1 box"); if (!verbose) ::testing::internal::GetCapturedStdout(); + TEST_FAILURE(".*ERROR: Illegal kim_interactions command.*", + lmp->input->one("kim_interactions Ar Ar");); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("clear"); + lmp->input->one("lattice fcc 4.4300"); + lmp->input->one("region box block 0 20 0 20 0 20"); + lmp->input->one("create_box 4 box"); + lmp->input->one("create_atoms 4 box"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + TEST_FAILURE(".*ERROR: Illegal kim_interactions command.*", lmp->input->one("kim_interactions Ar Ar");); @@ -172,6 +191,21 @@ TEST_F(KimCommandsTest, kim_interactions) TEST_FAILURE(".*ERROR: Species 'Ar' is not supported by this KIM Simulator Model.*", lmp->input->one("kim_interactions Ar");); + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("clear"); + lmp->input->one("kim_init Sim_LAMMPS_LJcut_AkersonElliott_Alchemy_PbAu metal"); + lmp->input->one("lattice fcc 4.08"); + lmp->input->one("region box block 0 10 0 10 0 10"); + lmp->input->one("create_box 1 box"); + lmp->input->one("create_atoms 1 box"); + lmp->input->one("kim_interactions Au"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + // ASSERT_EQ(lmp->output->var_kim_periodic, 1); + // TEST_FAILURE(".*ERROR: Incompatible units for KIM Simulator Model.*", + // lmp->input->one("kim_interactions Au");); + + if (!verbose) ::testing::internal::CaptureStdout(); lmp->input->one("clear"); lmp->input->one("kim_init LennardJones_Ar real"); @@ -302,6 +336,133 @@ TEST_F(KimCommandsTest, kim_property) } } +TEST_F(KimCommandsTest, kim_query) +{ + if (!LAMMPS::is_installed_pkg("KIM")) GTEST_SKIP(); + + TEST_FAILURE(".*ERROR: Illegal kim_query command.*", + lmp->input->one("kim_query");); + TEST_FAILURE(".*ERROR: Must use 'kim_init' before 'kim_query'.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic");); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("clear"); + lmp->input->one("kim_init LennardJones612_UniversalShifted__MO_959249795837_003 real"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + TEST_FAILURE(".*ERROR: Illegal kim_query command.\nThe keyword 'split' " + "must be followed by the name of the query function.*", + lmp->input->one("kim_query a0 split");); + + TEST_FAILURE(".*ERROR: Illegal kim_query command.\nThe 'list' keyword " + "can not be used after 'split'.*", + lmp->input->one("kim_query a0 split list");); + + TEST_FAILURE(".*ERROR: Illegal kim_query command.\nThe 'list' keyword " + "must be followed by \\('split' and\\) the name of the query " + "function.*", lmp->input->one("kim_query a0 list");); + + TEST_FAILURE(".*ERROR: Illegal 'model' key in kim_query command.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic " + "model=[MO_959249795837_003]");); + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `crystal` " + "to kim_query is wrong. The query format is the " + "keyword=\\[value\\], where value is always an array of one " + "or more comma-separated items.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic " + "crystal");); + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `" + "crystal=fcc` to kim_query is wrong. The query format is the " + "keyword=\\[value\\], where value is always an array of one " + "or more comma-separated items.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic " + "crystal=fcc");); + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `" + "crystal=\\[fcc` to kim_query is wrong. The query format is " + "the keyword=\\[value\\], where value is always an array of " + "one or more comma-separated items.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic " + "crystal=[fcc");); + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `" + "crystal=fcc\\]` to kim_query is wrong. The query format is " + "the keyword=\\[value\\], where value is always an array of " + "one or more comma-separated items.*", + lmp->input->one("kim_query a0 get_lattice_constant_cubic " + "crystal=fcc]");); + + std::string squery("kim_query a0 get_lattice_constant_cubic "); + squery += "crystal=[\"fcc\"] species=\"Al\",\"Ni\" units=[\"angstrom\"]"; + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `species=" + "\"Al\",\"Ni\"` to kim_query is wrong. The query format is " + "the keyword=\\[value\\], where value is always an array of " + "one or more comma-separated items.*", + lmp->input->one(squery);); + + squery = "kim_query a0 get_lattice_constant_cubic "; + squery += "crystal=[\"fcc\"] species=\"Al\",\"Ni\", units=[\"angstrom\"]"; + + TEST_FAILURE(".*ERROR: Illegal query format.\nInput argument of `species=" + "\"Al\",\"Ni\",` to kim_query is wrong. The query format is " + "the keyword=\\[value\\], where value is always an array of " + "one or more comma-separated items.*", + lmp->input->one(squery);); + + squery = "kim_query a0 get_lattice_constant_cubic crystal=[fcc] " + "species=[Al]"; + TEST_FAILURE(".*ERROR: OpenKIM query failed:.*", lmp->input->one(squery);); + + squery = "kim_query a0 get_lattice_constant_cubic crystal=[fcc] " + "units=[\"angstrom\"]"; + TEST_FAILURE(".*ERROR: OpenKIM query failed:.*", lmp->input->one(squery);); + + // if (!verbose) ::testing::internal::CaptureStdout(); + // lmp->input->one("clear"); + // lmp->input->one("kim_init EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000 metal"); + + // squery = "kim_query latconst split get_lattice_constant_hexagonal "; + // squery += "crystal=[\"hcp\"] species=[\"Zr\"] units=[\"angstrom\"]"; + // lmp->input->one(squery); + // if (!verbose) ::testing::internal::GetCapturedStdout(); + + // ASSERT_TRUE((std::string(lmp->input->variable->retrieve("latconst_1")) == + // std::string("3.234055244384789"))); + // ASSERT_TRUE((std::string(lmp->input->variable->retrieve("latconst_2")) == + // std::string("5.167650199630013"))); + + // if (!verbose) ::testing::internal::CaptureStdout(); + // lmp->input->one("clear"); + // lmp->input->one("kim_init EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000 metal"); + + // squery = "kim_query latconst list get_lattice_constant_hexagonal "; + // squery += "crystal=[hcp] species=[Zr] units=[angstrom]"; + // lmp->input->one(squery); + // if (!verbose) ::testing::internal::GetCapturedStdout(); + + // ASSERT_TRUE((std::string(lmp->input->variable->retrieve("latconst")) == + // std::string("3.234055244384789 5.167650199630013"))); + + // squery = "kim_query latconst list get_lattice_constant_hexagonal "; + // squery += "crystal=[bcc] species=[Zr] units=[angstrom]"; + // TEST_FAILURE(".*ERROR: OpenKIM query failed:.*", lmp->input->one(squery);); + + // if (!verbose) ::testing::internal::CaptureStdout(); + // lmp->input->one("clear"); + // lmp->input->one("kim_init EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005 metal"); + + // squery = "kim_query alpha get_linear_thermal_expansion_coefficient_cubic "; + // squery += "crystal=[fcc] species=[Al] units=[1/K] temperature=[293.15] "; + // squery += "temperature_units=[K]"; + // lmp->input->one(squery); + // if (!verbose) ::testing::internal::GetCapturedStdout(); + + // ASSERT_TRUE((std::string(lmp->input->variable->retrieve("alpha")) == + // std::string("1.654960564704273e-05"))); +} } // namespace LAMMPS_NS int main(int argc, char **argv) diff --git a/unittest/force-styles/tests/angle-gaussian.yaml b/unittest/force-styles/tests/angle-gaussian.yaml new file mode 100644 index 0000000000..0103811913 --- /dev/null +++ b/unittest/force-styles/tests/angle-gaussian.yaml @@ -0,0 +1,86 @@ +--- +lammps_version: 29 Oct 2020 +date_generated: Sat Nov 14 17:18:12 202 +epsilon: 2.5e-13 +prerequisites: ! | + atom full + angle gaussian +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +angle_style: gaussian +angle_coeff: ! | + 1 300.0 3 0.2189 8.66 88.1 0.5439 9.94 142.7 0.2307 17.14 167.0 + 2 300.0 3 0.0214 14.29 85.3 0.3934 18.22 118.1 0.9411 31.21 136.2 + 3 300.0 2 0.2189 8.66 88.1 0.5439 9.94 142.7 + 4 300.0 2 0.0214 14.29 85.3 0.3934 18.22 118.1 +equilibrium: 4 1.5376350710070041 1.4887658519511628 1.5376350710070041 1.4887658519511628 +extract: ! "" +natoms: 29 +init_energy: 57.7940091589437 +init_stress: ! |2- + 1.8540667929702014e-02 -2.1128225814185397e-02 2.5875578844833853e-03 1.9120388327532396e-02 4.4564216146325795e-03 9.8400107526492309e-03 +init_forces: ! |2 + 1 2.0043481356521267e-03 -2.0981767618283602e-04 -2.3209495543181176e-03 + 2 3.5529860865363239e-05 2.8767662992335252e-04 2.3132719898919820e-04 + 3 -3.8667546050125514e-04 2.0816982895064929e-03 3.3491343253347289e-03 + 4 -4.7471002298472755e-04 -1.0657600563862162e-03 2.4796171798629677e-04 + 5 -1.6406411064805123e-03 -1.9359456692002001e-03 -9.0394626833801086e-04 + 6 7.6979393607576755e-03 -8.1291748957998018e-03 -9.9817242122740776e-03 + 7 -4.2572919294439045e-03 3.8055081960951839e-03 -9.2914562787063215e-04 + 8 -5.6354813852788994e-03 -2.4077018853792594e-03 2.2075705642764603e-02 + 9 -8.5366182106296319e-04 7.7489027943554716e-04 4.3064306944232773e-06 + 10 5.1692145475614156e-03 1.1913856117509856e-02 -1.2167886610859078e-02 + 11 -1.1932769900266460e-03 -5.6695560883419894e-04 6.3912856495926295e-04 + 12 1.5380901789534356e-03 -4.4133578770737883e-03 -3.1448353955117270e-03 + 13 1.0633492902104076e-04 2.6203621306745674e-04 1.4238358783950607e-04 + 14 7.0133753297702093e-04 -2.4463789752656221e-04 -2.8963548243354467e-04 + 15 4.8989433563144716e-04 -2.8356800785119403e-04 8.3150095382329269e-04 + 16 -3.7177027531705075e-03 6.9883235686930205e-04 1.9103996085473906e-03 + 17 4.1675258752988943e-04 -5.6757850817313379e-04 3.0627512066648804e-04 + 18 2.4034769486436930e-04 2.5116260424105398e-03 -1.0192991295370015e-02 + 19 -2.9428816898762833e-03 -3.6662332265825250e-03 4.3652227193094872e-03 + 20 2.7025339950119140e-03 1.1546071841719852e-03 5.8277685760605271e-03 + 21 2.8475989038751108e-03 3.5647233204607149e-03 -1.0653298109293498e-02 + 22 -5.0707463777986093e-03 -3.4572321794689781e-03 3.7660889578173145e-03 + 23 2.2231474739234985e-03 -1.0749114099173709e-04 6.8872091514761841e-03 + 24 -2.2557311483102044e-03 8.7681606485922377e-03 -5.2202265270714686e-03 + 25 -1.8094194422620109e-03 -6.0622399916731113e-03 9.9144196300429460e-04 + 26 4.0651505905722153e-03 -2.7059206569191264e-03 4.2287845640671742e-03 + 27 -1.1419349153643011e-03 9.3743910881221686e-03 -3.1932287532239921e-03 + 28 -2.4969086656509735e-03 -5.6050861901068560e-03 1.3020567790072643e-04 + 29 3.6388435810152746e-03 -3.7693048980153125e-03 3.0630230753232656e-03 +run_energy: 57.7939328926191 +run_stress: ! |2- + 1.8333297775264226e-02 -2.1096486667701031e-02 2.7631888924368117e-03 1.9016013391593267e-02 4.3940373471828961e-03 9.6831700037143531e-03 +run_forces: ! |2 + 1 1.9869422162142472e-03 -1.9353557824978272e-04 -2.3006702071328641e-03 + 2 3.3476059893797837e-05 2.7319932557553025e-04 2.2063210015606759e-04 + 3 -3.4876540425126176e-04 2.0590710647431592e-03 3.3392846732858885e-03 + 4 -4.7868102566407460e-04 -1.0516647507108046e-03 2.6031681319326416e-04 + 5 -1.6499221858598965e-03 -1.9312689427161936e-03 -9.0479891915083840e-04 + 6 7.7122074651881124e-03 -8.1501276471632952e-03 -1.0025269020742612e-02 + 7 -4.2683877840090532e-03 3.8176212404067746e-03 -9.3444409464597659e-04 + 8 -5.6308797825345026e-03 -2.3959292864231662e-03 2.2111508090050601e-02 + 9 -8.5650859745461569e-04 7.7767395513096715e-04 5.1882091688871450e-06 + 10 5.1639019878229383e-03 1.1894316899845880e-02 -1.2166342881754515e-02 + 11 -1.1924185597197892e-03 -5.6615421744289282e-04 6.3860767740208757e-04 + 12 1.5195074473272755e-03 -4.4151667225361133e-03 -3.1180448718694286e-03 + 13 1.0363114665204553e-04 2.5647171532249322e-04 1.4526915428658157e-04 + 14 7.0368317634416993e-04 -2.2809763878071953e-04 -2.8897971051002481e-04 + 15 4.9492006640407434e-04 -2.7632678493888751e-04 8.0977392634859263e-04 + 16 -3.7088817774414895e-03 6.9716958266289306e-04 1.9018611181369511e-03 + 17 4.1617555108802267e-04 -5.6725221472584374e-04 3.0610794377733961e-04 + 18 2.4756430759883402e-04 2.4867826389144984e-03 -1.0052208969303231e-02 + 19 -2.8956245324249508e-03 -3.6148801680506949e-03 4.3154199446524145e-03 + 20 2.6480602248261167e-03 1.1280975291361967e-03 5.7367890246508169e-03 + 21 2.8270003409956433e-03 3.5141039792481541e-03 -1.0563471348354937e-02 + 22 -5.0172695982285236e-03 -3.4085472370788878e-03 3.7417759205389629e-03 + 23 2.1902692572328803e-03 -1.0555674216926607e-04 6.8216954278159728e-03 + 24 -2.2433452545004657e-03 8.7085710861042054e-03 -5.1795033574261800e-03 + 25 -1.7836848781861761e-03 -6.0195638750124483e-03 9.8492810012785017e-04 + 26 4.0270301326866421e-03 -2.6890072110917563e-03 4.1945752572983301e-03 + 27 -1.1493082457243627e-03 9.3604335737891321e-03 -3.1700077185484042e-03 + 28 -2.4838761558667495e-03 -5.5967979792576599e-03 1.2157136916697260e-04 + 29 3.6331844015911122e-03 -3.7636355945314713e-03 3.0484363493814316e-03 +... diff --git a/unittest/force-styles/tests/bond-gaussian.yaml b/unittest/force-styles/tests/bond-gaussian.yaml new file mode 100644 index 0000000000..92157d7155 --- /dev/null +++ b/unittest/force-styles/tests/bond-gaussian.yaml @@ -0,0 +1,87 @@ +--- +lammps_version: 29 Oct 2020 +date_generated: Sat Nov 14 16:49:01 202 +epsilon: 2.5e-13 +prerequisites: ! | + atom full + bond gaussian +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +bond_style: gaussian +bond_coeff: ! | + 1 300.0 1 0.0100 0.098 1.45 + 2 300.0 2 0.0128 0.375 1.37 0.0730 0.148 2.63 + 3 300.0 3 0.0003 0.657 1.61 0.0050 0.360 1.07 0.0048 0.157 1.31 + 4 300.0 1 0.0100 0.098 2.45 + 5 300.0 1 0.0100 0.098 2.85 +equilibrium: 5 1.45 1.37 1.61 2.45 2.85 +extract: ! "" +natoms: 29 +init_energy: 194.978038221632 +init_stress: ! |- + -2.4024989684355553e+01 -3.8521513996632500e+01 -1.0851224048428129e+01 1.2562604359180053e+01 1.3677283516797356e+01 4.3206731051245653e+00 +init_forces: ! |2 + 1 -1.7791337913398690e+00 -5.2745532425491986e+00 -1.9333096530222391e+00 + 2 7.8999913149794128e-293 6.5010830500033665e-293 -9.2980646648301405e-293 + 3 2.4197086198752562e+01 -1.2911571268065043e+01 -1.2153319969868038e+01 + 4 -3.5002110421521651e+00 9.8124800657318079e-01 -2.4834895420880554e+00 + 5 -8.7934593181833831e-01 -1.3513167937313169e+00 4.4900533574430685e+00 + 6 -1.9224405194016612e+01 1.9525383982308810e+01 1.1251608936919853e+01 + 7 2.6580140740726381e-131 -1.3633763941647238e-130 -6.8018769495047054e-130 + 8 1.4462104594211977e+00 -1.2568711136582216e+00 7.3991622652588918e-01 + 9 1.2099652614352605e-300 1.3032068217192395e-300 5.3545155818429412e-300 + 10 1.8282673669124623e+01 -6.7893037436650294e-01 1.0475143579619905e+01 + 11 -9.5181855408160265e-01 -2.3577388099405021e+00 -3.8685744266264179e+00 + 12 -1.1761121482537199e+01 -1.1840691118605761e+01 8.9587696830512531e+00 + 13 3.9348879648968196e+00 -1.5566010373601853e+00 -7.3956496855403397e-02 + 14 -1.5580348688551586e+00 3.1703943744370217e-01 -4.0404862787928506e+00 + 15 -1.0483110905921594e-01 4.0280962447539723e+00 1.4354708657826634e+00 + 16 -8.1019563183350432e+00 1.2376506087197068e+01 -1.2797826282089627e+01 + 17 -9.6845722000297944e-125 6.7536031200741501e-125 2.5693469616608658e-124 + 18 5.0042083741224387e-291 3.2014176819490257e-291 6.0624670892900674e-291 + 19 -5.0042167517970120e-291 -3.2014265949545701e-291 -6.0624614384187022e-291 + 20 8.3776745733654894e-297 8.9130055442585484e-297 -5.6508713648842736e-297 + 21 5.0373663727594610e-296 1.1676684296048456e-296 8.1823232295641435e-296 + 22 -5.1857245273845906e-296 -1.2567112623130275e-296 -8.1358238807042024e-296 + 23 1.4835815462512912e-297 8.9042832708182009e-298 -4.6499348859940937e-298 + 24 6.5124799547612842e-295 -1.0579059065054233e-295 5.4786730014873485e-295 + 25 -6.5176382072810523e-295 1.0492453069148130e-295 -5.4792561056911984e-295 + 26 5.1582525197680877e-298 8.6605995906103569e-298 5.8310420384964103e-299 + 27 -1.5677247388593395e-295 -1.8232011058192963e-295 -3.8038051984576450e-296 + 28 -3.2483754529644398e-299 1.3960035208884715e-299 -2.1978823514938368e-299 + 29 1.5680495764046360e-295 1.8230615054672073e-295 3.8060030808091389e-296 +run_energy: 194.96889016686 +run_stress: ! |- + -2.4084235648269384e+01 -3.8596877573973650e+01 -1.0971337511117875e+01 1.2627485208541385e+01 1.3589007837800324e+01 4.4518443361436777e+00 +run_forces: ! |2 + 1 -1.7800915383536471e+00 -5.2662174638478936e+00 -1.9311810441446928e+00 + 2 9.1200716389742962e-293 7.5205784896271243e-293 -1.0695855329374170e-292 + 3 2.4188774318819682e+01 -1.2910730800434983e+01 -1.2139174094227805e+01 + 4 -3.4905807721708837e+00 9.7423802985974728e-01 -2.4827066691937869e+00 + 5 -8.7826414385513407e-01 -1.3507945719900971e+00 4.4847167409249762e+00 + 6 -1.9198711248640532e+01 1.9501343007070176e+01 1.1259539605043198e+01 + 7 4.0781500460380220e-131 -2.0934766207882755e-130 -1.0411772151605081e-129 + 8 1.4035232720380466e+00 -1.2181526258990241e+00 7.2552718656771575e-01 + 9 1.4877356608185432e-300 1.5947265521745610e-300 6.5759628249586203e-300 + 10 1.8340705485218969e+01 -7.9602516938863732e-01 1.0533434146468263e+01 + 11 -9.4713695434855716e-01 -2.3455928036230933e+00 -3.8477133980837270e+00 + 12 -1.1753841378581289e+01 -1.1839528950721563e+01 8.9356024501072664e+00 + 13 3.9289793641831325e+00 -1.5460483921060724e+00 -7.3078087497547045e-02 + 14 -1.5515717239320088e+00 3.1019421574866657e-01 -4.0233193667488729e+00 + 15 -1.1312732638809736e-01 4.0290637402465492e+00 1.4439547691915919e+00 + 16 -8.1486573539896803e+00 1.2458251785086224e+01 -1.2885602238406578e+01 + 17 -8.5522515805489358e-125 5.9749160301406998e-125 2.2702237597406565e-124 + 18 2.5382954259673697e-291 1.6282298856292719e-291 3.0672317979786876e-291 + 19 -2.5383561239391082e-291 -1.6282944740463789e-291 -3.0671910793881731e-291 + 20 6.0697971738423079e-296 6.4588417107197222e-296 -4.0718590514496707e-296 + 21 3.1636825215784415e-296 7.4502521705718285e-297 5.0914419661316058e-296 + 22 -3.2413538119513539e-296 -7.9143971383319095e-297 -5.0672219270657353e-296 + 23 7.7671290372912634e-298 4.6414496776008138e-298 -2.4220039065870281e-298 + 24 1.1528889554480086e-295 -1.8584672369333140e-296 9.7061626349018667e-296 + 25 -1.1544439355951613e-295 1.8323577266329387e-296 -9.7079719071127095e-296 + 26 1.5549801471527681e-298 2.6109510300375245e-298 1.8092722108425850e-299 + 27 -1.0502291554946705e-295 -1.2226612584790533e-295 -2.5738911540368265e-296 + 28 -1.8342692926757559e-299 7.8715078988712594e-300 -1.2385711775450889e-299 + 29 1.0504125824239381e-295 1.2225825434000646e-295 2.5751297252143716e-296 +...