Merge branch 'nwchem' of github.com:lammps/lammps into nwchem

examples/README

@@ -98,6 +98,7 @@ nb3b:     use of nonbonded 3-body harmonic pair style
 neb:      nudged elastic band (NEB) calculation for barrier finding
 nemd:     non-equilibrium MD of 2d sheared system
 numdiff:  numerical difference computation of forces, virial, and Born matrix
+nwchem:   wrapping NWChem DFT functionality with LAMMPS
 obstacle: flow around two voids in a 2d channel
 peptide:  dynamics of a small solvated peptide chain (5-mer)
 peri:     Peridynamic model of cylinder impacted by indenter

examples/nwchem/LAMMPS/data.POSCAR_bcc

@@ -0,0 +1,12 @@
+LAMMPS data file from VASP file VASP/POSCAR_bcc
+
+2 atoms
+1 atom types
+0.0   3.1654062429393064 xlo xhi
+0.0   3.1654062429393064 ylo yhi
+0.0   3.1654062429393064 zlo zhi
+
+Atoms
+
+1 1                    0                    0                    0
+2 1   1.5827031214696532   1.5827031214696532   1.5827031214696532

examples/nwchem/LAMMPS/data.POSCAR_bcc_222

@@ -0,0 +1,29 @@
+LAMMPS data file from VASP file VASP/POSCAR_bcc_222
+
+18 atoms
+1 atom types
+0.0   8.9531208783704628 xlo xhi
+0.0   3.8768150619108339 ylo yhi
+0.0   45.534159244872015 zlo zhi
+  2.2382802195926157                    0                    0 xy xz yz
+
+Atoms
+
+1 1                    0                    0                   15
+2 1   6.7148406587778462   1.2922716873036106   15.913774073227765
+3 1   4.4765604391852314   2.5845433746072226   16.827548146455531
+4 1   11.191401097963077    3.876815061910833   17.741322219683298
+5 1   6.7148406587778471   1.2922716873036111   18.655096292911061
+6 1   4.4765604391852305   2.5845433746072222   19.568870366138828
+7 1   2.2382802195926152    3.876815061910833   20.482644439366592
+8 1   6.7148406587778453   1.2922716873036111   21.396418512594359
+9 1   4.4765604391852305   2.5845433746072222   22.310192585822122
+10 1   2.2382802195926148   3.8768150619108321   23.223966659049893
+11 1   6.7148406587778453   1.2922716873036102   24.137740732277656
+12 1   4.4765604391852305   2.5845433746072195   25.051514805505423
+13 1   2.2382802195926148   3.8768150619108321   25.965288878733183
+14 1   6.7148406587778444   1.2922716873036084   26.879062951960954
+15 1   4.4765604391852314   2.5845433746072204   27.792837025188717
+16 1   2.2382802195926139   3.8768150619108304   28.706611098416481
+17 1   6.7148406587778453   1.2922716873036102   29.620385171644248
+18 1   4.4765604391852269   2.5845433746072204   30.534159244872015

examples/nwchem/LAMMPS/data.POSCAR_diamond

@@ -0,0 +1,14 @@
+LAMMPS data file from VASP file VASP/POSCAR_diamond
+
+4 atoms
+1 atom types
+0.0   2.8200563642280549 xlo xhi
+0.0   2.8200563642280549 ylo yhi
+0.0   3.9881619569478763 zlo zhi
+
+Atoms
+
+1 1                    0                    0                    0
+2 1   1.4100281821140275                    0  0.99704048923696909
+3 1   1.4100281821140275   1.4100281821140275   1.9940809784739382
+4 1                    0   1.4100281821140275    2.991121467710907

examples/nwchem/LAMMPS/data.POSCAR_fcc_001

@@ -0,0 +1,28 @@
+LAMMPS data file from VASP file VASP/POSCAR_fcc_001
+
+18 atoms
+1 atom types
+0.0   4.4765604391852305 xlo xhi
+0.0   4.4765604391852305 ylo yhi
+0.0   113.81190612996821 zlo zhi
+
+Atoms
+
+1 1                    0                    0                   30
+2 1   2.2382802195926152   2.2382802195926152   33.165406242939305
+3 1                    0                    0    36.33081248587861
+4 1   2.2382802195926152   2.2382802195926152   39.496218728817922
+5 1                    0                    0   42.661624971757227
+6 1   2.2382802195926152   2.2382802195926152   45.827031214696532
+7 1                    0                    0   48.992437457635837
+8 1   2.2382802195926152   2.2382802195926152   52.157843700575143
+9 1                    0                    0   55.323249943514455
+10 1   2.2382802195926152   2.2382802195926152   58.488656186453753
+11 1                    0                    0   61.654062429393065
+12 1   2.2382802195926152   2.2382802195926152   64.819468672332363
+13 1                    0                    0   67.984874915271675
+14 1   2.2382802195926152   2.2382802195926152   71.150281158210987
+15 1                    0                    0   74.315687401150285
+16 1   2.2382802195926152   2.2382802195926152   77.481093644089597
+17 1                    0                    0   80.646499887028909
+18 1   2.2382802195926152   2.2382802195926152   83.811906129968207

examples/nwchem/LAMMPS/data.POSCAR_sc_001

@@ -0,0 +1,19 @@
+LAMMPS data file from VASP file VASP/POSCAR_sc_001
+
+9 atoms
+1 atom types
+0.0   2.5123845999742804 xlo xhi
+0.0   2.5123845999742804 ylo yhi
+0.0   50.099076799794247 zlo zhi
+
+Atoms
+
+1 1                    0                    0   15.000000000000002
+2 1                    0                    0   17.512384599974283
+3 1                    0                    0   20.024769199948562
+4 1                    0                    0   22.537153799922841
+5 1                    0                    0   25.049538399897123
+6 1                    0                    0   27.561922999871406
+7 1                    0                    0   30.074307599845682
+8 1                    0                    0   32.586692199819964
+9 1                    0                    0   35.099076799794247

examples/nwchem/VASP/POSCAR_bcc

@@ -0,0 +1,10 @@
+ W 
+ 1.0000000000000000
+     3.1654062429393064    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    3.1654062429393064    0.0000000000000000
+     0.0000000000000000    0.0000000000000000    3.1654062429393064
+   W
+   2
+Cartesian
+  0.0000000000000000  0.0000000000000000  0.0000000000000000
+  1.5827031214696532  1.5827031214696532  1.5827031214696532

examples/nwchem/VASP/POSCAR_bcc_222

@@ -0,0 +1,26 @@
+ W 
+ 1.0000000000000000
+     8.9531208783704628    0.0000000000000000    0.0000000000000000
+     2.2382802195926157    3.8768150619108339    0.0000000000000000
+     0.0000000000000000    0.0000000000000000   45.5341592448720149
+  W
+  18
+Cartesian
+  0.0000000000000000  0.0000000000000000 15.0000000000000000
+  6.7148406587778462  1.2922716873036106 15.9137740732277653
+  4.4765604391852314  2.5845433746072226 16.8275481464555305
+ 11.1914010979630767  3.8768150619108330 17.7413222196832976
+  6.7148406587778471  1.2922716873036111 18.6550962929110611
+  4.4765604391852305  2.5845433746072222 19.5688703661388281
+  2.2382802195926152  3.8768150619108330 20.4826444393665916
+  6.7148406587778453  1.2922716873036111 21.3964185125943587
+  4.4765604391852305  2.5845433746072222 22.3101925858221222
+  2.2382802195926148  3.8768150619108321 23.2239666590498928
+  6.7148406587778453  1.2922716873036102 24.1377407322776563
+  4.4765604391852305  2.5845433746072195 25.0515148055054233
+  2.2382802195926148  3.8768150619108321 25.9652888787331833
+  6.7148406587778444  1.2922716873036084 26.8790629519609539
+  4.4765604391852314  2.5845433746072204 27.7928370251887173
+  2.2382802195926139  3.8768150619108304 28.7066110984164808
+  6.7148406587778453  1.2922716873036102 29.6203851716442479
+  4.4765604391852269  2.5845433746072204 30.5341592448720149

examples/nwchem/VASP/POSCAR_diamond

@@ -0,0 +1,12 @@
+ W 
+ 1.0000000000000000
+     2.8200563642280549    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    2.8200563642280549    0.0000000000000000
+     0.0000000000000000    0.0000000000000000    3.9881619569478763
+   W
+   4
+Cartesian
+  0.0000000000000000  0.0000000000000000  0.0000000000000000
+  1.4100281821140275  0.0000000000000000  0.9970404892369691
+  1.4100281821140275  1.4100281821140275  1.9940809784739382
+  0.0000000000000000  1.4100281821140275  2.9911214677109070

examples/nwchem/VASP/POSCAR_fcc_001

@@ -0,0 +1,26 @@
+ W 
+ 2.0000000000000000
+     2.2382802195926152    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    2.2382802195926152    0.0000000000000000
+     0.0000000000000000    0.0000000000000000   56.9059530649841037
+  W
+  18
+Cartesian
+  0.0000000000000000  0.0000000000000000 15.0000000000000000
+  1.1191401097963076  1.1191401097963076 16.5827031214696525
+  0.0000000000000000  0.0000000000000000 18.1654062429393051
+  1.1191401097963076  1.1191401097963076 19.7481093644089611
+  0.0000000000000000  0.0000000000000000 21.3308124858786137
+  1.1191401097963076  1.1191401097963076 22.9135156073482662
+  0.0000000000000000  0.0000000000000000 24.4962187288179187
+  1.1191401097963076  1.1191401097963076 26.0789218502875713
+  0.0000000000000000  0.0000000000000000 27.6616249717572273
+  1.1191401097963076  1.1191401097963076 29.2443280932268763
+  0.0000000000000000  0.0000000000000000 30.8270312146965324
+  1.1191401097963076  1.1191401097963076 32.4097343361661814
+  0.0000000000000000  0.0000000000000000 33.9924374576358375
+  1.1191401097963076  1.1191401097963076 35.5751405791054935
+  0.0000000000000000  0.0000000000000000 37.1578437005751425
+  1.1191401097963076  1.1191401097963076 38.7405468220447986
+  0.0000000000000000  0.0000000000000000 40.3232499435144547
+  1.1191401097963076  1.1191401097963076 41.9059530649841037

examples/nwchem/VASP/POSCAR_sc_001

@@ -0,0 +1,17 @@
+ W 
+ 1.0000000000000000
+     2.5123845999742804    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    2.5123845999742804    0.0000000000000000
+     0.0000000000000000    0.0000000000000000   50.0990767997942470
+   W
+   9
+Cartesian
+  0.0000000000000000  0.0000000000000000 15.0000000000000018
+  0.0000000000000000  0.0000000000000000 17.5123845999742827
+  0.0000000000000000  0.0000000000000000 20.0247691999485617
+  0.0000000000000000  0.0000000000000000 22.5371537999228408
+  0.0000000000000000  0.0000000000000000 25.0495383998971235
+  0.0000000000000000  0.0000000000000000 27.5619229998714061
+  0.0000000000000000  0.0000000000000000 30.0743075998456817
+  0.0000000000000000  0.0000000000000000 32.5866921998199643
+  0.0000000000000000  0.0000000000000000 35.0990767997942470

examples/nwchem/data.bcc

@@ -0,0 +1,12 @@
+LAMMPS data file from VASP file POSCAR/POSCAR_bcc
+
+2 atoms
+1 atom types
+0.0 3.16540624294 xlo xhi
+0.0 3.16540624294 ylo yhi
+0.0 3.16540624294 zlo zhi
+
+Atoms
+
+1 1 0.0000000000000000 0.0000000000000000 0.0000000000000000
+2 1 1.5827031214696532 1.5827031214696532 1.5827031214696532

examples/nwchem/data.py

@@ -0,0 +1,376 @@
+# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
+# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+#
+# Copyright (2005) 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.
+
+# data tool
+
+oneline = "Read, write, manipulate LAMMPS data files"
+
+docstr = """
+d = data("data.poly")            read a LAMMPS data file, can be gzipped
+d = data()			 create an empty data file
+
+d.map(1,"id",3,"x")              assign names to atom columns (1-N)
+
+coeffs = d.get("Pair Coeffs")    extract info from data file section
+q = d.get("Atoms",4)
+
+  1 arg = all columns returned as 2d array of floats
+  2 args = Nth column returned as vector of floats
+
+d.reorder("Atoms",1,3,2,4,5)     reorder columns (1-N) in a data file section
+
+  1,3,2,4,5 = new order of previous columns, can delete columns this way
+
+d.title = "My LAMMPS data file"	 set title of the data file
+d.headers["atoms"] = 1500        set a header value
+d.sections["Bonds"] = lines      set a section to list of lines (with newlines)
+d.delete("bonds")		 delete a keyword or section of data file
+d.delete("Bonds")
+d.replace("Atoms",5,vec)      	 replace Nth column of section with vector
+d.newxyz(dmp,1000)		 replace xyz in Atoms with xyz of snapshot N
+
+  newxyz assumes id,x,y,z are defined in both data and dump files
+    also replaces ix,iy,iz if they are defined
+  
+index,time,flag = d.iterator(0/1)          loop over single data file snapshot
+time,box,atoms,bonds,tris,lines = d.viz(index)   return list of viz objects
+
+  iterator() and viz() are compatible with equivalent dump calls
+  iterator() called with arg = 0 first time, with arg = 1 on subsequent calls
+    index = timestep index within dump object (only 0 for data file)
+    time = timestep value (only 0 for data file)
+    flag = -1 when iteration is done, 1 otherwise
+  viz() returns info for specified timestep index (must be 0)
+    time = 0
+    box = [xlo,ylo,zlo,xhi,yhi,zhi]
+    atoms = id,type,x,y,z for each atom as 2d array
+    bonds = id,type,x1,y1,z1,x2,y2,z2,t1,t2 for each bond as 2d array
+      NULL if bonds do not exist
+    tris = NULL
+    lines = NULL
+    
+d.write("data.new")             write a LAMMPS data file
+"""
+
+# History
+#   8/05, Steve Plimpton (SNL): original version
+#   11/07, added triclinic box support
+
+# ToDo list
+
+# Variables
+#   title = 1st line of data file
+#   names = dictionary with atom attributes as keys, col #s as values 
+#   headers = dictionary with header name as key, value or tuple as values
+#   sections = dictionary with section name as key, array of lines as values
+#   nselect = 1 = # of snapshots
+
+# Imports and external programs
+
+from os import popen
+
+try: tmp = PIZZA_GUNZIP
+except: PIZZA_GUNZIP = "gunzip"
+
+# Class definition
+
+class data:
+
+  # --------------------------------------------------------------------
+
+  def __init__(self,*list):
+    self.nselect = 1
+    
+    if len(list) == 0:
+      self.title = "LAMMPS data file"
+      self.names = {}
+      self.headers = {}
+      self.sections = {}
+      return
+
+    file = list[0]
+    if file[-3:] == ".gz": f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r')
+    else: f = open(file)
+
+    self.title = f.readline()
+    self.names = {}
+    
+    headers = {}
+    while 1:
+      line = f.readline()
+      line = line.strip()
+      if len(line) == 0:
+        continue
+      found = 0
+      for keyword in hkeywords:
+        if line.find(keyword) >= 0:
+	  found = 1
+	  words = line.split()
+	  if keyword == "xlo xhi" or keyword == "ylo yhi" or \
+	    keyword == "zlo zhi":
+	    headers[keyword] = (float(words[0]),float(words[1]))
+	  elif keyword == "xy xz yz":
+	    headers[keyword] = \
+              (float(words[0]),float(words[1]),float(words[2]))
+          else:
+	    headers[keyword] = int(words[0])
+      if not found:
+        break
+
+    sections = {}
+    while 1:
+      found = 0
+      for pair in skeywords:
+        keyword,length = pair[0],pair[1]
+        if keyword == line:
+	  found = 1
+          if not headers.has_key(length):
+            raise StandardError, \
+                  "data section %s has no matching header value" % line
+	  f.readline()
+          list = []
+          for i in xrange(headers[length]): list.append(f.readline())
+          sections[keyword] = list
+      if not found:
+        raise StandardError,"invalid section %s in data file" % line
+      f.readline()
+      line = f.readline()
+      if not line:
+        break
+      line = line.strip()
+      
+    f.close()
+    self.headers = headers
+    self.sections = sections
+
+  # --------------------------------------------------------------------
+  # assign names to atom columns
+
+  def map(self,*pairs):
+    if len(pairs) % 2 != 0:
+      raise StandardError, "data map() requires pairs of mappings"
+    for i in range(0,len(pairs),2):
+      j = i + 1
+      self.names[pairs[j]] = pairs[i]-1
+
+  # --------------------------------------------------------------------
+  # extract info from data file fields
+
+  def get(self,*list):
+    if len(list) == 1:
+      field = list[0]
+      array = []
+      lines = self.sections[field]
+      for line in lines:
+        words = line.split()
+        values = map(float,words)
+        array.append(values)
+      return array
+    elif len(list) == 2:
+      field = list[0]
+      n = list[1] - 1
+      vec = []
+      lines = self.sections[field]
+      for line in lines:
+        words = line.split()
+        vec.append(float(words[n]))
+      return vec
+    else:
+      raise StandardError, "invalid arguments for data.get()"
+
+  # --------------------------------------------------------------------
+  # reorder columns in a data file field
+
+  def reorder(self,name,*order):
+    n = len(order)
+    natoms = len(self.sections[name])
+    oldlines = self.sections[name]
+    newlines = natoms*[""]
+    for index in order:
+      for i in xrange(len(newlines)):
+        words = oldlines[i].split()
+        newlines[i] += words[index-1] + " "
+    for i in xrange(len(newlines)):
+      newlines[i] += "\n"
+    self.sections[name] = newlines
+
+  # --------------------------------------------------------------------
+  # replace a column of named section with vector of values
+
+  def replace(self,name,icol,vector):
+    lines = self.sections[name]
+    newlines = []
+    j = icol - 1
+    for i in xrange(len(lines)):
+      line = lines[i]
+      words = line.split()
+      words[j] = str(vector[i])
+      newline = ' '.join(words) + '\n'
+      newlines.append(newline)
+    self.sections[name] = newlines
+
+  # --------------------------------------------------------------------
+  # replace x,y,z in Atoms with x,y,z values from snapshot ntime of dump object
+  # assumes id,x,y,z are defined in both data and dump files
+  # also replaces ix,iy,iz if they are defined
+
+  def newxyz(self,dm,ntime):
+    nsnap = dm.findtime(ntime)
+
+    dm.sort(ntime)
+    x,y,z = dm.vecs(ntime,"x","y","z")
+
+    self.replace("Atoms",self.names['x']+1,x)
+    self.replace("Atoms",self.names['y']+1,y)
+    self.replace("Atoms",self.names['z']+1,z)
+    
+    if dm.names.has_key("ix") and self.names.has_key("ix"):
+      ix,iy,iz = dm.vecs(ntime,"ix","iy","iz")
+      self.replace("Atoms",self.names['ix']+1,ix)
+      self.replace("Atoms",self.names['iy']+1,iy)
+      self.replace("Atoms",self.names['iz']+1,iz)
+      
+  # --------------------------------------------------------------------
+  # delete header value or section from data file
+
+  def delete(self,keyword):
+
+    if self.headers.has_key(keyword): del self.headers[keyword]
+    elif self.sections.has_key(keyword): del self.sections[keyword]
+    else: raise StandardError, "keyword not found in data object"
+
+  # --------------------------------------------------------------------
+  # write out a LAMMPS data file
+
+  def write(self,file):
+    f = open(file,"w")
+    print >>f,self.title
+    for keyword in hkeywords:
+      if self.headers.has_key(keyword):
+        if keyword == "xlo xhi" or keyword == "ylo yhi" or \
+               keyword == "zlo zhi":
+	  pair = self.headers[keyword]
+	  print >>f,pair[0],pair[1],keyword
+        elif keyword == "xy xz yz":
+	  triple = self.headers[keyword]
+	  print >>f,triple[0],triple[1],triple[2],keyword
+        else:
+	  print >>f,self.headers[keyword],keyword
+    for pair in skeywords:
+      keyword = pair[0]
+      if self.sections.has_key(keyword):
+        print >>f,"\n%s\n" % keyword
+        for line in self.sections[keyword]:
+	  print >>f,line,
+    f.close()
+
+  # --------------------------------------------------------------------
+  # iterator called from other tools
+
+  def iterator(self,flag):
+    if flag == 0: return 0,0,1
+    return 0,0,-1
+
+  # --------------------------------------------------------------------
+  # time query from other tools
+
+  def findtime(self,n):
+    if n == 0: return 0
+    raise StandardError, "no step %d exists" % (n)
+   
+  # --------------------------------------------------------------------
+  # return list of atoms and bonds to viz for data object
+
+  def viz(self,isnap):
+    if isnap: raise StandardError, "cannot call data.viz() with isnap != 0"
+    
+    id = self.names["id"]
+    type = self.names["type"]
+    x = self.names["x"]
+    y = self.names["y"]
+    z = self.names["z"]
+    
+    xlohi = self.headers["xlo xhi"]
+    ylohi = self.headers["ylo yhi"]
+    zlohi = self.headers["zlo zhi"]
+    box = [xlohi[0],ylohi[0],zlohi[0],xlohi[1],ylohi[1],zlohi[1]]
+
+    # create atom list needed by viz from id,type,x,y,z
+
+    atoms = []
+    atomlines = self.sections["Atoms"]
+    for line in atomlines:
+      words = line.split()
+      atoms.append([int(words[id]),int(words[type]),
+                    float(words[x]),float(words[y]),float(words[z])])
+
+    # create list of current bond coords from list of bonds
+    # assumes atoms are sorted so can lookup up the 2 atoms in each bond
+
+    bonds = []
+    if self.sections.has_key("Bonds"):
+      bondlines = self.sections["Bonds"]
+      for line in bondlines:
+        words = line.split()
+        bid,btype   = int(words[0]),int(words[1])
+        atom1,atom2 = int(words[2]),int(words[3])
+        atom1words  = atomlines[atom1-1].split()
+        atom2words  = atomlines[atom2-1].split()
+        bonds.append([bid,btype,
+                      float(atom1words[x]),float(atom1words[y]),
+                      float(atom1words[z]),
+                      float(atom2words[x]),float(atom2words[y]),
+                      float(atom2words[z]),
+		      float(atom1words[type]),float(atom2words[type])])
+      
+    tris = []
+    lines = []
+    return 0,box,atoms,bonds,tris,lines
+
+  # --------------------------------------------------------------------
+  # return box size
+
+  def maxbox(self):
+    xlohi = self.headers["xlo xhi"]
+    ylohi = self.headers["ylo yhi"]
+    zlohi = self.headers["zlo zhi"]
+    return [xlohi[0],ylohi[0],zlohi[0],xlohi[1],ylohi[1],zlohi[1]]
+
+  # --------------------------------------------------------------------
+  # return number of atom types
+
+  def maxtype(self):
+    return self.headers["atom types"]
+
+# --------------------------------------------------------------------
+# data file keywords, both header and main sections
+
+hkeywords = ["atoms","ellipsoids","lines","triangles","bodies",
+             "bonds","angles","dihedrals","impropers",
+	     "atom types","bond types","angle types","dihedral types",
+	     "improper types","xlo xhi","ylo yhi","zlo zhi","xy xz yz"]
+
+skeywords = [["Masses","atom types"],
+             ["Atoms","atoms"],["Ellipsoids","ellipsoids"],
+             ["Lines","lines"],["Triangles","triangles"],["Bodies","bodies"],
+             ["Bonds","bonds"],
+	     ["Angles","angles"],["Dihedrals","dihedrals"],
+	     ["Impropers","impropers"],["Velocities","atoms"],
+             ["Pair Coeffs","atom types"],
+	     ["Bond Coeffs","bond types"],["Angle Coeffs","angle types"],
+	     ["Dihedral Coeffs","dihedral types"],
+	     ["Improper Coeffs","improper types"],
+             ["BondBond Coeffs","angle types"],
+             ["BondAngle Coeffs","angle types"],
+             ["MiddleBondTorsion Coeffs","dihedral types"],
+             ["EndBondTorsion Coeffs","dihedral types"],
+             ["AngleTorsion Coeffs","dihedral types"],
+             ["AngleAngleTorsion Coeffs","dihedral types"],
+             ["BondBond13 Coeffs","dihedral types"],
+             ["AngleAngle Coeffs","improper types"],
+             ["Molecules","atoms"]]

examples/nwchem/data.water.dimer

@@ -0,0 +1,45 @@
+LAMMPS data file for water dimer in large box
+
+6 atoms
+4 bonds
+2 angles
+2 atom types
+1 bond types
+1 angle types
+-6.879301 6.879301 xlo xhi
+-6.879301 6.879301 ylo yhi
+-6.879301 6.879301 zlo zhi
+
+Masses
+
+1 15.99491
+2 1.008
+
+Bond Coeffs
+
+1 554.25 1.0
+
+Angle Coeffs
+
+1 47.744 109.4
+
+Atoms
+
+1 1 1 -0.8476   0.161560  -0.052912   0.033173
+2 1 2  0.4238   0.803054   0.369132  -0.511660
+3 1 2  0.4238  -0.325571  -0.669574  -0.488560
+4 2 1 -0.8476   0.021259   0.506771   2.831278
+5 2 2  0.4238  -0.721039   1.083100   2.758378
+6 2 2  0.4238   0.158220   0.181883   1.945696
+
+Bonds
+
+1 1 1 2
+2 1 1 3
+3 1 4 5
+4 1 4 6
+
+Angles
+
+1 1 2 1 3
+2 1 5 4 6

examples/nwchem/data.water.dimer.mm

@@ -0,0 +1,45 @@
+LAMMPS data file for water dimer in large box
+
+6 atoms
+4 bonds
+2 angles
+2 atom types
+1 bond types
+1 angle types
+-6.879301 6.879301 xlo xhi
+-6.879301 6.879301 ylo yhi
+-6.879301 6.879301 zlo zhi
+
+Masses
+
+1 15.99491
+2 1.008
+
+Bond Coeffs
+
+1 554.25 1.0
+
+Angle Coeffs
+
+1 47.744 109.4
+
+Atoms
+
+1 1 1 -0.8476   0.161560  -0.052912   0.033173
+2 1 2  0.4238   0.803054   0.369132  -0.511660
+3 1 2  0.4238  -0.325571  -0.669574  -0.488560
+4 2 1 -0.8476   0.021259   0.506771   2.831278
+5 2 2  0.4238  -0.721039   1.083100   2.758378
+6 2 2  0.4238   0.158220   0.181883   1.945696
+
+Bonds
+
+1 1 1 2
+2 1 1 3
+3 1 4 5
+4 1 4 6
+
+Angles
+
+1 1 2 1 3
+2 1 5 4 6

examples/nwchem/data.water.dimer.qmmm

@@ -0,0 +1,47 @@
+LAMMPS data file for water dimer in large box - for QMMM with different types
+
+6 atoms
+4 bonds
+2 angles
+4 atom types
+1 bond types
+1 angle types
+-6.879301 6.879301 xlo xhi
+-6.879301 6.879301 ylo yhi
+-6.879301 6.879301 zlo zhi
+
+Masses
+
+1 15.99491
+2 1.008
+3 15.99491
+4 1.008
+
+Bond Coeffs
+
+1 554.25 1.0
+
+Angle Coeffs
+
+1 47.744 109.4
+
+Atoms
+
+1 1 1 -0.8476   0.161560  -0.052912   0.033173
+2 1 2  0.4238   0.803054   0.369132  -0.511660
+3 1 2  0.4238  -0.325571  -0.669574  -0.488560
+4 2 3 -0.8476   0.021259   0.506771   2.831278
+5 2 4  0.4238  -0.721039   1.083100   2.758378
+6 2 4  0.4238   0.158220   0.181883   1.945696
+
+Bonds
+
+1 1 1 2
+2 1 1 3
+3 1 4 5
+4 1 4 6
+
+Angles
+
+1 1 2 1 3
+2 1 5 4 6

examples/nwchem/data.zeolite

@@ -0,0 +1,114 @@
+LAMMPS data file for SiO2 zeolite with one methane moleclue
+
+77 atoms
+4 atom types
+4 bonds
+1 bond types
+6 angles
+1 angle types
+-5.9266 5.9926 xlo xhi
+-5.9266 5.9926 ylo yhi
+-5.9266 5.9926 zlo zhi
+
+Masses
+
+1  28.0855
+2  15.99491
+3  12.0
+4  1.008
+
+Atoms
+
+1 0 1    1.910418    0.00000    4.38651    2.18123
+2 0 1    1.910418    0.00000   -4.38651    2.18123
+3 0 1    1.910418    0.00000    4.38651   -2.18123
+4 0 1    1.910418    0.00000   -4.38651   -2.18123
+5 0 1    1.910418    2.18123    0.00000    4.38651
+6 0 1    1.910418    2.18123    0.00000   -4.38651
+7 0 1    1.910418   -2.18123    0.00000    4.38651
+8 0 1    1.910418   -2.18123    0.00000   -4.38651
+9 0 1    1.910418    4.38651    2.18123    0.00000
+10 0 1   1.910418   -4.38651    2.18123    0.00000
+11 0 1   1.910418    4.38651   -2.18123    0.00000
+12 0 1   1.910418   -4.38651   -2.18123    0.00000
+13 0 1   1.910418    4.38651    0.00000   -2.18123
+14 0 1   1.910418   -4.38651    0.00000   -2.18123
+15 0 1   1.910418    4.38651    0.00000    2.18123
+16 0 1   1.910418   -4.38651    0.00000    2.18123
+17 0 1   1.910418    0.00000    2.18123   -4.38651
+18 0 1   1.910418    0.00000    2.18123    4.38651
+19 0 1   1.910418    0.00000   -2.18123   -4.38651
+20 0 1   1.910418    0.00000   -2.18123    4.38651
+21 0 1   1.910418    2.18123    4.38651    0.00000
+22 0 1   1.910418    2.18123   -4.38651    0.00000
+23 0 1   1.910418   -2.18123    4.38651    0.00000
+24 0 1   1.910418   -2.18123   -4.38651    0.00000
+25 0 2   -0.955209    0.00000   -5.92660    2.64860
+26 0 2   -0.955209    0.00000   -5.92660   -2.64860
+27 0 2   -0.955209    2.64860    0.00000   -5.92660
+28 0 2   -0.955209   -2.64860    0.00000   -5.92660
+29 0 2   -0.955209   -5.92660    2.64860    0.00000
+30 0 2   -0.955209   -5.92660   -2.64860    0.00000
+31 0 2   -0.955209   -5.92660    0.00000   -2.64860
+32 0 2   -0.955209   -5.92660    0.00000    2.64860
+33 0 2   -0.955209    0.00000    2.64860   -5.92660
+34 0 2   -0.955209    0.00000   -2.64860   -5.92660
+35 0 2   -0.955209    2.64860   -5.92660    0.00000
+36 0 2   -0.955209   -2.64860   -5.92660    0.00000
+37 0 2   -0.955209    0.00000    3.45272    3.45272
+38 0 2   -0.955209    0.00000   -3.45272    3.45272
+39 0 2   -0.955209    0.00000    3.45272   -3.45272
+40 0 2   -0.955209    0.00000   -3.45272   -3.45272
+41 0 2   -0.955209    3.45272    0.00000    3.45272
+42 0 2   -0.955209    3.45272    0.00000   -3.45272
+43 0 2   -0.955209   -3.45272    0.00000    3.45272
+44 0 2   -0.955209   -3.45272    0.00000   -3.45272
+45 0 2   -0.955209    3.45272    3.45272    0.00000
+46 0 2   -0.955209   -3.45272    3.45272    0.00000
+47 0 2   -0.955209    3.45272   -3.45272    0.00000
+48 0 2   -0.955209   -3.45272   -3.45272    0.00000
+49 0 2   -0.955209    1.28702    1.28702    4.12598
+50 0 2   -0.955209   -1.28702   -1.28702    4.12598
+51 0 2   -0.955209   -1.28702    1.28702   -4.12598
+52 0 2   -0.955209    1.28702   -1.28702   -4.12598
+53 0 2   -0.955209    4.12598    1.28702    1.28702
+54 0 2   -0.955209    4.12598   -1.28702   -1.28702
+55 0 2   -0.955209   -4.12598   -1.28702    1.28702
+56 0 2   -0.955209   -4.12598    1.28702   -1.28702
+57 0 2   -0.955209    1.28702    4.12598    1.28702
+58 0 2   -0.955209   -1.28702    4.12598   -1.28702
+59 0 2   -0.955209    1.28702   -4.12598   -1.28702
+60 0 2   -0.955209   -1.28702   -4.12598    1.28702
+61 0 2   -0.955209    1.28702    1.28702   -4.12598
+62 0 2   -0.955209   -1.28702   -1.28702   -4.12598
+63 0 2   -0.955209    1.28702   -1.28702    4.12598
+64 0 2   -0.955209   -1.28702    1.28702    4.12598
+65 0 2   -0.955209    1.28702    4.12598   -1.28702
+66 0 2   -0.955209   -1.28702    4.12598    1.28702
+67 0 2   -0.955209   -1.28702   -4.12598   -1.28702
+68 0 2   -0.955209    1.28702   -4.12598    1.28702
+69 0 2   -0.955209    4.12598    1.28702   -1.28702
+70 0 2   -0.955209    4.12598   -1.28702    1.28702
+71 0 2   -0.955209   -4.12598    1.28702    1.28702
+72 0 2   -0.955209   -4.12598   -1.28702   -1.28702
+73 1 3  -0.66        0.00000    0.00000    0.00000
+74 1 4   0.165       0.00000   -0.89000   -0.62930
+75 1 4   0.165       0.00000    0.89000   -0.62930
+76 1 4   0.165      -0.89000    0.00000    0.62930
+77 1 4   0.165       0.89000    0.00000    0.62930
+
+Bonds
+
+1 1 73 74
+2 1 73 75
+3 1 73 76
+4 1 73 77
+
+Angles
+
+1 1 74 73 75
+2 1 74 73 76
+3 1 74 73 77
+4 1 75 73 76
+5 1 75 73 77
+6 1 76 73 77

examples/nwchem/file.list

@@ -0,0 +1,5 @@
+data.POSCAR_bcc
+data.POSCAR_bcc_222
+data.POSCAR_diamond
+data.POSCAR_fcc_001
+data.POSCAR_sc_001

examples/nwchem/in.w.aimd

@@ -0,0 +1,26 @@
+# AIMD with NWChem
+
+units		metal
+atom_style	atomic
+atom_modify     map yes #sort 1000 0.5
+comm_modify     cutoff 2.0
+
+read_data       data.bcc
+
+mass		1 183.84
+velocity        all create 500.0 394893
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+
+fix             2 all nwchem template-w.nw w.aimd.nw log.pwdft.w.aimd W
+fix_modify      2 energy yes
+
+dump            1 all custom 1 dump.w.aimd id x y z vx vy vz fx fy fz
+
+timestep        0.001
+
+thermo          1
+run		10

examples/nwchem/in.w.many

@@ -0,0 +1,34 @@
+# many conformations with NWChem
+
+variable        datafiles file file.list
+
+label           loop
+
+clear
+
+units		metal
+atom_style	atomic
+atom_modify     map yes
+comm_modify     cutoff 2.0
+
+log             log.${datafiles}
+
+read_data       LAMMPS/${datafiles}
+
+mass		1 183.84
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix             2 all nwchem template-w.nw ${datafiles}.nw &
+                log.pwdft.${datafiles} W
+fix_modify      2 energy yes
+
+dump            1 all custom 1 dump.${datafiles} id x y z fx fy fz
+
+timestep        0.001
+
+run		0
+
+next            datafiles
+jump            SELF loop

examples/nwchem/in.w.single

@@ -0,0 +1,22 @@
+# single conformation with NWChem
+
+units		metal
+atom_style	atomic
+atom_modify     map yes
+comm_modify     cutoff 2.0
+
+read_data       data.bcc
+
+mass		1 183.84
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix             2 all nwchem template-w.nw w.single.nw log.pwdft.w.single W
+fix_modify      2 energy yes
+
+dump            1 all custom 1 dump.w.single id x y z fx fy fz
+
+timestep        0.001
+
+run		0

examples/nwchem/in.water.dimer.mm

@@ -0,0 +1,46 @@
+# MM for water dimer
+
+units		real
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.water.dimer.mm
+
+group           mm molecule 1
+group           qm molecule 2
+
+# pair style must define stand-alone short-range Coulombics
+
+pair_style      lj/cut/coul/cut 6.0
+pair_coeff      1 1 0.13506 3.166         
+pair_coeff      2 2 0.0 1.0         
+
+#velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+
+compute         1 all pair/local dist
+compute         2 all reduce max c_1
+
+variable        fxabs atom abs(fx)
+variable        fyabs atom abs(fy)
+variable        fzabs atom abs(fz)
+variable        qabs atom abs(q)
+compute         3 all reduce max v_fxabs v_fyabs v_fzabs v_qabs
+
+dump            1 all custom 1 dump.water.dimer.mm id x y z q fx fy fz
+dump_modify     1 sort id format float "%20.16g"
+
+timestep        1.0
+
+thermo_style    custom step cpu temp ke evdwl ecoul epair emol elong &
+                pe etotal press c_2 c_3[*]
+
+thermo          1
+
+run		10

examples/nwchem/in.water.dimer.qmmm

@@ -0,0 +1,69 @@
+# QMMM with NWChem for water dimer
+
+units		real
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.water.dimer.qmmm
+
+group           qm molecule 1
+group           mm molecule 2
+
+# remove bonds/angles in QM water molecule
+
+delete_bonds    qm multi remove special
+
+# pair style must define stand-alone short-range Coulombics
+# must specify mixing explicitly b/c hybrid/overlay
+# QM O,H = types 1,2
+# MM O,H = types 3,4
+# QM O,H atoms do not LJ interact with each other
+# only MM O atoms LJ interact with other b/c MM H is zero
+# MM/QM O do LJ interact with each other, same as pair of MM O atoms
+# MM O and QM H do LJ interact with each other with non-zero H epsilon = 0.044
+#   geometric mixing for epsilon, arithmetic for sigma
+#   this is to provide stability for QM H atoms
+
+# mixing only for MM-O/QM-O and MM-O/QM-H
+
+pair_style      hybrid/overlay lj/cut 6.0 coul/cut 6.0
+pair_coeff      1 1 lj/cut 0.0 3.165558         
+pair_coeff      2 2 lj/cut 0.0 0.7   
+pair_coeff      3 3 lj/cut 0.155394 3.165558         
+pair_coeff      4 4 lj/cut 0.0 0.7         
+pair_coeff      1 3 lj/cut 0.155394 3.165558         
+pair_coeff      2 3 lj/cut 0.08268818537130924 1.932779
+pair_coeff      * * coul/cut
+
+#velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+
+fix             2 qm nwchem template.water.nw water.dimer.nw &
+                log.pwdft.water.dimer O H O H
+fix_modify      2 energy yes
+
+# convert dump file forces to Hartree/Bohr for comparison to NWChem
+
+variable        fx atom fx/1185.8
+variable        fy atom fy/1185.8
+variable        fz atom fz/1185.8
+
+dump            1 all custom 1 dump.water.dimer.qmmm id x y z q v_fx v_fy v_fz
+dump_modify     1 sort id format float "%20.16g"
+
+# small timestep for QMMM
+
+timestep        0.1
+
+thermo_style    custom step cpu temp ke evdwl ecoul epair emol elong &
+                f_2 pe etotal press
+
+thermo          1
+
+run		20

examples/nwchem/in.zeolite.mm

@@ -0,0 +1,116 @@
+# MM for SiO2 zeolite with one methane molecule
+
+# CHIK potential
+# EPL, Carre, Horbach, Ispas, Kob, 82, 17001 (2008)
+# B = 1/rho
+
+#q Si = 1.910418
+#q O = -0.955209
+#A OO = 659.595398 eV
+#B OO = 2.590066 1/Ang
+#C OO = 26.836679 eV-Ang^6
+#A SiO = 27029.419922 eV
+#B SiO = 5.158606 1/Ang
+#C SiO = 148.099091 eV-Ang^6
+#A SiSi = 3150.462646 eV
+#B SiSi = 2.851451 1/Ang
+#C SiSi = 626.7519553 eV-Ang^6
+
+# LJ params for methane and O from Table 1
+# Bhatia and Nicholson, J Phys Chem C, 2012, 116, 2344-2355.
+
+#q C = -0.66
+#Q H = 0.165
+#sigma C = 0.34 nm
+#sigma H = 0.265 nm
+#sigma O = 0.28 nm
+#eps/kB C = 55.082 K = 0.004745993 eV
+#eps/kB H = 7.905 K = 0.000681113 eV
+#eps/kB O = 492.7 K = 0.0424522 eV
+
+# LJ params for silicon
+#e-Journal of Surf Sci and Nanotech, Inui and Iwasaki, 15, 40-49 (2017)
+
+#sigma Si = 3.826 Ang
+#eps Si = 17.4 meV = 0.0174 eV
+
+# C-H bond and methane angle params
+
+#OPLS C-H bond k = 29.40 ev/Ang^2
+#C-H bond r0 = 1.09 Angs
+#methane angles = 109.5 degrees
+#C-H angle k/kB = 2000 K/rad^2
+
+# conversions
+
+#1 eV = 11606 K
+#1 eV = 23.0609 kcal/mole
+#1 kcal/mole = 503.2761 K
+#1 kcal = 4.814 kJoule
+
+# -------------------------
+
+units		metal
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.zeolite
+
+group           mm type 1 2
+group           qm type 3 4
+
+# pair style must define stand-alone short-range Coulombics
+# arithmetic mixing
+
+pair_style      hybrid/overlay buck 6.5 lj/cut 6.5 coul/cut 6.5
+
+pair_coeff      1 1 buck 3150.462646 0.35032282 626.7519553
+pair_coeff      2 2 buck 659.595398 0.38609055 26.836679
+pair_coeff      1 2 buck 27029.419922 0.19385082 148.099091
+pair_coeff      1 2 buck 27029.419922 0.19385082 148.099091
+pair_coeff      1 3 lj/cut 0.009087 3.613
+pair_coeff      1 4 lj/cut 0.00344258 3.238
+pair_coeff      2 3 lj/cut 0.01419429 3.1
+pair_coeff      2 4 lj/cut 0.00537724 2.725
+pair_coeff      3 3 lj/cut 0.004746 3.4
+pair_coeff      4 4 lj/cut 0.00068111 2.65
+pair_coeff      * * coul/cut
+
+bond_style      harmonic
+bond_coeff      1 29.40 1.09
+
+angle_style      harmonic
+angle_coeff      1 0.172325 109.5
+
+#velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+# dynamic or frozen zeolite
+
+fix		1 all nve
+#fix		1 qm nve
+
+compute         1 all pair/local dist
+compute         2 all reduce max c_1
+
+variable        fxabs atom abs(fx)
+variable        fyabs atom abs(fy)
+variable        fzabs atom abs(fz)
+variable        qabs atom abs(q)
+compute         3 all reduce max v_fxabs v_fyabs v_fzabs v_qabs
+
+dump            1 all custom 1 dump.zeolite.mm id x y z q fx fy fz
+dump_modify     1 sort id format float "%20.16g"
+
+timestep        0.001
+
+thermo_style    custom step cpu temp ke evdwl ecoul epair emol elong &
+                pe etotal press c_2 c_3[*]
+
+thermo          1
+
+run		3

examples/nwchem/in.zeolite.qmmm

@@ -0,0 +1,83 @@
+# QMMM for SiO2 zeolite with one methane molecule
+
+units		metal
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.zeolite
+
+group           mm type 1 2
+group           qm type 3 4
+
+# pair style must define stand-alone short-range Coulombics
+# must specify mixing explicitly b/c hybrid/overlay
+# MM Si,O = types 1,2
+# QM C,H = types 3,4
+# MM Si,O atoms do not LJ interact with each other (just via Buckingham)
+# QM C,H atoms do not LJ interact with each other
+# MM Si,O and QM C,H do LJ interact with each other
+
+pair_style      hybrid/overlay buck 6.5 lj/cut 6.5 coul/cut 6.5
+
+pair_coeff      1 1 buck 3150.462646 0.35032282 626.7519553
+pair_coeff      2 2 buck 659.595398 0.38609055 26.836679
+pair_coeff      1 2 buck 27029.419922 0.19385082 148.099091
+
+#pair_coeff      1 3 lj/cut 0.009087 3.613
+#pair_coeff      1 4 lj/cut 0.00344258 3.238
+#pair_coeff      2 3 lj/cut 0.01419429 3.1
+#pair_coeff      2 4 lj/cut 0.0035857762359063315 1.932779  # same as water dimer
+pair_coeff      1 3 lj/cut 0.09087 3.613
+pair_coeff      1 4 lj/cut 0.0344258 3.238
+pair_coeff      2 3 lj/cut 0.1419429 3.1
+pair_coeff      2 4 lj/cut 0.035857762359063315 1.932779  # same as water dimer
+
+pair_coeff      3 3 lj/cut 0.0 3.4
+pair_coeff      4 4 lj/cut 0.0 2.65
+pair_coeff      * * coul/cut
+
+bond_style      harmonic
+bond_coeff      1 29.40 1.09
+
+angle_style     harmonic
+angle_coeff     1 0.172325 109.5
+
+# remove bonds/angles in QM methane molecule
+
+delete_bonds    qm multi remove special
+
+velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+# dynamic or frozen zeolite
+
+#fix		1 all nve
+fix		1 qm nve
+
+fix             2 qm nwchem template.methane.nw methane.nw &
+                log.pwdft.zeolite Si O C H
+fix_modify      2 energy yes
+
+# convert dump file forces to Hartree/Bohr for comparison to NWChem
+
+variable        fx atom fx/1185.8
+variable        fy atom fy/1185.8
+variable        fz atom fz/1185.8
+
+dump            1 all custom 1 dump.zeolite.qmmm id x y z q v_fx v_fy v_fz
+dump_modify     1 sort id format float "%20.16g"
+
+# small timestep for QMMM
+
+timestep        0.0001
+
+thermo_style    custom step cpu temp ke evdwl ecoul epair emol elong &
+                f_2 pe etotal press
+
+thermo          1
+
+run		20

examples/nwchem/log.lammps.zeolite.qmmm

@@ -0,0 +1,234 @@
+LAMMPS (4 May 2022)
+# QMMM for SiO2 zeolite with one methane molecule
+
+units		metal
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.zeolite
+Reading data file ...
+  orthogonal box = (-5.9266 -5.9266 -5.9266) to (5.9926 5.9926 5.9926)
+  2 by 4 by 4 MPI processor grid
+  reading atoms ...
+  77 atoms
+  scanning bonds ...
+  4 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  reading bonds ...
+  4 bonds
+  reading angles ...
+  6 angles
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     4 = max # of 1-2 neighbors
+     3 = max # of 1-3 neighbors
+     3 = max # of 1-4 neighbors
+     4 = max # of special neighbors
+  special bonds CPU = 0.002 seconds
+  read_data CPU = 0.021 seconds
+
+group           mm type 1 2
+72 atoms in group mm
+group           qm type 3 4
+5 atoms in group qm
+
+# pair style must define stand-alone short-range Coulombics
+# must specify mixing explicitly b/c hybrid/overlay
+# MM Si,O = types 1,2
+# QM C,H = types 3,4
+# MM Si,O atoms do not LJ interact with each other (just via Buckingham)
+# QM C,H atoms do not LJ interact with each other
+# MM Si,O and QM C,H do LJ interact with each other
+
+pair_style      hybrid/overlay buck 6.5 lj/cut 6.5 coul/cut 6.5
+
+pair_coeff      1 1 buck 3150.462646 0.35032282 626.7519553
+pair_coeff      2 2 buck 659.595398 0.38609055 26.836679
+pair_coeff      1 2 buck 27029.419922 0.19385082 148.099091
+pair_coeff      1 3 lj/cut 0.009087 3.613
+pair_coeff      1 4 lj/cut 0.00344258 3.238
+pair_coeff      2 3 lj/cut 0.01419429 3.1
+pair_coeff      2 4 lj/cut 0.0035857762359063315 1.932779  # same as water dimer
+pair_coeff      3 3 lj/cut 0.0 3.4
+pair_coeff      4 4 lj/cut 0.0 2.65
+pair_coeff      * * coul/cut
+
+bond_style      harmonic
+bond_coeff      1 29.40 1.09
+
+angle_style     harmonic
+angle_coeff     1 0.172325 109.5
+
+# remove bonds/angles in QM methane molecule
+
+delete_bonds    qm multi remove special
+System init for delete_bonds ...
+Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.5
+  ghost atom cutoff = 8.5
+  binsize = 4.25, bins = 3 3 3
+  3 neighbor lists, perpetual/occasional/extra = 3 0 0
+  (1) pair buck, perpetual, skip from (3)
+      attributes: half, newton on
+      pair build: skip
+      stencil: none
+      bin: none
+  (2) pair lj/cut, perpetual, skip from (3)
+      attributes: half, newton on
+      pair build: skip
+      stencil: none
+      bin: none
+  (3) pair coul/cut, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d
+      bin: standard
+Deleting bonds ...
+  0 total bonds, 0 turned on, 0 turned off
+  0 total angles, 0 turned on, 0 turned off
+  0 total dihedrals, 0 turned on, 0 turned off
+  0 total impropers, 0 turned on, 0 turned off
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     0 = max # of 1-2 neighbors
+     0 = max # of 1-3 neighbors
+     0 = max # of 1-4 neighbors
+     4 = max # of special neighbors
+  special bonds CPU = 0.001 seconds
+
+velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+# dynamic or frozen zeolite
+
+#fix		1 all nve
+fix		1 qm nve
+
+fix             2 qm nwchem template.methane.nw methane.nw                 log.pwdft.zeolite Si O C H
+fix_modify      2 energy yes
+
+# convert dump file forces to Hartree/Bohr for comparison to NWChem
+
+variable        fx atom fx/1185.8
+variable        fy atom fy/1185.8
+variable        fz atom fz/1185.8
+
+dump            1 all custom 1 dump.zeolite.qmmm id x y z q v_fx v_fy v_fz
+dump_modify     1 sort id format float "%20.16g"
+
+# small timestep for QMMM
+
+timestep        0.0001
+
+thermo_style    custom step cpu temp ke evdwl ecoul epair emol elong                 f_2 pe etotal press
+
+thermo          1
+
+run		20
+Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule
+Calling pspw_input() in NWChem ...
+  time = 4.968 seconds
+Calling pspw_qmmm() in NWChem ...
+  time = 6.503 seconds
+Per MPI rank memory allocation (min/avg/max) = 11.26 | 11.26 | 11.26 Mbytes
+   Step          CPU            Temp          KinEng         E_vdwl         E_coul         E_pair         E_mol          E_long          f_2           PotEng         TotEng         Press     
+         0   0              300            2.9471313     -102.3998      -1065.0369     -1167.4367      0              0             -198.9818      -1366.4185     -1363.4714      5517.5432    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.868 seconds
+         1   3.8686158      299.95086      2.9466486     -102.39975     -1065.1059     -1167.5057      0              0             -198.90989     -1366.4156     -1363.4689      5495.4255    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.839 seconds
+         2   7.7086679      299.89221      2.9460724     -102.39971     -1065.1824     -1167.5821      0              0             -198.83094     -1366.413      -1363.4669      5470.9047    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.828 seconds
+         3   11.53711       299.8247       2.9454092     -102.39967     -1065.2698     -1167.6695      0              0             -198.74163     -1366.4111     -1363.4657      5442.8579    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.874 seconds
+         4   15.411604      299.74946      2.94467       -102.39963     -1065.3685     -1167.7681      0              0             -198.64174     -1366.4099     -1363.4652      5411.2193    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.867 seconds
+         5   19.279369      299.66809      2.9438707     -102.39959     -1065.4787     -1167.8782      0              0             -198.53108     -1366.4093     -1363.4655      5375.9224    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.884 seconds
+         6   23.164177      299.58269      2.9430317     -102.39955     -1065.6006     -1168.0001      0              0             -198.4094      -1366.4095     -1363.4665      5336.8881    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.931 seconds
+         7   27.095936      299.49578      2.9421779     -102.39951     -1065.7347     -1168.1342      0              0             -198.27641     -1366.4106     -1363.4684      5294.0245    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.995 seconds
+         8   31.091727      299.41036      2.9413388     -102.39947     -1065.8812     -1168.2807      0              0             -198.13177     -1366.4125     -1363.4711      5247.225     
+Calling pspw_qmmm() in NWChem ...
+  time = 3.917 seconds
+         9   35.008957      299.32983      2.9405477     -102.39943     -1066.0408     -1168.4402      0              0             -197.9751      -1366.4153     -1363.4748      5196.3675    
+Calling pspw_qmmm() in NWChem ...
+  time = 3.938 seconds
+        10   38.948036      299.25799      2.9398419     -102.39939     -1066.2138     -1168.6132      0              0             -197.80594     -1366.4191     -1363.4793      5141.3097    
+Calling pspw_qmmm() in NWChem ...
+  time = 5.938 seconds
+        11   44.886417      300.13711      2.9484783     -102.3987      -1112.1634     -1214.5621      0              0             -163.28337     -1377.8455     -1374.897      -9344.0882    
+Calling pspw_qmmm() in NWChem ...
+  time = 4.213 seconds
+        12   49.100354      302.20838      2.9688259     -102.39866     -1113.0737     -1215.4724      0              0             -162.44687     -1377.9192     -1374.9504     -9618.386     
+Calling pspw_qmmm() in NWChem ...
+  time = 4.522 seconds
+        13   53.622594      304.73496      2.9936465     -102.39862     -1114.1542     -1216.5529      0              0             -161.45551     -1378.0084     -1375.0147     -9943.5638    
+Calling pspw_qmmm() in NWChem ...
+  time = 4.301 seconds
+        14   57.924103      307.82384      3.0239909     -102.39858     -1115.4521     -1217.8506      0              0             -160.26631     -1378.1169     -1375.093      -10333.782    
+Calling pspw_qmmm() in NWChem ...
+  time = 4.587 seconds
+        15   62.511538      311.60627      3.0611487     -102.39853     -1117.0209     -1219.4194      0              0             -158.83006     -1378.2495     -1375.1883     -10805.182    
+Calling pspw_qmmm() in NWChem ...
+  time = 6.373 seconds
+        16   68.884983      321.46762      3.1580243     -102.39783     -1190.4155     -1292.8133      0              0             -103.41532     -1396.2286     -1393.0706     -33890.787    
+Calling pspw_qmmm() in NWChem ...
+  time = 7.227 seconds
+        17   76.112329      354.58063      3.4833189     -102.39711     -1297.586      -1399.9831      0              0             -21.400411     -1421.3835     -1417.9002     -67484.92     
+Calling pspw_qmmm() in NWChem ...
+  time = 5.903 seconds
+        18   82.015607      426.41585      4.1890116     -102.39703     -1319.896      -1422.293       0              0             -1.3910438     -1423.6841     -1419.4951     -74076.209    
+Calling pspw_qmmm() in NWChem ...
+  time = 6.612 seconds
+        19   88.628069      543.11841      5.3354708     -102.39692     -1357.8053     -1460.2022      0              0              32.790891     -1427.4114     -1422.0759     -85309.358    
+Calling pspw_qmmm() in NWChem ...
+  time = 7.909 seconds
+        20   96.537191      722.72131      7.0998486     -102.39678     -1420.4861     -1522.8829      0              0              89.745396     -1433.1375     -1426.0376     -103965.42    
+Loop time of 96.5372 on 32 procs for 20 steps with 77 atoms
+
+Performance: 0.002 ns/day, 13407.950 hours/ns, 0.207 timesteps/s
+99.8% CPU use with 32 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.00030032 | 0.00063124 | 0.00098374 |   0.0 |  0.00
+Bond    | 4.4281e-05 | 7.7957e-05 | 0.00012963 |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0.0016224  | 0.0018103  | 0.0021141  |   0.3 |  0.00
+Output  | 0.0050295  | 0.0052005  | 0.0057519  |   0.2 |  0.01
+Modify  | 96.528     | 96.528     | 96.529     |   0.0 | 99.99
+Other   |            | 0.001124   |            |       |  0.00
+
+Nlocal:        2.40625 ave           7 max           0 min
+Histogram: 4 10 3 0 6 6 0 0 2 1
+Nghost:        313.344 ave         405 max         225 min
+Histogram: 7 1 0 0 12 4 0 0 0 8
+Neighs:        98.3125 ave         324 max           0 min
+Histogram: 5 10 2 5 5 2 0 0 2 1
+
+Total # of neighbors = 3146
+Ave neighs/atom = 40.857143
+Ave special neighs/atom = 0
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:01:48

examples/nwchem/log.lammps.zeolite.qmmm.1

@@ -0,0 +1,133 @@
+LAMMPS (4 May 2022)
+# QMMM for SiO2 zeolite with one methane molecule
+
+units		metal
+atom_style	full
+
+bond_style      harmonic
+angle_style     harmonic
+
+read_data       data.zeolite
+Reading data file ...
+  orthogonal box = (-5.9266 -5.9266 -5.9266) to (5.9926 5.9926 5.9926)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  77 atoms
+  scanning bonds ...
+  4 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  reading bonds ...
+  4 bonds
+  reading angles ...
+  6 angles
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     4 = max # of 1-2 neighbors
+     3 = max # of 1-3 neighbors
+     3 = max # of 1-4 neighbors
+     4 = max # of special neighbors
+  special bonds CPU = 0.000 seconds
+  read_data CPU = 0.007 seconds
+
+group           mm type 1 2
+72 atoms in group mm
+group           qm type 3 4
+5 atoms in group qm
+
+# pair style must define stand-alone short-range Coulombics
+# must specify mixing explicitly b/c hybrid/overlay
+# MM Si,O = types 1,2
+# QM C,H = types 3,4
+# MM Si,O atoms do not LJ interact with each other (just via Buckingham)
+# QM C,H atoms do not LJ interact with each other
+# MM Si,O and QM C,H do LJ interact with each other
+
+pair_style      hybrid/overlay buck 6.5 lj/cut 6.5 coul/cut 6.5
+
+pair_coeff      1 1 buck 3150.462646 0.35032282 626.7519553
+pair_coeff      2 2 buck 659.595398 0.38609055 26.836679
+pair_coeff      1 2 buck 27029.419922 0.19385082 148.099091
+
+#pair_coeff      1 3 lj/cut 0.009087 3.613
+#pair_coeff      1 4 lj/cut 0.00344258 3.238
+#pair_coeff      2 3 lj/cut 0.01419429 3.1
+#pair_coeff      2 4 lj/cut 0.0035857762359063315 1.932779  # same as water dimer
+pair_coeff      1 3 lj/cut 0.09087 3.613
+pair_coeff      1 4 lj/cut 0.0344258 3.238
+pair_coeff      2 3 lj/cut 0.1419429 3.1
+pair_coeff      2 4 lj/cut 0.035857762359063315 1.932779  # same as water dimer
+
+pair_coeff      3 3 lj/cut 0.0 3.4
+pair_coeff      4 4 lj/cut 0.0 2.65
+pair_coeff      * * coul/cut
+
+bond_style      harmonic
+bond_coeff      1 29.40 1.09
+
+angle_style     harmonic
+angle_coeff     1 0.172325 109.5
+
+# remove bonds/angles in QM methane molecule
+
+delete_bonds    qm multi remove special
+System init for delete_bonds ...
+Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.5
+  ghost atom cutoff = 8.5
+  binsize = 4.25, bins = 3 3 3
+  3 neighbor lists, perpetual/occasional/extra = 3 0 0
+  (1) pair buck, perpetual, skip from (3)
+      attributes: half, newton on
+      pair build: skip
+      stencil: none
+      bin: none
+  (2) pair lj/cut, perpetual, skip from (3)
+      attributes: half, newton on
+      pair build: skip
+      stencil: none
+      bin: none
+  (3) pair coul/cut, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d
+      bin: standard
+Deleting bonds ...
+  0 total bonds, 0 turned on, 0 turned off
+  0 total angles, 0 turned on, 0 turned off
+  0 total dihedrals, 0 turned on, 0 turned off
+  0 total impropers, 0 turned on, 0 turned off
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     0 = max # of 1-2 neighbors
+     0 = max # of 1-3 neighbors
+     0 = max # of 1-4 neighbors
+     4 = max # of special neighbors
+  special bonds CPU = 0.000 seconds
+
+velocity        all create 300.0 458732
+
+neighbor	1.0 bin
+neigh_modify	delay 0 every 1 check yes
+
+# dynamic or frozen zeolite
+
+#fix		1 all nve
+fix		1 qm nve
+
+fix             2 qm nwchem template.methane.nw methane.nw                 log.pwdft.zeolite Si O C H
+fix_modify      2 energy yes
+
+# convert dump file forces to Hartree/Bohr for comparison to NWChem
+
+variable        fx atom fx/1185.8
+variable        fy atom fy/1185.8
+variable        fz atom fz/1185.8
+
+dump            1 all custom 1 dump.zeolite.qmmm id x y z q v_fx v_fy v_fz
+dump_modify     1 sort id form

examples/nwchem/qcoul.py

@@ -0,0 +1,42 @@
+
+from math import sqrt
+
+lines = """
+1 1 1 -0.8476   0.13513   -0.05627    0.09445
+2 1 2  0.4238   0.79058    0.37197   -0.48187
+3 1 2  0.4238  -0.33549   -0.65246   -0.50563
+4 2 1 -0.8476   0.02105    0.48666    2.81228
+5 2 2  0.4238  -0.73085    1.07786    2.68447
+6 2 2  0.4238   0.21349    0.19973    1.90016
+"""
+
+qqrd2e = 332.06371
+
+q = []
+x = []
+
+lines = lines.split('\n')
+for line in lines[1:7]:
+  print line
+  id,imol,itype,qone,xone,yone,zone = line.split()
+  q.append(float(qone))
+  x.append((float(xone),float(yone),float(zone)))
+
+qp4 = 0.0
+qp5 = 0.0
+qp6 = 0.0
+
+for i in range(6):
+  for j in range(6):
+    if i == j: continue
+    dx = x[i][0] - x[j][0]
+    dy = x[i][1] - x[j][1]
+    dz = x[i][2] - x[j][2]
+    r = sqrt(dx*dx + dy*dy + dz*dz)
+    eng = qqrd2e * q[i]*q[j] / r
+    print "Eng of atoms %d-%d = %g",i+1,j+1,eng
+    if i+1 == 4: qp4 += eng/q[i]
+    if i+1 == 5: qp5 += eng/q[i]
+    if i+1 == 6: qp6 += eng/q[i]
+
+print "QPOT full: 4 %g 5 %g 6 %g" % (qp4,qp5,qp6)

examples/nwchem/template.methane.nw

@@ -0,0 +1,17 @@
+Title "LAMMPS wrapping of PWDFT"
+
+memory 1900 mb
+
+echo
+
+GEOMINSERT
+
+nwpw
+  xc pbe
+  cutoff 30.0
+  2d-hcurve
+  tolerances 1.0e-9 1.0-9
+  apc on
+end
+
+task pspw gradient

examples/nwchem/template.w.nw

@@ -0,0 +1,16 @@
+Title "LAMMPS wrapping of PWDFT"
+
+memory 1900 mb
+
+echo
+
+GEOMINSERT
+
+nwpw
+   2d-hcurve
+   initialize_wavefunction on
+   cutoff 10.0
+   xc pbe
+end
+
+task pspw gradient

examples/nwchem/template.water.nw

@@ -0,0 +1,17 @@
+Title "LAMMPS wrapping of PWDFT"
+
+memory 1900 mb
+
+echo
+
+GEOMINSERT
+
+nwpw
+  xc pbe
+  cutoff 30.0
+  2d-hcurve
+  tolerances 1.0e-9 1.0-9
+  apc on
+end
+
+task pspw gradient

examples/nwchem/vasp2lmp.py

@@ -0,0 +1,156 @@
+#!/bin/env python
+
+# convert VASP POSCAR file(s) to LAMMPS data file(s)
+
+# Syntax: python vasp2lmp.py vaspfile datafile
+#   if both args are files, vaspfile is POSCAR input, datafile is output
+#   if both args are dirs:
+#     each POSCAR file in vasp dir is converted to LAMMPS data file in data dir
+#     output filename = POSCAR file name prepended with "data."
+
+# POSCAR file settings this script recognizes
+#   comment on line 1 (ignored)
+#   scale factor > 0.0 on line 2 for box edge vectors and Cartesian coords
+#   3 box edge vectors on lines 3,4,5
+#   element names on line 6
+#   species counts on line 7
+#   optional Selective Dynamics line 8 (ignored)
+#   Cartesian or Direct line 8 (or 9 if Selective Dynamics)
+#   any T/F flags at end of atom coords lines are ignored
+#   any remaining lines are ignored
+
+# ---------------------------
+# error message
+
+def error(str=None):
+  if not str: print "Syntax: vasp2lmp.py POSCARfile datafile"
+  else: print str
+  sys.exit()
+  
+# ---------------------------
+# convert a VASP POSCAR file to a LAMMPS data file
+
+def convert(vaspfile,datafile):
+
+  # read VASP POSCAR file
+  
+  lines = open(vaspfile,"r").readlines()
+
+  comment = lines[0]
+  scale = float(lines[1])
+  if scale < 0.0: scale = 1.0
+  
+  avec = lines[2].split()
+  avec = [scale*float(one) for one in avec]
+  bvec = lines[3].split()
+  bvec = [scale*float(one) for one in bvec]
+  cvec = lines[4].split()
+  cvec = [scale*float(one) for one in cvec]
+
+  if avec[1] != 0.0 or avec[2] != 0.0 or bvec[2] != 0.0:
+    error("Triclinic VASP box does not meet LAMMPS triclinic restrictions")
+    
+  elements = lines[5].split()
+  counts = lines[6].split()
+  if len(elements) != len(counts):
+    error("POSCAR file elements and counts lines do not match")
+
+  ntypes = len(elements)
+  counts = [int(one) for one in counts]
+  natoms = sum(counts)
+  
+  label = lines[7].strip().lower()
+  firstatomline = 8
+
+  if label[0] == 's':
+    label = lines[8].strip().lower()
+    firstatomline = 9
+    
+  if label[0] == 'c' or label[0] == 'k': fractional = 0
+  elif label[0] == 'd': fractional = 1
+  else: error("POSCAR file Cartesian or Direct line is missing")
+
+  atomcoords = lines[firstatomline:firstatomline + natoms]
+
+  xyz = []
+  for line in atomcoords:
+    coords = line.split()
+    x = float(coords[0])
+    y = float(coords[1])
+    z = float(coords[2])
+    if fractional == 1:
+      xnew = x*avec[0] + y*bvec[0] + z*cvec[0]
+      ynew = x*avec[1] + y*bvec[1] + z*cvec[1]
+      znew = x*avec[2] + y*bvec[2] + z*cvec[2]
+      x,y,z = xnew,ynew,znew
+    else:
+      x *= scale
+      y *= scale
+      z *= scale
+    xyz.append((x,y,z))
+
+  # write LAMMPS data file
+
+  d = data()
+
+  d.title = "LAMMPS data file from VASP file %s\n" % vaspfile
+  d.headers["atoms"] = natoms
+  d.headers["atom types"] = ntypes
+  d.headers["xlo xhi"] = (0.0,"%20.17g" % avec[0])
+  d.headers["ylo yhi"] = (0.0,"%20.17g" % bvec[1])
+  d.headers["zlo zhi"] = (0.0,"%20.17g" % cvec[2])
+  if (bvec[0] != 0.0) or (cvec[0] != 0.0) or (cvec[1] != 0.0):
+    d.headers["xy xz yz"] = ("%20.17g" % bvec[0],
+                             "%20.17g" % cvec[0],
+                             "%20.17g" % cvec[1])
+
+  id = 0
+  itype = 0
+  previous = 0
+  
+  atomlines = []
+  for x,y,z in xyz:
+    if id == counts[itype] + previous:
+      previous += counts[itype]
+      itype += 1
+    line = "%d %d %20.17g %20.17g %20.17g\n" % (id+1,itype+1,x,y,z)
+    atomlines.append(line)
+    id += 1
+
+  d.sections["Atoms"] = atomlines
+
+  d.write(datafile)
+
+# ---------------------------
+# main program
+# ---------------------------
+
+import sys,os,glob
+from data import data
+
+if (len(sys.argv) != 3): error()
+
+vaspfile = sys.argv[1]
+datafile = sys.argv[2]
+
+vflag = os.path.isdir(vaspfile)
+dflag = os.path.isdir(datafile)
+
+if (vflag and not dflag) or (not vflag and dflag):
+  error("Both comand-line args must be files or both must be dirs")
+  sys.exit()
+  
+# if both are files, vaspfile is input, datafile is output
+
+if not vflag and not dflag:
+  print "Converting",vaspfile,"to",datafile
+  convert(vaspfile,datafile)
+
+# if both are dirs, then process all files in POSCAR dir
+
+if vflag and dflag:
+  vaspfiles = glob.glob("%s/*" % vaspfile)
+  for infile in vaspfiles:
+    outfile = ("%s/data." % datafile) + os.path.basename(infile)
+    print "Converting",infile,"to",outfile
+    convert(infile,outfile)

lib/README

@@ -45,6 +45,8 @@ mscg          hooks to the MSCG library, used by fix_mscg command
                 from Jacob Wagner and Greg Voth group (U Chicago)
 netcdf        hooks to a NetCDF library installed on your system
                 from Lars Pastewka (Karlsruhe Institute of Technology)
+nwchem        hooks to the NWChem PWDFT library for AIMD, QM/MM modeling
+                from Eric Bylaska (PNNL)
 plugin        settings to load styles into LAMMPS from plugins
                 from Axel Kohlmeyer (Temple U)
 poems         POEMS rigid-body integration package, POEMS package

lib/nwchem/Install.py

@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+
+"""
+Install.py tool to download, unpack, build, and link to the NWChem PWDFT
+library used to automate the steps described in the README file in this dir
+"""
+
+from __future__ import print_function
+import sys, os, subprocess, shutil, zipfile
+from argparse import ArgumentParser
+
+sys.path.append('..')
+from install_helpers import fullpath, geturl, checkmd5sum
+
+parser = ArgumentParser(prog='Install.py',
+                        description="LAMMPS library build wrapper script")
+
+# settings
+
+version = "PWDFT"
+url = "https://github.com/ebylaska/PWDFT/archive/master.zip"
+
+# known checksums for different PWDFT versions. used to validate the download.
+checksums = { \
+        'PWDFT' : '???' \
+        }
+
+# extra help message
+
+HELP = """
+Syntax from src dir: make lib-nwchem args="-b"
+                 or: make lib-nwchem args="-p /usr/local/PWDFT"
+                 or: make lib-nwchem args="-b -v PWDFT"
+
+Syntax from lib dir: python Install.py -b
+                 or: python Install.py -p /usr/local/PWDFT
+
+Example:
+
+make lib-nwchem args="-b"   # download/build in lib/nwchem/PWDFT
+make lib-nwchem args="-p $HOME/PWDFT" # use existing PWDFT installation in $HOME/PWDFT
+"""
+
+# parse and process arguments
+
+pgroup = parser.add_mutually_exclusive_group()
+pgroup.add_argument("-b", "--build", action="store_true",
+                    help="download and build the PWDFT library")
+pgroup.add_argument("-p", "--path",
+                    help="specify folder of existing PWDFT installation")
+parser.add_argument("-v", "--version", default=version,
+                    help="set version of PWDFT to download and build (default: %s)" % version)
+
+args = parser.parse_args()
+
+# print help message and exit, if neither build nor path options are given
+if not args.build and not args.path:
+  parser.print_help()
+  sys.exit(HELP)
+
+buildflag = args.build
+pathflag = args.path is not None
+PWDFTpath = args.path
+
+homepath = fullpath(".")
+homedir = os.path.join(homepath, version)
+
+if pathflag:
+    if not os.path.isdir(PWDFTpath):
+      sys.exit("PWDFT path %s does not exist" % PWDFTpath)
+    homedir = fullpath(PWDFTpath)
+
+# download and unpack PWDFT zipfile from GitHub
+
+if buildflag:
+  print("Downloading NWChem PWDFT ...")
+  PWDFTzip = os.path.join(homepath, version) + '.zip'
+  geturl(url, PWDFTzip)
+
+  # verify downloaded archive integrity via md5 checksum, if known.
+  #if version in checksums:
+  #  if not checkmd5sum(checksums[version], vorotar):
+  #    sys.exit("Checksum for Voro++ library does not match")
+
+  print("Unpacking PWDFT zipfile ...")
+  srcpath = os.path.join(homepath, version)
+  if os.path.exists(srcpath):
+    shutil.rmtree(srcpath)
+  if zipfile.is_zipfile(PWDFTzip):
+    zip = zipfile.ZipFile(PWDFTzip)
+    zip.extractall(path=homepath)
+    os.rename("PWDFT-master","PWDFT")
+    os.remove(PWDFTzip)
+  else:
+    sys.exit("File %s is not a supported archive" % vorotar)
+
+  #if os.path.basename(homedir) != version:
+  #  if os.path.exists(homedir):
+  #    shutil.rmtree(homedir)
+  #  os.rename(srcpath, homedir)
+
+# build PWDFT
+
+if buildflag:
+  print("Building PWDFT ...")
+  cmd = 'cd %s; mkdir build_library; cd build_library; cmake ../Nwpw -DMAKE_LIBRARY=true -DCMAKE_POSITION_INDEPENDENT_CODE=ON; make -j' % homedir
+  try:
+    txt = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
+    print(txt.decode('UTF-8'))
+  except subprocess.CalledProcessError as e:
+    print("Make failed with:\n %s" % e.output.decode('UTF-8'))
+    sys.exit(1)
+
+# create 2 links in lib/voronoi to Voro++ src dir
+
+print("Creating links to PWDFT lib file")
+if os.path.isfile("liblink") or os.path.islink("liblink"):
+  os.remove("liblink")
+os.symlink(os.path.join(homedir, 'build_library'), 'liblink')

lib/nwchem/README

@@ -0,0 +1,60 @@
+This directory contains links to the NWChem PWDFT library which is
+required to use the NWCHEM package and its fix nwchem command in a
+LAMMPS input script.
+
+The PWDFT library is available at GitHub URL and was
+developed by Eric Bylaska (PNNL).
+
+You can type "make lib-nwchem" from the src directory to see help on
+how to download and build this library via make commands, or you can
+do the same thing by typing "python Install.py" from within this
+directory, or you can do it manually by following the instructions
+below.
+
+-----------------
+
+Instructions:
+
+1.  Download PWDFT at GitHub URL ...
+    either as a tarball or via SVN, and unpack the
+    tarball either in this /lib/voronoi directory
+    or somewhere else on your system.
+
+2.  Compile PWDFT from within its home directory
+    % make
+
+3.  There is no need to install PWDFT on your system if you only wish
+    to use it from LAMMPS.  You can install it if you wish to use it
+    stand-alone or from other codes:
+    a) install under the default /usr/local
+       % sudo make install
+    b) install under a user-writeable location by first
+         changing the PREFIX variable in the config.mk file, then
+       % make install
+
+4.  Create a soft link in this dir (lib/nwchem)
+    to the PWDFT build_library directory.  E.g if you built PWDFT in this dir:
+      % ln -s PWDFT/build_library liblink
+    This link could instead be set to the lib
+    directories created by a build of PWDFT elsewher on your system, e.g.
+      % ln -s /user/sjplimp/nwchem/PWDFT/build_library liblink
+
+-----------------
+
+When these steps are complete you can build LAMMPS
+with the NWCHEM package installed:
+
+% cd lammps/src
+% make yes-nwchem
+% make g++ (or whatever target you wish)
+
+Note that if you download and unpack a new LAMMPS tarball, the
+"liblink" files will be lost and you will need to re-create them (step
+4).  If you built PWDFT in this directory (as opposed to somewhere
+else on your system) and did not install it somewhere else, you will
+also need to repeat steps 1,2,3.
+
+The Makefile.lammps file in this directory is there for compatibility
+with the way other libraries under the lib dir are linked with by
+LAMMPS.  However, PWDFT requires no auxiliary files or settings, so
+its variables are blank.

src/Makefile

@@ -110,6 +110,7 @@ PACKAGE = \
 	mpiio \
 	mscg \
 	netcdf \
+	nwchem \
 	openmp \
 	opt \
 	orient \

src/NWCHEM/Install.sh

@@ -0,0 +1,69 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+mode=$1
+
+# enforce using portable C locale
+LC_ALL=C
+export LC_ALL
+
+# arg1 = file, arg2 = file it depends on
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+# all package files with no dependencies
+
+for file in *.cpp *.h; do
+  test -f ${file} && action $file
+done
+
+# edit 2 Makefile.package files to include/exclude package info
+
+if (test $1 = 1) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*nwchem[^ \t]* //' ../Makefile.package
+    sed -i -e 's/[^ \t]*pwdft[^ \t]* //' ../Makefile.package
+    sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/nwchem/liblink -Wl,-rpath=/home/sjplimp/lammps/git/lib/nwchem/PWDFT/build_library |' ../Makefile.package
+    sed -i -e 's|^PKG_LIB =[ \t]*|&-lpwdft |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(nwchem_SYSINC) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(nwchem_SYSLIB) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(nwchem_SYSPATH) |' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*nwchem.*$/d' ../Makefile.package.settings
+    # multiline form needed for BSD sed on Macs
+    sed -i -e '4 i \
+include ..\/..\/lib\/nwchem\/Makefile.lammps
+' ../Makefile.package.settings
+  fi
+
+elif (test $1 = 0) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*nwchem[^ \t]* //' ../Makefile.package
+    sed -i -e 's/[^ \t]*nwchem[^ \t]* //' ../Makefile.package
+    sed -i -e 's/[^ \t]*pwdft[^ \t]* //' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*nwchem.*$/d' ../Makefile.package.settings
+  fi
+
+fi

src/NWCHEM/README

@@ -0,0 +1,23 @@
+The NWCHEM package enables LAMMPS to wrap the NWChem quantum code,
+specifically its PWDFT library for plane-wave electronic structure
+calculations.  With the fix nwchem command, either AIMD or QM/MM
+calculations can be done using quantum forces and energy computed by
+NWChem.
+
+The NWChem PWDFT library is available to GitHub URL.  It was developed
+by Eric Bylaska (PNNL), one of the NWChem developers.
+
+That library can be downloaded and built in lib/nwchem or elsewhere on
+your system, which must be done before building LAMMPS with this
+package.  Details of the download, build, and install process for
+PWDFT are given in the lib/nwchem/README file, and scripts are
+provided to help automate the process.  Also see the LAMMPS manual for
+general information on building LAMMPS with external libraries.  The
+settings in the Makefile.lammps file in lib/nwchem must be correct for
+LAMMPS to build correctly with this package installed.
+
+Once you have successfully built LAMMPS with this package and NWChem
+PWDFT, you can test it using one or more of the input files from the
+examples dir:
+
+mpirun -np 4 ./lmp_mpi < lammps/examples/nwchem/in.charges.pbc.aimd

src/NWCHEM/fix_nwchem.h

@@ -0,0 +1,162 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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
+// clang-format off
+FixStyle(nwchem,FixNWChem);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_NWCHEM_H
+#define LMP_FIX_NWCHEM_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixNWChem : public Fix {
+ public:
+  FixNWChem(class LAMMPS *, int, char **);
+  virtual ~FixNWChem();
+  int setmask() override;
+  void init() override;
+  void setup(int) override;
+  void setup_post_neighbor() override;
+  void setup_pre_force(int) override;
+  void post_neighbor() override;
+  void pre_force(int) override;
+  void post_force(int) override;
+  void min_setup(int) override;
+  void min_post_neighbor() override;
+  void min_pre_force(int) override;
+  void min_post_force(int) override;
+  int pack_forward_comm(int, int *, double *, int, int *) override;
+  void unpack_forward_comm(int, int, double *) override;
+  int pack_reverse_comm(int, int, double *) override;
+  void unpack_reverse_comm(int, int *, double *) override;
+  double compute_scalar() override;
+  double memory_usage() override;
+
+ protected:
+  char *nw_template;   // template file for NWChem input params
+  char *nw_input;      // generated NWChem input file w/ box and coords
+  char *nw_output;     // output from NWChem, NULL if none
+
+  int pbcflag;         // 1 if fully periodic, 0 if fully non-periodic
+  int mode;            // AIMD or QMMM
+  int qflag;           // 1 if per-atom charge defined, 0 if not
+  int qm_init;         // 1 if NWChem and qqm are initialized, 0 if not
+  char **elements;     // species name of each LAMMPS atom type
+
+  double qmenergy;     // QM energy
+       
+  // data for QMMM mode
+
+  int nqm;                   // # of QM atoms
+  tagint *qmIDs;             // IDs of QM atoms in ascending order
+  double **xqm,**fqm;        // QM coords and forces
+  double *qqm;               // QM charges
+  int *tqm;                  // QM atom types
+  double *qpotential;        // Coulomb potential
+  double **xqm_mine;         // same values for QM atoms I own
+  double *qqm_mine;
+  int *tqm_mine;
+  double *qpotential_mine;
+  int *qm2owned;             // index of local atom for each QM atom
+                             // index = -1 if this proc does not own
+  
+  double *ecoul;             // peratom Coulombic energy from LAMMPS
+  int ncoulmax;              // length of ecoul
+
+  // conversion factors between LAMMPS and NWChem units
+
+  double lmp2qm_distance,lmp2qm_energy,qm2lmp_force,qm2lmp_energy;
+
+  class Compute *c_pe;      // NOTE: not sure if need this
+  class Pair *pair_coul;    // ptr to instance of pair coul/long
+
+  // local methods
+
+  void nwchem_initialize();
+  void pre_force_qmmm(int);
+  void post_force_qmmm(int);
+  void post_force_aimd(int);
+
+  void set_qm2owned();
+  void set_qqm();
+  void set_tqm();
+  void set_xqm();
+
+  void dummy_pspw_input(MPI_Comm, std::string &) {}
+  int dummy_pspw_qmmm_minimizer(MPI_Comm, double *, double *, 
+                                double *, double *, double *);
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Must use units metal with fix nwchem command
+
+UNDOCUMENTED
+
+E: Fix nwchem currently runs only in serial
+
+UNDOCUMENTED
+
+E: LAMMPS is linked against incompatible NWChem library
+
+UNDOCUMENTED
+
+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: Fix nwchem does not yet support a LAMMPS calculation of a Coulomb potential
+
+UNDOCUMENTED
+
+E: Could not find fix nwchem compute ID
+
+UNDOCUMENTED
+
+E: Fix nwchem compute ID does not compute pe/atom
+
+UNDOCUMENTED
+
+E: Fix nwchem requires 3d problem
+
+UNDOCUMENTED
+
+E: Fix nwchem cannot compute Coulomb potential
+
+UNDOCUMENTED
+
+E: Fix nwchem requires 3d simulation
+
+UNDOCUMENTED
+
+W: Fix nwchem should come after all other integration fixes
+
+UNDOCUMENTED
+
+E: Internal NWChem problem
+
+UNDOCUMENTED
+
+*/

src/pair_coul_cut.cpp

@@ -106,6 +106,7 @@ void PairCoulCut::compute(int eflag, int vflag)
         f[i][0] += delx * fpair;
         f[i][1] += dely * fpair;
         f[i][2] += delz * fpair;
+
         if (newton_pair || j < nlocal) {
           f[j][0] -= delx * fpair;
           f[j][1] -= dely * fpair;

src/pair_lj_cut.cpp

@@ -119,6 +119,7 @@ void PairLJCut::compute(int eflag, int vflag)
         f[i][0] += delx * fpair;
         f[i][1] += dely * fpair;
         f[i][2] += delz * fpair;
+
         if (newton_pair || j < nlocal) {
           f[j][0] -= delx * fpair;
           f[j][1] -= dely * fpair;