diff --git a/src/Make.py b/src/Make.py index 2c8e479ac..29bcd930c 100755 --- a/src/Make.py +++ b/src/Make.py @@ -1,2067 +1,2067 @@ #!/usr/bin/env python2 # Make.py tool for managing packages and their auxiliary libs, # auto-editing machine Makefiles, and building LAMMPS # Syntax: Make.py -h (for help) # Notes: needs python 2.7 (not Python 3) import sys,os,commands,re,copy # switch abbrevs # switch classes = created class for each switch # lib classes = auxiliary package libs # build classes = build options with defaults # make classes = makefile options with no defaults # setargs = makefile settings # actionargs = allowed actions (also lib-dir and machine) abbrevs = "adhjmoprsv" switchclasses = ("actions","dir","help","jmake","makefile", "output","packages","redo","settings","verbose") libclasses = ("atc","awpmd","colvars","cuda","gpu", "meam","poems","qmmm","reax","voronoi") buildclasses = ("intel","kokkos") makeclasses = ("cc","mpi","fft","jpg","png") setargs = ("gzip","#gzip","ffmpeg","#ffmpeg","smallbig","bigbig","smallsmall") actionargs = ("lib-all","file","clean","exe") # ---------------------------------------------------------------- # functions # ---------------------------------------------------------------- # if flag = 1, print str and exit # if flag = 0, print str as warning and do not exit def error(str,flag=1): if flag: print "ERROR:",str sys.exit() else: print "WARNING:",str # store command-line args as sw = dict of key/value # key = switch word, value = list of following args # order = list of switches in order specified # enforce no switch more than once def parse_args(args): narg = len(args) sw = {} order = [] iarg = 0 while iarg < narg: if args[iarg][0] != '-': error("Arg %s is not a switch" % args[iarg]) switch = args[iarg][1:] if switch in sw: error("Duplicate switch %s" % args[iarg]) order.append(switch) first = iarg+1 last = first while last < narg and args[last][0] != '-': last += 1 sw[switch] = args[first:last] iarg = last return sw,order # convert info in switches dict back to a string, in switch_order def switch2str(switches,switch_order): txt = "" for switch in switch_order: if txt: txt += ' ' txt += "-%s" % switch txt += ' ' + ' '.join(switches[switch]) return txt # check if compiler works with ccflags on dummy one-line tmpauto.cpp file # return 1 if successful, else 0 # warn = 1 = print warning if not successful, warn = 0 = no warning # NOTE: unrecognized -override-limits can leave verride-limits file def compile_check(compiler,ccflags,warn): open("tmpauto.cpp",'w').write("int main(int, char **) {}") str = "%s %s -c tmpauto.cpp" % (compiler,ccflags) txt = commands.getoutput(str) flag = 1 if txt or not os.path.isfile("tmpauto.o"): flag = 0 if warn: print str if txt: print txt else: print "compile produced no output" os.remove("tmpauto.cpp") if os.path.isfile("tmpauto.o"): os.remove("tmpauto.o") return flag # check if linker works with linkflags on tmpauto.o file # return 1 if successful, else 0 # warn = 1 = print warning if not successful, warn = 0 = no warning def link_check(linker,linkflags,warn): open("tmpauto.cpp",'w').write("int main(int, char **) {}") str = "%s %s -o tmpauto tmpauto.cpp" % (linker,linkflags) txt = commands.getoutput(str) flag = 1 if txt or not os.path.isfile("tmpauto"): flag = 0 if warn: print str if txt: print txt else: print "link produced no output" os.remove("tmpauto.cpp") if os.path.isfile("tmpauto"): os.remove("tmpauto") return flag # ---------------------------------------------------------------- # switch classes, one per single-letter switch # ---------------------------------------------------------------- # actions class Actions: def __init__(self,list): self.inlist = list[:] def help(self): return """ -a action1 action2 ... possible actions = lib-all, lib-dir, file, clean, exe or machine machine is a Makefile.machine suffix actions can be specified in any order each action can appear only once lib-dir can appear multiple times for different dirs some actions depend on installed packages installed packages = currently installed + result of -p switch actions are invoked in this order, independent of specified order (1) lib-all or lib-dir = build auxiliary libraries lib-all builds all auxiliary libs needed by installed packages lib-dir builds a specific lib whether package installed or not dir is any dir in lib directory (atc, cuda, meam, etc) except linalg (2) file = create src/MAKE/MINE/Makefile.auto use -m switch for Makefile.machine to start from, else use existing Makefile.auto adds settings needed for installed accelerator packages existing Makefile.auto is NOT changed unless "file" action is specified (3) clean = invoke "make clean-auto" to insure full build useful if compiler flags have changed (4) exe or machine = build LAMMPS machine can be any existing Makefile.machine suffix machine is converted to "exe" action, as well as: "-m machine" is added if -m switch is not specified "-o machine" is added if -o switch is not specified if either "-m" or "-o" are specified, they are not overridden does not invoke any lib builds, since libs could be previously built exe always builds using src/MAKE/MINE/Makefile.auto if file action also specified, it creates Makefile.auto else if -m switch specified, existing Makefile.machine is copied to create Makefile.auto else Makefile.auto must already exist and is not changed produces src/lmp_auto, or error message if unsuccessful use -o switch to copy src/lmp_auto to new filename """ def check(self): if not self.inlist: error("-a args are invalid") alist = [] machine = 0 nlib = 0 for one in self.inlist: if one in alist: error("An action is duplicated") if one.startswith("lib-"): lib = one[4:] if lib != "all" and lib not in libclasses: error("Actions are invalid") alist.insert(nlib,one) nlib += 1 elif one == "file": if nlib == 0: alist.insert(0,"file") else: alist.insert(1,"file") elif one == "clean": if nlib == 0: alist.insert(0,"clean") elif "file" not in alist: alist.insert(1,"clean") else: alist.insert(2,"clean") elif one == "exe": if machine == 0: alist.append("exe") else: error("Actions are invalid") machine = 1 # one action can be unknown in case is a machine (checked in setup) elif machine == 0: alist.append(one) machine = 1 else: error("Actions are invalid") self.alist = alist # dedup list of actions concatenated from two lists # current self.inlist = specified -a switch + redo command -a switch # specified exe/machine action replaces redo exe/machine action # operates on and replaces self.inlist def dedup(self): alist = [] exemachine = 0 for one in self.inlist: if one == "exe" or (one not in actionargs and not one.startswith("lib-")): if exemachine: continue exemachine = 1 if one not in alist: alist.append(one) self.inlist = alist # if last action is unknown, assume machine and convert to exe # only done if action is a suffix for an existing Makefile.machine # return machine if conversion done, else None def setup(self): machine = self.alist[-1] if machine in actionargs or machine.startswith("lib-"): return None make = MakeReader(machine,2) self.alist[-1] = "exe" return machine # build one or more auxiliary package libraries def lib(self,suffix): if suffix != "all": print "building",suffix,"library ..." str = "%s.build()" % suffix exec(str) else: final = packages.final for one in packages.lib: if final[one]: if "user" in one: pkg = one[5:] else: pkg = one print "building",pkg,"library ..." str = "%s.build()" % pkg exec(str) # read Makefile.machine # if caller = "file", edit via switches # if caller = "exe", just read # write out new Makefile.auto def file(self,caller): # if caller = "file", create from mpi or read from makefile.machine or auto # if caller = "exe" and "file" action already invoked, read from auto # if caller = "exe" and no "file" action, read from makefile.machine or auto if caller == "file": if makefile and makefile.machine == "none": if cc and mpi: machine = "mpi" else: error("Cannot create makefile unless -cc and -mpi are used") elif makefile: machine = makefile.machine else: machine = "auto" elif caller == "exe" and "file" in self.alist: machine = "auto" elif caller == "exe" and "file" not in self.alist: if makefile and makefile.machine == "none": error("Cannot build with makefile = none") elif makefile: machine = makefile.machine else: machine = "auto" make = MakeReader(machine,1) # change makefile settings to user specifications precompiler = "" if caller == "file": # add compiler/linker and default CCFLAGS,LINKFLAGS # if cc.wrap, add wrapper setting for nvcc or mpi = ompi/mpich # if cc.wwrap, add 2nd wrapper setting for mpi = ompi/mpich # precompiler = env variable setting for OpenMPI wrapper compiler if cc: make.setvar("CC",cc.compiler) make.setvar("LINK",cc.compiler) if cc.wrap: abbrev = cc.abbrev if abbrev == "nvcc": make.addvar("CC","-ccbin=%s" % cc.wrap) make.addvar("LINK","-ccbin=%s" % cc.wrap) elif abbrev == "mpi": txt = commands.getoutput("mpicxx -show") if "-lmpich" in txt: make.addvar("CC","-cxx=%s" % cc.wrap) make.addvar("LINK","-cxx=%s" % cc.wrap) elif "-lmpi" in txt: make.addvar("OMPI_CXX",cc.wrap,"cc") precompiler = "env OMPI_CXX=%s " % cc.wrap else: error("Could not add MPI wrapper compiler, " + "did not recognize OpenMPI or MPICH") if cc.wwrap: txt = commands.getoutput("mpicxx -show") if "-lmpich" in txt: make.addvar("CC","-Xcompiler -cxx=%s" % cc.wwrap) make.addvar("LINK","-Xcompiler -cxx=%s" % cc.wwrap) elif "-lmpi" in txt: make.addvar("OMPI_CXX",cc.wwrap,"cc") precompiler = "env OMPI_CXX=%s " % cc.wwrap else: error("Could not add MPI wrapper compiler, " + "did not recognize OpenMPI or MPICH") make.setvar("CCFLAGS","-g") make.addvar("CCFLAGS","-O3") make.setvar("LINKFLAGS","-g") make.addvar("LINKFLAGS","-O") # add MPI settings if mpi: make.delvar("MPI_INC","*") make.delvar("MPI_PATH","*") make.delvar("MPI_LIB","*") if mpi.style == "mpi": make.addvar("MPI_INC","-DMPICH_SKIP_MPICXX") make.addvar("MPI_INC","-DOMPI_SKIP_MPICXX=1") elif mpi.style == "mpich": make.addvar("MPI_INC","-DMPICH_SKIP_MPICXX") make.addvar("MPI_INC","-DOMPI_SKIP_MPICXX=1") if mpi.dir: make.addvar("MPI_INC","-I%s/include" % mpi.dir) if mpi.dir: make.addvar("MPI_PATH","-L%s/lib" % mpi.dir) make.addvar("MPI_LIB","-lmpich") make.addvar("MPI_LIB","-lmpl") make.addvar("MPI_LIB","-lpthread") elif mpi.style == "ompi": make.addvar("MPI_INC","-DMPICH_SKIP_MPICXX") make.addvar("MPI_INC","-DOMPI_SKIP_MPICXX=1") if mpi.dir: make.addvar("MPI_INC","-I%s/include" % mpi.dir) if mpi.dir: make.addvar("MPI_PATH","-L%s/lib" % mpi.dir) make.addvar("MPI_LIB","-lmpi") make.addvar("MPI_LIB","-lmpi_cxx") elif mpi.style == "serial": make.addvar("MPI_INC","-I../STUBS") make.addvar("MPI_PATH","-L../STUBS") make.addvar("MPI_LIB","-lmpi_stubs") # add accelerator package CCFLAGS and LINKFLAGS and variables # pre = "" if compiler not nvcc, # else "-Xcompiler " to pass flag thru to wrapper compiler compiler = precompiler + ' '.join(make.getvar("CC")) linker = precompiler + ' '.join(make.getvar("LINK")) if "nvcc" in compiler: pre = "-Xcompiler " else: pre = "" final = packages.final if final["opt"]: if compile_check(compiler,pre + "-restrict",0): make.addvar("CCFLAGS",pre + "-restrict") if final["user-omp"]: if compile_check(compiler,pre + "-restrict",0): make.addvar("CCFLAGS",pre + "-restrict") #if "nvcc" not in compiler: if compile_check(compiler,pre + "-fopenmp",1): make.addvar("CCFLAGS",pre + "-fopenmp") make.addvar("LINKFLAGS",pre + "-fopenmp") if final["user-intel"]: if intel.mode == "cpu": if compile_check(compiler,pre + "-fopenmp",1): make.addvar("CCFLAGS",pre + "-fopenmp") make.addvar("LINKFLAGS",pre + "-fopenmp") make.addvar("CCFLAGS",pre + "-DLAMMPS_MEMALIGN=64") if compile_check(compiler,pre + "-restrict",1): make.addvar("CCFLAGS",pre + "-restrict") if compile_check(compiler,pre + "-xHost",1): make.addvar("CCFLAGS",pre + "-xHost") make.addvar("LINKFLAGS",pre + "-xHost") if compile_check(compiler,pre + "-fno-alias",1): make.addvar("CCFLAGS",pre + "-fno-alias") if compile_check(compiler,pre + "-ansi-alias",1): make.addvar("CCFLAGS",pre + "-ansi-alias") if compile_check(compiler,pre + "-override-limits",1): make.addvar("CCFLAGS",pre + "-override-limits") make.delvar("CCFLAGS",pre + "-DLMP_INTEL_OFFLOAD") make.delvar("LINKFLAGS",pre + "-offload") elif intel.mode == "phi": if compile_check(compiler,pre + "-fopenmp",1): make.addvar("CCFLAGS",pre + "-fopenmp") make.addvar("LINKFLAGS",pre + "-fopenmp") make.addvar("CCFLAGS",pre + "-DLAMMPS_MEMALIGN=64") if compile_check(compiler,pre + "-restrict",1): make.addvar("CCFLAGS",pre + "-restrict") if compile_check(compiler,pre + "-xHost",1): make.addvar("CCFLAGS",pre + "-xHost") make.addvar("CCFLAGS",pre + "-DLMP_INTEL_OFFLOAD") if compile_check(compiler,pre + "-fno-alias",1): make.addvar("CCFLAGS",pre + "-fno-alias") if compile_check(compiler,pre + "-ansi-alias",1): make.addvar("CCFLAGS",pre + "-ansi-alias") if compile_check(compiler,pre + "-override-limits",1): make.addvar("CCFLAGS",pre + "-override-limits") if compile_check(compiler,pre + '-offload-option,mic,compiler,' + '"-fp-model fast=2 -mGLOB_default_function_attrs=' + '\\"gather_scatter_loop_unroll=4\\""',1): make.addvar("CCFLAGS",pre + '-offload-option,mic,compiler,' + '"-fp-model fast=2 -mGLOB_default_function_attrs=' + '\\"gather_scatter_loop_unroll=4\\""') if link_check(linker,"-offload",1): make.addvar("LINKFLAGS",pre + "-offload") if final["kokkos"]: if kokkos.mode == "omp": make.addvar("OMP","yes","lmp") make.delvar("CUDA") make.delvar("MIC") elif kokkos.mode == "cuda": if "nvcc" not in compiler: error("Kokkos/cuda build appears to not be " + "using NVIDIA nvcc compiler",0) make.addvar("OMP","yes","lmp") make.addvar("CUDA","yes","lmp") make.delvar("MIC") if kokkos.archflag: make.delvar("CCFLAGS","-arch=sm_*") make.addvar("CCFLAGS","-arch=sm_%s" % kokkos.arch) elif kokkos.mode == "phi": make.addvar("OMP","yes","lmp") make.addvar("MIC","yes","lmp") make.delvar("CUDA") # add LMP settings if settings: list = settings.inlist for one in list: if one == "gzip": make.addvar("LMP_INC","-DLAMMPS_GZIP") elif one == "#gzip": make.delvar("LMP_INC","-DLAMMPS_GZIP") elif one == "ffmpeg": make.addvar("LMP_INC","-DLAMMPS_FFMPEG") elif one == "#ffmpeg": make.delvar("LMP_INC","-DLAMMPS_FFMPEG") elif one == "smallbig": make.delvar("LMP_INC","-DLAMMPS_BIGBIG") make.delvar("LMP_INC","-DLAMMPS_SMALLSMALL") elif one == "bigbig": make.delvar("LMP_INC","-DLAMMPS_SMALLBIG") make.delvar("LMP_INC","-DLAMMPS_SMALLSMALL") make.addvar("LMP_INC","-DLAMMPS_BIGBIG") elif one == "smallsmall": make.delvar("LMP_INC","-DLAMMPS_SMALLBIG") make.delvar("LMP_INC","-DLAMMPS_BIGBIG") make.addvar("LMP_INC","-DLAMMPS_SMALLSMALL") # add FFT, JPG, PNG settings if fft: make.delvar("FFT_INC","*") make.delvar("FFT_PATH","*") make.delvar("FFT_LIB","*") if fft.mode == "none": make.addvar("FFT_INC","-DFFT_NONE") else: make.addvar("FFT_INC","-DFFT_%s" % fft.mode.upper()) make.addvar("FFT_LIB",fft.lib) if fft.dir: make.addvar("FFT_INC","-I%s/include" % fft.dir) make.addvar("FFT_PATH","-L%s/lib" % fft.dir) else: if fft.incdir: make.addvar("FFT_INC","-I%s" % fft.incdir) if fft.libdir: make.addvar("FFT_PATH","-L%s" % fft.libdir) if jpg: if jpg.on == 0: make.delvar("LMP_INC","-DLAMMPS_JPEG") make.delvar("JPG_LIB","-ljpeg") else: make.addvar("LMP_INC","-DLAMMPS_JPEG") make.addvar("JPG_LIB","-ljpeg") if jpg.dir: make.addvar("JPG_INC","-I%s/include" % jpg.dir) make.addvar("JPG_PATH","-L%s/lib" % jpg.dir) else: if jpg.incdir: make.addvar("JPG_INC","-I%s" % jpg.incdir) if jpg.libdir: make.addvar("JPG_PATH","-L%s" % jpg.libdir) if png: if png.on == 0: make.delvar("LMP_INC","-DLAMMPS_PNG") make.delvar("JPG_LIB","-lpng") else: make.addvar("LMP_INC","-DLAMMPS_PNG") make.addvar("JPG_LIB","-lpng") if png.dir: make.addvar("JPG_INC","-I%s/include" % png.dir) make.addvar("JPG_PATH","-L%s/lib" % png.dir) else: if png.incdir: make.addvar("JPG_INC","-I%s" % png.incdir) if png.libdir: make.addvar("JPG_PATH","-L%s" % png.libdir) # set self.stubs if Makefile.auto uses STUBS lib in MPI settings if "-lmpi_stubs" in make.getvar("MPI_LIB"): self.stubs = 1 else: self.stubs = 0 # write out Makefile.auto # unless caller = "exe" and "file" action already invoked if caller == "file" or "file" not in self.alist: make.write("%s/MAKE/MINE/Makefile.auto" % dir.src,1) print "Created src/MAKE/MINE/Makefile.auto" # test full compile and link # unless caller = "file" and "exe" action will be invoked later if caller == "file" and "exe" in self.alist: return compiler = precompiler + ' '.join(make.getvar("CC")) ccflags = ' '.join(make.getvar("CCFLAGS")) linker = precompiler + ' '.join(make.getvar("LINK")) linkflags = ' '.join(make.getvar("LINKFLAGS")) if not compile_check(compiler,ccflags,1): error("Test of compilation failed") if not link_check(linker,linkflags,1): error("Test of link failed") # invoke "make clean-auto" to force clean before build def clean(self): str = "cd %s; make clean-auto" % dir.src commands.getoutput(str) if verbose: print "Performed make clean-auto" # build LAMMPS using Makefile.auto and -j setting # invoke self.file() first, to test makefile compile/link # delete existing lmp_auto, so can detect if build fails # build STUBS lib (if unbuilt) if Makefile.auto MPI settings need it def exe(self): self.file("exe") commands.getoutput("cd %s; rm -f lmp_auto" % dir.src) if self.stubs and not os.path.isfile("%s/STUBS/libmpi_stubs.a" % dir.src): print "building serial STUBS library ..." str = "cd %s/STUBS; make clean; make" % dir.src txt = commands.getoutput(str) if not os.path.isfile("%s/STUBS/libmpi_stubs.a" % dir.src): print txt error('Unsuccessful "make stubs"') print "Created src/STUBS/libmpi_stubs.a" if jmake: str = "cd %s; make -j %d auto" % (dir.src,jmake.n) else: str = "cd %s; make auto" % dir.src txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/lmp_auto" % dir.src): if not verbose: print txt error('Unsuccessful "make auto"') elif not output: print "Created src/lmp_auto" # dir switch class Dir: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] def help(self): return """ -d dir dir = LAMMPS home dir if -d not specified, working dir must be lammps/src """ def check(self): if self.inlist != None and len(self.inlist) != 1: error("-d args are invalid") # if inlist = None, check that cwd = lammps/src # store cwd and lammps dir # derive src,make,lib dirs from lammps dir # check that they all exist def setup(self): self.cwd = os.getcwd() if self.inlist == None: self.lammps = ".." else: self.lammps = self.inlist[0] self.lammps = os.path.realpath(self.lammps) self.src = self.lammps + "/src" self.make = self.lammps + "/src/MAKE" self.lib = self.lammps + "/lib" if not os.path.isdir(self.lammps): error("LAMMPS home dir is invalid") if not os.path.isdir(self.src): error("LAMMPS src dir is invalid") if not os.path.isdir(self.lib): error("LAMMPS lib dir is invalid") # help switch class Help: def __init__(self,list): pass def help(self): return """ Syntax: Make.py switch args ... switches can be listed in any order help switch: -h prints help and syntax for all other specified switches switch for actions: -a lib-all, lib-dir, clean, file, exe or machine list one or more actions, in any order - machine is a Makefile.machine suffix, must be last if used + machine is a Makefile.machine suffix one-letter switches: -d (dir), -j (jmake), -m (makefile), -o (output), -p (packages), -r (redo), -s (settings), -v (verbose) switches for libs: -atc, -awpmd, -colvars, -cuda -gpu, -meam, -poems, -qmmm, -reax, -voronoi switches for build and makefile options: -intel, -kokkos, -cc, -mpi, -fft, -jpg, -png """ # jmake switch class Jmake: def __init__(self,list): self.inlist = list[:] def help(self): return """ -j N use N procs for performing parallel make commands used when building a lib or LAMMPS itself if -j not specified, serial make commands run on single core """ def check(self): if len(self.inlist) != 1: error("-j args are invalid") if not self.inlist[0].isdigit(): error("-j args are invalid") n = int(self.inlist[0]) if n <= 0: error("-j args are invalid") self.n = n # makefile switch class Makefile: def __init__(self,list): self.inlist = list[:] def help(self): return """ -m machine use Makefile.machine under src/MAKE as starting point to create Makefile.auto if machine = "none", file action will create Makefile.auto from scratch must use -cc and -mpi switches to specify compiler and MPI if -m not specified, file/exe actions alter existing Makefile.auto """ def check(self): if len(self.inlist) != 1: error("-m args are invalid") self.machine = self.inlist[0] # output switch class Output: def __init__(self,list): self.inlist = list[:] def help(self): return """ -o machine copy final src/lmp_auto to lmp_machine in working dir if -o not specified, exe action only produces src/lmp_auto """ def check(self): if len(self.inlist) != 1: error("-o args are invalid") self.machine = self.inlist[0] # packages switch class Packages: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] def help(self): return """ -p = package1 package2 ... list of packages to install or uninstall in order specified operates on set of packages currently installed valid package names: and LAMMPS standard or user package (type "make package" to see list) prefix by yes/no to install/uninstall (see abbrevs) yes-molecule, yes-user-atc, no-molecule, no-user-atc can use LAMMPS categories (type "make package" to see list) all = all standard and user packages (also none = no-all) std (or standard) = all standard packages user = all user packages lib = all standard and user packages with auxiliary libs can abbreviate package names and yes/no omp = user-omp = yes-user-omp ^omp = ^user-omp = no-user-omp user = yes-user, ^user = no-user all = yes-all, ^all = none = no-all when action performed, list is processed in order, as if typed "make yes/no" for each if "orig" or "original" is last package in list, set of installed packages will be restored to original (current) list after "build" action is performed if -p not specified, currently installed packages are not changed """ def check(self): if self.inlist != None and not self.inlist: error("-p args are invalid") def setup(self): # extract package lists from src/Makefile make = MakeReader("%s/Makefile" % dir.src) std = make.getvar("PACKAGE") user = make.getvar("PACKUSER") lib = make.getvar("PACKLIB") all = std + user # plist = command line args expanded to yes-package or no-package plist = [] if self.inlist: for one in self.inlist: if one in std: plist.append("yes-%s" % one) elif one in user: plist.append("yes-%s" % one) elif "user-"+one in user: plist.append("yes-user-%s" % one) elif one == "std" or one == "standard" or one == "user" or \ one == "lib" or one == "all": plist.append("yes-%s" % one) elif one.startswith("yes-"): if one[4:] in std: plist.append("yes-%s" % one[4:]) elif one[4:] in user: plist.append("yes-%s" % one[4:]) elif "user-"+one[4:] in user: plist.append("yes-user-%s" % one[4:]) elif one == "yes-std" or one == "yes-standard" or \ one == "yes-user" or one == "yes-lib" or one == "yes-all": plist.append("yes-%s" % one[4:]) else: error("Invalid package name %s" % one) elif one.startswith("no-"): if one[3:] in std: plist.append("no-%s" % one[3:]) elif one[3:] in user: plist.append("no-%s" % one[3:]) elif "user-"+one[3:] in user: plist.append("no-user-%s" % one[3:]) elif one == "no-std" or one == "no-standard" or one == "no-user" or \ one == "no-lib" or one == "no-all": plist.append("no-%s" % one[3:]) else: error("Invalid package name %s" % one) elif one.startswith('^'): if one[1:] in std: plist.append("no-%s" % one[1:]) elif one[1:] in user: plist.append("no-%s" % one[1:]) elif "user-"+one[1:] in user: plist.append("no-user-%s" % one[1:]) elif one == "^std" or one == "^standard" or one == "^user" or \ one == "^lib" or one == "^all": plist.append("no-%s" % one[1:]) else: error("Invalid package name %s" % one) elif one == "none": plist.append("no-all") elif one == "orig": plist.append(one) else: error("Invalid package name %s" % one) if "orig" in plist and plist.index("orig") != len(plist)-1: error('-p orig arg must be last') if plist.count("orig") > 1: error('-p orig arg must be last') # original = dict of all packages # key = package name, value = 1 if currently installed, else 0 original = {} str = "cd %s; make ps" % dir.src output = commands.getoutput(str).split('\n') pattern = "Installed\s+(\w+): package (\S+)" for line in output: m = re.search(pattern,line) if not m: continue pkg = m.group(2).lower() if pkg not in all: error('Package list does not math "make ps" results') if m.group(1) == "NO": original[pkg] = 0 elif m.group(1) == "YES": original[pkg] = 1 # final = dict of all packages after plist applied to original # key = package name, value = 1 if installed, else 0 final = copy.deepcopy(original) for i,one in enumerate(plist): if "yes" in one: pkg = one[4:] yes = 1 else: pkg = one[3:] yes = 0 if pkg in all: final[pkg] = yes elif pkg == "std": for pkg in std: final[pkg] = yes elif pkg == "user": for pkg in user: final[pkg] = yes elif pkg == "lib": for pkg in lib: final[pkg] = yes elif pkg == "all": for pkg in all: final[pkg] = yes self.std = std self.user = user self.lib = lib self.all = all self.plist = plist self.original = original self.final = final # install packages in plist def install(self): if self.plist: print "Installing packages ..." for one in self.plist: if one == "orig": continue commands.getoutput("cd %s; make %s" % (dir.src,one)) if self.plist and verbose: txt = commands.getoutput("cd %s; make ps" % dir.src) print "Package status after installation:" print txt # restore packages to original list if requested # order of re-install should not matter matter b/c of Depend.sh def uninstall(self): if not self.plist or self.plist[-1] != "orig": return print "Restoring packages to original state ..." commands.getoutput("cd %s; make no-all" % dir.src) for one in self.all: if self.original[one]: commands.getoutput("cd %s; make yes-%s" % (dir.src,one)) if verbose: txt = commands.getoutput("cd %s; make ps" % dir.src) print "Restored package status:" print txt # redo switch class Redo: def __init__(self,list): self.inlist = list[:] def help(self): return """ -r file label1 label2 ... all args are optional invoke Make.py commands from a file other specified switches are merged with file commands (see below) redo file format: blank lines and lines starting with "#" are skipped other lines are treated as commands each command is a list of Make.py args, as if typed at command-line commands can have leading label, followed by ":" commands cannot contain a "-r" switch if no args, execute previous command, which is stored in src/Make.py.last if one arg, execute all commands from specified file unlabeled or labeled commands are all executed if multiple args, execute only matching labeled commands from file if other switches are specified, if file command does not have the switch, it is added if file command has the switch, the specified switch replaces it except if -a (action) switch is both specified and in the file command, two sets of actions are merged and duplicates removed if both switches have "exe or machine" action, the specified exe/machine overrides the file exe/machine """ def check(self): if len(self.inlist) == 0: self.dir = 1 self.file = "Make.py.last" self.labels = [] else: self.dir = 0 self.file = self.inlist[0] self.labels = self.inlist[1:] # read redo file # self.commands = list of commands to execute def setup(self): file = self.file if not os.path.isfile(file): error("Redo file %s does not exist" % file) lines = open(file,'r').readlines() cmdlines = [] for line in lines: line = line.strip() if not line or line[0] == '#' : continue cmdlines.append(line) # if no labels, add all file commands to command list # if labels, make a dict with key = label, value = command # and discard unlabeled commands dict = {} commands = [] for line in cmdlines: words = line.split() if "-r" in words: error("Redo command cannot contain -r switch") if words[0][-1] == ':': label = words[0][:-1] else: label = None if not self.labels: if label: commands.append(' '.join(words[1:])) else: commands.append(line) else: if not label: continue dict[label] = ' '.join(words[1:]) # extract labeled commands from dict and add to command list for label in self.labels: if label not in dict: error("Redo label not in redo file") commands.append(dict[label]) self.commands = commands # settings switch class Settings: def __init__(self,list): self.inlist = list[:] def help(self): return """ -s set1 set2 ... possible settings = gzip smallbig bigbig smallsmall add each setting as LAMMPS setting to created Makefile.auto if -s not specified, no settings are changed in Makefile.auto """ def check(self): if not self.inlist: error("-s args are invalid") for one in self.inlist: if one not in setargs: error("-s args are invalid") # verbose switch class Verbose: def __init__(self,list): self.inlist = list[:] def help(self): return """ -v (no arguments) produce verbose output as Make.py executes if -v not specified, minimal output is produced """ def check(self): if len(self.inlist): error("-v args are invalid") # ---------------------------------------------------------------- # lib classes, one per LAMMPS auxiliary lib # ---------------------------------------------------------------- # ATC lib class ATC: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "g++" self.lammpsflag = 0 def help(self): return """ -atc make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = g++) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-atc args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-atc args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-atc args are invalid") def build(self): libdir = dir.lib + "/atc" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libatc.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/atc library") else: print "Created lib/atc library" # AWPMD lib class AWPMD: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "mpicc" self.lammpsflag = 0 def help(self): return """ -awpmd make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = mpicc) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-awpmd args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-awpmd args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-awpmd args are invalid") def build(self): libdir = dir.lib + "/awpmd" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libawpmd.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/awpmd library") else: print "Created lib/awpmd library" # COLVARS lib class COLVARS: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "g++" self.lammpsflag = 0 def help(self): return """ -colvars make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = g++) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-colvars args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-colvars args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-colvars args are invalid") def build(self): libdir = dir.lib + "/colvars" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libcolvars.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/colvars library") else: print "Created lib/colvars library" # CUDA lib class CUDA: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.mode = "double" self.arch = "31" def help(self): return """ -cuda mode=double arch=31 all args are optional and can be in any order mode = double or mixed or single (def = double) arch = M (def = 31) M = 31 for Kepler M = 20 for CC2.0 (GF100/110, e.g. C2050,GTX580,GTX470) M = 21 for CC2.1 (GF104/114, e.g. GTX560, GTX460, GTX450) M = 13 for CC1.3 (GF200, e.g. C1060, GTX285) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-cuda args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-cuda args are invalid") if words[0] == "mode": self.mode = words[1] elif words[0] == "arch": self.arch = words[1] else: error("-cuda args are invalid") if self.mode != "double" and self.mode != "mixed" and \ self.mode != "single": error("-cuda args are invalid") if not self.arch.isdigit(): error("-cuda args are invalid") def build(self): libdir = dir.lib + "/cuda" commands.getoutput("cd %s; make clean" % libdir) if self.mode == "double": n = 2 elif self.mode == "mixed": n = 3 elif self.mode == "single": n = 1 if jmake: str = "cd %s; make -j %d precision=%d arch=%s" % \ (libdir,jmake.n,n,self.arch) else: str = str = "cd %s; make precision=%d arch=%s" % \ (libdir,n,self.arch) txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/liblammpscuda.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/cuda library") else: print "Created lib/cuda library" # GPU lib class GPU: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "linux.double" self.lammpsflag = self.modeflag = self.archflag = 0 def help(self): return """ -gpu make=suffix lammps=suffix2 mode=double arch=N all args are optional and can be in any order make = use Makefile.suffix (def = linux.double) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) mode = double or mixed or single (def = CUDA_PREC in makefile) arch = 31 (Kepler) or 21 (Fermi) (def = CUDA_ARCH in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-gpu args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-gpu args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 elif words[0] == "mode": self.mode = words[1] self.modeflag = 1 elif words[0] == "arch": self.arch = words[1] self.archflag = 1 else: error("-gpu args are invalid") if self.modeflag and (self.mode != "double" and self.mode != "mixed" and self.mode != "single"): error("-gpu args are invalid") if self.archflag and not self.arch.isdigit(): error("-gpu args are invalid") def build(self): libdir = dir.lib + "/gpu" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.modeflag: if self.mode == "double": make.setvar("CUDA_PRECISION","-D_DOUBLE_DOUBLE") elif self.mode == "mixed": make.setvar("CUDA_PRECISION","-D_SINGLE_DOUBLE") elif self.mode == "single": make.setvar("CUDA_PRECISION","-D_SINGLE_SINGLE") if self.archflag: make.setvar("CUDA_ARCH","-arch=sm_%s" % self.arch) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libgpu.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/gpu library") else: print "Created lib/gpu library" # MEAM lib class MEAM: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "gfortran" self.lammpsflag = 0 def help(self): return """ -meam make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = gfortran) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-meam args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-meam args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-meam args are invalid") def build(self): libdir = dir.lib + "/meam" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) # do not use -j for MEAM build, parallel build does not work str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libmeam.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/meam library") else: print "Created lib/meam library" # POEMS lib class POEMS: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "g++" self.lammpsflag = 0 def help(self): return """ -poems make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = g++) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-poems args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-poems args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-poems args are invalid") def build(self): libdir = dir.lib + "/poems" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libpoems.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/poems library") else: print "Created lib/poems library" # QMMM lib class QMMM: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "gfortran" self.lammpsflag = 0 def help(self): return """ -qmmm make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = gfortran) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-qmmm args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-qmmm args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-qmmm args are invalid") def build(self): libdir = dir.lib + "/qmmm" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libqmmm.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/qmmm library") else: print "Created lib/qmmm library" # REAX lib class REAX: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.make = "gfortran" self.lammpsflag = 0 def help(self): return """ -reax make=suffix lammps=suffix2 all args are optional and can be in any order make = use Makefile.suffix (def = gfortran) lammps = use Makefile.lammps.suffix2 (def = EXTRAMAKE in makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-reax args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-reax args are invalid") if words[0] == "make": self.make = words[1] elif words[0] == "lammps": self.lammps = words[1] self.lammpsflag = 1 else: error("-reax args are invalid") def build(self): libdir = dir.lib + "/reax" make = MakeReader("%s/Makefile.%s" % (libdir,self.make)) if self.lammpsflag: make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps) make.write("%s/Makefile.auto" % libdir) commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir) if jmake: str = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n) else: str = "cd %s; make -f Makefile.auto" % libdir txt = commands.getoutput(str) if verbose: print txt if not os.path.isfile("%s/libreax.a" % libdir) or \ not os.path.isfile("%s/Makefile.lammps" % libdir): if not verbose: print txt error("Unsuccessful build of lib/reax library") else: print "Created lib/reax library" # VORONOI lib class VORONOI: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.install = "" def help(self): return """ -voronoi install="-d dir -v version -g -b -i installdir -l incdir libdir" arg is optional, only needed if want to run install.py script install = args to use with lib/voronoi/install.py script must enclose in quotes since install.py args have switches install.py can download, build, install, setup links to the Voro++ library see lib/voronoi/README for details on Voro++ and using install.py """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-voronoi args are invalid") for one in self.inlist: words = one.split('=') if len(words) != 2: error("-voronoi args are invalid") if words[0] == "install": self.install = words[1] else: error("-voronoi args are invalid") def build(self): if not self.install: return libdir = dir.lib + "/voronoi" cmd = "cd %s; python install.py %s" % (libdir,self.install) txt = commands.getoutput(cmd) if verbose: print txt print "Created lib/voronoi library" # ---------------------------------------------------------------- # build classes for intel, kokkos build options # ---------------------------------------------------------------- # Intel class class Intel: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.mode = "cpu" def help(self): return """ -intel mode mode = cpu or phi (def = cpu) build Intel package for CPU or Xeon Phi """ def check(self): if self.inlist == None: return if len(self.inlist) != 1: error("-intel args are invalid") self.mode = self.inlist[0] if self.mode != "cpu" and self.mode != "phi": error("-intel args are invalid") # Kokkos class class Kokkos: def __init__(self,list): if list == None: self.inlist = None else: self.inlist = list[:] self.mode = "omp" self.archflag = 0 def help(self): return """ -kokkos mode arch=N mode is not optional, arch is optional mode = omp or cuda or phi (def = omp if -kokkos is not used) build Kokkos package for omp or cuda or phi arch = 31 (Kepler) or 21 (Fermi) (def = -arch setting in Makefile) """ def check(self): if self.inlist != None and len(self.inlist) == 0: error("-kokkos args are invalid") if self.inlist == None: return if len(self.inlist) < 1: error("-kokkos args are invalid") self.mode = self.inlist[0] if self.mode != "omp" and self.mode != "cuda" and self.mode != "phi": error("-kokkos args are invalid") for one in self.inlist[1:]: words = one.split('=') if len(words) != 2: error("-kokkos args are invalid") if words[0] == "arch": self.arch = words[1] self.archflag = 1 else: error("-kokkos args are invalid") # ---------------------------------------------------------------- # makefile classes for CC, MPI, JPG, PNG, FFT settings # ---------------------------------------------------------------- # Cc class class Cc: def __init__(self,list): self.inlist = list[:] self.compiler = self.abbrev = "" self.wrap = self.wabbrev = "" self.wwrap = "" def help(self): return """ -cc compiler wrap=wcompiler wwrap=wwcompiler change CC setting in makefile compiler is required, all other args are optional compiler = any string with g++ or icc or icpc or nvcc or mpi (or mpicxx, mpiCC, mpiicpc, etc) can be compiler name or full path to compiler mpi by itself is changed to mpicxx wcompiler = g++ or icc or icpc or mpi (or mpicxx, mpiCC, mpiicpc, etc) can only be used when compiler is a wrapper (mpi or nvcc) mpi and variants can only be used with compiler = nvcc mpi by itself is changed to mpicxx specify compiler for wrapper to use in CC setting wwcompiler = g++ or icc or icpc can only be used when wcompiler is itself a wrapper (mpi) specify compiler for wrapper of wrapper to use in CC setting """ def check(self): if len(self.inlist) < 1: error("-cc args are invalid") self.compiler = self.inlist[0] if self.compiler == "mpi": self.compiler = "mpicxx" self.abbrev = "mpi" elif self.compiler.startswith("mpi"): self.abbrev = "mpi" elif self.compiler == "g++" or self.compiler == "icc" or \ self.compiler == "icpc" or self.compiler == "nvcc": self.abbrev = self.compiler elif "mpi" in self.compiler: self.abbrev = "mpi" elif "g++" in self.compiler: self.abbrev = "g++" elif "icc" in self.compiler: self.abbrev = "icc" elif "icpc" in self.compiler: self.abbrev = "icpc" elif "nvcc" in self.compiler: self.abbrev = "nvcc" else: error("-cc args are invalid") for one in self.inlist[1:]: words = one.split('=') if len(words) != 2: error("-cc args are invalid") if words[0] == "wrap": if self.abbrev != "mpi" and self.abbrev != "nvcc": error("-cc compiler is not a wrapper") self.wrap = words[1] if self.wrap == "mpi": self.wrap = "mpicxx" self.wabbrev = "mpi" elif self.wrap.startswith("mpi"): self.wabbrev = "mpi" elif self.compiler == "g++" or self.compiler == "icc" or \ self.compiler == "icpc": self.wabbrev = self.wrap elif words[0] == "wwrap": self.wwrap = words[1] if self.wabbrev != "mpi": error("-cc wrap is not a wrapper") if self.wwrap != "g++" and self.wwrap != "icc" and self.wwrap != "icpc": error("-cc args are invalid") else: error("-cc args are invalid") # Mpi class class Mpi: def __init__(self,list): self.inlist = list[:] self.style = self.dir = "" def help(self): return """ -mpi style dir=path change MPI settings in makefile style is required, all other args are optional style = mpi or mpich or ompi or serial mpi = no MPI settings (assume compiler is MPI wrapper) mpich = use explicit settings for MPICH ompi = use explicit settings for OpenMPI serial = use settings for src/STUBS library dir = path for MPICH or OpenMPI directory add -I and -L settings for include and lib sub-dirs """ def check(self): if len(self.inlist) < 1: error("-mpi args are invalid") self.style = self.inlist[0] if self.style != "mpi" and self.style != "mpich" and \ self.style != "ompi" and self.style != "serial": error("-mpi args are invalid") for one in self.inlist[1:]: words = one.split('=') if len(words) != 2: error("-mpi args are invalid") if words[0] == "dir": self.dir = words[1] else: error("-mpi args are invalid") # Fft class class Fft: def __init__(self,list): self.inlist = list[:] self.dir = self.incdir = self.libdir = "" def help(self): return """ -fft mode lib=libname dir=homedir idir=incdir ldir=libdir change FFT settings in makefile mode is required, all other args are optional removes all current FFT variable settings mode = none or fftw or fftw3 of ... adds -DFFT_MODE setting lib = name of FFT library to link with (def is libname = mode) adds -lliblib setting, e.g. -llibfftw3 dir = home dir for include and library files (def = none) adds -Idir/include and -Ldir/lib settings if set, overrides idir and ldir args idir = dir for include file (def = none) adds -Iidir setting ldir = dir for library file (def = none) adds -Lldir setting """ def check(self): if not len(self.inlist): error("-fft args are invalid") self.mode = self.inlist[0] self.lib = "-l%s" % self.mode for one in self.inlist[1:]: words = one.split('=') if len(words) != 2: error("-fft args are invalid") if words[0] == "lib": self.lib = "-l%s" % words[1] elif words[0] == "dir": self.dir = words[1] elif words[0] == "idir": self.incdir = words[1] elif words[0] == "ldir": self.libdir = words[1] else: error("-fft args are invalid") # Jpg class class Jpg: def __init__(self,list): self.inlist = list[:] self.on = 1 self.dir = self.incdir = self.libdir = "" def help(self): return """ -jpg flag dir=homedir idir=incdir ldir=libdir change JPG settings in makefile all args are optional, flag must come first if specified flag = yes or no (def = yes) include or exclude JPEG support adds/removes -DLAMMPS_JPEG and -ljpeg settings dir = home dir for include and library files (def = none) adds -Idir/include and -Ldir/lib settings if set, overrides idir and ldir args idir = dir for include file (def = none) adds -Iidir setting ldir = dir for library file (def = none) adds -Lldir setting """ def check(self): for i,one in enumerate(self.inlist): if one == "no" and i == 0: self.on = 0 elif one == "yes" and i == 0: self.on = 1 else: words = one.split('=') if len(words) != 2: error("-jpeg args are invalid") if words[0] == "dir": self.dir = words[1] elif words[0] == "idir": self.incdir = words[1] elif words[0] == "ldir": self.libdir = words[1] else: error("-jpeg args are invalid") # Png class class Png: def __init__(self,list): self.inlist = list[:] self.on = 1 self.dir = self.incdir = self.libdir = "" def help(self): return """ -png flag dir=homedir idir=incdir ldir=libdir change PNG settings in makefile all args are optional, flag must come first if specified flag = yes or no (def = yes) include or exclude PNG support adds/removes -DLAMMPS_PNG and -lpng settings dir = home dir for include and library files (def = none) adds -Idir/include and -Ldir/lib settings if set, overrides idir and ldir args idir = dir for include file (def = none) adds -Iidir setting ldir = dir for library file (def = none) adds -Lldir setting """ def check(self): for i,one in enumerate(self.inlist): if one == "no" and i == 0: self.on = 0 elif one == "yes" and i == 0: self.on = 1 else: words = one.split('=') if len(words) != 2: error("-png args are invalid") if words[0] == "dir": self.dir = words[1] elif words[0] == "idir": self.incdir = words[1] elif words[0] == "ldir": self.libdir = words[1] else: error("-png args are invalid") # ---------------------------------------------------------------- # auxiliary classes # ---------------------------------------------------------------- # read, tweak, and write a Makefile class MakeReader: # read a makefile # flag = 0 if file is full path name # flag = 1,2 if file is suffix for any Makefile.machine under src/MAKE # look for this file in same order that src/Makefile does # if flag = 1, read the file # if flag = 2, just check if file exists def __init__(self,file,flag=0): if flag == 0: if not os.path.isfile(file): error("Makefile %s does not exist" % file) lines = open(file,'r').readlines() else: mfile = "%s/MAKE/MINE/Makefile.%s" % (dir.src,file) if not os.path.isfile(mfile): mfile = "%s/MAKE/Makefile.%s" % (dir.src,file) if not os.path.isfile(mfile): mfile = "%s/MAKE/OPTIONS/Makefile.%s" % (dir.src,file) if not os.path.isfile(mfile): mfile = "%s/MAKE/MACHINES/Makefile.%s" % (dir.src,file) if not os.path.isfile(mfile): error("Makefile.%s does not exist" % file) if flag == 1: lines = open(mfile,'r').readlines() else: return # scan lines of makefile # if not a variable line, just copy to newlines # if a variable line, concatenate any continuation lines # convert variable to var dict entry: key = name, value = list of words # discard any portion of value string with a comment char # varinfo = list of variable info: (name, name with whitespace for print) # add index into varinfo to newlines # ccindex = index of "CC =" line, to add OMPI var before it # lmpindex = index of "LAMMPS-specific settings" line to add KOKKOS vars before it var = {} varinfo = [] newlines = [] pattern = "(\S+\s+=\s+)(.*)" multiline = 0 self.ccindex = self.lmpindex = 0 for line in lines: line = line[:-1] if "CC =" in line: self.ccindex = len(newlines) if "LAMMPS-specific settings" in line: self.lmpindex = len(newlines) if multiline: if '#' in line: line = line[:line.find('#')] morevalues = line.split() values = values[:-1] + morevalues if values[-1] != '\\': var[name] = values multiline = 0 newlines.append(str(len(varinfo))) varinfo.append((name,namewhite)) continue varflag = 1 if len(line.strip()) == 0: varflag = 0 elif line.lstrip()[0] == '#': varflag = 0 else: m = re.match(pattern,line) if not m: varflag = 0 if varflag: namewhite = m.group(1) name = namewhite.split()[0] if name in var: error("Makefile variable %s appears more than once" % name) remainder = m.group(2) if '#' in remainder: remainder = remainder[:remainder.find('#')] values = remainder.split() if values and values[-1] == '\\': multiline = 1 else: var[name] = values newlines.append(str(len(varinfo))) varinfo.append((name,namewhite)) else: newlines.append(line) self.var = var self.varinfo = varinfo self.lines = newlines # return list of values associated with var # return None if var not defined def getvar(self,var): if var in self.var: return self.var[var] else: return None # set var to single value # if var not defined, error def setvar(self,var,value): if var not in self.var: error("Variable %s not in makefile" % var) self.var[var] = [value] # add value to var # do not add if value already defined by var # if var not defined, # create new variable using where # where="cc", line before "CC =" line, use ":=" # where="lmp", 2 lines before "LAMMPS-specific settings" line, use "=" def addvar(self,var,value,where=""): if var in self.var: if value not in self.var[var]: self.var[var].append(value) else: if not where: error("Variable %s with value %s is not in makefile" % (var,value)) if where == "cc": if not self.ccindex: error("No 'CC =' line in makefile to add variable") index = self.ccindex varwhite = "%s :=\t\t" % var elif where == "lmp": if not self.lmpindex: error("No 'LAMMPS-specific settings line' " + "in makefile to add variable") index = self.lmpindex - 2 varwhite = "%s =\t\t" % var self.var[var] = [value] varwhite = "%s =\t\t" % var self.lines.insert(index,str(len(self.varinfo))) self.varinfo.append((var,varwhite)) # if value = None, remove entire var # no need to update lines or varinfo, write() will ignore deleted vars # else remove value from var # value can have trailing '*' to remove wildcard match # if var or value not defined, ignore it def delvar(self,var,value=None): if var not in self.var: return if not value: del self.var[var] elif value and value[-1] != '*': if value not in self.var[var]: return self.var[var].remove(value) else: value = value[:-1] values = self.var[var] dellist = [] for i,one in enumerate(values): if one.startswith(value): dellist.append(i) while dellist: values.pop(dellist.pop()) self.var[var] = values # write stored makefile lines to file, using vars that may have been updated # do not write var if not in dict, since has been deleted # wrap var values into multiple lines if needed # file = 1 if this is Makefile.auto, change 1st line to use "auto" def write(self,file,flag=0): fp = open(file,'w') for i,line in enumerate(self.lines): if not line.isdigit(): if flag and i == 0: line = "# auto = makefile auto-generated by Make.py" print >>fp,line else: index = int(line) name = self.varinfo[index][0] txt = self.varinfo[index][1] if name not in self.var: continue values = self.var[name] print >>fp,"%s%s" % (txt,' '.join(values)) # ---------------------------------------------------------------- # main program # ---------------------------------------------------------------- # parse command-line args # switches dict: key = switch letter, value = list of args # switch_order = list of switches in order # will possibly be merged with redo file args below cmd_switches,cmd_switch_order = parse_args(sys.argv[1:]) if "v" in cmd_switches: print "Command-line parsing:" for switch in cmd_switch_order: print " %s: %s" % (switch,' '.join(cmd_switches[switch])) # check for redo switch, process redo file # redolist = list of commands to execute redoflag = 0 redolist = [] if 'r' in cmd_switches and 'h' not in cmd_switches: redoflag = 1 redo = Redo(cmd_switches['r']) redo.check() redo.setup() redolist = redo.commands redoindex = 0 del redo if not redolist: error("No commands to execute from redo file") # loop over Make.py commands # if no redo switch, loop once for command-line command # if redo, loop over one or more commands from redo file while 1: # if redo: # parse next command from redo file # use command-line switches to add/replace file command switches # do not add -r, since already processed # and don't want -r swtich to appear in Make.py.last file # if -a in both: concatenate, de-dup, # specified exe/machine action replaces file exe/machine action # print resulting new command # else just use command-line switches if redoflag: if redoindex == len(redolist): break args = redolist[redoindex].split() switches,switch_order = parse_args(args) redoindex += 1 for switch in cmd_switches: if switch == 'r': continue if switch == 'a' and switch in switches: tmp = Actions(cmd_switches[switch] + switches[switch]) tmp.dedup() switches[switch] = tmp.inlist continue if switch not in switches: switch_order.append(switch) switches[switch] = cmd_switches[switch] argstr = switch2str(switches,switch_order) print "Redo command: Make.py",argstr else: switches = cmd_switches switch_order = cmd_switch_order # initialize all class variables to None for one in switchclasses: exec("%s = None" % one) for one in libclasses: exec("%s = None" % one) for one in buildclasses: exec("%s = None" % one) for one in makeclasses: exec("%s = None" % one) # classes = dictionary of created classes # key = switch, value = class instance classes = {} for switch in switches: if len(switch) == 1 and switch in abbrevs: i = abbrevs.index(switch) capitalized = switchclasses[i][0].upper() + switchclasses[i][1:] txt = '%s = classes["%s"] = %s(switches["%s"])' % \ (switchclasses[i],switch,capitalized,switch) exec(txt) elif switch in libclasses: i = libclasses.index(switch) txt = '%s = classes["%s"] = %s(switches["%s"])' % \ (libclasses[i],switch,libclasses[i].upper(),switch) exec(txt) elif switch in buildclasses: i = buildclasses.index(switch) capitalized = buildclasses[i][0].upper() + buildclasses[i][1:] txt = '%s = classes["%s"] = %s(switches["%s"])' % \ (buildclasses[i],switch,capitalized,switch) exec(txt) elif switch in makeclasses: i = makeclasses.index(switch) capitalized = makeclasses[i][0].upper() + makeclasses[i][1:] txt = '%s = classes["%s"] = %s(switches["%s"])' % \ (makeclasses[i],switch,capitalized,switch) exec(txt) else: error("Unknown command-line switch -%s" % switch) # print help messages and exit if help or (actions and "-h" in actions.inlist) or not switches: if not help: help = Help(None) print help.help() for switch in switch_order: if switch == "h": continue print classes[switch].help()[1:] sys.exit() # create needed default classes if not specified with switch # dir and packages plus lib and build classes so defaults are set if not dir: dir = Dir(None) if not packages: packages = Packages(None) for one in libclasses: txt = "if not %s: %s = %s(None)" % (one,one,one.upper()) exec(txt) for one in buildclasses: capitalized = one[0].upper() + one[1:] txt = "if not %s: %s = %s(None)" % (one,one,capitalized) exec(txt) # error check on args for all classes for switch in classes: classes[switch].check() # prep for action # actions.setup() detects if last action = machine # if yes, induce addition of "-m" and "-o" switches dir.setup() packages.setup() if actions: machine = actions.setup() if machine: switches['a'][-1] = "exe" if 'm' not in switches: switches['m'] = [machine] switch_order.insert(-1,'m') makefile = classes['m'] = Makefile(switches['m']) makefile.check() if 'o' not in switches: switches['o'] = [machine] switch_order.insert(-1,'o') output = classes['o'] = Makefile(switches['o']) output.check() # perform actions packages.install() if actions: for action in actions.alist: print "Action %s ..." % action if action.startswith("lib-"): actions.lib(action[4:]) elif action == "file": actions.file("file") elif action == "clean": actions.clean() elif action == "exe": actions.exe() packages.uninstall() # create output file if requested and exe action performed if output and actions and "exe" in actions.alist: txt = "cp %s/lmp_auto %s/lmp_%s" % (dir.src,dir.cwd,output.machine) commands.getoutput(txt) print "Created lmp_%s in %s" % (output.machine,dir.cwd) # write current Make.py command to src/Make.py.last fp = open("%s/Make.py.last" % dir.src,'w') print >>fp,"# last invoked Make.py command" print >>fp,switch2str(switches,switch_order) fp.close() # if not redoflag, done if not redoflag: break diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index a03b394ee..0dd6c85d5 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -1,2658 +1,2667 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "stdlib.h" #include "string.h" #include "stdio.h" #include "fix_shake.h" #include "atom.h" #include "atom_vec.h" #include "molecule.h" #include "update.h" #include "respa.h" #include "modify.h" #include "domain.h" #include "force.h" #include "bond.h" #include "angle.h" #include "comm.h" #include "group.h" #include "fix_respa.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; // allocate space for static class variable FixShake *FixShake::fsptr; #define BIG 1.0e20 #define MASSDELTA 0.1 /* ---------------------------------------------------------------------- */ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); virial_flag = 1; create_attribute = 1; // error check molecular = atom->molecular; if (molecular == 0) error->all(FLERR,"Cannot use fix shake with non-molecular system"); // perform initial allocation of atom-based arrays // register with Atom class shake_flag = NULL; shake_atom = NULL; shake_type = NULL; xshake = NULL; grow_arrays(atom->nmax); atom->add_callback(0); // set comm size needed by this fix comm_forward = 3; // parse SHAKE args if (narg < 8) error->all(FLERR,"Illegal fix shake command"); tolerance = force->numeric(FLERR,arg[3]); max_iter = force->inumeric(FLERR,arg[4]); output_every = force->inumeric(FLERR,arg[5]); // parse SHAKE args for bond and angle types // will be used by find_clusters // store args for "b" "a" "t" as flags in (1:n) list for fast access // store args for "m" in list of length nmass for looping over // for "m" verify that atom masses have been set bond_flag = new int[atom->nbondtypes+1]; for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0; angle_flag = new int[atom->nangletypes+1]; for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0; type_flag = new int[atom->ntypes+1]; for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0; mass_list = new double[atom->ntypes]; nmass = 0; char mode = '\0'; int next = 6; while (next < narg) { if (strcmp(arg[next],"b") == 0) mode = 'b'; else if (strcmp(arg[next],"a") == 0) mode = 'a'; else if (strcmp(arg[next],"t") == 0) mode = 't'; else if (strcmp(arg[next],"m") == 0) { mode = 'm'; atom->check_mass(); // break if keyword that is not b,a,t,m } else if (isalpha(arg[next][0])) break; // read numeric args of b,a,t,m else if (mode == 'b') { int i = force->inumeric(FLERR,arg[next]); if (i < 1 || i > atom->nbondtypes) error->all(FLERR,"Invalid bond type index for fix shake"); bond_flag[i] = 1; } else if (mode == 'a') { int i = force->inumeric(FLERR,arg[next]); if (i < 1 || i > atom->nangletypes) error->all(FLERR,"Invalid angle type index for fix shake"); angle_flag[i] = 1; } else if (mode == 't') { int i = force->inumeric(FLERR,arg[next]); if (i < 1 || i > atom->ntypes) error->all(FLERR,"Invalid atom type index for fix shake"); type_flag[i] = 1; } else if (mode == 'm') { double massone = force->numeric(FLERR,arg[next]); if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake"); if (nmass == atom->ntypes) error->all(FLERR,"Too many masses for fix shake"); mass_list[nmass++] = massone; } else error->all(FLERR,"Illegal fix shake command"); next++; } // parse optional args onemols = NULL; int iarg = next; while (iarg < narg) { if (strcmp(arg[next],"mol") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command"); int imol = atom->find_molecule(arg[iarg+1]); if (imol == -1) error->all(FLERR,"Molecule template ID for fix shake does not exist"); if (atom->molecules[imol]->nset > 1 && comm->me == 0) error->warning(FLERR,"Molecule template for " "fix shake has multiple molecules"); onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; iarg += 2; } else error->all(FLERR,"Illegal fix shake command"); } // error check for Molecule template if (onemols) { for (int i = 0; i < nmol; i++) if (onemols[i]->shakeflag == 0) error->all(FLERR,"Fix shake molecule template must have shake info"); } // allocate bond and angle distance arrays, indexed from 1 to n bond_distance = new double[atom->nbondtypes+1]; angle_distance = new double[atom->nangletypes+1]; // allocate statistics arrays if (output_every) { int nb = atom->nbondtypes + 1; b_count = new int[nb]; b_count_all = new int[nb]; b_ave = new double[nb]; b_ave_all = new double[nb]; b_max = new double[nb]; b_max_all = new double[nb]; b_min = new double[nb]; b_min_all = new double[nb]; int na = atom->nangletypes + 1; a_count = new int[na]; a_count_all = new int[na]; a_ave = new double[na]; a_ave_all = new double[na]; a_max = new double[na]; a_max_all = new double[na]; a_min = new double[na]; a_min_all = new double[na]; } + // SHAKE vs RATTLE + + rattle = 0; + // identify all SHAKE clusters find_clusters(); // initialize list of SHAKE clusters to constrain maxlist = 0; list = NULL; } /* ---------------------------------------------------------------------- */ FixShake::~FixShake() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); // set bond_type and angle_type back to positive for SHAKE clusters // must set for all SHAKE bonds and angles stored by each atom int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { if (shake_flag[i] == 0) continue; else if (shake_flag[i] == 1) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1); angletype_findset(i,shake_atom[i][1],shake_atom[i][2],1); } else if (shake_flag[i] == 2) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1); } else if (shake_flag[i] == 3) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1); } else if (shake_flag[i] == 4) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],1); } } // delete locally stored arrays memory->destroy(shake_flag); memory->destroy(shake_atom); memory->destroy(shake_type); memory->destroy(xshake); delete [] bond_flag; delete [] angle_flag; delete [] type_flag; delete [] mass_list; delete [] bond_distance; delete [] angle_distance; if (output_every) { delete [] b_count; delete [] b_count_all; delete [] b_ave; delete [] b_ave_all; delete [] b_max; delete [] b_max_all; delete [] b_min; delete [] b_min_all; delete [] a_count; delete [] a_count_all; delete [] a_ave; delete [] a_ave_all; delete [] a_max; delete [] a_max_all; delete [] a_min; delete [] a_min_all; } memory->destroy(list); } /* ---------------------------------------------------------------------- */ int FixShake::setmask() { int mask = 0; mask |= PRE_NEIGHBOR; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; return mask; } /* ---------------------------------------------------------------------- set bond and angle distances this init must happen after force->bond and force->angle inits ------------------------------------------------------------------------- */ void FixShake::init() { int i,m,flag,flag_all,type1,type2,bond1_type,bond2_type; double rsq,angle; // error if more than one shake fix int count = 0; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"shake") == 0) count++; if (count > 1) error->all(FLERR,"More than one fix shake"); // cannot use with minimization since SHAKE turns off bonds // that should contribute to potential energy if (update->whichflag == 2) error->all(FLERR,"Fix shake cannot be used with minimization"); // error if npt,nph fix comes before shake fix for (i = 0; i < modify->nfix; i++) { if (strcmp(modify->fix[i]->style,"npt") == 0) break; if (strcmp(modify->fix[i]->style,"nph") == 0) break; } if (i < modify->nfix) { for (int j = i; j < modify->nfix; j++) if (strcmp(modify->fix[j]->style,"shake") == 0) error->all(FLERR,"Shake fix must come before NPT/NPH fix"); } // if rRESPA, find associated fix that must exist // could have changed locations in fix list since created // set ptrs to rRESPA variables if (strstr(update->integrate_style,"respa")) { for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"RESPA") == 0) ifix_respa = i; nlevels_respa = ((Respa *) update->integrate)->nlevels; loop_respa = ((Respa *) update->integrate)->loop; step_respa = ((Respa *) update->integrate)->step; } // set equilibrium bond distances if (force->bond == NULL) error->all(FLERR,"Bond potential must be defined for SHAKE"); for (i = 1; i <= atom->nbondtypes; i++) bond_distance[i] = force->bond->equilibrium_distance(i); // set equilibrium angle distances int nlocal = atom->nlocal; for (i = 1; i <= atom->nangletypes; i++) { if (angle_flag[i] == 0) continue; if (force->angle == NULL) error->all(FLERR,"Angle potential must be defined for SHAKE"); // scan all atoms for a SHAKE angle cluster // extract bond types for the 2 bonds in the cluster // bond types must be same in all clusters of this angle type, // else set error flag flag = 0; bond1_type = bond2_type = 0; for (m = 0; m < nlocal; m++) { if (shake_flag[m] != 1) continue; if (shake_type[m][2] != i) continue; type1 = MIN(shake_type[m][0],shake_type[m][1]); type2 = MAX(shake_type[m][0],shake_type[m][1]); if (bond1_type > 0) { if (type1 != bond1_type || type2 != bond2_type) { flag = 1; break; } } bond1_type = type1; bond2_type = type2; } // error check for any bond types that are not the same MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world); if (flag_all) error->all(FLERR,"Shake angles have different bond types"); // insure all procs have bond types MPI_Allreduce(&bond1_type,&flag_all,1,MPI_INT,MPI_MAX,world); bond1_type = flag_all; MPI_Allreduce(&bond2_type,&flag_all,1,MPI_INT,MPI_MAX,world); bond2_type = flag_all; // if bond types are 0, no SHAKE angles of this type exist // just skip this angle if (bond1_type == 0) { angle_distance[i] = 0.0; continue; } // compute the angle distance as a function of 2 bond distances angle = force->angle->equilibrium_angle(i); rsq = 2.0*bond_distance[bond1_type]*bond_distance[bond2_type] * (1.0-cos(angle)); angle_distance[i] = sqrt(rsq); } } /* ---------------------------------------------------------------------- SHAKE as pre-integrator constraint ------------------------------------------------------------------------- */ void FixShake::setup(int vflag) { pre_neighbor(); if (output_every) stats(); // setup SHAKE output bigint ntimestep = update->ntimestep; if (output_every) { next_output = ntimestep + output_every; if (ntimestep % output_every != 0) next_output = (ntimestep/output_every)*output_every + output_every; } else next_output = -1; // half timestep constraint on pre-step, full timestep thereafter if (strstr(update->integrate_style,"verlet")) { - dtv = update->dt; - dtfsq = 0.5 * update->dt * update->dt * force->ftm2v; - post_force(vflag); - dtfsq = update->dt * update->dt * force->ftm2v; + dtv = update->dt; + dtfsq = 0.5 * update->dt * update->dt * force->ftm2v; + FixShake::post_force(vflag); + if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v; + } else { - dtv = step_respa[0]; + dtv = step_respa[0]; dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v; dtf_inner = dtf_innerhalf; // apply correction to all rRESPA levels for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(vflag,ilevel,loop_respa[ilevel]-1); + FixShake::post_force_respa(vflag,ilevel,loop_respa[ilevel]-1); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } - - dtf_inner = step_respa[0] * force->ftm2v; + if (!rattle) dtf_inner = step_respa[0] * force->ftm2v; } } /* ---------------------------------------------------------------------- build list of SHAKE clusters to constrain if one or more atoms in cluster are on this proc, this proc lists the cluster exactly once ------------------------------------------------------------------------- */ void FixShake::pre_neighbor() { int atom1,atom2,atom3,atom4; // local copies of atom quantities // used by SHAKE until next re-neighboring x = atom->x; v = atom->v; f = atom->f; mass = atom->mass; rmass = atom->rmass; type = atom->type; nlocal = atom->nlocal; // extend size of SHAKE list if necessary if (nlocal > maxlist) { maxlist = nlocal; memory->destroy(list); memory->create(list,maxlist,"shake:list"); } // build list of SHAKE clusters I compute nlist = 0; for (int i = 0; i < nlocal; i++) if (shake_flag[i]) { if (shake_flag[i] == 2) { atom1 = atom->map(shake_atom[i][0]); atom2 = atom->map(shake_atom[i][1]); if (atom1 == -1 || atom2 == -1) { char str[128]; sprintf(str,"Shake atoms " TAGINT_FORMAT " " TAGINT_FORMAT " missing on proc %d at step " BIGINT_FORMAT, shake_atom[i][0],shake_atom[i][1],me,update->ntimestep); error->one(FLERR,str); } if (i <= atom1 && i <= atom2) list[nlist++] = i; } else if (shake_flag[i] % 2 == 1) { atom1 = atom->map(shake_atom[i][0]); atom2 = atom->map(shake_atom[i][1]); atom3 = atom->map(shake_atom[i][2]); if (atom1 == -1 || atom2 == -1 || atom3 == -1) { char str[128]; sprintf(str,"Shake atoms " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " missing on proc %d at step " BIGINT_FORMAT, shake_atom[i][0],shake_atom[i][1],shake_atom[i][2], me,update->ntimestep); error->one(FLERR,str); } if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i; } else { atom1 = atom->map(shake_atom[i][0]); atom2 = atom->map(shake_atom[i][1]); atom3 = atom->map(shake_atom[i][2]); atom4 = atom->map(shake_atom[i][3]); if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) { char str[128]; sprintf(str,"Shake atoms " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " missing on proc %d at step " BIGINT_FORMAT, shake_atom[i][0],shake_atom[i][1], shake_atom[i][2],shake_atom[i][3], me,update->ntimestep); error->one(FLERR,str); } if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4) list[nlist++] = i; } } } /* ---------------------------------------------------------------------- compute the force adjustment for SHAKE constraint ------------------------------------------------------------------------- */ void FixShake::post_force(int vflag) { if (update->ntimestep == next_output) stats(); // xshake = unconstrained move with current v,f // communicate results if necessary unconstrained_update(); if (nprocs > 1) comm->forward_comm_fix(this); // virial setup if (vflag) v_setup(vflag); else evflag = 0; // loop over clusters to add constraint forces int m; for (int i = 0; i < nlist; i++) { m = list[i]; if (shake_flag[m] == 2) shake(m); else if (shake_flag[m] == 3) shake3(m); else if (shake_flag[m] == 4) shake4(m); else shake3angle(m); } } /* ---------------------------------------------------------------------- enforce SHAKE constraints from rRESPA xshake prediction portion is different than Verlet ------------------------------------------------------------------------- */ void FixShake::post_force_respa(int vflag, int ilevel, int iloop) { // call stats only on outermost level if (ilevel == nlevels_respa-1 && update->ntimestep == next_output) stats(); // might be OK to skip enforcing SHAKE constraings // on last iteration of inner levels if pressure not requested // however, leads to slightly different trajectories //if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1 && !vflag) // return; // xshake = unconstrained move with current v,f as function of level // communicate results if necessary unconstrained_update_respa(ilevel); if (nprocs > 1) comm->forward_comm_fix(this); // virial setup only needed on last iteration of innermost level // and if pressure is requested // virial accumulation happens via evflag at last iteration of each level if (ilevel == 0 && iloop == loop_respa[ilevel]-1 && vflag) v_setup(vflag); if (iloop == loop_respa[ilevel]-1) evflag = 1; else evflag = 0; // loop over clusters to add constraint forces int m; for (int i = 0; i < nlist; i++) { m = list[i]; if (shake_flag[m] == 2) shake(m); else if (shake_flag[m] == 3) shake3(m); else if (shake_flag[m] == 4) shake4(m); else shake3angle(m); } } /* ---------------------------------------------------------------------- count # of degrees-of-freedom removed by SHAKE for atoms in igroup ------------------------------------------------------------------------- */ int FixShake::dof(int igroup) { int groupbit = group->bitmask[igroup]; int *mask = atom->mask; tagint *tag = atom->tag; int nlocal = atom->nlocal; // count dof in a cluster if and only if // the central atom is in group and atom i is the central atom int n = 0; for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; if (shake_flag[i] == 0) continue; if (shake_atom[i][0] != tag[i]) continue; if (shake_flag[i] == 1) n += 3; else if (shake_flag[i] == 2) n += 1; else if (shake_flag[i] == 3) n += 2; else if (shake_flag[i] == 4) n += 3; } int nall; MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world); return nall; } /* ---------------------------------------------------------------------- identify whether each atom is in a SHAKE cluster only include atoms in fix group and those bonds/angles specified in input test whether all clusters are valid set shake_flag, shake_atom, shake_type values set bond,angle types negative so will be ignored in neighbor lists ------------------------------------------------------------------------- */ void FixShake::find_clusters() { int i,j,m,n,imol,iatom; int flag,flag_all,nbuf,size; tagint tagprev; double massone; tagint *buf; - - if (me == 0 && screen) fprintf(screen,"Finding SHAKE clusters ...\n"); + + if (me == 0 && screen) { + if (!rattle) fprintf(screen,"Finding SHAKE clusters ...\n"); + else fprintf(screen,"Finding RATTLE clusters ...\n"); + } atommols = atom->avec->onemols; tagint *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; double *mass = atom->mass; double *rmass = atom->rmass; int **nspecial = atom->nspecial; tagint **special = atom->special; int *molindex = atom->molindex; int *molatom = atom->molatom; int nlocal = atom->nlocal; int angles_allow = atom->avec->angles_allow; // setup ring of procs int next = me + 1; int prev = me -1; if (next == nprocs) next = 0; if (prev < 0) prev = nprocs - 1; // ----------------------------------------------------- // allocate arrays for self (1d) and bond partners (2d) // max = max # of bond partners for owned atoms = 2nd dim of partner arrays // npartner[i] = # of bonds attached to atom i // nshake[i] = # of SHAKE bonds attached to atom i // partner_tag[i][] = global IDs of each partner // partner_mask[i][] = mask of each partner // partner_type[i][] = type of each partner // partner_massflag[i][] = 1 if partner meets mass criterion, 0 if not // partner_bondtype[i][] = type of bond attached to each partner // partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not // partner_nshake[i][] = nshake value for each partner // ----------------------------------------------------- int max = 0; if (molecular == 1) { for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]); } else { for (i = 0; i < nlocal; i++) { imol = molindex[i]; if (imol < 0) continue; iatom = molatom[i]; max = MAX(max,atommols[imol]->nspecial[iatom][0]); } } int *npartner; memory->create(npartner,nlocal,"shake:npartner"); memory->create(nshake,nlocal,"shake:nshake"); tagint **partner_tag; int **partner_mask,**partner_type,**partner_massflag; int **partner_bondtype,**partner_shake,**partner_nshake; memory->create(partner_tag,nlocal,max,"shake:partner_tag"); memory->create(partner_mask,nlocal,max,"shake:partner_mask"); memory->create(partner_type,nlocal,max,"shake:partner_type"); memory->create(partner_massflag,nlocal,max,"shake:partner_massflag"); memory->create(partner_bondtype,nlocal,max,"shake:partner_bondtype"); memory->create(partner_shake,nlocal,max,"shake:partner_shake"); memory->create(partner_nshake,nlocal,max,"shake:partner_nshake"); // ----------------------------------------------------- // set npartner and partner_tag from special arrays // ----------------------------------------------------- if (molecular == 1) { for (i = 0; i < nlocal; i++) { npartner[i] = nspecial[i][0]; for (j = 0; j < npartner[i]; j++) partner_tag[i][j] = special[i][j]; } } else { for (i = 0; i < nlocal; i++) { imol = molindex[i]; if (imol < 0) continue; iatom = molatom[i]; tagprev = tag[i] - iatom - 1; npartner[i] = atommols[imol]->nspecial[iatom][0]; for (j = 0; j < npartner[i]; j++) partner_tag[i][j] = atommols[imol]->special[iatom][j] + tagprev;; } } // ----------------------------------------------------- // set partner_mask, partner_type, partner_massflag, partner_bondtype // for bonded partners // requires communication for off-proc partners // ----------------------------------------------------- // fill in mask, type, massflag, bondtype if own bond partner // info to store in buf for each off-proc bond = nper = 6 // 2 atoms IDs in bond, space for mask, type, massflag, bondtype // nbufmax = largest buffer needed to hold info from any proc int nper = 6; nbuf = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < npartner[i]; j++) { partner_mask[i][j] = 0; partner_type[i][j] = 0; partner_massflag[i][j] = 0; partner_bondtype[i][j] = 0; m = atom->map(partner_tag[i][j]); if (m >= 0 && m < nlocal) { partner_mask[i][j] = mask[m]; partner_type[i][j] = type[m]; if (nmass) { if (rmass) massone = rmass[m]; else massone = mass[type[m]]; partner_massflag[i][j] = masscheck(massone); } n = bondtype_findset(i,tag[i],partner_tag[i][j],0); if (n) partner_bondtype[i][j] = n; else { n = bondtype_findset(m,tag[i],partner_tag[i][j],0); if (n) partner_bondtype[i][j] = n; } } else nbuf += nper; } } memory->create(buf,nbuf,"shake:buf"); // fill buffer with info size = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < npartner[i]; j++) { m = atom->map(partner_tag[i][j]); if (m < 0 || m >= nlocal) { buf[size] = tag[i]; buf[size+1] = partner_tag[i][j]; buf[size+2] = 0; buf[size+3] = 0; buf[size+4] = 0; n = bondtype_findset(i,tag[i],partner_tag[i][j],0); if (n) buf[size+5] = n; else buf[size+5] = 0; size += nper; } } } // cycle buffer around ring of procs back to self fsptr = this; comm->ring(size,sizeof(tagint),buf,1,ring_bonds,buf); // store partner info returned to me m = 0; while (m < size) { i = atom->map(buf[m]); for (j = 0; j < npartner[i]; j++) if (buf[m+1] == partner_tag[i][j]) break; partner_mask[i][j] = buf[m+2]; partner_type[i][j] = buf[m+3]; partner_massflag[i][j] = buf[m+4]; partner_bondtype[i][j] = buf[m+5]; m += nper; } memory->destroy(buf); // error check for unfilled partner info // if partner_type not set, is an error // partner_bondtype may not be set if special list is not consistent // with bondatom (e.g. due to delete_bonds command) // this is OK if one or both atoms are not in fix group, since // bond won't be SHAKEn anyway // else it's an error flag = 0; for (i = 0; i < nlocal; i++) for (j = 0; j < npartner[i]; j++) { if (partner_type[i][j] == 0) flag = 1; if (!(mask[i] & groupbit)) continue; if (!(partner_mask[i][j] & groupbit)) continue; if (partner_bondtype[i][j] == 0) flag = 1; } MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Did not find fix shake partner info"); // ----------------------------------------------------- // identify SHAKEable bonds // set nshake[i] = # of SHAKE bonds attached to atom i // set partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not // both atoms must be in group, bondtype must be > 0 // check if bondtype is in input bond_flag // check if type of either atom is in input type_flag // check if mass of either atom is in input mass_list // ----------------------------------------------------- int np; for (i = 0; i < nlocal; i++) { nshake[i] = 0; np = npartner[i]; for (j = 0; j < np; j++) { partner_shake[i][j] = 0; if (!(mask[i] & groupbit)) continue; if (!(partner_mask[i][j] & groupbit)) continue; if (partner_bondtype[i][j] <= 0) continue; if (bond_flag[partner_bondtype[i][j]]) { partner_shake[i][j] = 1; nshake[i]++; continue; } if (type_flag[type[i]] || type_flag[partner_type[i][j]]) { partner_shake[i][j] = 1; nshake[i]++; continue; } if (nmass) { if (partner_massflag[i][j]) { partner_shake[i][j] = 1; nshake[i]++; continue; } else { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if (masscheck(massone)) { partner_shake[i][j] = 1; nshake[i]++; continue; } } } } } // ----------------------------------------------------- // set partner_nshake for bonded partners // requires communication for off-proc partners // ----------------------------------------------------- // fill in partner_nshake if own bond partner // info to store in buf for each off-proc bond = // 2 atoms IDs in bond, space for nshake value // nbufmax = largest buffer needed to hold info from any proc nbuf = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < npartner[i]; j++) { m = atom->map(partner_tag[i][j]); if (m >= 0 && m < nlocal) partner_nshake[i][j] = nshake[m]; else nbuf += 3; } } memory->create(buf,nbuf,"shake:buf"); // fill buffer with info size = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < npartner[i]; j++) { m = atom->map(partner_tag[i][j]); if (m < 0 || m >= nlocal) { buf[size] = tag[i]; buf[size+1] = partner_tag[i][j]; size += 3; } } } // cycle buffer around ring of procs back to self fsptr = this; comm->ring(size,sizeof(tagint),buf,2,ring_nshake,buf); // store partner info returned to me m = 0; while (m < size) { i = atom->map(buf[m]); for (j = 0; j < npartner[i]; j++) if (buf[m+1] == partner_tag[i][j]) break; partner_nshake[i][j] = buf[m+2]; m += 3; } memory->destroy(buf); // ----------------------------------------------------- // error checks // no atom with nshake > 3 // no connected atoms which both have nshake > 1 // ----------------------------------------------------- flag = 0; for (i = 0; i < nlocal; i++) if (nshake[i] > 3) flag = 1; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Shake cluster of more than 4 atoms"); flag = 0; for (i = 0; i < nlocal; i++) { if (nshake[i] <= 1) continue; for (j = 0; j < npartner[i]; j++) if (partner_shake[i][j] && partner_nshake[i][j] > 1) flag = 1; } MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Shake clusters are connected"); // ----------------------------------------------------- // set SHAKE arrays that are stored with atoms & add angle constraints // zero shake arrays for all owned atoms // if I am central atom set shake_flag & shake_atom & shake_type // for 2-atom clusters, I am central atom if my atom ID < partner ID // for 3-atom clusters, test for angle constraint // angle will be stored by this atom if it exists // if angle type matches angle_flag, then it is angle-constrained // shake_flag[] = 0 if atom not in SHAKE cluster // 2,3,4 = size of bond-only cluster // 1 = 3-atom angle cluster // shake_atom[][] = global IDs of 2,3,4 atoms in cluster // central atom is 1st // for 2-atom cluster, lowest ID is 1st // shake_type[][] = bondtype of each bond in cluster // for 3-atom angle cluster, 3rd value is angletype // ----------------------------------------------------- for (i = 0; i < nlocal; i++) { shake_flag[i] = 0; shake_atom[i][0] = 0; shake_atom[i][1] = 0; shake_atom[i][2] = 0; shake_atom[i][3] = 0; shake_type[i][0] = 0; shake_type[i][1] = 0; shake_type[i][2] = 0; if (nshake[i] == 1) { for (j = 0; j < npartner[i]; j++) if (partner_shake[i][j]) break; if (partner_nshake[i][j] == 1 && tag[i] < partner_tag[i][j]) { shake_flag[i] = 2; shake_atom[i][0] = tag[i]; shake_atom[i][1] = partner_tag[i][j]; shake_type[i][0] = partner_bondtype[i][j]; } } if (nshake[i] > 1) { shake_flag[i] = 1; shake_atom[i][0] = tag[i]; for (j = 0; j < npartner[i]; j++) if (partner_shake[i][j]) { m = shake_flag[i]; shake_atom[i][m] = partner_tag[i][j]; shake_type[i][m-1] = partner_bondtype[i][j]; shake_flag[i]++; } } if (nshake[i] == 2 && angles_allow) { n = angletype_findset(i,shake_atom[i][1],shake_atom[i][2],0); if (n <= 0) continue; if (angle_flag[n]) { shake_flag[i] = 1; shake_type[i][2] = n; } } } // ----------------------------------------------------- // set shake_flag,shake_atom,shake_type for non-central atoms // requires communication for off-proc atoms // ----------------------------------------------------- // fill in shake arrays for each bond partner I own // info to store in buf for each off-proc bond = // all values from shake_flag, shake_atom, shake_type // nbufmax = largest buffer needed to hold info from any proc nbuf = 0; for (i = 0; i < nlocal; i++) { if (shake_flag[i] == 0) continue; for (j = 0; j < npartner[i]; j++) { if (partner_shake[i][j] == 0) continue; m = atom->map(partner_tag[i][j]); if (m >= 0 && m < nlocal) { shake_flag[m] = shake_flag[i]; shake_atom[m][0] = shake_atom[i][0]; shake_atom[m][1] = shake_atom[i][1]; shake_atom[m][2] = shake_atom[i][2]; shake_atom[m][3] = shake_atom[i][3]; shake_type[m][0] = shake_type[i][0]; shake_type[m][1] = shake_type[i][1]; shake_type[m][2] = shake_type[i][2]; } else nbuf += 9; } } memory->create(buf,nbuf,"shake:buf"); // fill buffer with info size = 0; for (i = 0; i < nlocal; i++) { if (shake_flag[i] == 0) continue; for (j = 0; j < npartner[i]; j++) { if (partner_shake[i][j] == 0) continue; m = atom->map(partner_tag[i][j]); if (m < 0 || m >= nlocal) { buf[size] = partner_tag[i][j]; buf[size+1] = shake_flag[i]; buf[size+2] = shake_atom[i][0]; buf[size+3] = shake_atom[i][1]; buf[size+4] = shake_atom[i][2]; buf[size+5] = shake_atom[i][3]; buf[size+6] = shake_type[i][0]; buf[size+7] = shake_type[i][1]; buf[size+8] = shake_type[i][2]; size += 9; } } } // cycle buffer around ring of procs back to self fsptr = this; comm->ring(size,sizeof(tagint),buf,3,ring_shake,NULL); memory->destroy(buf); // ----------------------------------------------------- // free local memory // ----------------------------------------------------- memory->destroy(npartner); memory->destroy(nshake); memory->destroy(partner_tag); memory->destroy(partner_mask); memory->destroy(partner_type); memory->destroy(partner_massflag); memory->destroy(partner_bondtype); memory->destroy(partner_shake); memory->destroy(partner_nshake); // ----------------------------------------------------- // set bond_type and angle_type negative for SHAKE clusters // must set for all SHAKE bonds and angles stored by each atom // ----------------------------------------------------- for (i = 0; i < nlocal; i++) { if (shake_flag[i] == 0) continue; else if (shake_flag[i] == 1) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1); angletype_findset(i,shake_atom[i][1],shake_atom[i][2],-1); } else if (shake_flag[i] == 2) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1); } else if (shake_flag[i] == 3) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1); } else if (shake_flag[i] == 4) { bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1); bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],-1); } } // ----------------------------------------------------- // print info on SHAKE clusters // ----------------------------------------------------- int count1,count2,count3,count4; count1 = count2 = count3 = count4 = 0; for (i = 0; i < nlocal; i++) { if (shake_flag[i] == 1) count1++; else if (shake_flag[i] == 2) count2++; else if (shake_flag[i] == 3) count3++; else if (shake_flag[i] == 4) count4++; } int tmp; tmp = count1; MPI_Allreduce(&tmp,&count1,1,MPI_INT,MPI_SUM,world); tmp = count2; MPI_Allreduce(&tmp,&count2,1,MPI_INT,MPI_SUM,world); tmp = count3; MPI_Allreduce(&tmp,&count3,1,MPI_INT,MPI_SUM,world); tmp = count4; MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world); if (me == 0) { if (screen) { fprintf(screen," %d = # of size 2 clusters\n",count2/2); fprintf(screen," %d = # of size 3 clusters\n",count3/3); fprintf(screen," %d = # of size 4 clusters\n",count4/4); fprintf(screen," %d = # of frozen angles\n",count1/3); } if (logfile) { fprintf(logfile," %d = # of size 2 clusters\n",count2/2); fprintf(logfile," %d = # of size 3 clusters\n",count3/3); fprintf(logfile," %d = # of size 4 clusters\n",count4/4); fprintf(logfile," %d = # of frozen angles\n",count1/3); } } } /* ---------------------------------------------------------------------- when receive buffer, scan bond partner IDs for atoms I own if I own partner: fill in mask and type and massflag search for bond with 1st atom and fill in bondtype ------------------------------------------------------------------------- */ void FixShake::ring_bonds(int ndatum, char *cbuf) { Atom *atom = fsptr->atom; double *rmass = atom->rmass; double *mass = atom->mass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; int nmass = fsptr->nmass; tagint *buf = (tagint *) cbuf; int m,n; double massone; for (int i = 0; i < ndatum; i += 6) { m = atom->map(buf[i+1]); if (m >= 0 && m < nlocal) { buf[i+2] = mask[m]; buf[i+3] = type[m]; if (nmass) { if (rmass) massone = rmass[m]; else massone = mass[type[m]]; buf[i+4] = fsptr->masscheck(massone); } if (buf[i+5] == 0) { n = fsptr->bondtype_findset(m,buf[i],buf[i+1],0); if (n) buf[i+5] = n; } } } } /* ---------------------------------------------------------------------- when receive buffer, scan bond partner IDs for atoms I own if I own partner, fill in nshake value ------------------------------------------------------------------------- */ void FixShake::ring_nshake(int ndatum, char *cbuf) { Atom *atom = fsptr->atom; int nlocal = atom->nlocal; int *nshake = fsptr->nshake; tagint *buf = (tagint *) cbuf; int m; for (int i = 0; i < ndatum; i += 3) { m = atom->map(buf[i+1]); if (m >= 0 && m < nlocal) buf[i+2] = nshake[m]; } } /* ---------------------------------------------------------------------- when receive buffer, scan bond partner IDs for atoms I own if I own partner, fill in nshake value ------------------------------------------------------------------------- */ void FixShake::ring_shake(int ndatum, char *cbuf) { Atom *atom = fsptr->atom; int nlocal = atom->nlocal; int *shake_flag = fsptr->shake_flag; tagint **shake_atom = fsptr->shake_atom; int **shake_type = fsptr->shake_type; tagint *buf = (tagint *) cbuf; int m; for (int i = 0; i < ndatum; i += 9) { m = atom->map(buf[i]); if (m >= 0 && m < nlocal) { shake_flag[m] = buf[i+1]; shake_atom[m][0] = buf[i+2]; shake_atom[m][1] = buf[i+3]; shake_atom[m][2] = buf[i+4]; shake_atom[m][3] = buf[i+5]; shake_type[m][0] = buf[i+6]; shake_type[m][1] = buf[i+7]; shake_type[m][2] = buf[i+8]; } } } /* ---------------------------------------------------------------------- check if massone is within MASSDELTA of any mass in mass_list return 1 if yes, 0 if not ------------------------------------------------------------------------- */ int FixShake::masscheck(double massone) { for (int i = 0; i < nmass; i++) if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1; return 0; } /* ---------------------------------------------------------------------- update the unconstrained position of each atom only for SHAKE clusters, else set to 0.0 assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well ------------------------------------------------------------------------- */ void FixShake::unconstrained_update() { double dtfmsq; if (rmass) { for (int i = 0; i < nlocal; i++) { if (shake_flag[i]) { dtfmsq = dtfsq / rmass[i]; xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; } } else { for (int i = 0; i < nlocal; i++) { if (shake_flag[i]) { dtfmsq = dtfsq / mass[type[i]]; xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; } } } /* ---------------------------------------------------------------------- update the unconstrained position of each atom in a rRESPA step only for SHAKE clusters, else set to 0.0 assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well ------------------------------------------------------------------------- */ void FixShake::unconstrained_update_respa(int ilevel) { // xshake = atom coords after next x update in innermost loop // depends on rRESPA level // for levels > 0 this includes more than one velocity update // xshake = predicted position from call to this routine at level N = // x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0) // also set dtfsq = dt0*dtN so that shake,shake3,etc can use it double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level; dtfsq = dtf_inner * step_respa[ilevel]; double invmass,dtfmsq; int jlevel; if (rmass) { for (int i = 0; i < nlocal; i++) { if (shake_flag[i]) { invmass = 1.0 / rmass[i]; dtfmsq = dtfsq * invmass; xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; for (jlevel = 0; jlevel < ilevel; jlevel++) { dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; } } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; } } else { for (int i = 0; i < nlocal; i++) { if (shake_flag[i]) { invmass = 1.0 / mass[type[i]]; dtfmsq = dtfsq * invmass; xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; for (jlevel = 0; jlevel < ilevel; jlevel++) { dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; } } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; } } } /* ---------------------------------------------------------------------- */ void FixShake::shake(int m) { int nlist,list[2]; double v[6]; double invmass0,invmass1; // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); int i1 = atom->map(shake_atom[m][1]); double bond1 = bond_distance[shake_type[m][0]]; // r01 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i0][0] - x[i1][0]; r01[1] = x[i0][1] - x[i1][1]; r01[2] = x[i0][2] - x[i1][2]; domain->minimum_image(r01); // s01 = distance vec after unconstrained update, with PBC double s01[3]; s01[0] = xshake[i0][0] - xshake[i1][0]; s01[1] = xshake[i0][1] - xshake[i1][1]; s01[2] = xshake[i0][2] - xshake[i1][2]; domain->minimum_image(s01); // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; // a,b,c = coeffs in quadratic equation for lamda if (rmass) { invmass0 = 1.0/rmass[i0]; invmass1 = 1.0/rmass[i1]; } else { invmass0 = 1.0/mass[type[i0]]; invmass1 = 1.0/mass[type[i1]]; } double a = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; double b = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); double c = s01sq - bond1*bond1; // error check double determ = b*b - 4.0*a*c; if (determ < 0.0) { error->warning(FLERR,"Shake determinant < 0.0",0); determ = 0.0; } // exact quadratic solution for lamda double lamda,lamda1,lamda2; lamda1 = (-b+sqrt(determ)) / (2.0*a); lamda2 = (-b-sqrt(determ)) / (2.0*a); if (fabs(lamda1) <= fabs(lamda2)) lamda = lamda1; else lamda = lamda2; // update forces if atom is owned by this processor lamda /= dtfsq; if (i0 < nlocal) { f[i0][0] += lamda*r01[0]; f[i0][1] += lamda*r01[1]; f[i0][2] += lamda*r01[2]; } if (i1 < nlocal) { f[i1][0] -= lamda*r01[0]; f[i1][1] -= lamda*r01[1]; f[i1][2] -= lamda*r01[2]; } if (evflag) { nlist = 0; if (i0 < nlocal) list[nlist++] = i0; if (i1 < nlocal) list[nlist++] = i1; v[0] = lamda*r01[0]*r01[0]; v[1] = lamda*r01[1]*r01[1]; v[2] = lamda*r01[2]*r01[2]; v[3] = lamda*r01[0]*r01[1]; v[4] = lamda*r01[0]*r01[2]; v[5] = lamda*r01[1]*r01[2]; v_tally(nlist,list,2.0,v); } } /* ---------------------------------------------------------------------- */ void FixShake::shake3(int m) { int nlist,list[3]; double v[6]; double invmass0,invmass1,invmass2; // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); int i1 = atom->map(shake_atom[m][1]); int i2 = atom->map(shake_atom[m][2]); double bond1 = bond_distance[shake_type[m][0]]; double bond2 = bond_distance[shake_type[m][1]]; // r01,r02 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i0][0] - x[i1][0]; r01[1] = x[i0][1] - x[i1][1]; r01[2] = x[i0][2] - x[i1][2]; domain->minimum_image(r01); double r02[3]; r02[0] = x[i0][0] - x[i2][0]; r02[1] = x[i0][1] - x[i2][1]; r02[2] = x[i0][2] - x[i2][2]; domain->minimum_image(r02); // s01,s02 = distance vec after unconstrained update, with PBC double s01[3]; s01[0] = xshake[i0][0] - xshake[i1][0]; s01[1] = xshake[i0][1] - xshake[i1][1]; s01[2] = xshake[i0][2] - xshake[i1][2]; domain->minimum_image(s01); double s02[3]; s02[0] = xshake[i0][0] - xshake[i2][0]; s02[1] = xshake[i0][1] - xshake[i2][1]; s02[2] = xshake[i0][2] - xshake[i2][2]; domain->minimum_image(s02); // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; // matrix coeffs and rhs for lamda equations if (rmass) { invmass0 = 1.0/rmass[i0]; invmass1 = 1.0/rmass[i1]; invmass2 = 1.0/rmass[i2]; } else { invmass0 = 1.0/mass[type[i0]]; invmass1 = 1.0/mass[type[i1]]; invmass2 = 1.0/mass[type[i2]]; } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); double a12 = 2.0 * invmass0 * (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); double a21 = 2.0 * invmass0 * (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); double a22 = 2.0 * (invmass0+invmass2) * (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); // inverse of matrix double determ = a11*a22 - a12*a21; if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); double determinv = 1.0/determ; double a11inv = a22*determinv; double a12inv = -a12*determinv; double a21inv = -a21*determinv; double a22inv = a11*determinv; // quadratic correction coeffs double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; double quad1_0202 = invmass0*invmass0 * r02sq; double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; double quad2_0101 = invmass0*invmass0 * r01sq; double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; // iterate until converged double lamda01 = 0.0; double lamda02 = 0.0; int niter = 0; int done = 0; double quad1,quad2,b1,b2,lamda01_new,lamda02_new; while (!done && niter < max_iter) { quad1 = quad1_0101 * lamda01*lamda01 + quad1_0202 * lamda02*lamda02 + quad1_0102 * lamda01*lamda02; quad2 = quad2_0101 * lamda01*lamda01 + quad2_0202 * lamda02*lamda02 + quad2_0102 * lamda01*lamda02; b1 = bond1*bond1 - s01sq - quad1; b2 = bond2*bond2 - s02sq - quad2; lamda01_new = a11inv*b1 + a12inv*b2; lamda02_new = a21inv*b1 + a22inv*b2; done = 1; if (fabs(lamda01_new-lamda01) > tolerance) done = 0; if (fabs(lamda02_new-lamda02) > tolerance) done = 0; lamda01 = lamda01_new; lamda02 = lamda02_new; niter++; } // update forces if atom is owned by this processor lamda01 = lamda01/dtfsq; lamda02 = lamda02/dtfsq; if (i0 < nlocal) { f[i0][0] += lamda01*r01[0] + lamda02*r02[0]; f[i0][1] += lamda01*r01[1] + lamda02*r02[1]; f[i0][2] += lamda01*r01[2] + lamda02*r02[2]; } if (i1 < nlocal) { f[i1][0] -= lamda01*r01[0]; f[i1][1] -= lamda01*r01[1]; f[i1][2] -= lamda01*r01[2]; } if (i2 < nlocal) { f[i2][0] -= lamda02*r02[0]; f[i2][1] -= lamda02*r02[1]; f[i2][2] -= lamda02*r02[2]; } if (evflag) { nlist = 0; if (i0 < nlocal) list[nlist++] = i0; if (i1 < nlocal) list[nlist++] = i1; if (i2 < nlocal) list[nlist++] = i2; v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0]; v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1]; v[2] = lamda01*r01[2]*r01[2] + lamda02*r02[2]*r02[2]; v[3] = lamda01*r01[0]*r01[1] + lamda02*r02[0]*r02[1]; v[4] = lamda01*r01[0]*r01[2] + lamda02*r02[0]*r02[2]; v[5] = lamda01*r01[1]*r01[2] + lamda02*r02[1]*r02[2]; v_tally(nlist,list,3.0,v); } } /* ---------------------------------------------------------------------- */ void FixShake::shake4(int m) { - int nlist,list[4]; + int nlist,list[4]; double v[6]; double invmass0,invmass1,invmass2,invmass3; // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); int i1 = atom->map(shake_atom[m][1]); int i2 = atom->map(shake_atom[m][2]); int i3 = atom->map(shake_atom[m][3]); double bond1 = bond_distance[shake_type[m][0]]; double bond2 = bond_distance[shake_type[m][1]]; double bond3 = bond_distance[shake_type[m][2]]; // r01,r02,r03 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i0][0] - x[i1][0]; r01[1] = x[i0][1] - x[i1][1]; r01[2] = x[i0][2] - x[i1][2]; domain->minimum_image(r01); double r02[3]; r02[0] = x[i0][0] - x[i2][0]; r02[1] = x[i0][1] - x[i2][1]; r02[2] = x[i0][2] - x[i2][2]; domain->minimum_image(r02); double r03[3]; r03[0] = x[i0][0] - x[i3][0]; r03[1] = x[i0][1] - x[i3][1]; r03[2] = x[i0][2] - x[i3][2]; domain->minimum_image(r03); // s01,s02,s03 = distance vec after unconstrained update, with PBC double s01[3]; s01[0] = xshake[i0][0] - xshake[i1][0]; s01[1] = xshake[i0][1] - xshake[i1][1]; s01[2] = xshake[i0][2] - xshake[i1][2]; domain->minimum_image(s01); double s02[3]; s02[0] = xshake[i0][0] - xshake[i2][0]; s02[1] = xshake[i0][1] - xshake[i2][1]; s02[2] = xshake[i0][2] - xshake[i2][2]; domain->minimum_image(s02); double s03[3]; s03[0] = xshake[i0][0] - xshake[i3][0]; s03[1] = xshake[i0][1] - xshake[i3][1]; s03[2] = xshake[i0][2] - xshake[i3][2]; domain->minimum_image(s03); // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; double r03sq = r03[0]*r03[0] + r03[1]*r03[1] + r03[2]*r03[2]; double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; double s03sq = s03[0]*s03[0] + s03[1]*s03[1] + s03[2]*s03[2]; // matrix coeffs and rhs for lamda equations if (rmass) { invmass0 = 1.0/rmass[i0]; invmass1 = 1.0/rmass[i1]; invmass2 = 1.0/rmass[i2]; invmass3 = 1.0/rmass[i3]; } else { invmass0 = 1.0/mass[type[i0]]; invmass1 = 1.0/mass[type[i1]]; invmass2 = 1.0/mass[type[i2]]; invmass3 = 1.0/mass[type[i3]]; } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); double a12 = 2.0 * invmass0 * (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); double a13 = 2.0 * invmass0 * (s01[0]*r03[0] + s01[1]*r03[1] + s01[2]*r03[2]); double a21 = 2.0 * invmass0 * (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); double a22 = 2.0 * (invmass0+invmass2) * (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); double a23 = 2.0 * invmass0 * (s02[0]*r03[0] + s02[1]*r03[1] + s02[2]*r03[2]); double a31 = 2.0 * invmass0 * (s03[0]*r01[0] + s03[1]*r01[1] + s03[2]*r01[2]); double a32 = 2.0 * invmass0 * (s03[0]*r02[0] + s03[1]*r02[1] + s03[2]*r02[2]); double a33 = 2.0 * (invmass0+invmass3) * (s03[0]*r03[0] + s03[1]*r03[1] + s03[2]*r03[2]); // inverse of matrix; double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31; if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); double determinv = 1.0/determ; double a11inv = determinv * (a22*a33 - a23*a32); double a12inv = -determinv * (a12*a33 - a13*a32); double a13inv = determinv * (a12*a23 - a13*a22); double a21inv = -determinv * (a21*a33 - a23*a31); double a22inv = determinv * (a11*a33 - a13*a31); double a23inv = -determinv * (a11*a23 - a13*a21); double a31inv = determinv * (a21*a32 - a22*a31); double a32inv = -determinv * (a11*a32 - a12*a31); double a33inv = determinv * (a11*a22 - a12*a21); // quadratic correction coeffs double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); double r0103 = (r01[0]*r03[0] + r01[1]*r03[1] + r01[2]*r03[2]); double r0203 = (r02[0]*r03[0] + r02[1]*r03[1] + r02[2]*r03[2]); double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; double quad1_0202 = invmass0*invmass0 * r02sq; double quad1_0303 = invmass0*invmass0 * r03sq; double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; double quad1_0103 = 2.0 * (invmass0+invmass1)*invmass0 * r0103; double quad1_0203 = 2.0 * invmass0*invmass0 * r0203; double quad2_0101 = invmass0*invmass0 * r01sq; double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; double quad2_0303 = invmass0*invmass0 * r03sq; double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; double quad2_0103 = 2.0 * invmass0*invmass0 * r0103; double quad2_0203 = 2.0 * (invmass0+invmass2)*invmass0 * r0203; double quad3_0101 = invmass0*invmass0 * r01sq; double quad3_0202 = invmass0*invmass0 * r02sq; double quad3_0303 = (invmass0+invmass3)*(invmass0+invmass3) * r03sq; double quad3_0102 = 2.0 * invmass0*invmass0 * r0102; double quad3_0103 = 2.0 * (invmass0+invmass3)*invmass0 * r0103; double quad3_0203 = 2.0 * (invmass0+invmass3)*invmass0 * r0203; // iterate until converged double lamda01 = 0.0; double lamda02 = 0.0; double lamda03 = 0.0; int niter = 0; int done = 0; double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda03_new; while (!done && niter < max_iter) { quad1 = quad1_0101 * lamda01*lamda01 + quad1_0202 * lamda02*lamda02 + quad1_0303 * lamda03*lamda03 + quad1_0102 * lamda01*lamda02 + quad1_0103 * lamda01*lamda03 + quad1_0203 * lamda02*lamda03; quad2 = quad2_0101 * lamda01*lamda01 + quad2_0202 * lamda02*lamda02 + quad2_0303 * lamda03*lamda03 + quad2_0102 * lamda01*lamda02 + quad2_0103 * lamda01*lamda03 + quad2_0203 * lamda02*lamda03; quad3 = quad3_0101 * lamda01*lamda01 + quad3_0202 * lamda02*lamda02 + quad3_0303 * lamda03*lamda03 + quad3_0102 * lamda01*lamda02 + quad3_0103 * lamda01*lamda03 + quad3_0203 * lamda02*lamda03; b1 = bond1*bond1 - s01sq - quad1; b2 = bond2*bond2 - s02sq - quad2; b3 = bond3*bond3 - s03sq - quad3; lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3; lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3; lamda03_new = a31inv*b1 + a32inv*b2 + a33inv*b3; done = 1; if (fabs(lamda01_new-lamda01) > tolerance) done = 0; if (fabs(lamda02_new-lamda02) > tolerance) done = 0; if (fabs(lamda03_new-lamda03) > tolerance) done = 0; lamda01 = lamda01_new; lamda02 = lamda02_new; lamda03 = lamda03_new; niter++; } // update forces if atom is owned by this processor lamda01 = lamda01/dtfsq; lamda02 = lamda02/dtfsq; lamda03 = lamda03/dtfsq; if (i0 < nlocal) { f[i0][0] += lamda01*r01[0] + lamda02*r02[0] + lamda03*r03[0]; f[i0][1] += lamda01*r01[1] + lamda02*r02[1] + lamda03*r03[1]; f[i0][2] += lamda01*r01[2] + lamda02*r02[2] + lamda03*r03[2]; } if (i1 < nlocal) { f[i1][0] -= lamda01*r01[0]; f[i1][1] -= lamda01*r01[1]; f[i1][2] -= lamda01*r01[2]; } if (i2 < nlocal) { f[i2][0] -= lamda02*r02[0]; f[i2][1] -= lamda02*r02[1]; f[i2][2] -= lamda02*r02[2]; } if (i3 < nlocal) { f[i3][0] -= lamda03*r03[0]; f[i3][1] -= lamda03*r03[1]; f[i3][2] -= lamda03*r03[2]; } if (evflag) { nlist = 0; if (i0 < nlocal) list[nlist++] = i0; if (i1 < nlocal) list[nlist++] = i1; if (i2 < nlocal) list[nlist++] = i2; if (i3 < nlocal) list[nlist++] = i3; v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0]; v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1]; v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda03*r03[2]*r03[2]; v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda03*r03[0]*r03[1]; v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda03*r03[0]*r03[2]; v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda03*r03[1]*r03[2]; v_tally(nlist,list,4.0,v); } } /* ---------------------------------------------------------------------- */ void FixShake::shake3angle(int m) { int nlist,list[3]; double v[6]; double invmass0,invmass1,invmass2; // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); int i1 = atom->map(shake_atom[m][1]); int i2 = atom->map(shake_atom[m][2]); double bond1 = bond_distance[shake_type[m][0]]; double bond2 = bond_distance[shake_type[m][1]]; double bond12 = angle_distance[shake_type[m][2]]; // r01,r02,r12 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i0][0] - x[i1][0]; r01[1] = x[i0][1] - x[i1][1]; r01[2] = x[i0][2] - x[i1][2]; domain->minimum_image(r01); double r02[3]; r02[0] = x[i0][0] - x[i2][0]; r02[1] = x[i0][1] - x[i2][1]; r02[2] = x[i0][2] - x[i2][2]; domain->minimum_image(r02); double r12[3]; r12[0] = x[i1][0] - x[i2][0]; r12[1] = x[i1][1] - x[i2][1]; r12[2] = x[i1][2] - x[i2][2]; domain->minimum_image(r12); // s01,s02,s12 = distance vec after unconstrained update, with PBC double s01[3]; s01[0] = xshake[i0][0] - xshake[i1][0]; s01[1] = xshake[i0][1] - xshake[i1][1]; s01[2] = xshake[i0][2] - xshake[i1][2]; domain->minimum_image(s01); double s02[3]; s02[0] = xshake[i0][0] - xshake[i2][0]; s02[1] = xshake[i0][1] - xshake[i2][1]; s02[2] = xshake[i0][2] - xshake[i2][2]; domain->minimum_image(s02); double s12[3]; s12[0] = xshake[i1][0] - xshake[i2][0]; s12[1] = xshake[i1][1] - xshake[i2][1]; s12[2] = xshake[i1][2] - xshake[i2][2]; domain->minimum_image(s12); // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; double r12sq = r12[0]*r12[0] + r12[1]*r12[1] + r12[2]*r12[2]; double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; double s12sq = s12[0]*s12[0] + s12[1]*s12[1] + s12[2]*s12[2]; // matrix coeffs and rhs for lamda equations if (rmass) { invmass0 = 1.0/rmass[i0]; invmass1 = 1.0/rmass[i1]; invmass2 = 1.0/rmass[i2]; } else { invmass0 = 1.0/mass[type[i0]]; invmass1 = 1.0/mass[type[i1]]; invmass2 = 1.0/mass[type[i2]]; } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); double a12 = 2.0 * invmass0 * (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); double a13 = - 2.0 * invmass1 * (s01[0]*r12[0] + s01[1]*r12[1] + s01[2]*r12[2]); double a21 = 2.0 * invmass0 * (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); double a22 = 2.0 * (invmass0+invmass2) * (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); double a23 = 2.0 * invmass2 * (s02[0]*r12[0] + s02[1]*r12[1] + s02[2]*r12[2]); double a31 = - 2.0 * invmass1 * (s12[0]*r01[0] + s12[1]*r01[1] + s12[2]*r01[2]); double a32 = 2.0 * invmass2 * (s12[0]*r02[0] + s12[1]*r02[1] + s12[2]*r02[2]); double a33 = 2.0 * (invmass1+invmass2) * (s12[0]*r12[0] + s12[1]*r12[1] + s12[2]*r12[2]); // inverse of matrix double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31; if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); double determinv = 1.0/determ; double a11inv = determinv * (a22*a33 - a23*a32); double a12inv = -determinv * (a12*a33 - a13*a32); double a13inv = determinv * (a12*a23 - a13*a22); double a21inv = -determinv * (a21*a33 - a23*a31); double a22inv = determinv * (a11*a33 - a13*a31); double a23inv = -determinv * (a11*a23 - a13*a21); double a31inv = determinv * (a21*a32 - a22*a31); double a32inv = -determinv * (a11*a32 - a12*a31); double a33inv = determinv * (a11*a22 - a12*a21); // quadratic correction coeffs double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); double r0112 = (r01[0]*r12[0] + r01[1]*r12[1] + r01[2]*r12[2]); double r0212 = (r02[0]*r12[0] + r02[1]*r12[1] + r02[2]*r12[2]); double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; double quad1_0202 = invmass0*invmass0 * r02sq; double quad1_1212 = invmass1*invmass1 * r12sq; double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; double quad1_0112 = - 2.0 * (invmass0+invmass1)*invmass1 * r0112; double quad1_0212 = - 2.0 * invmass0*invmass1 * r0212; double quad2_0101 = invmass0*invmass0 * r01sq; double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; double quad2_1212 = invmass2*invmass2 * r12sq; double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; double quad2_0112 = 2.0 * invmass0*invmass2 * r0112; double quad2_0212 = 2.0 * (invmass0+invmass2)*invmass2 * r0212; double quad3_0101 = invmass1*invmass1 * r01sq; double quad3_0202 = invmass2*invmass2 * r02sq; double quad3_1212 = (invmass1+invmass2)*(invmass1+invmass2) * r12sq; double quad3_0102 = - 2.0 * invmass1*invmass2 * r0102; double quad3_0112 = - 2.0 * (invmass1+invmass2)*invmass1 * r0112; double quad3_0212 = 2.0 * (invmass1+invmass2)*invmass2 * r0212; // iterate until converged double lamda01 = 0.0; double lamda02 = 0.0; double lamda12 = 0.0; int niter = 0; int done = 0; double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda12_new; while (!done && niter < max_iter) { quad1 = quad1_0101 * lamda01*lamda01 + quad1_0202 * lamda02*lamda02 + quad1_1212 * lamda12*lamda12 + quad1_0102 * lamda01*lamda02 + quad1_0112 * lamda01*lamda12 + quad1_0212 * lamda02*lamda12; quad2 = quad2_0101 * lamda01*lamda01 + quad2_0202 * lamda02*lamda02 + quad2_1212 * lamda12*lamda12 + quad2_0102 * lamda01*lamda02 + quad2_0112 * lamda01*lamda12 + quad2_0212 * lamda02*lamda12; quad3 = quad3_0101 * lamda01*lamda01 + quad3_0202 * lamda02*lamda02 + quad3_1212 * lamda12*lamda12 + quad3_0102 * lamda01*lamda02 + quad3_0112 * lamda01*lamda12 + quad3_0212 * lamda02*lamda12; b1 = bond1*bond1 - s01sq - quad1; b2 = bond2*bond2 - s02sq - quad2; b3 = bond12*bond12 - s12sq - quad3; lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3; lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3; lamda12_new = a31inv*b1 + a32inv*b2 + a33inv*b3; done = 1; if (fabs(lamda01_new-lamda01) > tolerance) done = 0; if (fabs(lamda02_new-lamda02) > tolerance) done = 0; if (fabs(lamda12_new-lamda12) > tolerance) done = 0; lamda01 = lamda01_new; lamda02 = lamda02_new; lamda12 = lamda12_new; niter++; } // update forces if atom is owned by this processor lamda01 = lamda01/dtfsq; lamda02 = lamda02/dtfsq; lamda12 = lamda12/dtfsq; if (i0 < nlocal) { f[i0][0] += lamda01*r01[0] + lamda02*r02[0]; f[i0][1] += lamda01*r01[1] + lamda02*r02[1]; f[i0][2] += lamda01*r01[2] + lamda02*r02[2]; } if (i1 < nlocal) { f[i1][0] -= lamda01*r01[0] - lamda12*r12[0]; f[i1][1] -= lamda01*r01[1] - lamda12*r12[1]; f[i1][2] -= lamda01*r01[2] - lamda12*r12[2]; } if (i2 < nlocal) { f[i2][0] -= lamda02*r02[0] + lamda12*r12[0]; f[i2][1] -= lamda02*r02[1] + lamda12*r12[1]; f[i2][2] -= lamda02*r02[2] + lamda12*r12[2]; } if (evflag) { nlist = 0; if (i0 < nlocal) list[nlist++] = i0; if (i1 < nlocal) list[nlist++] = i1; if (i2 < nlocal) list[nlist++] = i2; v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0]; v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1]; v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda12*r12[2]*r12[2]; v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda12*r12[0]*r12[1]; v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda12*r12[0]*r12[2]; v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda12*r12[1]*r12[2]; v_tally(nlist,list,3.0,v); } } /* ---------------------------------------------------------------------- print-out bond & angle statistics ------------------------------------------------------------------------- */ void FixShake::stats() { int i,j,m,n,iatom,jatom,katom; double delx,dely,delz; double r,r1,r2,r3,angle; // zero out accumulators int nb = atom->nbondtypes + 1; int na = atom->nangletypes + 1; for (i = 0; i < nb; i++) { b_count[i] = 0; b_ave[i] = b_max[i] = 0.0; b_min[i] = BIG; } for (i = 0; i < na; i++) { a_count[i] = 0; a_ave[i] = a_max[i] = 0.0; a_min[i] = BIG; } // log stats for each bond & angle // OK to double count since are just averaging double **x = atom->x; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) { if (shake_flag[i] == 0) continue; // bond stats n = shake_flag[i]; if (n == 1) n = 3; iatom = atom->map(shake_atom[i][0]); for (j = 1; j < n; j++) { jatom = atom->map(shake_atom[i][j]); delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); m = shake_type[i][j-1]; b_count[m]++; b_ave[m] += r; b_max[m] = MAX(b_max[m],r); b_min[m] = MIN(b_min[m],r); } // angle stats if (shake_flag[i] == 1) { iatom = atom->map(shake_atom[i][0]); jatom = atom->map(shake_atom[i][1]); katom = atom->map(shake_atom[i][2]); delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; domain->minimum_image(delx,dely,delz); r1 = sqrt(delx*delx + dely*dely + delz*delz); delx = x[iatom][0] - x[katom][0]; dely = x[iatom][1] - x[katom][1]; delz = x[iatom][2] - x[katom][2]; domain->minimum_image(delx,dely,delz); r2 = sqrt(delx*delx + dely*dely + delz*delz); delx = x[jatom][0] - x[katom][0]; dely = x[jatom][1] - x[katom][1]; delz = x[jatom][2] - x[katom][2]; domain->minimum_image(delx,dely,delz); r3 = sqrt(delx*delx + dely*dely + delz*delz); angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); angle *= 180.0/MY_PI; m = shake_type[i][2]; a_count[m]++; a_ave[m] += angle; a_max[m] = MAX(a_max[m],angle); a_min[m] = MIN(a_min[m],angle); } } // sum across all procs MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world); MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world); MPI_Allreduce(a_count,a_count_all,na,MPI_INT,MPI_SUM,world); MPI_Allreduce(a_ave,a_ave_all,na,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(a_max,a_max_all,na,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(a_min,a_min_all,na,MPI_DOUBLE,MPI_MIN,world); // print stats only for non-zero counts if (me == 0) { if (screen) { fprintf(screen, "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n", update->ntimestep); for (i = 1; i < nb; i++) if (b_count_all[i]) fprintf(screen," %d %g %g %d\n",i, b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i], b_count_all[i]); for (i = 1; i < na; i++) if (a_count_all[i]) fprintf(screen," %d %g %g\n",i, a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]); } if (logfile) { fprintf(logfile, "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n", update->ntimestep); for (i = 0; i < nb; i++) if (b_count_all[i]) fprintf(logfile," %d %g %g\n",i, b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i]); for (i = 0; i < na; i++) if (a_count_all[i]) fprintf(logfile," %d %g %g\n",i, a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]); } } // next timestep for stats next_output += output_every; } /* ---------------------------------------------------------------------- find a bond between global atom IDs n1 and n2 stored with local atom i if find it: if setflag = 0, return bond type if setflag = -1/1, set bond type to negative/positive and return 0 if do not find it, return 0 ------------------------------------------------------------------------- */ int FixShake::bondtype_findset(int i, tagint n1, tagint n2, int setflag) { int m,nbonds; int *btype; if (molecular == 1) { tagint *tag = atom->tag; tagint **bond_atom = atom->bond_atom; nbonds = atom->num_bond[i]; for (m = 0; m < nbonds; m++) { if (n1 == tag[i] && n2 == bond_atom[i][m]) break; if (n1 == bond_atom[i][m] && n2 == tag[i]) break; } } else { int imol = atom->molindex[i]; int iatom = atom->molatom[i]; tagint *tag = atom->tag; tagint tagprev = tag[i] - iatom - 1; tagint *batom = atommols[imol]->bond_atom[iatom]; btype = atommols[imol]->bond_type[iatom]; nbonds = atommols[imol]->num_bond[iatom]; for (m = 0; m < nbonds; m++) { if (n1 == tag[i] && n2 == batom[m]+tagprev) break; if (n1 == batom[m]+tagprev && n2 == tag[i]) break; } } if (m < nbonds) { if (setflag == 0) { if (molecular == 1) return atom->bond_type[i][m]; else return btype[m]; } if (molecular == 1) { if ((setflag < 0 && atom->bond_type[i][m] > 0) || (setflag > 0 && atom->bond_type[i][m] < 0)) atom->bond_type[i][m] = -atom->bond_type[i][m]; } else { if ((setflag < 0 && btype[m] > 0) || (setflag > 0 && btype[m] < 0)) btype[m] = -btype[m]; } } return 0; } /* ---------------------------------------------------------------------- find an angle with global end atom IDs n1 and n2 stored with local atom i if find it: if setflag = 0, return angle type if setflag = -1/1, set angle type to negative/positive and return 0 if do not find it, return 0 ------------------------------------------------------------------------- */ int FixShake::angletype_findset(int i, tagint n1, tagint n2, int setflag) { int m,nangles; int *atype; if (molecular == 1) { tagint **angle_atom1 = atom->angle_atom1; tagint **angle_atom3 = atom->angle_atom3; nangles = atom->num_angle[i]; for (m = 0; m < nangles; m++) { if (n1 == angle_atom1[i][m] && n2 == angle_atom3[i][m]) break; if (n1 == angle_atom3[i][m] && n2 == angle_atom1[i][m]) break; } } else { int imol = atom->molindex[i]; int iatom = atom->molatom[i]; tagint *tag = atom->tag; tagint tagprev = tag[i] - iatom - 1; tagint *aatom1 = atommols[imol]->angle_atom1[iatom]; tagint *aatom3 = atommols[imol]->angle_atom3[iatom]; atype = atommols[imol]->angle_type[iatom]; nangles = atommols[imol]->num_angle[iatom]; for (m = 0; m < nangles; m++) { if (n1 == aatom1[m]+tagprev && n2 == aatom3[m]+tagprev) break; if (n1 == aatom3[m]+tagprev && n2 == aatom1[m]+tagprev) break; } } if (m < nangles) { if (setflag == 0) { if (molecular == 1) return atom->angle_type[i][m]; else return atype[m]; } if (molecular == 1) { if ((setflag < 0 && atom->angle_type[i][m] > 0) || (setflag > 0 && atom->angle_type[i][m] < 0)) atom->angle_type[i][m] = -atom->angle_type[i][m]; } else { if ((setflag < 0 && atype[m] > 0) || (setflag > 0 && atype[m] < 0)) atype[m] = -atype[m]; } } return 0; } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixShake::memory_usage() { int nmax = atom->nmax; double bytes = nmax * sizeof(int); bytes += nmax*4 * sizeof(int); bytes += nmax*3 * sizeof(int); bytes += nmax*3 * sizeof(double); bytes += maxvatom*6 * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- allocate local atom-based arrays ------------------------------------------------------------------------- */ void FixShake::grow_arrays(int nmax) { memory->grow(shake_flag,nmax,"shake:shake_flag"); memory->grow(shake_atom,nmax,4,"shake:shake_atom"); memory->grow(shake_type,nmax,3,"shake:shake_type"); memory->destroy(xshake); memory->create(xshake,nmax,3,"shake:xshake"); } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ void FixShake::copy_arrays(int i, int j, int delflag) { int flag = shake_flag[j] = shake_flag[i]; if (flag == 1) { shake_atom[j][0] = shake_atom[i][0]; shake_atom[j][1] = shake_atom[i][1]; shake_atom[j][2] = shake_atom[i][2]; shake_type[j][0] = shake_type[i][0]; shake_type[j][1] = shake_type[i][1]; shake_type[j][2] = shake_type[i][2]; } else if (flag == 2) { shake_atom[j][0] = shake_atom[i][0]; shake_atom[j][1] = shake_atom[i][1]; shake_type[j][0] = shake_type[i][0]; } else if (flag == 3) { shake_atom[j][0] = shake_atom[i][0]; shake_atom[j][1] = shake_atom[i][1]; shake_atom[j][2] = shake_atom[i][2]; shake_type[j][0] = shake_type[i][0]; shake_type[j][1] = shake_type[i][1]; } else if (flag == 4) { shake_atom[j][0] = shake_atom[i][0]; shake_atom[j][1] = shake_atom[i][1]; shake_atom[j][2] = shake_atom[i][2]; shake_atom[j][3] = shake_atom[i][3]; shake_type[j][0] = shake_type[i][0]; shake_type[j][1] = shake_type[i][1]; shake_type[j][2] = shake_type[i][2]; } } /* ---------------------------------------------------------------------- initialize one atom's array values, called when atom is created ------------------------------------------------------------------------- */ void FixShake::set_arrays(int i) { shake_flag[i] = 0; } /* ---------------------------------------------------------------------- update one atom's array values called when molecule is created from fix gcmc ------------------------------------------------------------------------- */ void FixShake::update_arrays(int i, int atom_offset) { int flag = shake_flag[i]; if (flag == 1) { shake_atom[i][0] += atom_offset; shake_atom[i][1] += atom_offset; shake_atom[i][2] += atom_offset; } else if (flag == 2) { shake_atom[i][0] += atom_offset; shake_atom[i][1] += atom_offset; } else if (flag == 3) { shake_atom[i][0] += atom_offset; shake_atom[i][1] += atom_offset; shake_atom[i][2] += atom_offset; } else if (flag == 4) { shake_atom[i][0] += atom_offset; shake_atom[i][1] += atom_offset; shake_atom[i][2] += atom_offset; shake_atom[i][3] += atom_offset; } } /* ---------------------------------------------------------------------- initialize a molecule inserted by another fix, e.g. deposit or pour called when molecule is created nlocalprev = # of atoms on this proc before molecule inserted tagprev = atom ID previous to new atoms in the molecule xgeom,vcm,quat ignored ------------------------------------------------------------------------- */ void FixShake::set_molecule(int nlocalprev, tagint tagprev, int imol, double *xgeom, double *vcm, double *quat) { int m,flag; int nlocal = atom->nlocal; if (nlocalprev == nlocal) return; tagint *tag = atom->tag; tagint **mol_shake_atom = onemols[imol]->shake_atom; int **mol_shake_type = onemols[imol]->shake_type; for (int i = nlocalprev; i < nlocal; i++) { m = tag[i] - tagprev-1; flag = shake_flag[i] = onemols[imol]->shake_flag[m]; if (flag == 1) { shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; shake_type[i][0] = mol_shake_type[m][0]; shake_type[i][1] = mol_shake_type[m][1]; shake_type[i][2] = mol_shake_type[m][2]; } else if (flag == 2) { shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; shake_type[i][0] = mol_shake_type[m][0]; } else if (flag == 3) { shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; shake_type[i][0] = mol_shake_type[m][0]; shake_type[i][1] = mol_shake_type[m][1]; } else if (flag == 4) { shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; shake_atom[i][3] = mol_shake_atom[m][3] + tagprev; shake_type[i][0] = mol_shake_type[m][0]; shake_type[i][1] = mol_shake_type[m][1]; shake_type[i][2] = mol_shake_type[m][2]; } } } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ int FixShake::pack_exchange(int i, double *buf) { int m = 0; buf[m++] = shake_flag[i]; int flag = shake_flag[i]; if (flag == 1) { buf[m++] = shake_atom[i][0]; buf[m++] = shake_atom[i][1]; buf[m++] = shake_atom[i][2]; buf[m++] = shake_type[i][0]; buf[m++] = shake_type[i][1]; buf[m++] = shake_type[i][2]; } else if (flag == 2) { buf[m++] = shake_atom[i][0]; buf[m++] = shake_atom[i][1]; buf[m++] = shake_type[i][0]; } else if (flag == 3) { buf[m++] = shake_atom[i][0]; buf[m++] = shake_atom[i][1]; buf[m++] = shake_atom[i][2]; buf[m++] = shake_type[i][0]; buf[m++] = shake_type[i][1]; } else if (flag == 4) { buf[m++] = shake_atom[i][0]; buf[m++] = shake_atom[i][1]; buf[m++] = shake_atom[i][2]; buf[m++] = shake_atom[i][3]; buf[m++] = shake_type[i][0]; buf[m++] = shake_type[i][1]; buf[m++] = shake_type[i][2]; } return m; } /* ---------------------------------------------------------------------- unpack values in local atom-based arrays from exchange with another proc ------------------------------------------------------------------------- */ int FixShake::unpack_exchange(int nlocal, double *buf) { int m = 0; int flag = shake_flag[nlocal] = static_cast (buf[m++]); if (flag == 1) { shake_atom[nlocal][0] = static_cast (buf[m++]); shake_atom[nlocal][1] = static_cast (buf[m++]); shake_atom[nlocal][2] = static_cast (buf[m++]); shake_type[nlocal][0] = static_cast (buf[m++]); shake_type[nlocal][1] = static_cast (buf[m++]); shake_type[nlocal][2] = static_cast (buf[m++]); } else if (flag == 2) { shake_atom[nlocal][0] = static_cast (buf[m++]); shake_atom[nlocal][1] = static_cast (buf[m++]); shake_type[nlocal][0] = static_cast (buf[m++]); } else if (flag == 3) { shake_atom[nlocal][0] = static_cast (buf[m++]); shake_atom[nlocal][1] = static_cast (buf[m++]); shake_atom[nlocal][2] = static_cast (buf[m++]); shake_type[nlocal][0] = static_cast (buf[m++]); shake_type[nlocal][1] = static_cast (buf[m++]); } else if (flag == 4) { shake_atom[nlocal][0] = static_cast (buf[m++]); shake_atom[nlocal][1] = static_cast (buf[m++]); shake_atom[nlocal][2] = static_cast (buf[m++]); shake_atom[nlocal][3] = static_cast (buf[m++]); shake_type[nlocal][0] = static_cast (buf[m++]); shake_type[nlocal][1] = static_cast (buf[m++]); shake_type[nlocal][2] = static_cast (buf[m++]); } return m; } /* ---------------------------------------------------------------------- */ int FixShake::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = xshake[j][0]; buf[m++] = xshake[j][1]; buf[m++] = xshake[j][2]; } } else { if (domain->triclinic == 0) { dx = pbc[0]*domain->xprd; dy = pbc[1]*domain->yprd; dz = pbc[2]*domain->zprd; } else { dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; dz = pbc[2]*domain->zprd; } for (i = 0; i < n; i++) { j = list[i]; buf[m++] = xshake[j][0] + dx; buf[m++] = xshake[j][1] + dy; buf[m++] = xshake[j][2] + dz; } } return m; } /* ---------------------------------------------------------------------- */ void FixShake::unpack_forward_comm(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { xshake[i][0] = buf[m++]; xshake[i][1] = buf[m++]; xshake[i][2] = buf[m++]; } } /* ---------------------------------------------------------------------- */ void FixShake::reset_dt() { if (strstr(update->integrate_style,"verlet")) { dtv = update->dt; - dtfsq = update->dt * update->dt * force->ftm2v; + if (rattle) dtfsq = 0.5 * update->dt * update->dt * force->ftm2v; + else dtfsq = update->dt * update->dt * force->ftm2v; } else { dtv = step_respa[0]; dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v; - dtf_inner = step_respa[0] * force->ftm2v; + if (rattle) dtf_inner = dtf_innerhalf; + else dtf_inner = step_respa[0] * force->ftm2v; } } /* ---------------------------------------------------------------------- extract Molecule ptr ------------------------------------------------------------------------- */ void *FixShake::extract(const char *str, int &dim) { dim = 0; if (strcmp(str,"onemol") == 0) { return onemols; } return NULL; } diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 8ba856aa2..74e09d12d 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -1,256 +1,257 @@ /* -*- 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(shake,FixShake) #else #ifndef LMP_FIX_SHAKE_H #define LMP_FIX_SHAKE_H #include "fix.h" namespace LAMMPS_NS { class FixShake : public Fix { public: FixShake(class LAMMPS *, int, char **); - ~FixShake(); - int setmask(); - void init(); + virtual ~FixShake(); + virtual int setmask(); + virtual void init(); void setup(int); void pre_neighbor(); - void post_force(int); - void post_force_respa(int, int, int); + virtual void post_force(int); + virtual void post_force_respa(int, int, int); - double memory_usage(); - void grow_arrays(int); - void copy_arrays(int, int, int); + virtual double memory_usage(); + virtual void grow_arrays(int); + virtual void copy_arrays(int, int, int); void set_arrays(int); - void update_arrays(int, int); + virtual void update_arrays(int, int); void set_molecule(int, tagint, int, double *, double *, double *); - int pack_exchange(int, double *); - int unpack_exchange(int, double *); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); + virtual int pack_exchange(int, double *); + virtual int unpack_exchange(int, double *); + virtual int pack_forward_comm(int, int *, double *, int, int *); + virtual void unpack_forward_comm(int, int, double *); int dof(int); - void reset_dt(); + virtual void reset_dt(); void *extract(const char *, int &); - private: + protected: int me,nprocs; + int rattle; // 0 = SHAKE, 1 = RATTLE double tolerance; // SHAKE tolerance int max_iter; // max # of SHAKE iterations int output_every; // SHAKE stat output every so often bigint next_output; // timestep for next output // settings from input command int *bond_flag,*angle_flag; // bond/angle types to constrain int *type_flag; // constrain bonds to these types double *mass_list; // constrain bonds to these masses int nmass; // # of masses in mass_list int molecular; // copy of atom->molecular double *bond_distance,*angle_distance; // constraint distances int ifix_respa; // rRESPA fix needed by SHAKE int nlevels_respa; // copies of needed rRESPA variables int *loop_respa; double *step_respa; double **x,**v,**f; // local ptrs to atom class quantities double *mass,*rmass; int *type; int nlocal; // atom-based arrays int *shake_flag; // 0 if atom not in SHAKE cluster // 1 = size 3 angle cluster // 2,3,4 = size of bond-only cluster tagint **shake_atom; // global IDs of atoms in cluster // central atom is 1st // lowest global ID is 1st for size 2 int **shake_type; // bondtype of each bond in cluster // for angle cluster, 3rd value // is angletype double **xshake; // unconstrained atom coords int *nshake; // count double dtv,dtfsq; // timesteps for trial move double dtf_inner,dtf_innerhalf; // timesteps for rRESPA trial move int *list; // list of clusters to SHAKE int nlist,maxlist; // size and max-size of list // stat quantities int *b_count,*b_count_all; // counts for each bond type double *b_ave,*b_max,*b_min; // ave/max/min dist for each bond type double *b_ave_all,*b_max_all,*b_min_all; // MPI summing arrays int *a_count,*a_count_all; // ditto for angle types double *a_ave,*a_max,*a_min; double *a_ave_all,*a_max_all,*a_min_all; class Molecule **atommols; // atom style template pointer class Molecule **onemols; // molecule added on-the-fly int nmol; void find_clusters(); int masscheck(double); void unconstrained_update(); void unconstrained_update_respa(int); void shake(int); void shake3(int); void shake4(int); void shake3angle(int); void stats(); int bondtype_findset(int, tagint, tagint, int); int angletype_findset(int, tagint, tagint, int); // static variable for ring communication callback to access class data // callback functions for ring communication static FixShake *fsptr; static void ring_bonds(int, char *); static void ring_nshake(int, char *); static void ring_shake(int, char *); }; } #endif #endif /* ERROR/WARNING messages: E: Cannot use fix shake with non-molecular system Your choice of atom style does not have bonds. 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. E: Invalid bond type index for fix shake Self-explanatory. Check the fix shake command in the input script. E: Invalid angle type index for fix shake Self-explanatory. E: Invalid atom type index for fix shake Atom types must range from 1 to Ntypes inclusive. E: Invalid atom mass for fix shake Mass specified in fix shake command must be > 0.0. E: Too many masses for fix shake The fix shake command cannot list more masses than there are atom types. E: Molecule template ID for fix shake does not exist Self-explanatory. W: Molecule template for fix shake has multiple molecules The fix shake command will only recoginze molecules of a single type, i.e. the first molecule in the template. E: Fix shake molecule template must have shake info The defined molecule does not specify SHAKE information. E: More than one fix shake Only one fix shake can be defined. E: Fix shake cannot be used with minimization Cannot use fix shake while doing an energy minimization since it turns off bonds that should contribute to the energy. E: Shake fix must come before NPT/NPH fix NPT fix must be defined in input script after SHAKE fix, else the SHAKE fix contribution to the pressure virial is incorrect. E: Bond potential must be defined for SHAKE Cannot use fix shake unless bond potential is defined. E: Angle potential must be defined for SHAKE When shaking angles, an angle_style potential must be used. E: Shake angles have different bond types All 3-atom angle-constrained SHAKE clusters specified by the fix shake command that are the same angle type, must also have the same bond types for the 2 bonds in the angle. E: Shake atoms %d %d missing on proc %d at step %ld The 2 atoms in a single shake cluster specified by the fix shake command are not all accessible to a processor. This probably means an atom has moved too far. E: Shake atoms %d %d %d missing on proc %d at step %ld The 3 atoms in a single shake cluster specified by the fix shake command are not all accessible to a processor. This probably means an atom has moved too far. E: Shake atoms %d %d %d %d missing on proc %d at step %ld The 4 atoms in a single shake cluster specified by the fix shake command are not all accessible to a processor. This probably means an atom has moved too far. E: Did not find fix shake partner info Could not find bond partners implied by fix shake command. This error can be triggered if the delete_bonds command was used before fix shake, and it removed bonds without resetting the 1-2, 1-3, 1-4 weighting list via the special keyword. E: Shake cluster of more than 4 atoms A single cluster specified by the fix shake command can have no more than 4 atoms. E: Shake clusters are connected A single cluster specified by the fix shake command must have a single central atom with up to 3 other atoms bonded to it. W: Shake determinant < 0.0 The determinant of the quadratic equation being solved for a single cluster specified by the fix shake command is numerically suspect. LAMMPS will set it to 0.0 and continue. E: Shake determinant = 0.0 The determinant of the matrix being solved for a single cluster specified by the fix shake command is numerically invalid. */ diff --git a/src/fix_respa.cpp b/src/fix_respa.cpp index d346712b9..03c227fd9 100644 --- a/src/fix_respa.cpp +++ b/src/fix_respa.cpp @@ -1,121 +1,120 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "stdlib.h" #include "fix_respa.h" #include "atom.h" -#include "memory.h" -#include "error.h" #include "force.h" +#include "memory.h" using namespace LAMMPS_NS; using namespace FixConst; /* ---------------------------------------------------------------------- */ FixRespa::FixRespa(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { // nlevels = # of rRESPA levels nlevels = force->inumeric(FLERR,arg[3]); // perform initial allocation of atom-based arrays // register with Atom class f_level = NULL; grow_arrays(atom->nmax); atom->add_callback(0); } /* ---------------------------------------------------------------------- */ FixRespa::~FixRespa() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); // delete locally stored arrays memory->destroy(f_level); } /* ---------------------------------------------------------------------- */ int FixRespa::setmask() { return 0; } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixRespa::memory_usage() { double bytes = atom->nmax*nlevels*3 * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- allocate local atom-based arrays ------------------------------------------------------------------------- */ void FixRespa::grow_arrays(int nmax) { memory->grow(f_level,nmax,nlevels,3,"fix_respa:f_level"); } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ void FixRespa::copy_arrays(int i, int j, int delflag) { for (int k = 0; k < nlevels; k++) { f_level[j][k][0] = f_level[i][k][0]; f_level[j][k][1] = f_level[i][k][1]; f_level[j][k][2] = f_level[i][k][2]; } } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ int FixRespa::pack_exchange(int i, double *buf) { int m = 0; for (int k = 0; k < nlevels; k++) { buf[m++] = f_level[i][k][0]; buf[m++] = f_level[i][k][1]; buf[m++] = f_level[i][k][2]; } return m; } /* ---------------------------------------------------------------------- unpack values in local atom-based arrays from exchange with another proc ------------------------------------------------------------------------- */ int FixRespa::unpack_exchange(int nlocal, double *buf) { int m = 0; for (int k = 0; k < nlevels; k++) { f_level[nlocal][k][0] = buf[m++]; f_level[nlocal][k][1] = buf[m++]; f_level[nlocal][k][2] = buf[m++]; } return m; } diff --git a/src/fix_respa.h b/src/fix_respa.h index f01d9339a..d65587cdc 100644 --- a/src/fix_respa.h +++ b/src/fix_respa.h @@ -1,55 +1,56 @@ /* -*- 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(RESPA,FixRespa) #else #ifndef LMP_FIX_RESPA_H #define LMP_FIX_RESPA_H #include "fix.h" namespace LAMMPS_NS { class FixRespa : public Fix { friend class Respa; friend class FixShake; - friend class FixShake2; + friend class FixRattle; public: FixRespa(class LAMMPS *, int, char **); ~FixRespa(); int setmask(); void init() {} double memory_usage(); void grow_arrays(int); void copy_arrays(int, int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); private: int nlevels; double ***f_level; // force at each rRESPA level }; } #endif #endif + /* ERROR/WARNING messages: */