diff --git a/.github/workflows/coverity.yml b/.github/workflows/coverity.yml index c0c3e3f89a..2691b9e895 100644 --- a/.github/workflows/coverity.yml +++ b/.github/workflows/coverity.yml @@ -67,7 +67,6 @@ jobs: -D PKG_MANIFOLD=on \ -D PKG_MDI=on \ -D PKG_MGPT=on \ - -D PKG_ML-PACE=on \ -D PKG_ML-RANN=on \ -D PKG_MOLFILE=on \ -D PKG_NETCDF=on \ diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml new file mode 100644 index 0000000000..d7dc602476 --- /dev/null +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -0,0 +1,53 @@ +# GitHub action to build LAMMPS-GUI as a flatpak bundle +name: "Build LAMMPS-GUI as flatpak bundle" + +on: + push: + branches: + - develop + + workflow_dispatch: + +concurrency: + group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{github.event_name == 'pull_request'}} + +jobs: + build: + name: LAMMPS-GUI flatpak build + if: ${{ github.repository == 'lammps/lammps' }} + runs-on: ubuntu-latest + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + with: + fetch-depth: 2 + + - name: Install extra packages + run: | + sudo apt-get update + sudo apt-get install -y ccache \ + libeigen3-dev \ + libcurl4-openssl-dev \ + mold \ + ninja-build \ + python3-dev \ + flatpak \ + flatpak-builder + + - name: Set up access to flatpak repo + run: flatpak --user remote-add --if-not-exists flathub https://dl.flathub.org/repo/flathub.flatpakrepo + + - name: Build flatpak + run: | + mkdir flatpack-state + sed -i -e 's/branch:.*/branch: develop/' tools/lammps-gui/org.lammps.lammps-gui.yml + flatpak-builder --force-clean --verbose --repo=flatpak-repo \ + --install-deps-from=flathub --state-dir=flatpak-state \ + --user --ccache --default-branch=${{ github.ref_name }} \ + flatpak-build tools/lammps-gui/org.lammps.lammps-gui.yml + flatpak build-bundle --runtime-repo=https://flathub.org/repo/flathub.flatpakrepo \ + --verbose flatpak-repo LAMMPS-Linux-x86_64-GUI.flatpak \ + org.lammps.lammps-gui ${{ github.ref_name }} + flatpak install -y -v --user LAMMPS-Linux-x86_64-GUI.flatpak diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 7c88fb58ed..ff0d69e316 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -605,6 +605,16 @@ foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE endif() endforeach() +# settings for misc packages and styles +if(PKG_MISC) + option(LAMMPS_ASYNC_IMD "Asynchronous IMD processing" OFF) + mark_as_advanced(LAMMPS_ASYNC_IMD) + if(LAMMPS_ASYNC_IMD) + target_compile_definitions(lammps PRIVATE -DLAMMPS_ASYNC_IMD) + message(STATUS "Using IMD in asynchronous mode") + endif() +endif() + # optionally enable building script wrappers using swig option(WITH_SWIG "Build scripting language wrappers with SWIG" OFF) if(WITH_SWIG) diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index 1f13525f98..dab2267ee8 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -48,6 +48,7 @@ This is the list of packages that may require additional steps. * :ref:`LEPTON ` * :ref:`MACHDYN ` * :ref:`MDI ` + * :ref:`MISC ` * :ref:`ML-HDNNP ` * :ref:`ML-IAP ` * :ref:`ML-PACE ` @@ -2031,7 +2032,7 @@ TBB and MKL. .. _mdi: MDI package ------------------------------ +----------- .. tabs:: @@ -2058,6 +2059,37 @@ MDI package ---------- +.. _misc: + +MISC package +------------ + +The :doc:`fix imd ` style in this package can be run either +synchronously (communication with IMD clients is done in the main +process) or asynchronously (the fix spawns a separate thread that can +communicate with IMD clients concurrently to the LAMMPS execution). + +.. tabs:: + + .. tab:: CMake build + + .. code-block:: bash + + -D LAMMPS_ASYNC_IMD=value # Run IMD server asynchronously + # value = no (default) or yes + + .. tab:: Traditional make + + To enable asynchronous mode the ``-DLAMMPS_ASYNC_IMD`` define + needs to be added to the ``LMP_INC`` variable in the + ``Makefile.machine`` you are using. For example: + + .. code-block:: make + + LMP_INC = -DLAMMPS_ASYNC_IMD -DLAMMPS_MEMALIGN=64 + +---------- + .. _molfile: MOLFILE package diff --git a/doc/src/Build_package.rst b/doc/src/Build_package.rst index 8b2da8ad30..c4c4889806 100644 --- a/doc/src/Build_package.rst +++ b/doc/src/Build_package.rst @@ -49,6 +49,7 @@ packages: * :ref:`LEPTON ` * :ref:`MACHDYN ` * :ref:`MDI ` + * :ref:`MISC ` * :ref:`ML-HDNNP ` * :ref:`ML-IAP ` * :ref:`ML-PACE ` diff --git a/doc/src/Commands_removed.rst b/doc/src/Commands_removed.rst index 777be3bf16..be775f4c19 100644 --- a/doc/src/Commands_removed.rst +++ b/doc/src/Commands_removed.rst @@ -1,6 +1,10 @@ Removed commands and packages ============================= +.. contents:: + +------ + This page lists LAMMPS commands and packages that have been removed from the distribution and provides suggestions for alternatives or replacements. LAMMPS has special dummy styles implemented, that will @@ -8,47 +12,60 @@ stop LAMMPS and print a suitable error message in most cases, when a style/command is used that has been removed or will replace the command with the direct alternative (if available) and print a warning. -restart2data tool ------------------ +LAMMPS shell +------------ -.. versionchanged:: 23Nov2013 +.. versionchanged:: 29Aug2024 -The functionality of the restart2data tool has been folded into the -LAMMPS executable directly instead of having a separate tool. A -combination of the commands :doc:`read_restart ` and -:doc:`write_data ` can be used to the same effect. For -added convenience this conversion can also be triggered by -:doc:`command-line flags ` +The LAMMPS shell has been removed from the LAMMPS distribution. Users +are encouraged to use the :ref:`LAMMPS-GUI ` tool instead. -Fix ave/spatial and fix ave/spatial/sphere ------------------------------------------- +i-PI tool +--------- -.. deprecated:: 11Dec2015 +.. versionchanged:: 27Jun2024 -The fixes ave/spatial and ave/spatial/sphere have been removed from LAMMPS -since they were superseded by the more general and extensible "chunk -infrastructure". Here the system is partitioned in one of many possible -ways through the :doc:`compute chunk/atom ` command -and then averaging is done using :doc:`fix ave/chunk `. -Please refer to the :doc:`chunk HOWTO ` section for an overview. +The i-PI tool has been removed from the LAMMPS distribution. Instead, +instructions to install i-PI from PyPI via pip are provided. -Box command ------------ +USER-REAXC package +------------------ -.. deprecated:: 22Dec2022 +.. deprecated:: 7Feb2024 -The *box* command has been removed and the LAMMPS code changed so it won't -be needed. If present, LAMMPS will ignore the command and print a warning. +The USER-REAXC package has been renamed to :ref:`REAXFF `. +In the process also the pair style and related fixes were renamed to use +the "reaxff" string instead of "reax/c". For a while LAMMPS was maintaining +backward compatibility by providing aliases for the styles. These have +been removed, so using "reaxff" is now *required*. -Reset_ids, reset_atom_ids, reset_mol_ids commands -------------------------------------------------- +MPIIO package +------------- -.. deprecated:: 22Dec2022 +.. deprecated:: 21Nov2023 -The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have -been folded into the :doc:`reset_atoms ` command. If -present, LAMMPS will replace the commands accordingly and print a -warning. +The MPIIO package has been removed from LAMMPS since it was unmaintained +for many years and thus not updated to incorporate required changes that +had been applied to the corresponding non-MPIIO commands. As a +consequence the MPIIO commands had become unreliable and sometimes +crashing LAMMPS or corrupting data. Similar functionality is available +through the :ref:`ADIOS package ` and the :ref:`NETCDF +package `. Also, the :doc:`dump_modify nfile or dump_modify +fileper ` keywords may be used for an efficient way of +writing out dump files when running on large numbers of processors. +Similarly, the "nfile" and "fileper" keywords exist for restarts: +see :doc:`restart `, :doc:`read_restart `, +:doc:`write_restart `. + +MSCG package +------------ + +.. deprecated:: 21Nov2023 + +The MSCG package has been removed from LAMMPS since it was unmaintained +for many years and instead superseded by the `OpenMSCG software +`_ of the Voth group at the +University of Chicago, which can be used independent from LAMMPS. LATTE package ------------- @@ -64,18 +81,6 @@ packages, including LATTE. See the ``examples/QUANTUM`` dir and the with LATTE as a plugin library (similar to the way fix latte worked), as well as on a different set of MPI processors. -MEAM package ------------- - -The MEAM package in Fortran has been replaced by a C++ implementation. -The code in the :ref:`MEAM package ` is a translation of the -Fortran code of MEAM into C++, which removes several restrictions -(e.g. there can be multiple instances in hybrid pair styles) and allows -for some optimizations leading to better performance. The pair style -:doc:`meam ` has the exact same syntax. For a transition -period the C++ version of MEAM was called USER-MEAMC so it could -coexist with the Fortran version. - Minimize style fire/old ----------------------- @@ -97,38 +102,38 @@ The same functionality is available through :doc:`bond style mesocnt ` and :doc:`angle style mesocnt `. -MPIIO package -------------- +Box command +----------- -.. deprecated:: 21Nov2023 +.. deprecated:: 22Dec2022 -The MPIIO package has been removed from LAMMPS since it was unmaintained -for many years and thus not updated to incorporate required changes that -had been applied to the corresponding non-MPIIO commands. As a -consequence the MPIIO commands had become unreliable and sometimes -crashing LAMMPS or corrupting data. Similar functionality is available -through the :ref:`ADIOS package ` and the :ref:`NETCDF -package `. Also, the :doc:`dump_modify nfile or dump_modify -fileper ` keywords may be used for an efficient way of -writing out dump files when running on large numbers of processors. -Similarly, the "nfile" and "fileper" keywords exist for restarts: -see :doc:`restart `, :doc:`read_restart `, -:doc:`write_restart `. +The *box* command has been removed and the LAMMPS code changed so it won't +be needed. If present, LAMMPS will ignore the command and print a warning. +Reset_ids, reset_atom_ids, reset_mol_ids commands +------------------------------------------------- -MSCG package ------------- +.. deprecated:: 22Dec2022 -.. deprecated:: 21Nov2023 +The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have +been folded into the :doc:`reset_atoms ` command. If +present, LAMMPS will replace the commands accordingly and print a +warning. -The MSCG package has been removed from LAMMPS since it was unmaintained -for many years and instead superseded by the `OpenMSCG software -`_ of the Voth group at the -University of Chicago, which can be used independent from LAMMPS. +MESSAGE package +--------------- + +.. deprecated:: 4May2022 + +The MESSAGE package has been removed since it was superseded by the +:ref:`MDI package `. MDI implements the same functionality +and in a more general way with direct support for more applications. REAX package ------------ +.. deprecated:: 4Jan2019 + The REAX package has been removed since it was superseded by the :ref:`REAXFF package `. The REAXFF package has been tested to yield equivalent results to the REAX package, offers better @@ -138,20 +143,25 @@ syntax compatible with the removed reax pair style, so input files will have to be adapted. The REAXFF package was originally called USER-REAXC. -USER-REAXC package ------------------- +MEAM package +------------ -.. deprecated:: 7Feb2024 +.. deprecated:: 4Jan2019 -The USER-REAXC package has been renamed to :ref:`REAXFF `. -In the process also the pair style and related fixes were renamed to use -the "reaxff" string instead of "reax/c". For a while LAMMPS was maintaining -backward compatibility by providing aliases for the styles. These have -been removed, so using "reaxff" is now *required*. +The MEAM package in Fortran has been replaced by a C++ implementation. +The code in the :ref:`MEAM package ` is a translation of the +Fortran code of MEAM into C++, which removes several restrictions +(e.g. there can be multiple instances in hybrid pair styles) and allows +for some optimizations leading to better performance. The pair style +:doc:`meam ` has the exact same syntax. For a transition +period the C++ version of MEAM was called USER-MEAMC so it could +coexist with the Fortran version. USER-CUDA package ----------------- +.. deprecated:: 31May2016 + The USER-CUDA package had been removed, since it had been unmaintained for a long time and had known bugs and problems. Significant parts of the design were transferred to the @@ -160,19 +170,27 @@ performance characteristics on NVIDIA GPUs. Both, the KOKKOS and the :ref:`GPU package ` are maintained and allow running LAMMPS with GPU acceleration. -i-PI tool ---------- +Fix ave/spatial and fix ave/spatial/sphere +------------------------------------------ -.. versionchanged:: 27Jun2024 +.. deprecated:: 11Dec2015 -The i-PI tool has been removed from the LAMMPS distribution. Instead, -instructions to install i-PI from PyPI via pip are provided. +The fixes ave/spatial and ave/spatial/sphere have been removed from LAMMPS +since they were superseded by the more general and extensible "chunk +infrastructure". Here the system is partitioned in one of many possible +ways through the :doc:`compute chunk/atom ` command +and then averaging is done using :doc:`fix ave/chunk `. +Please refer to the :doc:`chunk HOWTO ` section for an overview. -LAMMPS shell ------------- +restart2data tool +----------------- -.. versionchanged:: 29Aug2024 +.. deprecated:: 23Nov2013 -The LAMMPS shell has been removed from the LAMMPS distribution. Users -are encouraged to use the :ref:`LAMMPS-GUI ` tool instead. +The functionality of the restart2data tool has been folded into the +LAMMPS executable directly instead of having a separate tool. A +combination of the commands :doc:`read_restart ` and +:doc:`write_data ` can be used to the same effect. For +added convenience this conversion can also be triggered by +:doc:`command-line flags ` diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index d50c7bdcf3..bc641a237f 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -321,6 +321,8 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype set_string_variable: subroutine :f set_internal_variable: :f:subr:`set_internal_variable` :ftype set_internal_variable: subroutine + :f eval: :f:func:`eval` + :ftype eval: function :f gather_atoms: :f:subr:`gather_atoms` :ftype gather_atoms: subroutine :f gather_atoms_concat: :f:subr:`gather_atoms_concat` diff --git a/doc/src/Howto_lammps_gui.rst b/doc/src/Howto_lammps_gui.rst index 08be9fb00c..63cf859c57 100644 --- a/doc/src/Howto_lammps_gui.rst +++ b/doc/src/Howto_lammps_gui.rst @@ -64,13 +64,18 @@ simple LAMMPS simulations. It is very suitable for tutorials on LAMMPS since you only need to learn how to use a single program for most tasks and thus time can be saved and people can focus on learning LAMMPS. The tutorials at https://lammpstutorials.github.io/ are specifically -updated for use with LAMMPS-GUI. +updated for use with LAMMPS-GUI and can their tutorial materials can +be downloaded and loaded directly from the GUI. Another design goal is to keep the barrier low when replacing part of the functionality of LAMMPS-GUI with external tools. That said, LAMMPS-GUI has some unique functionality that is not found elsewhere: - auto-adapting to features available in the integrated LAMMPS library +- auto-completion for LAMMPS commands and options +- context-sensitive online help +- start and stop of simulations via mouse or keyboard +- monitoring of simulation progress - interactive visualization using the :doc:`dump image ` command with the option to copy-paste the resulting settings - automatic slide show generation from dump image out at runtime diff --git a/doc/src/Library_execute.rst b/doc/src/Library_execute.rst index 44b601ba4c..8560fbb148 100644 --- a/doc/src/Library_execute.rst +++ b/doc/src/Library_execute.rst @@ -7,6 +7,7 @@ This section documents the following functions: - :cpp:func:`lammps_command` - :cpp:func:`lammps_commands_list` - :cpp:func:`lammps_commands_string` +- :cpp:func:`lammps_expand` -------------------- @@ -79,3 +80,8 @@ Below is a short example using some of these functions. .. doxygenfunction:: lammps_commands_string :project: progguide +----------------------- + +.. doxygenfunction:: lammps_expand + :project: progguide + diff --git a/doc/src/Library_objects.rst b/doc/src/Library_objects.rst index 7c0ca824d7..53edfce9e6 100644 --- a/doc/src/Library_objects.rst +++ b/doc/src/Library_objects.rst @@ -1,5 +1,5 @@ -Compute, fixes, variables -========================= +Computes, fixes, variables +========================== This section documents accessing or modifying data stored by computes, fixes, or variables in LAMMPS using the following functions: @@ -12,6 +12,7 @@ fixes, or variables in LAMMPS using the following functions: - :cpp:func:`lammps_set_string_variable` - :cpp:func:`lammps_set_internal_variable` - :cpp:func:`lammps_variable_info` +- :cpp:func:`lammps_eval` ----------------------- @@ -55,6 +56,11 @@ fixes, or variables in LAMMPS using the following functions: ----------------------- +.. doxygenfunction:: lammps_eval + :project: progguide + +----------------------- + .. doxygenenum:: _LMP_DATATYPE_CONST .. doxygenenum:: _LMP_STYLE_CONST diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index 2974262c81..a4d7d9d387 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -26,6 +26,29 @@ Syntax *nowait* arg = *on* or *off* off = LAMMPS waits to be connected to an IMD client before continuing (default) on = LAMMPS listens for an IMD client, but continues with the run + *version* arg = *2* or *3* + 2 = use IMD protocol version 2 (default) + 3 = use IMD protocol version 3. + + The following keywords are only supported for IMD protocol version 3. + + .. parsed-literal:: + + *time* arg = *on* or *off* + off = simulation time is not transmitted (default) + on = simulation time is transmitted. + *box* arg = *on* or *off* + off = simulation box data is not transmitted (default) + on = simulation box data is transmitted. + *coordinates* arg = *on* or *off* + off = atomic coordinates are not transmitted (default) + on = atomic coordinates are transmitted. + *velocities* arg = *on* or *off* + off = atomic velocities are not transmitted (default) + on = atomic velocities are transmitted. + *forces* arg = *on* or *off* + off = atomic forces are not transmitted (default) + on = atomic forces are transmitted. Examples """""""" @@ -40,16 +63,19 @@ Description This fix implements the "Interactive MD" (IMD) protocol which allows realtime visualization and manipulation of MD simulations through the -IMD protocol, as initially implemented in VMD and NAMD. Specifically -it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_, so that it can monitor the progress of the -simulation and interactively apply forces to selected atoms. +IMD protocol, as initially implemented in VMD and NAMD. Specifically it +allows LAMMPS to connect an IMD client, for example the `VMD +visualization program `_ (currently only supports IMDv2) or the +`Python IMDClient `_ (supports both IMDv2 and IMDv3), so +that it can monitor the progress of the simulation and interactively +apply forces to selected atoms. -If LAMMPS is compiled with the pre-processor flag -DLAMMPS_ASYNC_IMD -then fix imd will use POSIX threads to spawn a IMD communication -thread on MPI rank 0 in order to offload data reading and writing -from the main execution thread and potentially lower the inferred -latencies for slow communication links. This feature has only been -tested under linux. +If LAMMPS is compiled with the pre-processor flag +:ref:`-DLAMMPS_ASYNC_IMD ` then fix imd will use POSIX threads to +spawn an IMD communication thread on MPI rank 0 in order to offload data +exchange with the IMD client from the main execution thread and +potentially lower the inferred latencies for slow communication +links. This feature has only been tested under linux. The source code for this fix includes code developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced @@ -94,6 +120,15 @@ with different units or as a measure to tweak the forces generated by the manipulation of the IMD client, this option allows to make adjustments. +.. versionadded:: TBD + +In `IMDv3 `_, the IMD protocol has been extended to allow for +the transmission of simulation time, box dimensions, atomic coordinates, +velocities, and forces. The *version* keyword allows to select the +version of the protocol to be used. The *time*, *box*, *coordinates*, +*velocities*, and *forces* keywords allow to select which data is +transmitted to the IMD client. The default is to transmit all data. + To connect VMD to a listening LAMMPS simulation on the same machine with fix imd enabled, one needs to start VMD and load a coordinate or topology file that matches the fix group. When the VMD command @@ -129,6 +164,10 @@ screen output is active. .. _imdvmd: https://www.ks.uiuc.edu/Research/vmd/imd/ +.. _IMDClient: https://github.com/Becksteinlab/imdclient/tree/main/imdclient + +.. _IMDv3: https://imdclient.readthedocs.io/en/latest/protocol_v3.html + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -147,14 +186,14 @@ This fix is part of the MISC package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. -When used in combination with VMD, a topology or coordinate file has -to be loaded, which matches (in number and ordering of atoms) the -group the fix is applied to. The fix internally sorts atom IDs by -ascending integer value; in VMD (and thus the IMD protocol) those will -be assigned 0-based consecutive index numbers. +When used in combination with VMD, a topology or coordinate file has to +be loaded, which matches (in number and ordering of atoms) the group the +fix is applied to. The fix internally sorts atom IDs by ascending +integer value; in VMD (and thus the IMD protocol) those will be assigned +0-based consecutive index numbers. When using multiple active IMD connections at the same time, each -needs to use a different port number. +fix instance needs to use a different port number. Related commands """""""""""""""" diff --git a/examples/COUPLE/plugin/liblammpsplugin.c b/examples/COUPLE/plugin/liblammpsplugin.c index 5b1308a5cc..99e38df32b 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.c +++ b/examples/COUPLE/plugin/liblammpsplugin.c @@ -117,6 +117,7 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib) ADDSYM(set_string_variable); ADDSYM(set_internal_variable); ADDSYM(variable_info); + ADDSYM(eval); ADDSYM(gather_atoms); ADDSYM(gather_atoms_concat); diff --git a/examples/COUPLE/plugin/liblammpsplugin.h b/examples/COUPLE/plugin/liblammpsplugin.h index b2df98d630..dd1c9628f5 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.h +++ b/examples/COUPLE/plugin/liblammpsplugin.h @@ -163,6 +163,7 @@ struct _liblammpsplugin { int (*set_string_variable)(void *, const char *, const char *); int (*set_internal_variable)(void *, const char *, double); int (*variable_info)(void *, int, char *, int); + double (*eval)(void *, const char *); void (*gather_atoms)(void *, const char *, int, int, void *); void (*gather_atoms_concat)(void *, const char *, int, int, void *); diff --git a/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py b/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py new file mode 100644 index 0000000000..176e15e811 --- /dev/null +++ b/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py @@ -0,0 +1,13 @@ +""" +For use with 'in.deca-ala-solv_imd_v3'. + +Tested with imdclient v0.1.4 and MDAnalysis v2.8.0 +""" +from imdclient.IMD import IMDReader +import MDAnalysis as mda + +u = mda.Universe('data.deca-ala-solv', "imd://localhost:5678", topology_format='DATA') + +for ts in u.trajectory: + print(ts.time) + print(ts.velocities) \ No newline at end of file diff --git a/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 b/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 new file mode 100644 index 0000000000..a7e8244cdc --- /dev/null +++ b/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 @@ -0,0 +1,31 @@ +# +units real +neighbor 2.5 bin +neigh_modify delay 1 every 1 + +atom_style full +bond_style harmonic +angle_style charmm +dihedral_style charmm +improper_style harmonic + +pair_style lj/charmm/coul/long 8 10 +pair_modify mix arithmetic +special_bonds charmm +read_data data.deca-ala-solv + + +group peptide id <= 103 +fix rigidh all shake 1e-6 100 1000 t 1 2 3 4 5 a 23 + +thermo 100 +thermo_style multi +timestep 2.0 +kspace_style pppm 1e-5 + +fix ensemble all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 0.2 + +# IMD setup. Client code available in 'deca-ala-solv_imd_v3.py' +fix comm all imd 5678 unwrap on trate 10 version 3 time on box on coordinates on velocities on forces off + +run 5000000 \ No newline at end of file diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 1cc6992c27..552b3dfad3 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -126,6 +126,7 @@ MODULE LIBLAMMPS PROCEDURE :: set_variable => lmp_set_variable PROCEDURE :: set_string_variable => lmp_set_string_variable PROCEDURE :: set_internal_variable => lmp_set_internal_variable + PROCEDURE :: eval => lmp_eval PROCEDURE, PRIVATE :: lmp_gather_atoms_int PROCEDURE, PRIVATE :: lmp_gather_atoms_double GENERIC :: gather_atoms => lmp_gather_atoms_int, & @@ -618,7 +619,14 @@ MODULE LIBLAMMPS INTEGER(c_int) :: lammps_set_internal_variable END FUNCTION lammps_set_internal_variable - SUBROUTINE lammps_gather_atoms(handle, name, type, count, data) BIND(C) + FUNCTION lammps_eval(handle, expr) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, expr + REAL(c_double) :: lammps_eval + END FUNCTION lammps_eval + + SUBROUTINE lammps_gather_atoms(handle, name, TYPE, count, DATA) BIND(C) IMPORT :: c_int, c_ptr IMPLICIT NONE TYPE(c_ptr), VALUE :: handle, name, data @@ -1812,7 +1820,7 @@ CONTAINS SUBROUTINE lmp_set_internal_variable(self, name, val) CLASS(lammps), INTENT(IN) :: self CHARACTER(LEN=*), INTENT(IN) :: name - REAL(KIND=c_double), INTENT(IN) :: val + REAL(c_double), INTENT(IN) :: val INTEGER :: err TYPE(c_ptr) :: Cname @@ -1826,6 +1834,18 @@ CONTAINS END IF END SUBROUTINE lmp_set_internal_variable + ! equivalent function to lammps_eval + FUNCTION lmp_eval(self, expr) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: expr + REAL(c_double) :: lmp_eval + TYPE(c_ptr) :: Cexpr + + Cexpr = f2c_string(expr) + lmp_eval = lammps_eval(self%handle, Cexpr) + CALL lammps_free(Cexpr) + END FUNCTION lmp_eval + ! equivalent function to lammps_gather_atoms (for integers) SUBROUTINE lmp_gather_atoms_int(self, name, count, data) CLASS(lammps), INTENT(IN) :: self diff --git a/python/lammps/core.py b/python/lammps/core.py index c7be0c84b5..7e70b5392a 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -103,7 +103,7 @@ class command_wrapper(object): This method is where the Python 'magic' happens. If a method is not defined by the class command_wrapper, it assumes it is a LAMMPS command. It takes all the arguments, concatinates them to a single string, and executes it using - :py:meth:`lammps.command()`. + :py:meth:`lammps.command`. Starting with Python 3.6 it also supports keyword arguments. key=value is transformed into 'key value'. Note, since these have come last in the @@ -426,6 +426,9 @@ class lammps(object): self.lib.lammps_compute_addstep.argtype = [c_void_p, self.c_bigint] self.lib.lammps_compute_addstep_all.argtype = [c_void_p, self.c_bigint] + self.lib.lammps_eval.argtypes = [c_void_p, c_char_p] + self.lib.lammps_eval.restype = c_double + self.lib.lammps_fix_external_get_force.argtypes = [c_void_p, c_char_p] self.lib.lammps_fix_external_get_force.restype = POINTER(POINTER(c_double)) @@ -1693,6 +1696,30 @@ class lammps(object): # ------------------------------------------------------------------------- + def eval(self, expr): + """ Evaluate a LAMMPS immediate variable expression + + .. versionadded:: TBD + + This function is a wrapper around the function :cpp:func:`lammps_eval` + of the C library interface. It evaluates and expression like in + immediate variables and returns the value. + + :param expr: immediate variable expression + :type name: string + :return: the result of the evaluation + :rtype: c_double + """ + + if expr: newexpr = expr.encode() + else: return None + + with ExceptionCheck(self): + return self.lib.lammps_eval(self.lmp, newexpr) + return None + + # ------------------------------------------------------------------------- + # return vector of atom properties gathered across procs # 3 variants to match src/library.cpp # name = atom property recognized by LAMMPS in atom->extract() diff --git a/python/lammps/ipython/wrapper.py b/python/lammps/ipython/wrapper.py index 729c0d62bf..6f6669f55a 100644 --- a/python/lammps/ipython/wrapper.py +++ b/python/lammps/ipython/wrapper.py @@ -17,7 +17,7 @@ ################################################################################ class wrapper(object): - """lammps API IPython Wrapper + """ lammps API IPython Wrapper This is a wrapper class that provides additional methods on top of an existing :py:class:`lammps` instance. It provides additional methods diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index 0683e63752..ba6807641d 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -63,6 +63,8 @@ void ComputeReaxFFAtomKokkos::init() template void ComputeReaxFFAtomKokkos::compute_bonds() { + invoked_bonds = update->ntimestep; + if (atom->nmax > nmax) { memory->destroy(array_atom); nmax = atom->nmax; diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 41ce97ee65..0dbafe49cb 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -43,7 +43,7 @@ negotiate an appropriate license for such distribution." ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing authors: Axel Kohlmeyer (Temple U), Lawson Woods (Arizona State U) IMD API, hash, and socket code written by: John E. Stone, Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) ------------------------------------------------------------------------- */ @@ -80,6 +80,8 @@ negotiate an appropriate license for such distribution." using namespace LAMMPS_NS; using namespace FixConst; +/* ---------------------------------------------------------------------- */ + /* re-usable integer hash table code with static linkage. */ /** hash table top level data structure */ @@ -358,7 +360,6 @@ typedef struct { } IMDheader; #define IMDHEADERSIZE 8 -#define IMDVERSION 2 typedef enum IMDType_t { IMD_DISCONNECT, /**< close IMD connection, leaving sim running */ @@ -370,7 +371,15 @@ typedef enum IMDType_t { IMD_MDCOMM, /**< MDComm style force data */ IMD_PAUSE, /**< pause the running simulation */ IMD_TRATE, /**< set IMD update transmission rate */ - IMD_IOERROR /**< indicate an I/O error */ + IMD_IOERROR, /**< indicate an I/O error */ + /* IMDv3 only headers */ + IMD_SESSIONINFO, + IMD_RESUME, + IMD_TIME, + IMD_BOX, + IMD_VELOCITIES, + IMD_FORCES, + IMD_WAIT, } IMDType; /**< IMD command message type enumerations */ typedef struct { @@ -386,8 +395,20 @@ typedef struct { float Eimpr; /**< Improper energy, Kcal/mol */ } IMDEnergies; /**< IMD simulation energy report structure */ +/* IMDv3 only */ +typedef struct IMDSessionInfo { + bool time; + bool box; + bool coords; + bool unwrap; + bool velocities; + bool forces; + bool energies; +} IMDSessionInfo; + /** Send control messages - these consist of a header with no subsequent data */ -static int imd_handshake(void *); /**< check endianness, version compat */ +static int imd_handshake_v2(void *); /**< check endianness, version compat */ +static int imd_handshake_v3(void *, IMDSessionInfo *); /** Receive header and data */ static IMDType imd_recv_header(void *, int32 *); /** Receive MDComm-style forces, units are Kcal/mol/angstrom */ @@ -426,33 +447,46 @@ static void imdsock_destroy(void *); * The implementation follows at the end of the file. * ***************************************************************/ -/* struct for packed data communication of coordinates and forces. */ +/* struct for packed data communication of coordinates, velocities, and forces. */ struct commdata { tagint tag; float x,y,z; }; +MPI_Datatype MPI_CommData; + /*************************************************************** - * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID imd [unwrap (on|off)] [fscale ] + * create class and parse arguments in LAMMPS script. ***************************************************************/ + FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + Fix(lmp, narg, arg), localsock(nullptr), clientsock(nullptr), coord_data(nullptr), + vel_data(nullptr), force_data(nullptr), idmap(nullptr), rev_idmap(nullptr), + recv_force_buf(nullptr), imdsinfo(nullptr), msgdata(nullptr) { - if (narg < 4) - error->all(FLERR,"Illegal fix imd command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "fix imd", error); imd_port = utils::inumeric(FLERR,arg[3],false,lmp); - if (imd_port < 1024) - error->all(FLERR,"Illegal fix imd parameter: port < 1024"); + if (imd_port < 1024) error->all(FLERR,"Illegal fix imd parameter: port < 1024"); /* default values for optional flags */ + + imd_version = 2; unwrap_flag = 0; nowait_flag = 0; connect_msg = 1; imd_fscale = 1.0; imd_trate = 1; + // IMDv3-only flags. + // Aren't stored as class attributes since they're converted into IMDSessionInfo + + int time_flag = 1; + int box_flag = 1; + int coord_flag = 1; + int vel_flag = 1; + int force_flag = 1; + /* parse optional arguments */ int iarg = 4; while (iarg+1 < narg) { @@ -464,36 +498,94 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : imd_fscale = utils::numeric(FLERR,arg[iarg+1],false,lmp); } else if (0 == strcmp(arg[iarg], "trate")) { imd_trate = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + } else if (0 == strcmp(arg[iarg], "version")) { + imd_version = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + } else if (0 == strcmp(arg[iarg], "time")) { + time_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "box")) { + box_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "coordinates")) { + coord_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "velocities")) { + vel_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "forces")) { + force_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); } else { - error->all(FLERR,"Unknown fix imd parameter"); + error->all(FLERR,"Unknown fix imd parameter: {}", arg[iarg]); } ++iarg; ++iarg; } /* sanity check on parameters */ if (imd_trate < 1) - error->all(FLERR,"Illegal fix imd parameter. trate < 1."); + error->all(FLERR,"Illegal fix imd trate parameter. trate < 1."); + + if (imd_version != 2 && imd_version != 3) + error->all(FLERR,"Illegal fix imd version parameter: version != 2 or 3."); + + imdsinfo = new IMDSessionInfo; + + /* In IMDv2 in LAMMPS, only coordinates are sent*/ + if (imd_version == 2) { + imdsinfo->time = false; + imdsinfo->box = false; + imdsinfo->coords = true; + imdsinfo->unwrap = unwrap_flag; + imdsinfo->velocities = false; + imdsinfo->forces = false; + imdsinfo->energies = false; + } + + if (imd_version == 3) { + imdsinfo->time = time_flag; + imdsinfo->box = box_flag; + imdsinfo->coords = coord_flag; + imdsinfo->unwrap = unwrap_flag; + imdsinfo->velocities = vel_flag; + imdsinfo->forces = force_flag; + imdsinfo->energies = false; + } bigint n = group->count(igroup); if (n > MAXSMALLINT) error->all(FLERR,"Too many atoms for fix imd"); num_coords = static_cast (n); + // Define MPI dtype for passing x/v/f data + MPI_Type_contiguous(4, MPI_FLOAT, &MPI_CommData); + MPI_Type_commit(&MPI_CommData); + MPI_Comm_rank(world,&me); /* initialize various imd state variables. */ - clientsock = nullptr; - localsock = nullptr; nlevels_respa = 0; imd_inactive = 0; imd_terminate = 0; imd_forces = 0; - force_buf = nullptr; maxbuf = 0; - msgdata = nullptr; - msglen = 0; - comm_buf = nullptr; - idmap = nullptr; - rev_idmap = nullptr; + + if (imd_version == 3) { + msglen = 0; + if (imdsinfo->time) { + msglen += 24+IMDHEADERSIZE; + } + if (imdsinfo->box) { + msglen += 9*4+IMDHEADERSIZE; + } + if (imdsinfo->coords) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->velocities) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->forces) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + msgdata = new char[msglen]; + } + else { + msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; + msgdata = new char[msglen]; + } if (me == 0) { /* set up incoming socket on MPI rank 0. */ @@ -512,7 +604,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : if (imd_terminate) error->all(FLERR,"LAMMPS Terminated on error in IMD."); - /* storage required to communicate a single coordinate or force. */ + /* storage required to communicate a single coordinate, velocity, or force. */ size_one = sizeof(struct commdata); #if defined(LAMMPS_ASYNC_IMD) @@ -561,32 +653,39 @@ FixIMD::~FixIMD() pthread_cond_destroy(&write_cond); } #endif - auto hashtable = (taginthash_t *)idmap; - memory->destroy(comm_buf); - memory->destroy(force_buf); + memory->destroy(coord_data); + memory->destroy(vel_data); + memory->destroy(force_data); + + delete[] msgdata; + memory->destroy(recv_force_buf); taginthash_destroy(hashtable); delete hashtable; free(rev_idmap); + delete imdsinfo; // close sockets imdsock_shutdown(clientsock); imdsock_destroy(clientsock); imdsock_shutdown(localsock); imdsock_destroy(localsock); - clientsock=nullptr; - localsock=nullptr; + clientsock = nullptr; + localsock = nullptr; } /* ---------------------------------------------------------------------- */ + int FixIMD::setmask() { int mask = 0; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; + mask |= END_OF_STEP; return mask; } /* ---------------------------------------------------------------------- */ + void FixIMD::init() { if (utils::strmatch(update->integrate_style,"^respa")) @@ -611,6 +710,7 @@ int FixIMD::reconnect() } else { fprintf(screen,"Waiting for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); } + fflush(screen); } connect_msg = 0; clientsock = nullptr; @@ -637,7 +737,8 @@ int FixIMD::reconnect() return 0; } else { /* check endianness and IMD protocol version. */ - if (imd_handshake(clientsock)) { + if ((imd_version == 2 && imd_handshake_v2(clientsock)) || + (imd_version == 3 && imd_handshake_v3(clientsock, imdsinfo))) { if (screen) fprintf(screen, "IMD handshake error. Dropping connection.\n"); imdsock_destroy(clientsock); @@ -666,6 +767,17 @@ int FixIMD::reconnect() * buffers and collect tag/id maps. */ void FixIMD::setup(int) { + if (imd_version == 2) { + setup_v2(); + } + else { + setup_v3(); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixIMD::setup_v2() { /* nme: number of atoms in group on this MPI task * nmax: max number of atoms in group across all MPI tasks * nlocal: all local atoms @@ -680,9 +792,9 @@ void FixIMD::setup(int) if (mask[i] & groupbit) ++nme; MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); - memory->destroy(comm_buf); + memory->destroy(coord_data); maxbuf = nmax*size_one; - comm_buf = (void *) memory->smalloc(maxbuf,"imd:comm_buf"); + coord_data = (void *) memory->smalloc(maxbuf,"imd:coord_data"); connect_msg = 1; reconnect(); @@ -697,7 +809,7 @@ void FixIMD::setup(int) idmap = (void *)hashtable; int tmp, ndata; - auto buf = static_cast(comm_buf); + auto buf = static_cast(coord_data); if (me == 0) { MPI_Status status; @@ -714,7 +826,7 @@ void FixIMD::setup(int) /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); @@ -747,11 +859,125 @@ void FixIMD::setup(int) } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); - MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); } } +/* ---------------------------------------------------------------------- */ + +void FixIMD::setup_v3() +{ + /* nme: number of atoms in group on this MPI task + * nmax: max number of atoms in group across all MPI tasks + * nlocal: all local atoms + */ + int i,j; + int nmax,nme,nlocal; + int *mask = atom->mask; + tagint *tag = atom->tag; + nlocal = atom->nlocal; + nme=0; + for (i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); + + maxbuf = nmax*size_one; + + if (imdsinfo->coords) { + memory->destroy(coord_data); + coord_data = (void *) memory->smalloc(maxbuf,"imd:coord_data"); + } + if (imdsinfo->velocities) { + memory->destroy(vel_data); + vel_data = (void *) memory->smalloc(maxbuf,"imd:vel_data"); + } + if (imdsinfo->forces) { + memory->destroy(force_data); + force_data = (void *) memory->smalloc(maxbuf,"imd:force_data"); + } + + connect_msg = 1; + reconnect(); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on error in setting up IMD connection."); + + /* initialize and build hashtable. */ + auto hashtable=new taginthash_t; + taginthash_init(hashtable, num_coords); + idmap = (void *)hashtable; + + int tmp, ndata; + + struct commdata *buf = nullptr; + if (imdsinfo->coords) { + buf = static_cast(coord_data); + } + else if (imdsinfo->velocities) { + buf = static_cast(vel_data); + } + else if (imdsinfo->forces) { + buf = static_cast(force_data); + } + + if (me == 0) { + if (buf == nullptr) { + return; + } + MPI_Status status; + MPI_Request request; + auto taglist = new tagint[num_coords]; + int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ + + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + taglist[numtag] = tag[i]; + ++numtag; + } + } + + /* loop over procs to receive remote data */ + for (i=1; i < comm->nprocs; ++i) { + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Send(&tmp, 0, MPI_INT, i, 0, world); + MPI_Wait(&request, &status); + MPI_Get_count(&status, MPI_BYTE, &ndata); + ndata /= size_one; + + for (j=0; j < ndata; ++j) { + taglist[numtag] = buf[j].tag; + ++numtag; + } + } + + /* sort list of tags by value to have consistently the + * same list when running in parallel and build hash table. */ + id_sort(taglist, 0, num_coords-1); + for (i=0; i < num_coords; ++i) { + taginthash_insert(hashtable, taglist[i], i); + } + delete[] taglist; + + /* generate reverse index-to-tag map for communicating + * IMD forces back to the proper atoms */ + rev_idmap=taginthash_keys(hashtable); + } else { + nme=0; + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[nme].tag = tag[i]; + ++nme; + } + } + /* blocking receive to wait until it is our turn to send data. */ + MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); + } + + } /* worker threads for asynchronous i/o */ #if defined(LAMMPS_ASYNC_IMD) /* c bindings wrapper */ @@ -775,10 +1001,12 @@ void FixIMD::ioworker() pthread_exit(nullptr); } else if (buf_has_data > 0) { /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + fprintf(screen,"Asynchronous I/O thread terminated on error in sending IMDFrame.\n"); + pthread_mutex_unlock(&write_mutex); + pthread_exit(nullptr); } - delete[] msgdata; + /* Don't delete msgdata- allow memory to persist. */ buf_has_data=0; pthread_mutex_unlock(&write_mutex); } else { @@ -795,6 +1023,45 @@ void FixIMD::ioworker() * Send coodinates, energies, and add IMD forces to atoms. */ void FixIMD::post_force(int /*vflag*/) { + fflush(screen); + if (imd_version == 2) { + handle_step_v2(); + } + else if (imd_version == 3) { + handle_client_input_v3(); + } + + } + +/* ---------------------------------------------------------------------- */ + +void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) +{ + /* only process IMD on the outmost RESPA level. */ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixIMD::end_of_step() +{ + if (imd_version == 3 && update->ntimestep % imd_trate == 0) { + handle_output_v3(); + } +} + +/* ---------------------------------------------------------------------- */ +/* local memory usage. approximately. */ + +double FixIMD::memory_usage() +{ + return static_cast(num_coords+maxbuf+imd_forces)*size_one; +} + +/* ---------------------------------------------------------------------- */ + +void FixIMD::handle_step_v2() { + /* check for reconnect */ if (imd_inactive) { reconnect(); @@ -825,22 +1092,12 @@ void FixIMD::post_force(int /*vflag*/) switch(msg) { - case IMD_GO: - if (screen) - fprintf(screen, "Ignoring unexpected IMD_GO message.\n"); - break; - - case IMD_IOERROR: - if (screen) - fprintf(screen, "IMD connection lost.\n"); - /* fallthrough */ - case IMD_DISCONNECT: { /* disconnect from client. wait for new connection. */ imd_paused = 0; imd_forces = 0; - memory->destroy(force_buf); - force_buf = nullptr; + memory->destroy(recv_force_buf); + recv_force_buf = nullptr; imdsock_destroy(clientsock); clientsock = nullptr; if (screen) @@ -866,9 +1123,9 @@ void FixIMD::post_force(int /*vflag*/) case IMD_PAUSE: /* pause the running simulation. wait for second IMD_PAUSE to continue. */ if (imd_paused) { - if (screen) - fprintf(screen, "Continuing run on IMD client request.\n"); - imd_paused = 0; + if (screen) + fprintf(screen, "Continuing run on IMD client request.\n"); + imd_paused = 0; } else { if (screen) fprintf(screen, "Pausing run on IMD client request.\n"); @@ -884,31 +1141,18 @@ void FixIMD::post_force(int /*vflag*/) fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); break; - case IMD_ENERGIES: { - IMDEnergies dummy_energies; - imd_recv_energies(clientsock, &dummy_energies); - break; - } - - case IMD_FCOORDS: { - auto dummy_coords = new float[3*length]; - imd_recv_fcoords(clientsock, length, dummy_coords); - delete[] dummy_coords; - break; - } - case IMD_MDCOMM: { auto imd_tags = new int32[length]; auto imd_fdat = new float[3*length]; imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); if (imd_forces < length) { /* grow holding space for forces, if needed. */ - memory->destroy(force_buf); - force_buf = (void *) memory->smalloc((bigint)length*size_one, - "imd:force_buf"); + memory->destroy(recv_force_buf); + recv_force_buf = (void *) memory->smalloc((bigint)length*size_one, + "imd:recv_force_buf"); } imd_forces = length; - buf = static_cast(force_buf); + buf = static_cast(recv_force_buf); /* compare data to hash table */ for (int ii=0; ii < length; ++ii) { @@ -943,12 +1187,12 @@ void FixIMD::post_force(int /*vflag*/) /* check if we need to readjust the forces comm buffer on the receiving nodes. */ if (me != 0) { if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ - if (force_buf != nullptr) - memory->sfree(force_buf); - force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:force_buf"); + if (recv_force_buf != nullptr) + memory->sfree(recv_force_buf); + recv_force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:recv_force_buf"); } } - MPI_Bcast(force_buf, imd_forces*size_one, MPI_BYTE, 0, world); + MPI_Bcast(recv_force_buf, imd_forces*size_one, MPI_BYTE, 0, world); } /* Check if we need to communicate coordinates to the client. @@ -961,7 +1205,7 @@ void FixIMD::post_force(int /*vflag*/) if (update->ntimestep % imd_trate) { if (imd_forces > 0) { double **f = atom->f; - buf = static_cast(force_buf); + buf = static_cast(recv_force_buf); /* XXX. this is in principle O(N**2) == not good. * however we assume for now that the number of atoms @@ -989,27 +1233,25 @@ void FixIMD::post_force(int /*vflag*/) MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); if (nmax*size_one > maxbuf) { - memory->destroy(comm_buf); + memory->destroy(coord_data); maxbuf = nmax*size_one; - comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); + coord_data = memory->smalloc(maxbuf,"imd:coord_data"); } int tmp, ndata; - buf = static_cast(comm_buf); + buf = static_cast(coord_data); if (me == 0) { MPI_Status status; MPI_Request request; /* collect data into new array. we bypass the IMD API to save * us one extra copy of the data. */ - msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; - msgdata = new char[msglen]; imd_fill_header((IMDheader *)msgdata, IMD_FCOORDS, num_coords); /* array pointer, to the offset where we receive the coordinates. */ auto recvcoord = (float *) (msgdata+IMDHEADERSIZE); /* add local data */ - if (unwrap_flag) { + if (imdsinfo->unwrap) { double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; @@ -1049,10 +1291,9 @@ void FixIMD::post_force(int /*vflag*/) } } } - /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); @@ -1067,7 +1308,6 @@ void FixIMD::post_force(int /*vflag*/) } } } - /* done collecting frame data now communicate with IMD client. */ #if defined(LAMMPS_ASYNC_IMD) @@ -1078,16 +1318,15 @@ void FixIMD::post_force(int /*vflag*/) pthread_mutex_unlock(&write_mutex); #else /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + error->all(FLERR, "LAMMPS terminated on error in sending IMDFrame"); } - delete[] msgdata; #endif } else { /* copy coordinate data into communication buffer */ nme = 0; - if (unwrap_flag) { + if (imdsinfo->unwrap) { double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; @@ -1128,23 +1367,442 @@ void FixIMD::post_force(int /*vflag*/) } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); - MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); } - } - -/* ---------------------------------------------------------------------- */ -void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) -{ - /* only process IMD on the outmost RESPA level. */ - if (ilevel == nlevels_respa-1) post_force(vflag); } /* ---------------------------------------------------------------------- */ -/* local memory usage. approximately. */ -double FixIMD::memory_usage() -{ - return static_cast(num_coords+maxbuf+imd_forces)*size_one; + +void FixIMD::handle_client_input_v3() { + // IMDV3 + + /* check for reconnect */ + if (imd_inactive) { + reconnect(); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on error in setting up IMD connection."); + if (imd_inactive) + return; /* IMD client has detached and not yet come back. do nothing. */ + } + + struct commdata *buf; + int nlocal = atom->nlocal; + int *mask = atom->mask; + tagint *tag = atom->tag; + double **f = atom->f; + + if (me == 0) { + /* process all pending incoming data. */ + int imd_paused=0; + while ((imdsock_selread(clientsock, 0) > 0) || imd_paused) { + /* if something requested to turn off IMD while paused get out */ + if (imd_inactive) break; + + int32 length; + int msg = imd_recv_header(clientsock, &length); + + switch(msg) { + + case IMD_DISCONNECT: { + /* disconnect from client. wait for new connection. */ + imd_paused = 0; + imd_forces = 0; + memory->destroy(recv_force_buf); + recv_force_buf = nullptr; + imdsock_destroy(clientsock); + clientsock = nullptr; + if (screen) + fprintf(screen, "IMD client detached. LAMMPS run continues.\n"); + + connect_msg = 1; + reconnect(); + if (imd_terminate) imd_inactive = 1; + break; + } + + case IMD_KILL: + /* stop the simulation job and shutdown IMD */ + if (screen) + fprintf(screen, "IMD client requested termination of run.\n"); + imd_inactive = 1; + imd_terminate = 1; + imd_paused = 0; + imdsock_destroy(clientsock); + clientsock = nullptr; + break; + + case IMD_PAUSE: + /* pause the running simulation. wait for second IMD_PAUSE to continue. */ + if (!imd_paused) { + if (screen) + fprintf(screen, "Pausing run on IMD client request.\n"); + imd_paused = 1; + } else { + // Pause in IMDv3 is idempotent + continue; + } + break; + + case IMD_RESUME: + /* resume the running simulation. */ + if (imd_paused) { + if (screen) + fprintf(screen, "Continuing run on IMD client request.\n"); + imd_paused = 0; + } else { + // Resume in IMDv3 is idempotent + continue; + } + break; + + case IMD_TRATE: + /* change the IMD transmission data rate */ + if (length > 0) + imd_trate = length; + if (screen) + fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); + break; + + case IMD_MDCOMM: { + auto imd_tags = new int32[length]; + auto imd_fdat = new float[3*length]; + imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); + + if (imd_forces < length) { /* grow holding space for forces, if needed. */ + memory->destroy(recv_force_buf); + recv_force_buf = (void *) memory->smalloc((bigint)length*size_one, + "imd:recv_force_buf"); + } + imd_forces = length; + buf = static_cast(recv_force_buf); + + /* compare data to hash table */ + for (int ii=0; ii < length; ++ii) { + buf[ii].tag = rev_idmap[imd_tags[ii]]; + buf[ii].x = imd_fdat[3*ii]; + buf[ii].y = imd_fdat[3*ii+1]; + buf[ii].z = imd_fdat[3*ii+2]; + } + delete[] imd_tags; + delete[] imd_fdat; + break; + } + case IMD_WAIT: { + /* Change IMD waiting behavior mid-session */ + if (length) { + nowait_flag = 0; + } + else { + nowait_flag = 1; + } + break; + } + + default: + if (screen) + fprintf(screen, "Unhandled incoming IMD message #%d. length=%d\n", msg, length); + break; + } + } + } + + /* update all tasks with current settings. */ + int old_imd_forces = imd_forces; + MPI_Bcast(&imd_trate, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_forces, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on IMD request."); + + if (imd_forces > 0) { + /* check if we need to readjust the forces comm buffer on the receiving nodes. */ + if (me != 0) { + if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ + if (recv_force_buf != nullptr) + memory->sfree(recv_force_buf); + recv_force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:recv_force_buf"); + } + } + MPI_Bcast(recv_force_buf, imd_forces*size_one, MPI_BYTE, 0, world); + } + + /* Check if we need to communicate coordinates to the client. + * Tuning imd_trate allows to keep the overhead for IMD low + * at the expense of a more jumpy display. Rather than using + * end_of_step() we do everything here in one go. + * + * If we don't communicate, only check if we have forces + * stored away and apply them. */ + if (imd_forces > 0) { + buf = static_cast(recv_force_buf); + + /* XXX. this is in principle O(N**2) == not good. + * however we assume for now that the number of atoms + * that we manipulate via IMD will be small compared + * to the total system size, so we don't hurt too much. */ + for (int j=0; j < imd_forces; ++j) { + for (int i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + if (buf[j].tag == tag[i]) { + f[i][0] += imd_fscale*buf[j].x; + f[i][1] += imd_fscale*buf[j].y; + f[i][2] += imd_fscale*buf[j].z; + } + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixIMD::handle_output_v3() { + + tagint *tag = atom->tag; + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + imageint *image = atom->image; + int nlocal = atom->nlocal; + int *mask = atom->mask; + struct commdata *buf; + + // Only main process should use: + float *global_coords = nullptr; + float *global_vel = nullptr; + float *global_force = nullptr; + + // Prepare offsets in outgoing buffer + // Fill what we can wihtout collecting from other processes + if (me == 0) { + int offset = 0; + if (imdsinfo->time) { + imd_fill_header((IMDheader *)msgdata, IMD_TIME, 1); + double dt = update->dt; + + double currtime = update->atime + ((update->ntimestep - update->atimestep) * update->dt); + long long currstep = update->ntimestep; + char *time = (msgdata+IMDHEADERSIZE); + + memcpy(time, &dt, sizeof(dt)); + memcpy(time+sizeof(dt), &currtime, sizeof(currtime)); + memcpy(time+sizeof(dt)+sizeof(currtime), &currstep, sizeof(currstep)); + offset += IMDHEADERSIZE+sizeof(dt)+sizeof(currtime)+sizeof(currstep); + } + if (imdsinfo->box) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_BOX, 1); + // Get triclinic box vectors + float *box = (float *)(msgdata+offset+IMDHEADERSIZE); + box[0] = domain->h[0]; + box[1] = 0.0; + box[2] = 0.0; + box[3] = domain->h[5]; + box[4] = domain->h[1]; + box[5] = 0.0; + box[6] = domain->h[4]; + box[7] = domain->h[3]; + box[8] = domain->h[2]; + + offset += (9*4)+IMDHEADERSIZE; + + } + if (imdsinfo->coords) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_FCOORDS, num_coords); + global_coords = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->velocities) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_VELOCITIES, num_coords); + global_vel = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->forces) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_FORCES, num_coords); + global_force = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + } + + int ntotal, nmax, nme=0; + for (int i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + // Atoms per proc + int *recvcounts = nullptr; + // Displacements in recv for each proc + int *displs = nullptr; + + + if (me == 0) { + recvcounts = new int[comm->nprocs]; + displs = new int[comm->nprocs]; + } + + MPI_Gather(&nme, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, world); + if (me == 0) { + displs[0] = 0; + ntotal = recvcounts[0]; + for (int i = 1; i < comm->nprocs; ++i) { + displs[i] = displs[i - 1] + recvcounts[i - 1]; + ntotal += recvcounts[i]; + } + } + + if (imdsinfo->coords) { + commdata *recvcoord = nullptr; + memory->destroy(coord_data); + coord_data = memory->smalloc(nme*size_one, "imd:coord_data"); + buf = static_cast(coord_data); + int idx = 0; + if (imdsinfo->unwrap) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xy = domain->xy; + double xz = domain->xz; + double yz = domain->yz; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + int ix = (image[i] & IMGMASK) - IMGMAX; + int iy = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + int iz = (image[i] >> IMG2BITS) - IMGMAX; + + if (domain->triclinic) { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + ix * xprd + iy * xy + iz * xz; + buf[idx].y = x[i][1]; + iy * yprd + iz * yz; + buf[idx].z = x[i][2]; + iz * zprd; + } + else { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + ix * xprd; + buf[idx].y = x[i][1]; + iy * yprd; + buf[idx].z = x[i][2]; + iz * zprd; + } + ++idx; + } + } + } + else { + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + buf[idx].y = x[i][1]; + buf[idx].z = x[i][2]; + ++idx; + } + } + } + if (me == 0) { + recvcoord = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvcoord, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_coords + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvcoord[idx].tag); + if (t != 3*HASH_FAIL) { + global_coords[t] = recvcoord[idx].x; + global_coords[t + 1] = recvcoord[idx].y; + global_coords[t + 2] = recvcoord[idx].z; + } + } + } + } + } + if (imdsinfo->velocities) { + commdata *recvvel = nullptr; + memory->destroy(vel_data); + vel_data = memory->smalloc(nme*size_one, "imd:vel_data"); + buf = static_cast(vel_data); + int idx = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = v[i][0]; + buf[idx].y = v[i][1]; + buf[idx].z = v[i][2]; + ++idx; + } + } + if (me == 0) { + recvvel = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvvel, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_vels + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvvel[idx].tag); + if (t != 3*HASH_FAIL) { + global_vel[t] = recvvel[idx].x; + global_vel[t + 1] = recvvel[idx].y; + global_vel[t + 2] = recvvel[idx].z; + } + } + } + } + } + if (imdsinfo->forces) { + commdata *recvforce = nullptr; + memory->destroy(force_data); + force_data = memory->smalloc(nme*size_one, "imd:force_data"); + buf = static_cast(force_data); + int idx = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = f[i][0]; + buf[idx].y = f[i][1]; + buf[idx].z = f[i][2]; + ++idx; + } + } + if (me == 0) { + recvforce = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvforce, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_coords + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvforce[idx].tag); + if (t != 3*HASH_FAIL) { + global_force[t] = recvforce[idx].x; + global_force[t + 1] = recvforce[idx].y; + global_force[t + 2] = recvforce[idx].z; + } + } + } + } + } + +/* done collecting frame data now communicate with IMD client. */ + +#if defined(LAMMPS_ASYNC_IMD) + /* wake up i/o worker thread and release lock on i/o buffer + * we can go back to our MD and let the i/o thread do the rest */ + buf_has_data=1; + pthread_cond_signal(&write_cond); + pthread_mutex_unlock(&write_mutex); +#else + /* send IMDFrame data, blocking until client accepts */ + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + error->all(FLERR, "LAMMPS terminated on error in sending IMDFrame"); + } +#endif } /* End of FixIMD class implementation. */ @@ -1178,6 +1836,7 @@ int imdsock_init() { #endif } +/* ---------------------------------------------------------------------- */ void * imdsock_create() { imdsocket * s; @@ -1196,6 +1855,8 @@ void * imdsock_create() { return (void *) s; } +/* ---------------------------------------------------------------------- */ + int imdsock_bind(void * v, int port) { auto s = (imdsocket *) v; auto *addr = &(s->addr); @@ -1207,11 +1868,15 @@ int imdsock_bind(void * v, int port) { return bind(s->sd, (struct sockaddr *) addr, s->addrlen); } +/* ---------------------------------------------------------------------- */ + int imdsock_listen(void * v) { auto s = (imdsocket *) v; return listen(s->sd, 5); } +/* ---------------------------------------------------------------------- */ + void *imdsock_accept(void * v) { int rc; imdsocket *new_s = nullptr, *s = (imdsocket *) v; @@ -1241,6 +1906,8 @@ void *imdsock_accept(void * v) { return (void *)new_s; } +/* ---------------------------------------------------------------------- */ + int imdsock_write(void * v, const void *buf, int len) { auto s = (imdsocket *) v; #if defined(_MSC_VER) || defined(__MINGW32__) @@ -1250,6 +1917,8 @@ int imdsock_write(void * v, const void *buf, int len) { #endif } +/* ---------------------------------------------------------------------- */ + int imdsock_read(void * v, void *buf, int len) { auto s = (imdsocket *) v; #if defined(_MSC_VER) || defined(__MINGW32__) @@ -1260,6 +1929,8 @@ int imdsock_read(void * v, void *buf, int len) { } +/* ---------------------------------------------------------------------- */ + void imdsock_shutdown(void *v) { auto s = (imdsocket *) v; if (s == nullptr) @@ -1272,6 +1943,8 @@ void imdsock_shutdown(void *v) { #endif } +/* ---------------------------------------------------------------------- */ + void imdsock_destroy(void * v) { auto s = (imdsocket *) v; if (s == nullptr) @@ -1285,6 +1958,8 @@ void imdsock_destroy(void * v) { free(s); } +/* ---------------------------------------------------------------------- */ + int imdsock_selread(void *v, int sec) { auto s = (imdsocket *)v; fd_set rfd; @@ -1304,6 +1979,8 @@ int imdsock_selread(void *v, int sec) { } +/* ---------------------------------------------------------------------- */ + int imdsock_selwrite(void *v, int sec) { auto s = (imdsocket *)v; fd_set wfd; @@ -1339,6 +2016,8 @@ typedef union { } b; } netint; +/* ---------------------------------------------------------------------- */ + static int32 imd_htonl(int32 h) { netint n; n.b.highest = h >> 24; @@ -1348,22 +2027,30 @@ static int32 imd_htonl(int32 h) { return n.i; } +/* ---------------------------------------------------------------------- */ + static int32 imd_ntohl(int32 n) { netint u; u.i = n; return (u.b.highest << 24 | u.b.high << 16 | u.b.low << 8 | u.b.lowest); } +/* ---------------------------------------------------------------------- */ + static void imd_fill_header(IMDheader *header, IMDType type, int32 length) { header->type = imd_htonl((int32)type); header->length = imd_htonl(length); } +/* ---------------------------------------------------------------------- */ + static void swap_header(IMDheader *header) { header->type = imd_ntohl(header->type); header->length= imd_ntohl(header->length); } +/* ---------------------------------------------------------------------- */ + static int32 imd_readn(void *s, char *ptr, int32 n) { int32 nleft; int32 nread; @@ -1383,6 +2070,8 @@ static int32 imd_readn(void *s, char *ptr, int32 n) { return n-nleft; } +/* ---------------------------------------------------------------------- */ + static int32 imd_writen(void *s, const char *ptr, int32 n) { int32 nleft; int32 nwritten; @@ -1401,13 +2090,39 @@ static int32 imd_writen(void *s, const char *ptr, int32 n) { return n; } -int imd_handshake(void *s) { +/* ---------------------------------------------------------------------- */ + +int imd_handshake_v2(void *s) { IMDheader header; imd_fill_header(&header, IMD_HANDSHAKE, 1); - header.length = IMDVERSION; /* Not byteswapped! */ + header.length = 2; /* Not byteswapped! */ return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); } +/* ---------------------------------------------------------------------- */ + +int imd_handshake_v3(void *s, IMDSessionInfo *imdsinfo) { + IMDheader header; + imd_fill_header(&header, IMD_HANDSHAKE, 1); + header.length = 3; /* Not byteswapped so client can determine native endinaness */ + + if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE) return -1; + + imd_fill_header(&header, IMD_SESSIONINFO, 7); + unsigned char body[7] = {0}; + body[0] = imdsinfo->time; + body[1] = imdsinfo->energies; + body[2] = imdsinfo->box; + body[3] = imdsinfo->coords; + body[4] = !imdsinfo->unwrap; + body[5] = imdsinfo->velocities; + body[6] = imdsinfo->forces; + + if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE || + imd_writen(s, (char *)&body, 7) != 7) return -1; + return 0; +} + /* The IMD receive functions */ IMDType imd_recv_header(void *s, int32 *length) { @@ -1419,17 +2134,23 @@ IMDType imd_recv_header(void *s, int32 *length) { return IMDType(header.type); } +/* ---------------------------------------------------------------------- */ + int imd_recv_mdcomm(void *s, int32 n, int32 *indices, float *forces) { if (imd_readn(s, (char *)indices, 4*n) != 4*n) return 1; if (imd_readn(s, (char *)forces, 12*n) != 12*n) return 1; return 0; } +/* ---------------------------------------------------------------------- */ + int imd_recv_energies(void *s, IMDEnergies *energies) { return (imd_readn(s, (char *)energies, sizeof(IMDEnergies)) != sizeof(IMDEnergies)); } +/* ---------------------------------------------------------------------- */ + int imd_recv_fcoords(void *s, int32 n, float *coords) { return (imd_readn(s, (char *)coords, 12*n) != 12*n); } diff --git a/src/MISC/fix_imd.h b/src/MISC/fix_imd.h index 2312540e66..3d35376868 100644 --- a/src/MISC/fix_imd.h +++ b/src/MISC/fix_imd.h @@ -56,6 +56,9 @@ FixStyle(imd,FixIMD); #include #endif +/* IMDv3 session information */ +struct IMDSessionInfo; + /* prototype for c wrapper that calls the real worker */ extern "C" void *fix_imd_ioworker(void *); @@ -69,8 +72,12 @@ class FixIMD : public Fix { void init() override; void setup(int) override; void post_force(int) override; + void end_of_step() override; void post_force_respa(int, int, int) override; double memory_usage() override; + // Fix nevery at 1, use trate to skip in 'end_of_step` + int nevery = 1; + int imd_version; // version of the IMD protocol to be used. protected: int imd_port; @@ -80,13 +87,15 @@ class FixIMD : public Fix { int num_coords; // total number of atoms controlled by this fix int size_one; // bytes per atom in communication buffer. int maxbuf; // size of atom communication buffer. - void *comm_buf; // communication buffer + void *coord_data; // communication buffer for coordinates + void *vel_data; // communication buffer for velocities + void *force_data; // communication buffer for forces void *idmap; // hash for mapping atom indices to consistent order. tagint *rev_idmap; // list of the hash keys for reverse mapping. - int imd_forces; // number of forces communicated via IMD. - void *force_buf; // force data buffer - double imd_fscale; // scale factor for forces. in case VMD's units are off. + int imd_forces; // number of forces communicated via IMD. + void *recv_force_buf; // force data buffer + double imd_fscale; // scale factor for forces. in case VMD's units are off. int imd_inactive; // true if IMD connection stopped. int imd_terminate; // true if IMD requests termination of run. @@ -96,12 +105,22 @@ class FixIMD : public Fix { int nowait_flag; // true if LAMMPS should not wait with the execution for VMD. int connect_msg; // flag to indicate whether a "listen for connection message" is needed. + /* IMDv3-only */ + IMDSessionInfo *imdsinfo; // session information for IMDv3 + int me; // my MPI rank in this "world". int nlevels_respa; // flag to determine respa levels. int msglen; char *msgdata; + private: + void setup_v2(); + void setup_v3(); + void handle_step_v2(); + void handle_client_input_v3(); + void handle_output_v3(); + #if defined(LAMMPS_ASYNC_IMD) int buf_has_data; // flag to indicate to the i/o thread what to do. pthread_mutex_t write_mutex; // mutex for sending coordinates to i/o thread diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp index cd331b5f30..bbfb4d8a0b 100644 --- a/src/QEQ/fix_qeq.cpp +++ b/src/QEQ/fix_qeq.cpp @@ -41,6 +41,9 @@ using namespace LAMMPS_NS; using namespace FixConst; static constexpr double QSUMSMALL = 0.00001; +static constexpr int MIN_CAP = 50; +static constexpr double SAFE_ZONE = 1.2; +static constexpr bigint MIN_NBRS = 100; namespace { class qeq_parser_error : public std::exception { diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h index 5d0e05c7c5..0124a6ee8a 100644 --- a/src/QEQ/fix_qeq.h +++ b/src/QEQ/fix_qeq.h @@ -16,12 +16,6 @@ #include "fix.h" -#define EV_TO_KCAL_PER_MOL 14.4 -#define DANGER_ZONE 0.90 -#define MIN_CAP 50 -#define SAFE_ZONE 1.2 -#define MIN_NBRS 100 - namespace LAMMPS_NS { class FixQEq : public Fix { @@ -36,6 +30,7 @@ class FixQEq : public Fix { void min_pre_force(int) override; double compute_scalar() override; + static constexpr double DANGER_ZONE = 0.90; // derived child classes must provide these functions diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp index b4827d8694..8d1ec7f87e 100644 --- a/src/QEQ/fix_qeq_shielded.cpp +++ b/src/QEQ/fix_qeq_shielded.cpp @@ -33,6 +33,8 @@ using namespace LAMMPS_NS; +static constexpr double EV_TO_KCAL_PER_MOL = 14.4; + /* ---------------------------------------------------------------------- */ FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg) diff --git a/src/finish.cpp b/src/finish.cpp index c774f51e58..3b824f2ac4 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -121,9 +121,9 @@ void Finish::end(int flag) MPI_Allreduce(&cpu_loop,&tmp,1,MPI_DOUBLE,MPI_SUM,world); cpu_loop = tmp/nprocs; if (time_loop > 0.0) cpu_loop = cpu_loop/time_loop*100.0; + output->thermo->footer(); if (me == 0) { - output->thermo->footer(); int ntasks = nprocs * nthreads; utils::logmesg(lmp,"Loop time of {:.6g} on {} procs for {} steps with {} atoms\n\n", time_loop,ntasks,update->nsteps,atom->natoms); diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index d6630bcbf0..c9d09fa8b7 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -208,14 +208,14 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : FixDeposit::~FixDeposit() { delete random; - delete [] molfrac; - delete [] idrigid; - delete [] idshake; - delete [] idregion; - delete [] vstr; - delete [] xstr; - delete [] ystr; - delete [] zstr; + delete[] molfrac; + delete[] idrigid; + delete[] idshake; + delete[] idregion; + delete[] vstr; + delete[] xstr; + delete[] ystr; + delete[] zstr; memory->destroy(coords); memory->destroy(imageflags); } @@ -737,7 +737,7 @@ void FixDeposit::options(int narg, char **arg) mode = MOLECULE; onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; - delete [] molfrac; + delete[] molfrac; molfrac = new double[nmol]; molfrac[0] = 1.0/nmol; for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol; @@ -755,13 +755,13 @@ void FixDeposit::options(int narg, char **arg) iarg += nmol+1; } else if (strcmp(arg[iarg],"rigid") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - delete [] idrigid; + delete[] idrigid; idrigid = utils::strdup(arg[iarg+1]); rigidflag = 1; iarg += 2; } else if (strcmp(arg[iarg],"shake") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - delete [] idshake; + delete[] idshake; idshake = utils::strdup(arg[iarg+1]); shakeflag = 1; iarg += 2; diff --git a/src/library.cpp b/src/library.cpp index ce0842ec2e..71243611a2 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -516,6 +516,9 @@ treated as part of the command and will **not** start a second command. The function returns the expanded string in a new string buffer that must be freed with :cpp:func:`lammps_free` after use to avoid a memory leak. +*See also* + :cpp:func:`lammps_eval` + \endverbatim * * \param handle pointer to a previously created LAMMPS instance @@ -850,7 +853,11 @@ This differs from :cpp:func:`lammps_get_thermo` in that it does **not** trigger an evaluation. Instead it provides direct access to a read-only location of the last thermo output data and the corresponding keyword strings. How to handle the return value depends on the value of the *what* -argument string. +argument string. When accessing the data from a concurrent thread while +LAMMPS is running, the cache needs to be locked first and then unlocked +after the data is obtained, so that the data is not corrupted while +reading in case LAMMPS wants to update it at the same time. Outside +of a run, the lock/unlock calls have no effect. .. note:: @@ -899,6 +906,14 @@ argument string. - actual field data for column - pointer to int, int64_t or double - yes + * - lock + - acquires lock to thermo data cache + - NULL pointer + - no + * - unlock + - releases lock to thermo data cache + - NULL pointer + - no \endverbatim * @@ -954,8 +969,14 @@ void *lammps_last_thermo(void *handle, const char *what, int index) } else if (field.type == multitype::LAMMPS_DOUBLE) { val = (void *) &field.data.d; } - + } else if (strcmp(what, "lock") == 0) { + th->lock_cache(); + val = nullptr; + } else if (strcmp(what, "unlock") == 0) { + th->unlock_cache(); + val = nullptr; } else val = nullptr; + } END_CAPTURE return val; @@ -2906,6 +2927,41 @@ int lammps_variable_info(void *handle, int idx, char *buffer, int buf_size) { /* ---------------------------------------------------------------------- */ +/** Evaluate an immediate variable expression + * +\verbatim embed:rst + +.. versionadded:: TBD + +This function takes a string with an expression that can be used +for :doc:`equal style variables `, evaluates it and returns +the resulting (scalar) value as a floating point number. + +*See also* + :cpp:func:`lammps_expand` + +\endverbatim + + * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. + * \param expr string with expression + * \return result from expression */ + +double lammps_eval(void *handle, const char *expr) +{ + auto lmp = (LAMMPS *) handle; + double result = 0.0; + + BEGIN_CAPTURE + { + result = lmp->input->variable->compute_equal(expr); + } + END_CAPTURE + + return result; +} + +/* ---------------------------------------------------------------------- */ + /** Clear whether a compute has been invoked. * \verbatim embed:rst diff --git a/src/library.h b/src/library.h index b1213b67c8..54600ed422 100644 --- a/src/library.h +++ b/src/library.h @@ -189,6 +189,7 @@ int lammps_set_variable(void *handle, const char *name, const char *str); int lammps_set_string_variable(void *handle, const char *name, const char *str); int lammps_set_internal_variable(void *handle, const char *name, double value); int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +double lammps_eval(void *handle, const char *expr); void lammps_compute_clearstep(void * handle); #if defined(LAMMPS_SMALLSMALL) diff --git a/src/thermo.cpp b/src/thermo.cpp index b2e4f7f934..35b5016118 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -100,8 +100,8 @@ static char fmtbuf[512]; /* ---------------------------------------------------------------------- */ Thermo::Thermo(LAMMPS *_lmp, int narg, char **arg) : - Pointers(_lmp), style(nullptr), vtype(nullptr), field2index(nullptr), argindex1(nullptr), - argindex2(nullptr), temperature(nullptr), pressure(nullptr), pe(nullptr) + Pointers(_lmp), style(nullptr), vtype(nullptr), cache_mutex(nullptr), field2index(nullptr), + argindex1(nullptr), argindex2(nullptr), temperature(nullptr), pressure(nullptr), pe(nullptr) { style = utils::strdup(arg[0]); @@ -208,6 +208,7 @@ void Thermo::init() ValueTokenizer *format_line = nullptr; if (format_line_user.size()) format_line = new ValueTokenizer(format_line_user); + lock_cache(); field_data.clear(); field_data.resize(nfield); std::string format_this, format_line_user_def; @@ -277,6 +278,7 @@ void Thermo::init() format[i] += fmt::format("{:<8} = {} ", keyword[i], format_this); } } + unlock_cache(); // chop off trailing blank or add closing bracket if needed and then add newline if (lineflag == ONELINE) @@ -320,6 +322,9 @@ void Thermo::init() if (index_press_scalar >= 0) pressure = computes[index_press_scalar]; if (index_press_vector >= 0) pressure = computes[index_press_vector]; if (index_pe >= 0) pe = computes[index_pe]; + + // create mutex to protect access to cached thermo data + cache_mutex = new std::mutex; } /* ---------------------------------------------------------------------- @@ -366,9 +371,17 @@ void Thermo::header() /* ---------------------------------------------------------------------- */ +// called at the end of a run from Finish class + void Thermo::footer() { - if (lineflag == YAMLLINE) utils::logmesg(lmp, "...\n"); + if (comm->me == 0) { + if (lineflag == YAMLLINE) utils::logmesg(lmp, "...\n"); + } + + // no more locking for cached thermo data access needed + delete cache_mutex; + cache_mutex = nullptr; } /* ---------------------------------------------------------------------- */ @@ -422,6 +435,7 @@ void Thermo::compute(int flag) } // add each thermo value to line with its specific format + lock_cache(); field_data.clear(); field_data.resize(nfield); @@ -441,6 +455,7 @@ void Thermo::compute(int flag) field_data[ifield] = bivalue; } } + unlock_cache(); // print line to screen and logfile @@ -579,7 +594,8 @@ void Thermo::modify_params(int narg, char **arg) if (iarg + 2 > narg) error->all(FLERR, "Illegal thermo_modify command"); triclinic_general = utils::logical(FLERR, arg[iarg + 1], false, lmp); if (triclinic_general && !domain->triclinic_general) - error->all(FLERR,"Thermo_modify triclinic/general cannot be used " + error->all(FLERR, + "Thermo_modify triclinic/general cannot be used " "if simulation box is not general triclinic"); iarg += 2; @@ -1566,6 +1582,26 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) return 0; } +/* ---------------------------------------------------------------------- */ + +// lock cache for current thermo data + +void Thermo::lock_cache() +{ + // no locking outside of a run + if (!cache_mutex) return; + cache_mutex->lock(); +} + +// unlock cache for current thermo data + +void Thermo::unlock_cache() +{ + // no locking outside of a run + if (!cache_mutex) return; + cache_mutex->unlock(); +} + /* ---------------------------------------------------------------------- extraction of Compute, Fix, Variable results compute/fix are normalized by atoms if returning extensive value diff --git a/src/thermo.h b/src/thermo.h index 52f69e7609..f5466e70cd 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -16,6 +16,7 @@ #include "pointers.h" #include +#include namespace LAMMPS_NS { @@ -43,6 +44,8 @@ class Thermo : protected Pointers { int evaluate_keyword(const std::string &, double *); // for accessing cached thermo and related data + void lock_cache(); + void unlock_cache(); const int *get_line() const { return &nline; } const char *get_image_fname() const { return image_fname.c_str(); } @@ -82,6 +85,9 @@ class Thermo : protected Pointers { int nline; std::string image_fname; + // mutex for locking the cache + std::mutex *cache_mutex; + // data used by routines that compute single values int ivalue; // integer value to print diff --git a/src/variable.cpp b/src/variable.cpp index 135882fde9..031709166b 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -457,15 +457,21 @@ void Variable::set(int narg, char **arg) // data = 2 values, 1st is string to eval, 2nd is filled on retrieval } else if (strcmp(arg[1],"equal") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable equal", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != EQUAL) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); replaceflag = 1; } else { if (nvar == maxvar) grow(); @@ -474,7 +480,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); data[nvar][1] = new char[VALUELENGTH]; strcpy(data[nvar][1],"(undefined)"); } @@ -485,15 +491,21 @@ void Variable::set(int narg, char **arg) // data = 1 value, string to eval } else if (strcmp(arg[1],"atom") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable atom", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != ATOM) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); replaceflag = 1; } else { if (nvar == maxvar) grow(); @@ -502,7 +514,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); } // VECTOR @@ -513,16 +525,22 @@ void Variable::set(int narg, char **arg) // immediately store it as N-length vector and set dynamic flag to 0 } else if (strcmp(arg[1],"vector") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable atom", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != VECTOR) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; delete[] data[ivar][1]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); if (data[ivar][0][0] != '[') vecs[ivar].dynamic = 1; else { @@ -539,7 +557,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); if (data[nvar][0][0] != '[') { vecs[nvar].dynamic = 1; data[nvar][1] = nullptr; diff --git a/tools/lammps-gui/lammps-gui.appdata.xml b/tools/lammps-gui/lammps-gui.appdata.xml index abf0a4bf45..c019ab1ce3 100644 --- a/tools/lammps-gui/lammps-gui.appdata.xml +++ b/tools/lammps-gui/lammps-gui.appdata.xml @@ -61,6 +61,8 @@ Highlight warnings and error messages in Output window Make Tutorial wizards more compact Include download and compilation of WHAM software from Alan Grossfield + Add dialog to run WHAM directly from LAMMPS-GUI + Use mutex to avoid corruption of thermo data diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index 8749e24821..59e8017fc0 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -1015,6 +1015,7 @@ void LammpsGui::logupdate() void *ptr = lammps.last_thermo("setup", 0); if (ptr && *(int *)ptr) return; + lammps.last_thermo("lock", 0); ptr = lammps.last_thermo("num", 0); if (ptr) { int ncols = *(int *)ptr; @@ -1066,6 +1067,7 @@ void LammpsGui::logupdate() chartwindow->add_data(step, data, i); } } + lammps.last_thermo("unlock", 0); } // update list of available image file names diff --git a/tools/lammps-gui/update-wham-2.0.11.patch b/tools/lammps-gui/update-wham-2.0.11.patch index 3b6e0e372a..b884c25336 100644 --- a/tools/lammps-gui/update-wham-2.0.11.patch +++ b/tools/lammps-gui/update-wham-2.0.11.patch @@ -55,8 +55,79 @@ index 0000000..b4f0fe6 + TYPE DOC + PERMISSIONS OWNER_READ GROUP_READ WORLD_READ +) +diff --git a/nr/locate.c b/nr/locate.c +index 9f92dc0..f3bf294 100644 +--- a/nr/locate.c ++++ b/nr/locate.c +@@ -11,7 +11,7 @@ void locate(double xx[], int n, double x, int *j) + ascnd=(xx[n] > xx[0]); // I think this makes it zero based + while (ju-jl > 1) { + jm=(ju+jl) >> 1; +- if (x > xx[jm] == ascnd) ++ if ((x > xx[jm]) == ascnd) + jl=jm; + else + ju=jm; +diff --git a/wham-2d/histogram.c b/wham-2d/histogram.c +index 1bd1329..b5d1c01 100644 +--- a/wham-2d/histogram.c ++++ b/wham-2d/histogram.c +@@ -78,14 +78,14 @@ return hp; + } + + /* Get a value from a histogram structure +- * Given i and j, the indices into the global histogram, return ++ * Given i and j, the indices into the global histogram, return + * the correct histogram value. Since we only store range of histogram +- * indices containing the nonzero values, we have to check the index value ++ * indices containing the nonzero values, we have to check the index value + * against the range in the structure. + */ + double get_histval(struct histogram *hist, int i, int j) + { +-if ( (i < hist->first_x) || (i > hist->last_x) ++if ( (i < hist->first_x) || (i > hist->last_x) + ||(j < hist->first_y) || (j > hist->last_y)) + { + return 0.0; +@@ -239,7 +239,7 @@ return error; + + + // Calculate the free energy, setting the minimum value to 0 +-void calc_free(double **free, double **prob, double kT, ++void calc_free(double **free, double **prob, double kT, + int use_mask, int **mask) + { + int i,j; +@@ -257,7 +257,14 @@ for (i=0; i 0.0) ++ { ++ free[i][j] = -kT * log(prob[i][j]); ++ } ++ else ++ { ++ free[i][j] = 0.0; ++ } + if (free[i][j] < min) + { + min = free[i][j]; +@@ -321,8 +328,6 @@ return 0.5*(springx*(dx*dx) + springy*(dy*dy)); + + void calc_coor(int i, int j, double *coor) + { +-coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); +-coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); ++coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); ++coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); + } +- +- diff --git a/wham-2d/wham-2d.c b/wham-2d/wham-2d.c -index fb6e059..2c5594f 100644 +index fb6e059..a6b8483 100644 --- a/wham-2d/wham-2d.c +++ b/wham-2d/wham-2d.c @@ -25,7 +25,7 @@ @@ -77,6 +148,15 @@ index fb6e059..2c5594f 100644 int main(int argc, char *argv[]) { +@@ -57,7 +57,7 @@ double sum; + int iteration; + int max_iteration = 100000; + int numpad; +-int **mask; ++int **mask = NULL; + int use_mask; + + cpu1 = ((double) clock())/CLOCKS_PER_SEC; @@ -76,6 +76,61 @@ for (i=0; i 0.0) ++ { ++ free[i] = -kT * log(prob[i]); ++ } ++ else ++ { ++ free[i] = 0.0; ++ } ++ ++ if (free[i] < min) + { + min = free[i]; + bin_min = i; +@@ -252,5 +260,5 @@ return 0.5*dx*dx*spring; + + double calc_coor(int i) + { +-return HIST_MIN + BIN_WIDTH*((double)i+0.5); ++return HIST_MIN + BIN_WIDTH*((double)i+0.5); + } diff --git a/wham/wham.c b/wham/wham.c -index 487871b..1496eed 100644 +index 487871b..edb8125 100644 --- a/wham/wham.c +++ b/wham/wham.c +@@ -1,4 +1,4 @@ +-/* ++/* + * WHAM -- Weighted Histogram Analysis Method + * + * Reference 1: Computer Physics Communications, 91(1995) 275-282, Benoit Roux @@ -21,7 +21,7 @@ #include "wham.h" @@ -184,6 +458,15 @@ index 487871b..1496eed 100644 int main(int argc, char *argv[]) { +@@ -41,7 +42,7 @@ int first; + int bin_min; + int have_energy; + char *freefile; +-FILE *METAFILE, *FREEFILE; ++FILE *METAFILE, *FREEFILE; + struct hist_group *hist_group; + struct histogram *hp; + double coor; @@ -82,6 +83,61 @@ for (i=0; inum_windows); + // Figure out if we have trajectories at different temperatures. + // Missing temperatures are set to -1 in read_metadata, and + // since we require that either all trajectories specify a temperature +-// or all trajectories are assumed to be at the wham temperature, we only ++// or all trajectories are assumed to be at the wham temperature, we only + // have to check one of them + if (hist_group->kT[0] > 0) + { +@@ -257,7 +313,7 @@ if (hist_group->kT[0] > 0) + else + { + have_energy = 0; +- for (i=0; i< hist_group->num_windows; i++) ++ for (i=0; i< hist_group->num_windows; i++) + { + hist_group->kT[i] = kT; + } +@@ -269,7 +325,7 @@ if (!final_f) + { + printf("couldn't allocate space for final_f: %s\n", strerror(errno)); + exit(errno); +- } ++ } + + free(HISTOGRAM); + +@@ -305,7 +361,8 @@ while (! is_converged(hist_group) || first ) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); + } + printf("\n"); + +@@ -319,7 +376,7 @@ while (! is_converged(hist_group) || first ) + } + } + // Cheesy bailout if we're going on too long +- if (iteration >= max_iteration) ++ if (iteration >= max_iteration) + { + printf("Too many iterations: %d\n", iteration); + break; +@@ -383,11 +440,11 @@ for (i=0; i< num_mc_runs; i++) + //printf("Faking %d: %d %d\n", i,j,hp->num_points); + num_used = hp->last - hp->first + 1; + mk_new_hist(hp->cum, hp->data, num_used, hp->num_mc_samples, &idum); +- ++ + hist_group->F_old[j] = 0.0; + hist_group->F[j] = 0.0; + } +- ++ + // perform WHAM iterations on the fake data sets + iteration = 0; + first = 1; +@@ -403,7 +460,7 @@ for (i=0; i< num_mc_runs; i++) + printf("Too many iterations: %d\n", iteration); + break; + } +- } ++ } + printf("#MC trial %d: %d iterations\n", i, iteration); + printf("#PMF values\n"); + // accumulate the average and stdev of the resulting probabilities +@@ -419,18 +476,19 @@ for (i=0; i< num_mc_runs; i++) + for (j=0; j < NUM_BINS; j++) + { + pdf = -kT*log(prob[j]); +- ++ + ave_p[j] += prob[j]; + ave_pdf[j] += pdf; + ave_p2[j] += prob[j] * prob[j]; + ave_pdf2[j] += pdf*pdf; + } +- for (j=0; jnum_windows;j++) ++ for (j=0; jnum_windows;j++) + { + ave_F[j] += hist_group->F[j] - hist_group->F[0]; +- ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; ++ ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; + } +- } ++ } ++ if (num_mc_runs == 0) num_mc_runs = 1; + for (i=0; i < NUM_BINS; i++) + { + ave_p[i] /= (double)num_mc_runs; +@@ -457,12 +515,12 @@ if (!FREEFILE) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], +- prob[i], ave_p2[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], prob[i], ave_p2[i]); + } + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ printf("%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + + exit(errno); +@@ -470,38 +528,37 @@ if (!FREEFILE) + else + { + // write out header +- fprintf(FREEFILE, "#Coor\t\tFree\t+/-\t\tProb\t\t+/-\n"); ++ fprintf(FREEFILE, "#Coor\t\tFree\t\t+/-\t\tProb\t\t+/-\n"); + // write out the leading padded values + for (i=-numpad; i<0; i++) + { + coor = calc_coor(i); +- fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], +- ave_pdf2[NUM_BINS+i], +- final_prob[NUM_BINS+i], +- ave_p2[NUM_BINS+i]); ++ if (prob[NUM_BINS+1] != 0.0) ++ fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], ++ ave_pdf2[NUM_BINS+i], final_prob[NUM_BINS+i], ave_p2[NUM_BINS+i]); + } + // write out the center values + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"#%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ fprintf(FREEFILE,"#%d\t\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + } + +@@ -515,7 +572,7 @@ exit(0); + /* + * Perform a single WHAM iteration + */ +-void wham_iteration(struct hist_group* hist_group, double *prob, ++void wham_iteration(struct hist_group* hist_group, double *prob, + int have_energy) + { + int i,j; +@@ -535,9 +592,9 @@ for (i=0; ihists[j]),i); + bias = calc_bias(hist_group,j,coor); + bf = exp((hist_group->F_old[j] - bias) / hist_group->kT[j]); +- /* If we have energies, then the histograms contain boltzmann +- * factors. If not, they contain integer counts. Accordingly, +- * we either use a partition function to normalize, or simply the ++ /* If we have energies, then the histograms contain boltzmann ++ * factors. If not, they contain integer counts. Accordingly, ++ * we either use a partition function to normalize, or simply the + * number of points. + */ + if (have_energy) diff --git a/wham/wham.h b/wham/wham.h index aacc1e8..7d509f2 100644 --- a/wham/wham.h diff --git a/tools/swig/lammps.i b/tools/swig/lammps.i index 119edc4367..bc9f46f407 100644 --- a/tools/swig/lammps.i +++ b/tools/swig/lammps.i @@ -141,6 +141,8 @@ extern int lammps_extract_variable_datatype(void *handle, const char *name); extern int lammps_set_variable(void *, const char *, const char *); extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); +extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +extern double lammps_eval(void *handle, const char *expr); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); @@ -332,6 +334,8 @@ extern int lammps_extract_variable_datatype(void *handle, const char *name); extern int lammps_set_variable(void *, const char *, const char *); extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); +extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +extern double lammps_eval(void *handle, const char *expr); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); diff --git a/unittest/c-library/test_library_objects.cpp b/unittest/c-library/test_library_objects.cpp index a175952dfd..4a51381e64 100644 --- a/unittest/c-library/test_library_objects.cpp +++ b/unittest/c-library/test_library_objects.cpp @@ -233,3 +233,34 @@ TEST_F(LibraryObjects, expand) EXPECT_THAT((char *)ptr, StrEq("'xx_$(4+5)_$(PI) ${one}-${two}-${three}'")); lammps_free(ptr); } + +TEST_F(LibraryObjects, eval) +{ + lammps_get_last_error_message(lmp, nullptr, 1); + ::testing::internal::CaptureStdout(); + lammps_commands_string(lmp, "region box1 block 0 10 0 5 -0.5 0.5\n" + "lattice fcc 0.8\n" + "create_box 1 box1\n" + "create_atoms 1 box\n" + "mass * 1.0\n" + "pair_style lj/cut 4.0\n" + "pair_coeff * * 1.0 1.0\n" + "variable t equal 15.0\n" + "velocity all create 1.5 532656\n" + "fix 1 all nve\n" + "run 0 post no\n"); + lammps_command(lmp, "variable one index 1 2 3 4"); + lammps_command(lmp, "variable two equal 2"); + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + ASSERT_EQ(lammps_has_error(lmp), 0); + + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "4+5"), 9.0); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "v_one / 2.0"), 0.5); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "count(all)"), 36.0); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "pe"), -3.9848867644689534); + EXPECT_EQ(lammps_has_error(lmp), 0); +} diff --git a/unittest/commands/test_variables.cpp b/unittest/commands/test_variables.cpp index 466dac17cc..0d8d28426c 100644 --- a/unittest/commands/test_variables.cpp +++ b/unittest/commands/test_variables.cpp @@ -140,8 +140,8 @@ TEST_F(VariableTest, CreateDelete) command("variable ten1 universe 1 2 3 4"); command("variable ten2 uloop 4"); command("variable ten3 uloop 4 pad"); - command("variable ten4 vector [0,1,2,3,5,7,11]"); - command("variable ten5 vector [0.5,1.25]"); + command("variable ten4 vector [0,1, 2,3, 5,7,11]"); + command("variable ten5 vector [ 0.5, 1.25 ]"); command("variable dummy index 0"); command("variable file equal is_file(MYFILE)"); command("variable iswin equal is_os(^Windows)"); @@ -323,13 +323,13 @@ TEST_F(VariableTest, Expressions) BEGIN_HIDE_OUTPUT(); command("variable one index 1"); command("variable two equal 2"); - command("variable three equal v_one+v_two"); + command("variable three equal v_one + v_two"); command("variable four equal PI"); command("variable five equal version"); command("variable six equal XXX"); command("variable seven equal -v_one"); command("variable eight equal v_three-0.5"); - command("variable nine equal v_two*(v_one+v_three)"); + command("variable nine equal v_two * (v_one+v_three)"); command("variable ten equal (1.0/v_two)^2"); command("variable eleven equal v_three%2"); command("variable twelve equal 1==2"); @@ -341,7 +341,7 @@ TEST_F(VariableTest, Expressions) command("variable ten8 equal 1|^0"); command("variable ten9 equal v_one-v_ten9"); command("variable ten10 internal 100.0"); - command("variable ten11 equal (1!=1)+(2<1)+(2<=1)+(1>2)+(1>=2)+(1&&0)+(0||0)+(1|^1)+10^0"); + command("variable ten11 equal (1 != 1)+(2 < 1)+(2<=1)+(1>2)+(1>=2)+(1&&0)+(0||0)+(1|^1)+10^0"); command("variable ten12 equal yes+no+on+off+true+false"); command("variable err1 equal v_one/v_ten7"); command("variable err2 equal v_one%v_ten7"); @@ -350,7 +350,7 @@ TEST_F(VariableTest, Expressions) command("variable vec2 vector v_vec1*0.5"); command("variable vec3 equal v_vec2[3]"); command("variable vec4 vector '[1, 5, 2.5, -10, -5, 20, 120, 4, 3, 3]'"); - command("variable sort vector sort(v_vec4)"); + command("variable sort vector sort(v_vec4 )"); command("variable rsrt vector rsort(v_vec4)"); command("variable max2 equal sort(v_vec4)[2]"); command("variable rmax equal rsort(v_vec4)[1]"); diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py index 2066514318..06dcacd412 100644 --- a/unittest/python/python-commands.py +++ b/unittest/python/python-commands.py @@ -538,6 +538,27 @@ create_atoms 1 single & expanded = self.lmp.expand("'xx_$(4+5)_$(PI) ${one}-${two}-${three}'") self.assertEqual(expanded, "'xx_$(4+5)_$(PI) ${one}-${two}-${three}'") + def test_eval(self): + self.lmp.commands_string( + """region box1 block 0 10 0 5 -0.5 0.5 + lattice fcc 0.8 + create_box 1 box1 + create_atoms 1 box + mass * 1.0 + pair_style lj/cut 4.0 + pair_coeff * * 1.0 1.0 + variable t equal 15.0 + velocity all create 1.5 532656 + fix 1 all nve + run 0 post no""") + self.lmp.command("variable one index 1 2 3 4") + self.lmp.command("variable two equal 2") + + self.assertEqual(self.lmp.eval("4+5"), 9.0) + self.assertEqual(self.lmp.eval("v_one / 2.0"), 0.5) + self.assertEqual(self.lmp.eval("count(all)"), 36.0) + self.assertEqual(self.lmp.eval("pe"), -3.9848867644689534) + def test_get_thermo(self): self.lmp.command("units lj") self.lmp.command("atom_style atomic")