new examples

contrib/python/README

@@ -5,7 +5,6 @@ nwgeom.py - look in the file for (a little) more documentation
    Defines the procedures
 
    - coords = geom_get_coords(name)
-
    - geom_set_coords(name,coords)
      Get/set cartesian coordinates for named geometry
 
@@ -22,6 +21,31 @@ nwgeom.py - look in the file for (a little) more documentation
      the NWChem input string and calling the specified
      task(theory), returning results from each step in results.
 
+ts_search.nw
+     Example use of minimize1d from nwgeom.py
+
+scanexp.nw
+     Example use of scan_input from nwgeom.py
+
+hcn.nw
+     Scans the HCN -> CNH reaction path displaying bond lengths
+     and the energy as function of angle using Gnuplot.  You
+     need Python and the GNUplot and Numeric extensions
+     (www.python.org)
+
+nh3.nw
+     Scans the NH3 inversion path displaying bond lengths
+     and the energy as function of angle using Gnuplot.  You
+     need Python and the GNUplot and Numeric extensions
+     (www.python.org)
+     
+isdoesmic.nw
+     Uses python to compute a reaction energy as a function
+     of basis set.
+
+
+
+          
    
 
   

contrib/python/hcn.nw

@@ -0,0 +1,89 @@
+start
+
+# To run this you need the numeric and gnuplot
+# modules ... http://www.python.org
+
+
+# Use GNUplot to display the HCN <-> CNH reaction path
+# as it is computed.
+
+basis
+  c library 3-21g
+  n library 3-21g
+  h library 3-21g
+end
+
+print none
+
+scf; print none; end
+driver; print low; end
+
+python
+  import Gnuplot, time, signal, os
+  from math import *
+  from nwgeom import *
+
+  geometry = '''
+    geometry noprint
+      zmatrix
+        c 
+        n 1 cn
+        x 1 1. 2 90.
+        h 1 ch 3 90. 2 hcn
+      variables
+        cn  %f     # hcn=1.137 cnh=1.160
+        ch  %f     # hcn=1.050 cnh=2.143
+      constants
+        hcn %f     # hcn=180.  cnh=0.
+      end
+    end
+  '''
+
+  signal.signal(signal.SIGCHLD, signal.SIG_DFL)
+  g = Gnuplot.Gnuplot()
+  g.xlabel('HCN angle')
+  g.ylabel('Energy')
+  g.title('HCN --- CNH isomerization barrier - Energy vs. HCN angle')
+  g('set data style linespoints')
+  g('set xrange [0:360]')
+  b = Gnuplot.Gnuplot()
+  b.xlabel('HCN angle')
+  b.ylabel('Bond')
+  b.title('HCN --- CNH isomerization barrier - Bond-lengths vs. HCN angle')
+  b('set data style linespoints')
+  b('set xrange [0:360]')
+
+  cn = 1.137
+  ch = 1.050
+  edata = []
+  cndata = []
+  chdata = []
+  nhdata = []
+  for i in range(-10,11):
+    try:
+      os.remove('/tmp/hcn.drv.hess')
+    except:
+      pass
+    angle = 18.0*i
+    input_parse(geometry % (cn, ch, angle))
+    (energy, gradient) = task_optimize('scf')
+    cn = bond_length(1,2)
+    nh = bond_length(2,4)
+    ch = bond_length(1,4)
+    print ' angle=%6.1f => cn=%5.3f ch=%5.3f nh=%5.3f energy=%10.6f ' % \
+                (angle,cn,ch,nh,energy)
+    angle = angle + 180
+    edata = edata + [[angle,energy]]
+    cndata = cndata + [[angle,cn]]
+    chdata = chdata + [[angle,ch]]
+    nhdata = nhdata + [[angle,nh]]
+    g.plot(edata)
+    b.plot(Gnuplot.Data(cndata,title='CN'),\
+           Gnuplot.Data(chdata,title='CH'),\
+           Gnuplot.Data(nhdata,title='NH'))
+  time.sleep(90)
+end
+
+task python
+
+  

contrib/python/isodesmic.nw

@@ -0,0 +1,50 @@
+start iso
+
+mp2; freeze atomic; end
+
+print none
+
+python
+  energies = {}
+
+  c2h4 = '''
+    geometry noprint
+      symmetry d2h
+      C 0 0     0.672
+      H 0 0.935 1.238
+    end
+  '''
+
+  ch4 = '''
+    geometry noprint
+      symmetry td 
+      C 0     0     0
+      H 0.634 0.634 0.634
+    end
+  '''
+
+  h2 = '''
+    geometry noprint
+      H 0 0  0.378
+      H 0 0 -0.378
+    end
+  '''
+  
+  def energy(basis, geometry):
+    input_parse('''
+      basis spherical noprint
+        c library %s ; h library %s 
+      end
+    ''' % (basis, basis))
+    input_parse(geometry)
+    return task_energy('mp2')
+ 
+  for basis in ('sto-3g', '6-31g', '6-31g*', 'cc-pvdz', 'cc-pvtz'):
+     energies[basis] = 2*energy(basis, ch4) - \
+                       2*energy(basis, h2)  - \
+                         energy(basis, c2h4)
+     if (ga_nodeid() == 0): print basis, ' %8.6f' % energies[basis]
+end 
+
+
+task python

contrib/python/nh3.nw

@@ -0,0 +1,105 @@
+start
+
+basis
+  n library 3-21g
+  h library 3-21g
+end
+
+driver
+  cvgopt 0.0001
+end
+
+print none
+
+python
+  import Gnuplot, time, signal
+  from math import *
+
+  geometry = '''
+    geometry noprint
+      zmatrix
+        x
+        n 1 1.
+        h 2 r 1 a
+        h 2 r 1 a 3  120.
+        h 2 r 1 a 3 -120.
+        variables
+          r  1.0
+        constants
+          a %f
+      end
+    end
+  '''
+
+  def geom_get_coords(name):
+    try:
+      actualname = rtdb_get(name)
+    except NWChemError:
+      actualname = name
+    coords = rtdb_get('geometry:' + actualname + ':coords')
+    units = rtdb_get('geometry:'+actualname+':user units')
+    if (units == 'a.u.'):
+      factor = 1.0
+    elif (units == 'angstroms'):
+      factor = rtdb_get('geometry:'+actualname+':angstrom_to_au')
+    else:
+      raise NWChemError,'unknown units'
+    i = 0
+    while (i < len(coords)):
+      coords[i] = coords[i] / factor
+      i = i + 1
+    return coords
+
+  def get_bond_length():
+    coords = geom_get_coords('geometry')
+    x = coords[6]-coords[3]
+    y = coords[7]-coords[4]
+    z = coords[8]-coords[5]
+    return sqrt(x*x + y*y + z*z)
+
+  signal.signal(signal.SIGCHLD, signal.SIG_DFL)
+  g = Gnuplot.Gnuplot()
+  g.xlabel('Out-of-plane angle')
+  g.ylabel('Energy')
+  g.title('Inversion of NH3 - Energy vs. angle')
+  g('set data style linespoints')
+  g('set xrange [-25:25]')
+
+  gr = Gnuplot.Gnuplot()
+  gr.xlabel('Out-of-plane angle')
+  gr.ylabel('Bond')
+  gr.title('Inversion of NH3 - Bond-length vs. angle')
+  gr('set data style linespoints')
+  gr('set xrange [-25:25]')
+
+  # Generate points in a visually interesting order
+  points = range(-24,30,6) + range(21,-27,-6)
+  for angle in range(-23,25,3):
+    points = points + [angle+0.5]
+  for angle in range(22,-26,-3):
+    points = points + [angle+1.25]
+    points = points + [angle-0.25]
+
+  energies = []
+  bonds = []
+  for angle in points:
+    input_parse(geometry % (angle+90))
+    (energy,gradient) = task_optimize('scf')
+    r = get_bond_length()
+    print ' angle=%6.2f  bond=%6.3f energy=%10.6f ' % (angle,r,energy)
+    energies = energies + [[angle,energy]]
+    bonds = bonds + [[angle,r]]
+    energies.sort()
+    bonds.sort()
+    if (len(energies)> 1 and ga_nodeid() == 0): 
+      g.plot(energies)
+      gr.plot(bonds)
+
+  print ' Done!'
+
+  time.sleep(90)  # time to look at the plot
+end
+
+task python
+
+

contrib/python/scanexp.nw

@@ -0,0 +1,23 @@
+start scanexp
+
+# Scan the d exponent on oxygen through a range
+
+geometry units au noprint
+ O 0       0        0
+ H 0       1.430   -1.107
+ H 0      -1.430   -1.107
+end
+
+print none
+
+python
+  from nwgeom import *
+  basis = 'basis noprint; H library 3-21g; O library 3-21g; O d; %f 1.0; end'
+
+  results = scan_input(basis, [0.5], [0.6], 20, 'scf', task_energy)
+
+  for (p,e) in results:
+     print ' Exponent = %7.4f  Energy = %10.6f ' % (p[0], e)
+end
+
+task python

contrib/python/ts_search.nw

@@ -0,0 +1,75 @@
+start
+
+title; Hydroxylamine <-> Ammonia-n-oxide transition state
+
+# Search for the hydroxylamine <-> ammonia-n-oxide
+# transition state.  Have already done the optimization for
+# the two products and know that the major coordinate is
+# 
+#               hydroxylamine     ammonia-n-oxide
+# hon = 3-2-1       103.3              26.96
+# oh  = 2-3           0.967             2.109
+#
+# The python program does a crude 1d search for the
+# saddle point and starts the transition state
+# search from that point.
+
+
+print none
+
+geometry autosym  # Hydroxylamine minimum
+  n  -0.239    -0.678     0.0
+  o   0.237     0.710     0.0
+  h  -0.579     1.226     0.0
+  h   0.179    -1.084     0.822
+  h   0.179    -1.084    -0.822
+end
+
+basis
+  n library 3-21g
+  h library 3-21g
+  o library 3-21g
+end
+
+scf; thresh 1e-6; end
+
+python
+  from nwgeom import *
+
+  def geometry(alpha):
+     hon = 103.3*(1.0-alpha) + 26.96*alpha
+     oh  = 0.967*(1.0-alpha) + 2.109*alpha
+     input_parse('''
+       geometry noprint adjust
+         zcoord
+           bond  1 3 nh
+           bond  1 2 no
+           bond  3 2   %f oh
+           angle 3 2 1 %f hon constant
+         end
+       end
+     ''' % (oh, hon))
+
+  def energy(alpha):
+     geometry(alpha)
+     (energy,gradient) = task_optimize('scf')
+     return -energy    # -ve sign since we seek the maximum
+
+  (alpha, e) = minimize1d(energy, 0.5, 0.8, 0.02, 10)
+
+end
+
+task python
+
+print low
+
+geometry adjust  # Remove constraints but use same coord system
+  zcoord
+    bond  1 3 nh
+    bond  1 2 no
+    bond  3 2 oh
+    angle 3 2 1 hon
+  end
+end
+
+task scf saddle