From b21479e20a94e27986282a6e97ca1503abd8a9c4 Mon Sep 17 00:00:00 2001 From: Taylor Barnes Date: Mon, 14 May 2018 14:03:28 -0400 Subject: [PATCH] Added initial MolSSI Driver interface --- cmake/CMakeLists.txt | 6 +- cmake/Modules/Packages/USER-MDI.cmake | 66 ++ lib/mdi/.gitignore | 69 ++ lib/mdi/Install.py | 210 ++++++ lib/mdi/Makefile.g++ | 18 + lib/mdi/Makefile.gcc | 18 + lib/mdi/Makefile.icc | 18 + lib/mdi/Makefile.lammps.empty | 5 + lib/mdi/Makefile.mpi | 18 + src/.gitignore | 4 + src/Makefile | 9 +- src/USER-MDI/Install.sh | 68 ++ src/USER-MDI/fix_mdi_engine.cpp | 926 ++++++++++++++++++++++++++ src/USER-MDI/fix_mdi_engine.h | 137 ++++ src/USER-MDI/mdi_engine.cpp | 356 ++++++++++ src/USER-MDI/mdi_engine.h | 52 ++ src/main.cpp | 24 +- src/mdi_interface.h | 34 + 18 files changed, 2027 insertions(+), 11 deletions(-) create mode 100644 cmake/Modules/Packages/USER-MDI.cmake create mode 100644 lib/mdi/.gitignore create mode 100644 lib/mdi/Install.py create mode 100644 lib/mdi/Makefile.g++ create mode 100644 lib/mdi/Makefile.gcc create mode 100644 lib/mdi/Makefile.icc create mode 100644 lib/mdi/Makefile.lammps.empty create mode 100644 lib/mdi/Makefile.mpi create mode 100755 src/USER-MDI/Install.sh create mode 100644 src/USER-MDI/fix_mdi_engine.cpp create mode 100644 src/USER-MDI/fix_mdi_engine.h create mode 100644 src/USER-MDI/mdi_engine.cpp create mode 100644 src/USER-MDI/mdi_engine.h create mode 100644 src/mdi_interface.h diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index b3dcd36913..b1cbf33c41 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -121,7 +121,7 @@ set(STANDARD_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE PLUGIN QEQ REPLICA RIGID SHOCK SPIN SNAP SRD KIM PYTHON MSCG MPIIO VORONOI USER-ADIOS USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD - USER-LB USER-MANIFOLD USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF + USER-LB USER-MANIFOLD USER-MDI USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REACTION USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF USER-PACE USER-BROWNIAN) @@ -324,8 +324,8 @@ else() set(CUDA_REQUEST_PIC) endif() -foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM - USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE) +foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MDI USER-MOLFILE USER-NETCDF USER-PLUMED + USER-QMMM USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE) if(PKG_${PKG_WITH_INCL}) include(Packages/${PKG_WITH_INCL}) endif() diff --git a/cmake/Modules/Packages/USER-MDI.cmake b/cmake/Modules/Packages/USER-MDI.cmake new file mode 100644 index 0000000000..b2030cc29a --- /dev/null +++ b/cmake/Modules/Packages/USER-MDI.cmake @@ -0,0 +1,66 @@ +find_package(mdi) +if(${mdi_FOUND}) + set(DOWNLOAD_MDI_DEFAULT OFF) +else() + set(DOWNLOAD_MDI_DEFAULT ON) +endif() +option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an already installed one" ${DOWNLOAD_MDI_DEFAULT}) +if(DOWNLOAD_MDI) + message(STATUS "MDI download requested - we will build our own") + set(mdi_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.2.7.tar.gz" CACHE STRING "URL for MDI tarball") + set(mdi_MD5 "2f3177b30ccdbd6ae28ea3bdd5fed0db" CACHE STRING "MD5 checksum for MDI tarball") + mark_as_advanced(mdi_URL) + mark_as_advanced(mdi_MD5) + + set(LAMMPS_LIB_MDI_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/mdi) + + include(ExternalProject) + message(STATUS "Building mdi.") + ExternalProject_Add(mdi_external + URL ${mdi_URL} + URL_MD5 ${mdi_MD5} + UPDATE_COMMAND "" + CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${LAMMPS_LIB_MDI_BIN_DIR} + -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} + -DCMAKE_INSTALL_LIBDIR=${CMAKE_INSTALL_LIBDIR} + -DCMAKE_INSTALL_INCLUDEDIR=${CMAKE_INSTALL_INCLUDEDIR} + -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} + -DENABLE_OPENMP=${ENABLE_OPENMP} + -DENABLE_XHOST=${ENABLE_XHOST} + -DBUILD_FPIC=${BUILD_FPIC} + -DENABLE_GENERIC=${ENABLE_GENERIC} + -DLIBC_INTERJECT=${LIBC_INTERJECT} + -Dlanguage=C + CMAKE_CACHE_ARGS -DCMAKE_C_FLAGS:STRING=${CMAKE_C_FLAGS} + -DCMAKE_CXX_FLAGS:STRING=${CMAKE_CXX_FLAGS} + -DTargetOpenMP_FIND_COMPONENTS:STRING=C;CXX) + + # Link the lammps library against MDI + target_include_directories(lammps PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/include/mdi) + target_link_directories(lammps PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/lib/mdi) + target_link_libraries(lammps PRIVATE mdi) + + # Link the lammps executable against MDI + target_include_directories(lmp PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/include/mdi) + target_link_directories(lmp PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/lib/mdi) + target_link_libraries(lmp PRIVATE mdi) + + add_definitions(-DLMP_USER_MDI=1) + +else() + + find_package(mdi) + if(NOT mdi_FOUND) + message(FATAL_ERROR "MDI library not found. Help CMake to find it by setting mdi_LIBRARY and mdi_INCLUDE_DIR, or set DOWNLOAD_MDI=ON to download it") + endif() + + # Link the lammps library against MDI + target_include_directories(lammps PRIVATE ${mdi_INCLUDE_DIR}) + target_link_libraries(lammps PRIVATE ${mdi_LIBRARIES}) + + # Link the lammps executable against MDI + target_include_directories(lmp PRIVATE ${mdi_INCLUDE_DIR}) + target_link_libraries(lmp PRIVATE ${mdi_LIBRARIES}) + +endif() diff --git a/lib/mdi/.gitignore b/lib/mdi/.gitignore new file mode 100644 index 0000000000..7fcc712431 --- /dev/null +++ b/lib/mdi/.gitignore @@ -0,0 +1,69 @@ +# Prerequisites +*.d + +# Object files +*.o +*.ko +*.obj +*.elf + +# Linker output +*.ilk +*.map +*.exp + +# Precompiled Headers +*.gch +*.pch + +# Libraries +*.lib +*.a +*.la +*.lo + +# Shared objects (inc. Windows DLLs) +*.dll +*.so +*.so.* +*.dylib + +# Executables +*.exe +*.out +*.app +*.i*86 +*.x86_64 +*.hex + +# Debug files +*.dSYM/ +*.su +*.idb +*.pdb + +# Kernel Module Compile Results +*.mod* +*.cmd +.tmp_versions/ +modules.order +Module.symvers +Mkfile.old +dkms.conf + +# Emacs temporary files +*~ + +# MacOS +.DS_Store + +# pip +*.egg-info + +# MDI files +MDI_Library/ +build/ +libmdi.a +mdi.h +includelink +liblink \ No newline at end of file diff --git a/lib/mdi/Install.py b/lib/mdi/Install.py new file mode 100644 index 0000000000..b2dba5b400 --- /dev/null +++ b/lib/mdi/Install.py @@ -0,0 +1,210 @@ +#!/usr/bin/env python + +# install.py tool to do a generic build of a library +# soft linked to by many of the lib/Install.py files +# used to automate the steps described in the corresponding lib/README + +from __future__ import print_function +import sys,os,subprocess +import glob + +sys.path.append('..') +from install_helpers import checkmd5sum + +# help message + +help = """ +Syntax from src dir: make lib-libname args="-m machine -e suffix" +Syntax from lib dir: python Install.py -m machine -e suffix + +libname = name of lib dir (e.g. atc, h5md, meam, poems, etc) +specify -m and optionally -e, order does not matter + + -m = peform a clean followed by "make -f Makefile.machine" + machine = suffix of a lib/Makefile.* file + -e = set EXTRAMAKE variable in Makefile.machine to Makefile.lammps.suffix + does not alter existing Makefile.machine + +Examples: + +make lib-poems args="-m serial" # build POEMS lib with same settings as in the serial Makefile in src +make lib-colvars args="-m mpi" # build USER-COLVARS lib with same settings as in the mpi Makefile in src +make lib-meam args="-m ifort" # build MEAM lib with custom Makefile.ifort (using Intel Fortran) +""" + +# settings + +version = "1.2.7" +url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version + +# known checksums for different MDI versions. used to validate the download. +checksums = { \ + '1.2.7' : '2f3177b30ccdbd6ae28ea3bdd5fed0db' \ + } + +# print error message or help + +def error(str=None): + if not str: print(help) + else: print("ERROR",str) + sys.exit() + + # expand to full path name +# process leading '~' or relative path + +def fullpath(path): + return os.path.abspath(os.path.expanduser(path)) + +def which(program): + def is_exe(fpath): + return os.path.isfile(fpath) and os.access(fpath, os.X_OK) + + fpath, fname = os.path.split(program) + if fpath: + if is_exe(program): + return program + else: + for path in os.environ["PATH"].split(os.pathsep): + path = path.strip('"') + exe_file = os.path.join(path, program) + if is_exe(exe_file): + return exe_file + + return None + +def geturl(url,fname): + success = False + + if which('curl') != None: + cmd = 'curl -L -o "%s" %s' % (fname,url) + try: + subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) + success = True + except subprocess.CalledProcessError as e: + print("Calling curl failed with: %s" % e.output.decode('UTF-8')) + + if not success and which('wget') != None: + cmd = 'wget -O "%s" %s' % (fname,url) + try: + subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) + success = True + except subprocess.CalledProcessError as e: + print("Calling wget failed with: %s" % e.output.decode('UTF-8')) + + if not success: + error("Failed to download source code with 'curl' or 'wget'") + return + +# parse args + +args = sys.argv[1:] +nargs = len(args) +if nargs == 0: error() + +machine = None +extraflag = 0 + +iarg = 0 +while iarg < nargs: + if args[iarg] == "-m": + if iarg+2 > nargs: error() + machine = args[iarg+1] + iarg += 2 + elif args[iarg] == "-e": + if iarg+2 > nargs: error() + extraflag = 1 + suffix = args[iarg+1] + iarg += 2 + else: error() + +# set lib from working dir + +cwd = os.getcwd() +lib = os.path.basename(cwd) + +# download and unpack MDI_Library tarball + +homepath = "." +homedir = "%s/MDI_Library" % homepath + +print("Downloading MDI_Library ...") +mditar = "%s/v%s.tar.gz" % (homepath,version) +geturl(url, mditar) + +# verify downloaded archive integrity via md5 checksum, if known. +if version in checksums: + if not checkmd5sum(checksums[version], mditar): + sys.exit("Checksum for MDI library does not match") + +print("Unpacking MDI_Library tarball ...") +if os.path.exists("%s/v%s" % (homepath,version)): + cmd = 'rm -rf "%s/v%s"' % (homepath,version) + subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) +cmd = 'cd "%s"; tar -xzvf v%s.tar.gz' % (homepath,version) +subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) +os.remove("%s/v%s.tar.gz" % (homepath,version)) +if os.path.basename(homedir) != version: + if os.path.exists(homedir): + cmd = 'rm -rf "%s"' % homedir + subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) + os.rename("%s/MDI_Library-%s" % (homepath,version),homedir) + +# create Makefile.auto as copy of Makefile.machine +# reset EXTRAMAKE if requested + +if not os.path.exists("Makefile.%s" % machine): + error("lib/%s/Makefile.%s does not exist" % (lib,machine)) + +lines = open("Makefile.%s" % machine,'r').readlines() +fp = open("Makefile.auto",'w') + +has_extramake = False +for line in lines: + words = line.split() + if len(words) == 3 and words[0] == "EXTRAMAKE" and words[1] == '=': + has_extramake = True + if extraflag: + line = line.replace(words[2],"Makefile.lammps.%s" % suffix) + fp.write(line) + +fp.close() + +# make the library via Makefile.auto optionally with parallel make + +try: + import multiprocessing + n_cpus = multiprocessing.cpu_count() +except: + n_cpus = 1 + +print("Building lib%s.so ..." % lib) +cmd = "make -f Makefile.auto clean; make -f Makefile.auto -j%d" % n_cpus +txt = subprocess.check_output(cmd,shell=True,stderr=subprocess.STDOUT) +print(txt.decode('UTF-8')) + +# create 2 links in lib/mdi to MDI Library src dir + +print("Creating links to MDI Library include and lib files") +if os.path.isfile("includelink") or os.path.islink("includelink"): + os.remove("includelink") +if os.path.isfile("liblink") or os.path.islink("liblink"): + os.remove("liblink") +os.symlink(os.path.join(homedir, 'MDI_Library'), 'includelink') +os.symlink(os.path.join(homepath, 'build', 'MDI_Library'), 'liblink') + +# Append the -rpath option to Makefile.lammps + +dir_path = os.path.dirname(os.path.realpath(__file__)) +rpath_option = "-Wl,-rpath=" + str(dir_path) + "/liblink" +makefile_lammps = open(str(dir_path) + "/Makefile.lammps", "a") +makefile_lammps.write(str(rpath_option) + "\n") +makefile_lammps.close() + + +shared_files = glob.glob( os.path.join( homepath, "liblink", "lib%s.so*" % lib) ) +if len(shared_files) > 0: + print("Build was successful") +else: + error("Build of lib/%s/lib%s.so was NOT successful" % (lib,lib)) +if has_extramake and not os.path.exists("Makefile.lammps"): + print("lib/%s/Makefile.lammps was NOT created" % lib) diff --git a/lib/mdi/Makefile.g++ b/lib/mdi/Makefile.g++ new file mode 100644 index 0000000000..15cc9c3a3b --- /dev/null +++ b/lib/mdi/Makefile.g++ @@ -0,0 +1,18 @@ +SHELL = /bin/sh + +# which file will be copied to Makefile.lammps + +EXTRAMAKE = Makefile.lammps.empty + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + mkdir -p build + cd build; cmake -Dlibtype=SHARED -Dlanguage=CXX -D CMAKE_C_COMPILER=gcc -D CMAKE_CXX_COMPILER=g++ ../MDI_Library; make + @cp $(EXTRAMAKE) Makefile.lammps + +# ------ CLEAN ------ + +clean: + -rm *.o *.h $(LIB) + -rm -r build diff --git a/lib/mdi/Makefile.gcc b/lib/mdi/Makefile.gcc new file mode 100644 index 0000000000..ef8a96f5ad --- /dev/null +++ b/lib/mdi/Makefile.gcc @@ -0,0 +1,18 @@ +SHELL = /bin/sh + +# which file will be copied to Makefile.lammps + +EXTRAMAKE = Makefile.lammps.empty + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + mkdir -p build + cd build; cmake -Dlibtype=SHARED -Dlanguage=C -D CMAKE_C_COMPILER=gcc ../MDI_Library; make + @cp $(EXTRAMAKE) Makefile.lammps + +# ------ CLEAN ------ + +clean: + -rm *.o *.h $(LIB) + -rm -r build diff --git a/lib/mdi/Makefile.icc b/lib/mdi/Makefile.icc new file mode 100644 index 0000000000..d70ff28f97 --- /dev/null +++ b/lib/mdi/Makefile.icc @@ -0,0 +1,18 @@ +SHELL = /bin/sh + +# which file will be copied to Makefile.lammps + +EXTRAMAKE = Makefile.lammps.empty + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + mkdir -p build + cd build; cmake -Dlibtype=SHARED -Dlanguage=C -D CMAKE_C_COMPILER=icc -D CMAKE_CXX_COMPILER=icc ../MDI_Library; make + @cp $(EXTRAMAKE) Makefile.lammps + +# ------ CLEAN ------ + +clean: + -rm *.o *.h $(LIB) + -rm -r build diff --git a/lib/mdi/Makefile.lammps.empty b/lib/mdi/Makefile.lammps.empty new file mode 100644 index 0000000000..3aabe162d2 --- /dev/null +++ b/lib/mdi/Makefile.lammps.empty @@ -0,0 +1,5 @@ +# Settings that the LAMMPS build will import when this package library is used + +mdi_SYSINC = +mdi_SYSLIB = +mdi_SYSPATH = \ No newline at end of file diff --git a/lib/mdi/Makefile.mpi b/lib/mdi/Makefile.mpi new file mode 100644 index 0000000000..cff60f5e0b --- /dev/null +++ b/lib/mdi/Makefile.mpi @@ -0,0 +1,18 @@ +SHELL = /bin/sh + +# which file will be copied to Makefile.lammps + +EXTRAMAKE = Makefile.lammps.empty + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + mkdir -p build + cd build; cmake -Dlibtype=SHARED -Dlanguage=CXX -D CMAKE_C_COMPILER=mpicc -D CMAKE_CXX_COMPILER=mpicxx ../MDI_Library; make + @cp $(EXTRAMAKE) Makefile.lammps + +# ------ CLEAN ------ + +clean: + -rm *.o *.h $(LIB) + -rm -r build diff --git a/src/.gitignore b/src/.gitignore index 20a0bde94e..bad9d70360 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -627,6 +627,8 @@ /fix_lb_viscous.h /fix_load_report.cpp /fix_load_report.h +/fix_mdi_engine.cpp +/fix_mdi_engine.h /fix_meso.cpp /fix_meso.h /fix_meso_move.cpp @@ -839,6 +841,8 @@ /lj_sdk_common.h /math_complex.h /math_vector.h +/mdi_engine.cpp +/mdi_engine.h /message.cpp /message.h /mgpt_*.cpp diff --git a/src/Makefile b/src/Makefile index 6523870d75..a0de708070 100644 --- a/src/Makefile +++ b/src/Makefile @@ -54,20 +54,21 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \ PACKUSER = user-adios user-atc user-awpmd user-brownian user-bocs user-cgdna \ user-cgsdk user-colvars user-diffraction user-dpd user-drude \ user-eff user-fep user-h5md user-intel user-lb user-manifold \ - user-meamc user-mesodpd user-mesont user-mgpt user-misc \ + user-mdi user-meamc user-mesodpd user-mesont user-mgpt user-misc \ user-mofff \user-molfile user-netcdf user-omp user-phonon \ user-pace user-plumed user-ptm user-qmmm user-qtb user-quip \ user-reaction user-reaxc user-scafacos user-smd user-smtbq \ user-sdpd user-sph user-tally user-uef user-vtk user-yaff PACKLIB = compress gpu kim kokkos latte message mpiio mscg poems python voronoi \ - user-adios user-atc user-awpmd user-colvars user-h5md user-lb user-mesont \ - user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \ + user-adios user-atc user-awpmd user-colvars user-h5md user-lb user-mdi \ + user-mesont user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \ user-scafacos user-smd user-vtk PACKSYS = compress mpiio python user-lb -PACKINT = gpu kokkos message poems user-atc user-awpmd user-colvars user-mesont +PACKINT = gpu kokkos message poems user-atc user-awpmd user-colvars user-mesont \ + user-mdi PACKEXT = kim latte mscg voronoi \ user-adios user-h5md user-molfile user-netcdf user-pace user-plumed \ diff --git a/src/USER-MDI/Install.sh b/src/USER-MDI/Install.sh new file mode 100755 index 0000000000..7a1f0f6b9c --- /dev/null +++ b/src/USER-MDI/Install.sh @@ -0,0 +1,68 @@ +# Install/unInstall package files in LAMMPS +# mode = 0/1/2 for uninstall/install/update + +mode=$1 + +# enforce using portable C locale +LC_ALL=C +export LC_ALL + +# arg1 = file, arg2 = file it depends on + +action () { + if (test $mode = 0) then + rm -f ../$1 + elif (! cmp -s $1 ../$1) then + if (test -z "$2" || test -e ../$2) then + cp $1 .. + if (test $mode = 2) then + echo " updating src/$1" + fi + fi + elif (test -n "$2") then + if (test ! -e ../$2) then + rm -f ../$1 + fi + fi +} + +# all package files with no dependencies + +for file in *.cpp *.h; do + test -f ${file} && action $file +done + +# edit 2 Makefile.package files to include/exclude package info + +if (test $1 = 1) then + + if (test -e ../Makefile.package) then + sed -i -e 's/[^ \t]*mdi[^ \t]* //g' ../Makefile.package + sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/mdi/includelink -DLMP_USER_MDI |' ../Makefile.package + sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/mdi/liblink |' ../Makefile.package + sed -i -e 's|^PKG_LIB =[ \t]*|&-lmdi |' ../Makefile.package + sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(mdi_SYSINC) |' ../Makefile.package + sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(mdi_SYSLIB) |' ../Makefile.package + sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(mdi_SYSPATH) |' ../Makefile.package + fi + + if (test -e ../Makefile.package.settings) then + sed -i -e '/^include.*mdi.*$/d' ../Makefile.package.settings + # multiline form needed for BSD sed on Macs + sed -i -e '4 i \ +include ..\/..\/lib\/mdi\/Makefile.lammps +' ../Makefile.package.settings + + fi + +elif (test $1 = 0) then + + if (test -e ../Makefile.package) then + sed -i -e 's/[^ \t]*mdi[^ \t]* //g' ../Makefile.package + fi + + if (test -e ../Makefile.package.settings) then + sed -i -e '/^include.*mdi.*$/d' ../Makefile.package.settings + fi + +fi diff --git a/src/USER-MDI/fix_mdi_engine.cpp b/src/USER-MDI/fix_mdi_engine.cpp new file mode 100644 index 0000000000..8c0e2b2072 --- /dev/null +++ b/src/USER-MDI/fix_mdi_engine.cpp @@ -0,0 +1,926 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Taylor Barnes (MolSSI) + MolSSI Driver Interface (MDI) support for LAMMPS +------------------------------------------------------------------------- */ + +#include "fix_mdi_engine.h" + +#include "atom.h" +#include "comm.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "irregular.h" +#include "mdi.h" +#include "memory.h" +#include "modify.h" +#include "update.h" +#include "verlet.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/****************************************************************************/ + + +/*************************************************************** + * create class and parse arguments in LAMMPS script. Syntax: + * fix ID group-ID mdi_engine [couple ] + ***************************************************************/ +FixMDIEngine::FixMDIEngine(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + id_pe(NULL), pe(NULL), + id_ke(NULL), ke(NULL) +{ + + if (narg > 3) + error->all(FLERR,"Illegal fix mdi command"); + + // allocate arrays + memory->create(add_force,3*atom->natoms,"mdi:add_force"); + for (int i=0; i< 3*atom->natoms; i++) { + add_force[i] = 0.0; + } + + // create a new compute pe style + // id = fix-ID + pe, compute group = all + + master = (comm->me==0) ? 1 : 0; + + // create instance of the Irregular class + irregular = new Irregular(lmp); + + most_recent_init = 0; + exit_flag = false; + local_exit_flag = false; + target_command = new char[MDI_COMMAND_LENGTH+1]; + command = new char[MDI_COMMAND_LENGTH+1]; + current_node = new char[MDI_COMMAND_LENGTH]; + target_node = new char[MDI_COMMAND_LENGTH]; + strncpy(target_node, "\0", MDI_COMMAND_LENGTH); + strncpy(current_node, "@DEFAULT", MDI_COMMAND_LENGTH); + + // create and add a compute for PE calculation + int n_pe = strlen(id) + 4; + id_pe = new char[n_pe]; + strcpy(id_pe,id); + strcat(id_pe,"_pe"); + char **newarg_pe = new char*[3]; + newarg_pe[0] = id_pe; + newarg_pe[1] = (char *) "all"; + newarg_pe[2] = (char *) "pe"; + modify->add_compute(3,newarg_pe); + delete [] newarg_pe; + + // create and add a compute for KE calculation + int n_ke = strlen(id) + 4; + id_ke = new char[n_ke]; + strcpy(id_ke,id); + strcat(id_ke,"_ke"); + char **newarg_ke = new char*[3]; + newarg_ke[0] = id_ke; + newarg_ke[1] = (char *) "all"; + newarg_ke[2] = (char *) "ke"; + modify->add_compute(3,newarg_ke); + delete [] newarg_ke; + + // accept a communicator to the driver + int ierr; + if (master) { + MDI_Accept_Communicator(&driver_socket); + if (driver_socket <= 0) + error->all(FLERR,"Unable to connect to driver"); + } else driver_socket=0; + +} + +/********************************* + * Clean up on deleting the fix. * + *********************************/ +FixMDIEngine::~FixMDIEngine() +{ + modify->delete_compute(id_pe); + modify->delete_compute(id_ke); + delete irregular; + delete [] id_pe; + delete [] id_ke; + delete [] target_command; + delete [] command; +} + +/* ---------------------------------------------------------------------- */ +int FixMDIEngine::setmask() +{ + int mask = 0; + + // MD masks + mask |= POST_INTEGRATE; + mask |= PRE_REVERSE; + mask |= POST_FORCE; + + // Minimizer masks + mask |= MIN_PRE_FORCE; + mask |= MIN_POST_FORCE; + + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::exchange_forces() +{ + double **f = atom->f; + const int * const mask = atom->mask; + const int nlocal = atom->nlocal; + + // add the forces from the driver + for (int i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + f[i][0] += add_force[3*(atom->tag[i]-1)+0]; + f[i][1] += add_force[3*(atom->tag[i]-1)+1]; + f[i][2] += add_force[3*(atom->tag[i]-1)+2]; + } + } + +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::init() +{ + // Confirm that the required computes are available + int icompute_pe = modify->find_compute(id_pe); + if (icompute_pe < 0) + error->all(FLERR,"Potential energy ID for fix mdi does not exist"); + int icompute_ke = modify->find_compute(id_ke); + if (icompute_pe < 0) + error->all(FLERR,"Kinetic energy ID for fix mdi does not exist"); + + pe = modify->compute[icompute_pe]; + ke = modify->compute[icompute_ke]; + + return; + +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::setup(int vflag) +{ + //compute the potential energy + potential_energy = pe->compute_scalar(); + kinetic_energy = ke->compute_scalar(); + + // trigger potential energy computation on next timestep + pe->addstep(update->ntimestep+1); + ke->addstep(update->ntimestep+1); + + if ( most_recent_init == 1 ) { // md + // @PRE-FORCES + engine_mode("@PRE-FORCES"); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::min_setup(int vflag) +{ + potential_energy = pe->compute_scalar(); + kinetic_energy = ke->compute_scalar(); + + // trigger potential energy computation on next timestep + pe->addstep(update->ntimestep+1); + ke->addstep(update->ntimestep+1); + + // @FORCES + engine_mode("@FORCES"); +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::post_integrate() +{ + engine_mode("@COORDS"); +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::pre_reverse(int eflag, int vflag) +{ + // calculate the energy + potential_energy = pe->compute_scalar(); + kinetic_energy = ke->compute_scalar(); + + // @PRE-FORCES + engine_mode("@PRE-FORCES"); + + // trigger potential energy computation on next timestep + pe->addstep(update->ntimestep+1); + ke->addstep(update->ntimestep+1); +} + + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::min_pre_force(int vflag) +{ + // @COORDS + engine_mode("@COORDS"); +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::min_post_force(int vflag) +{ + // calculate the energy + potential_energy = pe->compute_scalar(); + kinetic_energy = ke->compute_scalar(); + + // @FORCES + engine_mode("@FORCES"); + + // trigger potential energy computation on next timestep + pe->addstep(update->ntimestep+1); + ke->addstep(update->ntimestep+1); +} + +/* ---------------------------------------------------------------------- */ + +void FixMDIEngine::post_force(int vflag) +{ + if ( most_recent_init == 1 ) { // md + // @FORCES + engine_mode("@FORCES"); + } + else if ( most_recent_init == 2 ) { // optg + // @FORCES + engine_mode("@FORCES"); + } +} + +/* ---------------------------------------------------------------------- */ + +char *FixMDIEngine::engine_mode(const char *node) +{ + /* + if (screen) + fprintf(screen,"MDI ENGINE MODE: %i\n",node); + if (logfile) + fprintf(logfile,"MDI ENGINE MODE: %i\n",node); + */ + + // flag to indicate whether the engine should continue listening for commands at this node + strncpy(current_node, node, MDI_COMMAND_LENGTH); + if ( strcmp(target_node,"\0") != 0 and strcmp(target_node, current_node) != 0 ) { + local_exit_flag = true; + } + + /* ----------------------------------------------------------------- */ + // Answer commands from the driver + /* ----------------------------------------------------------------- */ + + while (not exit_flag and not local_exit_flag) { + + if (master) { + // read the next command from the driver + ierr = MDI_Recv_Command(command, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to receive command from driver"); + command[MDI_COMMAND_LENGTH]=0; + } + // broadcast the command to the other tasks + MPI_Bcast(command,MDI_COMMAND_LENGTH,MPI_CHAR,0,world); + + /* + if (screen) + fprintf(screen,"MDI command: %s\n",command); + if (logfile) + fprintf(logfile,"MDI command: %s:\n",command); + */ + + // confirm that this command is supported at this node + int command_exists = 0; + ierr = MDI_Check_Command_Exists(current_node, command, MDI_COMM_NULL, &command_exists); + if (ierr != 0) + error->all(FLERR,"Unable to check whether the current command is supported"); + if ( command_exists != 1 ) + error->all(FLERR,"Received a command that is not supported at the current node"); + + if (strcmp(command,"STATUS ") == 0 ) { + // send the calculation status to the driver + if (master) { + ierr = MDI_Send_Command("READY", driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to return status to driver"); + } + } + else if (strcmp(command,">NATOMS") == 0 ) { + // receive the number of atoms from the driver + if (master) { + ierr = MDI_Recv((char*) &atom->natoms, 1, MDI_INT, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to receive number of atoms from driver"); + } + MPI_Bcast(&atom->natoms,1,MPI_INT,0,world); + } + else if (strcmp(command,"natoms; + ierr = MDI_Send((char*) &mdi_natoms, 1, MDI_INT64_T, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send number of atoms to driver"); + } + } + else if (strcmp(command,"ntypes, 1, MDI_INT, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send number of atom types to driver"); + } + } + else if (strcmp(command,"COORDS") == 0 ) { + // receive the coordinate information + receive_coordinates(error); + } + else if (strcmp(command,"FORCES") == 0 ) { + // receive the forces from the driver + receive_forces(error); + } + else if (strcmp(command,"+FORCES") == 0 ) { + // receive additional forces from the driver + // these are added prior to SHAKE or other post-processing + add_forces(error); + } + else if (strcmp(command,"@INIT_MD") == 0 ) { + if ( most_recent_init != 0 ) { + error->all(FLERR,"MDI is already performing a simulation"); + } + + // initialize a new MD simulation + most_recent_init = 1; + local_exit_flag = true; + } + else if (strcmp(command,"@INIT_OPTG") == 0 ) { + if ( most_recent_init != 0 ) { + error->all(FLERR,"MDI is already performing a simulation"); + } + + // initialize a new geometry optimization + most_recent_init = 2; + local_exit_flag = true; + //optg_init(error); + } + else if (strcmp(command,"@") == 0 ) { + strncpy(target_node, "\0", MDI_COMMAND_LENGTH); + local_exit_flag = true; + } + else if (strcmp(command,"<@") == 0 ) { + if (master) { + ierr = MDI_Send(current_node, MDI_NAME_LENGTH, MDI_CHAR, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send node to driver"); + } + } + else if (strcmp(command,"etol = std::numeric_limits::max(); + update->ftol = std::numeric_limits::max(); + + // set the maximum number of force evaluations to 0 + update->max_eval = 0; + } + } + else if (strcmp(command,"EXIT") == 0 ) { + // exit the driver code + exit_flag = true; + + // are we in the middle of a geometry optimization? + if ( most_recent_init == 2 ) { + // ensure that the energy and force tolerances are met + update->etol = std::numeric_limits::max(); + update->ftol = std::numeric_limits::max(); + + // set the maximum number of force evaluations to 0 + update->max_eval = 0; + } + } + else { + // the command is not supported + error->all(FLERR,"Unknown command from driver"); + } + + // check if the target node is something other than the current node + if ( strcmp(target_node,"\0") != 0 and strcmp(target_node, current_node) != 0 ) { + local_exit_flag = true; + } + + } + + // a local exit has completed, so turn off the local exit flag + local_exit_flag = false; + + return command; + +} + + +void FixMDIEngine::receive_coordinates(Error* error) +{ + double posconv; + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + posconv=force->angstrom/angstrom_to_bohr; + + // create a buffer to hold the coordinates + double *buffer; + buffer = new double[3*atom->natoms]; + + if (master) { + ierr = MDI_Recv((char*) buffer, 3*atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to receive coordinates from driver"); + } + MPI_Bcast(buffer,3*atom->natoms,MPI_DOUBLE,0,world); + + // pick local atoms from the buffer + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + x[i][0]=buffer[3*(atom->tag[i]-1)+0]*posconv; + x[i][1]=buffer[3*(atom->tag[i]-1)+1]*posconv; + x[i][2]=buffer[3*(atom->tag[i]-1)+2]*posconv; + } + + // ensure atoms are in current box & update box via shrink-wrap + // has to be be done before invoking Irregular::migrate_atoms() + // since it requires atoms be inside simulation box + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + if (domain->triclinic) domain->lamda2x(atom->nlocal); + + // move atoms to new processors via irregular() + // only needed if migrate_check() says an atom moves to far + if (domain->triclinic) domain->x2lamda(atom->nlocal); + if (irregular->migrate_check()) irregular->migrate_atoms(); + if (domain->triclinic) domain->lamda2x(atom->nlocal); + + delete [] buffer; +} + + +void FixMDIEngine::send_coordinates(Error* error) +{ + double posconv; + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + posconv=force->angstrom/angstrom_to_bohr; + + double *coords; + double *coords_reduced; + + coords = new double[3*atom->natoms]; + coords_reduced = new double[3*atom->natoms]; + + // zero the coordinates array + for (int icoord = 0; icoord < 3*atom->natoms; icoord++) { + coords[icoord] = 0.0; + } + + // pick local atoms from the buffer + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + coords[3*(atom->tag[i]-1)+0] = x[i][0]/posconv; + coords[3*(atom->tag[i]-1)+1] = x[i][1]/posconv; + coords[3*(atom->tag[i]-1)+2] = x[i][2]/posconv; + } + + MPI_Reduce(coords, coords_reduced, 3*atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world); + + if (master) { + ierr = MDI_Send((char*) coords_reduced, 3*atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send coordinates to driver"); + } + + delete [] coords; + delete [] coords_reduced; +} + + +void FixMDIEngine::send_charges(Error* error) +{ + double *charges; + double *charges_reduced; + + charges = new double[atom->natoms]; + charges_reduced = new double[atom->natoms]; + + // zero the charges array + for (int icharge = 0; icharge < atom->natoms; icharge++) { + charges[icharge] = 0.0; + } + + // pick local atoms from the buffer + double *charge = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + charges[atom->tag[i]-1] = charge[i]; + } + + MPI_Reduce(charges, charges_reduced, atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world); + + if (master) { + ierr = MDI_Send((char*) charges_reduced, atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send charges to driver"); + } + + delete [] charges; + delete [] charges_reduced; +} + + +void FixMDIEngine::send_energy(Error* error) +{ + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double pe = potential_energy; + double ke = kinetic_energy; + double total_energy; + double *send_energy = &total_energy; + + // convert the energy to atomic units + pe *= kelvin_to_hartree/force->boltz; + ke *= kelvin_to_hartree/force->boltz; + total_energy = pe + ke; + + if (master) { + ierr = MDI_Send((char*) send_energy, 1, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send potential energy to driver"); + } +} + + +void FixMDIEngine::send_pe(Error* error) +{ + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double pe = potential_energy; + double *send_energy = &pe; + + // convert the energy to atomic units + pe *= kelvin_to_hartree/force->boltz; + + if (master) { + ierr = MDI_Send((char*) send_energy, 1, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send potential energy to driver"); + } +} + + +void FixMDIEngine::send_ke(Error* error) +{ + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double ke = kinetic_energy; + double *send_energy = &ke; + + // convert the energy to atomic units + ke *= kelvin_to_hartree/force->boltz; + + if (master) { + ierr = MDI_Send((char*) send_energy, 1, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send potential energy to driver"); + } +} + + +void FixMDIEngine::send_types(Error* error) +{ + int * const type = atom->type; + + if (master) { + ierr = MDI_Send((char*) type, atom->natoms, MDI_INT, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send atom types to driver"); + } +} + + +void FixMDIEngine::send_labels(Error* error) +{ + char *labels = new char[atom->natoms * MDI_LABEL_LENGTH]; + memset(labels, ' ', atom->natoms * MDI_LABEL_LENGTH); + + for (int iatom=0; iatom < atom->natoms; iatom++) { + std::string label = std::to_string( atom->type[iatom] ); + int label_len = std::min( int(label.length()), MDI_LABEL_LENGTH ); + strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len); + } + + if (master) { + ierr = MDI_Send( labels, atom->natoms * MDI_LABEL_LENGTH, MDI_CHAR, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send atom types to driver"); + } + + delete [] labels; +} + + +void FixMDIEngine::send_masses(Error* error) +{ + double * const mass = atom->mass; + int * const type = atom->type; + + if (master) { + double *mass_by_atom = new double[atom->natoms]; + for (int iatom=0; iatom < atom->natoms; iatom++) { + mass_by_atom[ atom->tag[iatom] - 1 ] = mass[ type[iatom] ]; + } + ierr = MDI_Send((char*) mass_by_atom, atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send atom masses to driver"); + delete [] mass_by_atom; + } + +} + + +void FixMDIEngine::send_forces(Error* error) +{ + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double potconv, posconv, forceconv; + potconv=kelvin_to_hartree/force->boltz; + posconv=force->angstrom/angstrom_to_bohr; + forceconv=potconv*posconv; + + double *forces; + double *forces_reduced; + double *x_buf; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + forces = new double[3*atom->natoms]; + forces_reduced = new double[3*atom->natoms]; + x_buf = new double[3*nlocal]; + + // zero the forces array + for (int iforce = 0; iforce < 3*atom->natoms; iforce++) { + forces[iforce] = 0.0; + } + + // if not at a node, calculate the forces + if ( strcmp(current_node, "@DEFAULT") == 0 ) { + // certain fixes, such as shake, move the coordinates + // to ensure that the coordinates do not change, store a copy + double **x = atom->x; + for (int i = 0; i < nlocal; i++) { + x_buf[3*i+0] = x[i][0]; + x_buf[3*i+1] = x[i][1]; + x_buf[3*i+2] = x[i][2]; + } + + // calculate the forces + update->whichflag = 1; // 1 for dynamics + update->nsteps = 1; + lmp->init(); + update->integrate->setup_minimal(1); + } + + // pick local atoms from the buffer + double **f = atom->f; + for (int i = 0; i < nlocal; i++) { + forces[3*(atom->tag[i]-1)+0] = f[i][0]*forceconv; + forces[3*(atom->tag[i]-1)+1] = f[i][1]*forceconv; + forces[3*(atom->tag[i]-1)+2] = f[i][2]*forceconv; + } + + // reduce the forces onto rank 0 + MPI_Reduce(forces, forces_reduced, 3*atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world); + + // send the forces through MDI + if (master) { + ierr = MDI_Send((char*) forces_reduced, 3*atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send atom forces to driver"); + } + + if ( strcmp(current_node, "@DEFAULT") == 0 ) { + // restore the original set of coordinates + double **x_new = atom->x; + for (int i = 0; i < nlocal; i++) { + x_new[i][0] = x_buf[3*i+0]; + x_new[i][1] = x_buf[3*i+1]; + x_new[i][2] = x_buf[3*i+2]; + } + } + + delete [] forces; + delete [] forces_reduced; + delete [] x_buf; + +} + + +void FixMDIEngine::receive_forces(Error* error) +{ + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double potconv, posconv, forceconv; + potconv=kelvin_to_hartree/force->boltz; + posconv=force->angstrom/angstrom_to_bohr; + forceconv=potconv*posconv; + + double *forces; + forces = new double[3*atom->natoms]; + + if (master) { + ierr = MDI_Recv((char*) forces, 3*atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to receive atom forces to driver"); + } + MPI_Bcast(forces,3*atom->natoms,MPI_DOUBLE,0,world); + + // pick local atoms from the buffer + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + f[i][0] = forces[3*(atom->tag[i]-1)+0]/forceconv; + f[i][1] = forces[3*(atom->tag[i]-1)+1]/forceconv; + f[i][2] = forces[3*(atom->tag[i]-1)+2]/forceconv; + } + + delete [] forces; +} + + +void FixMDIEngine::add_forces(Error* error) +{ + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + double kelvin_to_hartree; + MDI_Conversion_Factor("kelvin_energy", "hartree", &kelvin_to_hartree); + + double potconv, posconv, forceconv; + potconv=kelvin_to_hartree/force->boltz; + posconv=force->angstrom * angstrom_to_bohr; + forceconv=potconv*posconv; + + double *forces; + forces = new double[3*atom->natoms]; + + if (master) { + ierr = MDI_Recv((char*) forces, 3*atom->natoms, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to receive atom forces to driver"); + } + MPI_Bcast(forces,3*atom->natoms,MPI_DOUBLE,0,world); + + // pick local atoms from the buffer + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + f[i][0] += forces[3*(atom->tag[i]-1)+0]/forceconv; + f[i][1] += forces[3*(atom->tag[i]-1)+1]/forceconv; + f[i][2] += forces[3*(atom->tag[i]-1)+2]/forceconv; + } + + delete [] forces; +} + + +void FixMDIEngine::send_cell(Error* error) +{ + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + + double celldata[9]; + + celldata[0] = domain->boxhi[0] - domain->boxlo[0]; + celldata[1] = 0.0; + celldata[2] = 0.0; + celldata[3] = domain->xy; + celldata[4] = domain->boxhi[1] - domain->boxlo[1]; + celldata[5] = 0.0; + celldata[6] = domain->xz; + celldata[7] = domain->yz; + celldata[8] = domain->boxhi[2] - domain->boxlo[2]; + /* + celldata[9 ] = domain->boxlo[0]; + celldata[10] = domain->boxlo[1]; + celldata[11] = domain->boxlo[2]; + */ + + // convert the units to bohr + double unit_conv = force->angstrom * angstrom_to_bohr; + for (int icell=0; icell < 9; icell++) { + celldata[icell] *= unit_conv; + } + + if (master) { + ierr = MDI_Send((char*) celldata, 9, MDI_DOUBLE, driver_socket); + if (ierr != 0) + error->all(FLERR,"Unable to send cell dimensions to driver"); + } +} diff --git a/src/USER-MDI/fix_mdi_engine.h b/src/USER-MDI/fix_mdi_engine.h new file mode 100644 index 0000000000..fcb8f6c62a --- /dev/null +++ b/src/USER-MDI/fix_mdi_engine.h @@ -0,0 +1,137 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + + +#ifdef FIX_CLASS + +FixStyle(mdi/engine,FixMDIEngine) + +#else + +#ifndef LMP_FIX_MDI_ENGINE_H +#define LMP_FIX_MDI_ENGINE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixMDIEngine : public Fix { + public: + FixMDIEngine(class LAMMPS *, int, char **); + ~FixMDIEngine(); + int setmask(); + void init(); + + char *engine_mode(const char *node); + + // receive and update forces + void setup(int); + void min_setup(int); + void post_integrate(); + void post_force(int); + void pre_reverse(int, int); + void min_pre_force(int); //@COORDS + void min_post_force(int); //@FORCES + + double *add_force; // stores forces added using +FORCE command + double potential_energy; // stores potential energy + double kinetic_energy; // stores kinetic energy + + // current command + char *command; + + protected: + void exchange_forces(); + + private: + int master, ierr; + int driver_socket; + int most_recent_init; // which MDI init command was most recently received? + // 0 - none + // 1 - MD + // 2 - OPTG + bool exit_flag; + bool local_exit_flag; + char *current_node; + char *target_node; // is the code supposed to advance to a particular node? + // 0 - none + // 1 - @COORDS (before pre-force calculation) + // 2 - @PRE-FORCES (before final force calculation) + // 3 - @FORCES (before time integration) + // -1 - after MD_INIT command + // -2 - after MD_INIT command followed by @PRE-FORCES (actually @INIT_OPTG?) + + // command to be executed at the target node + char *target_command; + + char *id_pe; + char *id_ke; + class Irregular *irregular; + class Minimize *minimizer; + class Compute *pe; + class Compute *ke; + void send_types(Error *); + void send_labels(Error *); + void send_masses(Error *); + void receive_coordinates(Error *); + void send_coordinates(Error *); + void send_charges(Error *); + void send_energy(Error *); + void send_forces(Error *); + void send_pe(Error *); + void send_ke(Error *); + void add_forces(Error *); + void receive_forces(Error *); + void send_cell(Error *); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. + +E: Potential energy ID for fix mdi does not exist + +Self-explanatory. + +E: Cannot use MDI command without atom IDs + +Self-explanatory. + +E: MDI command requires consecutive atom IDs + +Self-explanatory. + +E: Unable to connect to driver + +Self-explanatory. + +E: Unable to ... driver + +Self-explanatory. + +E: Unknown command from driver + +The driver sent a command that is not supported by the LAMMPS +interface. In some cases this might be because a nonsensical +command was sent (i.e. "SCF"). In other cases, the LAMMPS +interface might benefit from being expanded. + +*/ diff --git a/src/USER-MDI/mdi_engine.cpp b/src/USER-MDI/mdi_engine.cpp new file mode 100644 index 0000000000..f3403329af --- /dev/null +++ b/src/USER-MDI/mdi_engine.cpp @@ -0,0 +1,356 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Taylor Barnes (MolSSI) + MolSSI Driver Interface (MDI) support for LAMMPS +------------------------------------------------------------------------- */ + +#include "mdi_engine.h" + +#include "atom.h" +#include "error.h" +#include "fix_mdi_engine.h" +#include "force.h" +#include "mdi.h" +#include "min.h" +#include "minimize.h" +#include "modify.h" +#include "output.h" +#include "timer.h" +#include "update.h" +#include "verlet.h" + +#include +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +CommandMDIEngine::CommandMDIEngine(LAMMPS *lmp) : Command(lmp) { + return; +} + +CommandMDIEngine::~CommandMDIEngine() { + return; +} + +/* ---------------------------------------------------------------------- */ + +void CommandMDIEngine::command(int narg, char **arg) +{ + + // register the default node + MDI_Register_Node("@DEFAULT"); + MDI_Register_Command("@DEFAULT", "<@"); + MDI_Register_Command("@DEFAULT", "COORDS"); + MDI_Register_Command("@DEFAULT", "@INIT_MD"); + MDI_Register_Command("@DEFAULT", "@INIT_OPTG"); + MDI_Register_Command("@DEFAULT", "EXIT"); + + // register the MD initialization node + MDI_Register_Node("@INIT_MD"); + MDI_Register_Command("@INIT_MD", "<@"); + MDI_Register_Command("@INIT_MD", "COORDS"); + MDI_Register_Command("@INIT_MD", ">FORCES"); + MDI_Register_Command("@INIT_MD", "@"); + MDI_Register_Command("@INIT_MD", "@COORDS"); + MDI_Register_Command("@INIT_MD", "@DEFAULT"); + MDI_Register_Command("@INIT_MD", "@FORCES"); + MDI_Register_Command("@INIT_MD", "@PRE-FORCES"); + MDI_Register_Command("@INIT_MD", "EXIT"); + + // register the OPTG initialization node + MDI_Register_Node("@INIT_OPTG"); + MDI_Register_Command("@INIT_OPTG", "<@"); + MDI_Register_Command("@INIT_OPTG", "COORDS"); + MDI_Register_Command("@INIT_OPTG", ">FORCES"); + MDI_Register_Command("@INIT_OPTG", "@"); + MDI_Register_Command("@INIT_OPTG", "@COORDS"); + MDI_Register_Command("@INIT_OPTG", "@DEFAULT"); + MDI_Register_Command("@INIT_OPTG", "@FORCES"); + MDI_Register_Command("@INIT_OPTG", "EXIT"); + + // register the pre-forces node + MDI_Register_Node("@PRE-FORCES"); + MDI_Register_Command("@PRE-FORCES", "<@"); + MDI_Register_Command("@PRE-FORCES", "COORDS"); + MDI_Register_Command("@PRE-FORCES", ">FORCES"); + MDI_Register_Command("@PRE-FORCES", "@"); + MDI_Register_Command("@PRE-FORCES", "@COORDS"); + MDI_Register_Command("@PRE-FORCES", "@DEFAULT"); + MDI_Register_Command("@PRE-FORCES", "@FORCES"); + MDI_Register_Command("@PRE-FORCES", "@PRE-FORCES"); + MDI_Register_Command("@PRE-FORCES", "EXIT"); + + // register the forces node + MDI_Register_Node("@FORCES"); + MDI_Register_Callback("@FORCES", ">FORCES"); + MDI_Register_Command("@FORCES", "<@"); + MDI_Register_Command("@FORCES", "COORDS"); + MDI_Register_Command("@FORCES", ">FORCES"); + MDI_Register_Command("@FORCES", "@"); + MDI_Register_Command("@FORCES", "@COORDS"); + MDI_Register_Command("@FORCES", "@DEFAULT"); + MDI_Register_Command("@FORCES", "@FORCES"); + MDI_Register_Command("@FORCES", "@PRE-FORCES"); + MDI_Register_Command("@FORCES", "EXIT"); + + // register the coordinates node + MDI_Register_Node("@COORDS"); + MDI_Register_Command("@COORDS", "<@"); + MDI_Register_Command("@COORDS", "COORDS"); + MDI_Register_Command("@COORDS", ">FORCES"); + MDI_Register_Command("@COORDS", "@"); + MDI_Register_Command("@COORDS", "@COORDS"); + MDI_Register_Command("@COORDS", "@DEFAULT"); + MDI_Register_Command("@COORDS", "@FORCES"); + MDI_Register_Command("@COORDS", "@PRE-FORCES"); + MDI_Register_Command("@COORDS", "EXIT"); + + + // identify the mdi_engine fix + int ifix = modify->find_fix_by_style("mdi/engine"); + if (ifix < 0) error->all(FLERR,"The mdi_engine command requires the mdi/engine fix"); + mdi_fix = static_cast(modify->fix[ifix]); + + /* format for MDI Engine command: + * mdi_engine + */ + if (narg > 0) error->all(FLERR,"Illegal MDI command"); + + if (atom->tag_enable == 0) + error->all(FLERR,"Cannot use MDI command without atom IDs"); + + if (atom->tag_consecutive() == 0) + error->all(FLERR,"MDI command requires consecutive atom IDs"); + + // begin engine_mode + char *command = NULL; + while ( true ) { + // listen for MDI commands at the default command + // the response to most MDI commands is handled here + command = mdi_fix->engine_mode("@DEFAULT"); + + // MDI commands that involve large-scale program flow are handled here + if (strcmp(command,"@INIT_MD") == 0 ) { + // enter MD control loop + int received_exit = mdi_md(); + if ( received_exit == 1 ) { + return; + } + } + if (strcmp(command,"@INIT_OPTG") == 0 ) { + // enter minimizer control loop + int received_exit = mdi_optg(); + if ( received_exit == 1 ) { + return; + } + } + else if (strcmp(command,"EXIT") == 0 ) { + return; + } + else { + error->all(FLERR,fmt::format("MDI node exited with invalid command: {}",command)); + } + } + + return; +} + + + +int CommandMDIEngine::mdi_md() +{ + // initialize an MD simulation + update->whichflag = 1; // 1 for dynamics + timer->init_timeout(); + update->nsteps = 1; + update->ntimestep = 0; + update->firststep = update->ntimestep; + update->laststep = update->ntimestep + update->nsteps; + update->beginstep = update->firststep; + update->endstep = update->laststep; + lmp->init(); + + // the MD simulation is now at the @INIT_MD node + char *command = NULL; + command = mdi_fix->engine_mode("@INIT_MD"); + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + // continue the MD simulation + update->integrate->setup(1); + + // the MD simulation is now at the @FORCES node + command = mdi_fix->engine_mode("@FORCES"); + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + // do MD iterations until told to exit + while ( true ) { + + // run an MD timestep + update->whichflag = 1; // 1 for dynamics + timer->init_timeout(); + update->nsteps += 1; + update->laststep += 1; + update->endstep = update->laststep; + output->next = update->ntimestep + 1; + update->integrate->run(1); + + // get the most recent command the MDI engine received + command = mdi_fix->command; + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + } + +} + + + +int CommandMDIEngine::mdi_optg() +{ + char *command = NULL; + + // create instance of the Minimizer class + Minimize *minimizer = new Minimize(lmp); + + // initialize the minimizer in a way that ensures optimization will continue until MDI exits + update->etol = std::numeric_limits::min(); + update->ftol = std::numeric_limits::min(); + update->nsteps = std::numeric_limits::max(); + update->max_eval = std::numeric_limits::max(); + + update->whichflag = 2; // 2 for minimization + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + update->nsteps; + lmp->init(); + + command = mdi_fix->engine_mode("@INIT_OPTG"); + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + update->minimize->setup(); + command = mdi_fix->command; + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + update->minimize->iterate(update->nsteps); + command = mdi_fix->command; + + if (strcmp(command,"@DEFAULT") == 0 ) { + // return, and flag for @DEFAULT node + return 0; + } + else if (strcmp(command,"EXIT") == 0 ) { + // return, and flag for global exit + return 1; + } + + error->all(FLERR,fmt::format("MDI reached end of OPTG simulation with invalid command: {}",command)); + return 0; +} diff --git a/src/USER-MDI/mdi_engine.h b/src/USER-MDI/mdi_engine.h new file mode 100644 index 0000000000..8e74e05e6a --- /dev/null +++ b/src/USER-MDI/mdi_engine.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMMAND_CLASS + +CommandStyle(mdi_engine,CommandMDIEngine) + +#else + +#ifndef LMP_COMMAND_MDI_ENGINE_H +#define LMP_COMMAND_MDI_ENGINE_H + +#include "command.h" + +namespace LAMMPS_NS { + + class CommandMDIEngine : public Command { + public: + CommandMDIEngine(class LAMMPS *); + virtual ~CommandMDIEngine(); + void command(int, char **); + int mdi_md(); + int mdi_optg(); + +private: + class FixMDIEngine *mdi_fix; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/src/main.cpp b/src/main.cpp index 26f12683b9..41f9cfc825 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,6 +13,7 @@ #include "lammps.h" #include "input.h" +#include "mdi_interface.h" #include #include @@ -35,6 +36,21 @@ int main(int argc, char **argv) { MPI_Init(&argc,&argv); + // initialize MDI + if ( MDI_Init(&argc,&argv) ) + MPI_Abort(MPI_COMM_WORLD, 1); + + // get the MPI communicator that spans all ranks running LAMMPS + // when using MDI, this may be a subset of MPI_COMM_WORLD + MPI_Comm lammps_comm = MPI_COMM_WORLD; + int mdi_flag; + if ( MDI_Initialized(&mdi_flag) ) + MPI_Abort(MPI_COMM_WORLD, 1); + if ( mdi_flag ) { + if ( MDI_MPI_get_world_comm( &lammps_comm ) ) + MPI_Abort(MPI_COMM_WORLD, 1); + } + // enable trapping selected floating point exceptions. // this uses GNU extensions and is only tested on Linux // therefore we make it depend on -D_GNU_SOURCE, too. @@ -49,13 +65,13 @@ int main(int argc, char **argv) #ifdef LAMMPS_EXCEPTIONS try { - LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); + LAMMPS *lammps = new LAMMPS(argc,argv,lammps_comm); lammps->input->file(); delete lammps; } catch(LAMMPSAbortException &ae) { MPI_Abort(ae.universe, 1); } catch(LAMMPSException &e) { - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(lammps_comm); MPI_Finalize(); exit(1); } catch(fmt::format_error &fe) { @@ -65,7 +81,7 @@ int main(int argc, char **argv) } #else try { - LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); + LAMMPS *lammps = new LAMMPS(argc,argv,lammps_comm); lammps->input->file(); delete lammps; } catch(fmt::format_error &fe) { @@ -74,6 +90,6 @@ int main(int argc, char **argv) exit(1); } #endif - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(lammps_comm); MPI_Finalize(); } diff --git a/src/mdi_interface.h b/src/mdi_interface.h new file mode 100644 index 0000000000..d76758253f --- /dev/null +++ b/src/mdi_interface.h @@ -0,0 +1,34 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_MDI_INTERFACE_H +#define LMP_MDI_INTERFACE_H + +#if defined(LMP_USER_MDI) + +// true interface to MDI +// used when MDI is installed + +#include "mdi.h" + +#else + +// dummy interface to MDI +// needed for compiling when MDI is not installed + + int MDI_Init(int* argc, char ***argv) {return 0;}; + int MDI_Initialized(int* flag) {return 0;}; + int MDI_MPI_get_world_comm(void* world_comm) {return 0;}; + +#endif +#endif