Added initial MolSSI Driver interface

This commit is contained in:
Taylor Barnes
2018-05-14 14:03:28 -04:00
committed by taylor-a-barnes
parent 058d56cf5c
commit b21479e20a
18 changed files with 2027 additions and 11 deletions

View File

@ -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 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-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-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-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-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) 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) set(CUDA_REQUEST_PIC)
endif() endif()
foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MDI USER-MOLFILE USER-NETCDF USER-PLUMED
USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE) USER-QMMM USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE)
if(PKG_${PKG_WITH_INCL}) if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL}) include(Packages/${PKG_WITH_INCL})
endif() endif()

View File

@ -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()

69
lib/mdi/.gitignore vendored Normal file
View File

@ -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

210
lib/mdi/Install.py Normal file
View File

@ -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)

18
lib/mdi/Makefile.g++ Normal file
View File

@ -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

18
lib/mdi/Makefile.gcc Normal file
View File

@ -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

18
lib/mdi/Makefile.icc Normal file
View File

@ -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

View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
mdi_SYSINC =
mdi_SYSLIB =
mdi_SYSPATH =

18
lib/mdi/Makefile.mpi Normal file
View File

@ -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

4
src/.gitignore vendored
View File

@ -627,6 +627,8 @@
/fix_lb_viscous.h /fix_lb_viscous.h
/fix_load_report.cpp /fix_load_report.cpp
/fix_load_report.h /fix_load_report.h
/fix_mdi_engine.cpp
/fix_mdi_engine.h
/fix_meso.cpp /fix_meso.cpp
/fix_meso.h /fix_meso.h
/fix_meso_move.cpp /fix_meso_move.cpp
@ -839,6 +841,8 @@
/lj_sdk_common.h /lj_sdk_common.h
/math_complex.h /math_complex.h
/math_vector.h /math_vector.h
/mdi_engine.cpp
/mdi_engine.h
/message.cpp /message.cpp
/message.h /message.h
/mgpt_*.cpp /mgpt_*.cpp

View File

@ -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 \ PACKUSER = user-adios user-atc user-awpmd user-brownian user-bocs user-cgdna \
user-cgsdk user-colvars user-diffraction user-dpd user-drude \ user-cgsdk user-colvars user-diffraction user-dpd user-drude \
user-eff user-fep user-h5md user-intel user-lb user-manifold \ 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-mofff \user-molfile user-netcdf user-omp user-phonon \
user-pace user-plumed user-ptm user-qmmm user-qtb user-quip \ user-pace user-plumed user-ptm user-qmmm user-qtb user-quip \
user-reaction user-reaxc user-scafacos user-smd user-smtbq \ user-reaction user-reaxc user-scafacos user-smd user-smtbq \
user-sdpd user-sph user-tally user-uef user-vtk user-yaff user-sdpd user-sph user-tally user-uef user-vtk user-yaff
PACKLIB = compress gpu kim kokkos latte message mpiio mscg poems python voronoi \ 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-adios user-atc user-awpmd user-colvars user-h5md user-lb user-mdi \
user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \ user-mesont user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \
user-scafacos user-smd user-vtk user-scafacos user-smd user-vtk
PACKSYS = compress mpiio python user-lb 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 \ PACKEXT = kim latte mscg voronoi \
user-adios user-h5md user-molfile user-netcdf user-pace user-plumed \ user-adios user-h5md user-molfile user-netcdf user-pace user-plumed \

68
src/USER-MDI/Install.sh Executable file
View File

@ -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

View File

@ -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 <limits>
#include <string.h>
using namespace LAMMPS_NS;
using namespace FixConst;
/****************************************************************************/
/***************************************************************
* create class and parse arguments in LAMMPS script. Syntax:
* fix ID group-ID mdi_engine [couple <group-ID>]
***************************************************************/
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") == 0 ) {
// send the number of atoms to the driver
if (master) {
int64_t mdi_natoms = atom->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") == 0 ) {
// send the number of atom types to the driver
if (master) {
ierr = MDI_Send((char*) &atom->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,"<TYPES") == 0 ) {
// send the atom types
send_types(error);
}
else if (strcmp(command,"<LABELS") == 0 ) {
// send the atom labels
send_labels(error);
}
else if (strcmp(command,"<MASSES") == 0 ) {
// send the atom types
send_masses(error);
}
else if (strcmp(command,"<CELL") == 0 ) {
// send the cell dimensions to the driver
send_cell(error);
}
else if (strcmp(command,">COORDS") == 0 ) {
// receive the coordinate information
receive_coordinates(error);
}
else if (strcmp(command,"<COORDS") == 0 ) {
// send the coordinate information
send_coordinates(error);
}
else if (strcmp(command,"<CHARGES") == 0 ) {
// send the charges
send_charges(error);
}
else if (strcmp(command,"<ENERGY") == 0 ) {
// send the total energy to the driver
send_energy(error);
}
else if (strcmp(command,"<FORCES") == 0 ) {
// send the forces to the driver
send_forces(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,"<KE") == 0 ) {
// send the kinetic energy to the driver
send_ke(error);
}
else if (strcmp(command,"<PE") == 0 ) {
// send the potential energy to the driver
send_pe(error);
}
else if (strcmp(command,"@COORDS") == 0 ) {
strncpy(target_node, "@COORDS", MDI_COMMAND_LENGTH);
local_exit_flag = true;
}
else if (strcmp(command,"@PRE-FORCES") == 0 ) {
strncpy(target_node, "@PRE-FORCES", MDI_COMMAND_LENGTH);
local_exit_flag = true;
}
else if (strcmp(command,"@FORCES") == 0 ) {
strncpy(target_node, "@FORCES", MDI_COMMAND_LENGTH);
local_exit_flag = true;
}
else if (strcmp(command,"EXIT_SIM") == 0 ) {
most_recent_init = 0;
local_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<double>::max();
update->ftol = std::numeric_limits<double>::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<double>::max();
update->ftol = std::numeric_limits<double>::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");
}
}

View File

@ -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.
*/

356
src/USER-MDI/mdi_engine.cpp Normal file
View File

@ -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 <string.h>
#include <limits>
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", "<CELL");
MDI_Register_Command("@DEFAULT", "<CHARGES");
MDI_Register_Command("@DEFAULT", "<COORDS");
MDI_Register_Command("@DEFAULT", "<LABELS");
MDI_Register_Command("@DEFAULT", "<NATOMS");
MDI_Register_Command("@DEFAULT", "<MASSES");
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", "<CELL");
MDI_Register_Command("@INIT_MD", "<CHARGES");
MDI_Register_Command("@INIT_MD", "<COORDS");
MDI_Register_Command("@INIT_MD", "<ENERGY");
MDI_Register_Command("@INIT_MD", "<FORCES");
MDI_Register_Command("@INIT_MD", "<KE");
MDI_Register_Command("@INIT_MD", "<LABELS");
MDI_Register_Command("@INIT_MD", "<MASSES");
MDI_Register_Command("@INIT_MD", "<NATOMS");
MDI_Register_Command("@INIT_MD", "<PE");
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", "<CELL");
MDI_Register_Command("@INIT_OPTG", "<CHARGES");
MDI_Register_Command("@INIT_OPTG", "<COORDS");
MDI_Register_Command("@INIT_OPTG", "<ENERGY");
MDI_Register_Command("@INIT_OPTG", "<FORCES");
MDI_Register_Command("@INIT_OPTG", "<KE");
MDI_Register_Command("@INIT_OPTG", "<LABELS");
MDI_Register_Command("@INIT_OPTG", "<MASSES");
MDI_Register_Command("@INIT_OPTG", "<NATOMS");
MDI_Register_Command("@INIT_OPTG", "<PE");
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", "<CELL");
MDI_Register_Command("@PRE-FORCES", "<CHARGES");
MDI_Register_Command("@PRE-FORCES", "<COORDS");
MDI_Register_Command("@PRE-FORCES", "<ENERGY");
MDI_Register_Command("@PRE-FORCES", "<FORCES");
MDI_Register_Command("@PRE-FORCES", "<KE");
MDI_Register_Command("@PRE-FORCES", "<LABELS");
MDI_Register_Command("@PRE-FORCES", "<MASSES");
MDI_Register_Command("@PRE-FORCES", "<NATOMS");
MDI_Register_Command("@PRE-FORCES", "<PE");
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", "<CELL");
MDI_Register_Command("@FORCES", "<CHARGES");
MDI_Register_Command("@FORCES", "<COORDS");
MDI_Register_Command("@FORCES", "<ENERGY");
MDI_Register_Command("@FORCES", "<FORCES");
MDI_Register_Command("@FORCES", "<KE");
MDI_Register_Command("@FORCES", "<LABELS");
MDI_Register_Command("@FORCES", "<MASSES");
MDI_Register_Command("@FORCES", "<NATOMS");
MDI_Register_Command("@FORCES", "<PE");
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", "<CELL");
MDI_Register_Command("@COORDS", "<CHARGES");
MDI_Register_Command("@COORDS", "<COORDS");
MDI_Register_Command("@COORDS", "<ENERGY");
MDI_Register_Command("@COORDS", "<FORCES");
MDI_Register_Command("@COORDS", "<KE");
MDI_Register_Command("@COORDS", "<LABELS");
MDI_Register_Command("@COORDS", "<MASSES");
MDI_Register_Command("@COORDS", "<NATOMS");
MDI_Register_Command("@COORDS", "<PE");
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<FixMDIEngine*>(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<double>::min();
update->ftol = std::numeric_limits<double>::min();
update->nsteps = std::numeric_limits<int>::max();
update->max_eval = std::numeric_limits<int>::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;
}

52
src/USER-MDI/mdi_engine.h Normal file
View File

@ -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.
*/

View File

@ -13,6 +13,7 @@
#include "lammps.h" #include "lammps.h"
#include "input.h" #include "input.h"
#include "mdi_interface.h"
#include <mpi.h> #include <mpi.h>
#include <cstdlib> #include <cstdlib>
@ -35,6 +36,21 @@ int main(int argc, char **argv)
{ {
MPI_Init(&argc,&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. // enable trapping selected floating point exceptions.
// this uses GNU extensions and is only tested on Linux // this uses GNU extensions and is only tested on Linux
// therefore we make it depend on -D_GNU_SOURCE, too. // therefore we make it depend on -D_GNU_SOURCE, too.
@ -49,13 +65,13 @@ int main(int argc, char **argv)
#ifdef LAMMPS_EXCEPTIONS #ifdef LAMMPS_EXCEPTIONS
try { try {
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); LAMMPS *lammps = new LAMMPS(argc,argv,lammps_comm);
lammps->input->file(); lammps->input->file();
delete lammps; delete lammps;
} catch(LAMMPSAbortException &ae) { } catch(LAMMPSAbortException &ae) {
MPI_Abort(ae.universe, 1); MPI_Abort(ae.universe, 1);
} catch(LAMMPSException &e) { } catch(LAMMPSException &e) {
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(lammps_comm);
MPI_Finalize(); MPI_Finalize();
exit(1); exit(1);
} catch(fmt::format_error &fe) { } catch(fmt::format_error &fe) {
@ -65,7 +81,7 @@ int main(int argc, char **argv)
} }
#else #else
try { try {
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); LAMMPS *lammps = new LAMMPS(argc,argv,lammps_comm);
lammps->input->file(); lammps->input->file();
delete lammps; delete lammps;
} catch(fmt::format_error &fe) { } catch(fmt::format_error &fe) {
@ -74,6 +90,6 @@ int main(int argc, char **argv)
exit(1); exit(1);
} }
#endif #endif
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(lammps_comm);
MPI_Finalize(); MPI_Finalize();
} }

34
src/mdi_interface.h Normal file
View File

@ -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