mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-27 05:35:37 -04:00
97 lines
2.1 KiB
Text
97 lines
2.1 KiB
Text
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
|
|
|
|
dft; sym off;adapt off; print none; xc hfexch ;end
|
|
driver; maxiter 99; trust 0.15;print low; end
|
|
|
|
python
|
|
from __future__ import print_function
|
|
"""
|
|
import Gnuplot
|
|
"""
|
|
import 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(5,13):
|
|
try:
|
|
os.remove('/tmp/hcn.drv.hess')
|
|
except:
|
|
pass
|
|
angle = 18.0*i
|
|
input_parse(geometry % (cn, ch, angle))
|
|
(energy, gradient) = task_optimize('dft')
|
|
cn = bond_length(1,2)
|
|
nh = bond_length(2,3)
|
|
ch = bond_length(1,3)
|
|
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'))
|
|
"""
|
|
print(' Done!')
|
|
time.sleep(20)
|
|
end
|
|
|
|
task python
|
|
|
|
|