mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-27 13:45:27 -04:00
92 lines
1.9 KiB
Text
92 lines
1.9 KiB
Text
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
|
|
|
|
|