start

basis
  n library 3-21g
  h library 3-21g
end

driver
  cvgopt 0.0001
end

print none

python
  from __future__ import print_function
  """
  import Gnuplot
  """
  import time, signal
  from math import *
  from nwgeom 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 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 style data 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 style data linespoints')
  gr('set xrange [-25:25]')
  """
  # Generate points in a visually interesting order
  points = list(range(-24,30,6)) + list(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