-
Notifications
You must be signed in to change notification settings - Fork 4
/
gxdmopt.sh
executable file
·145 lines (125 loc) · 3.46 KB
/
gxdmopt.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#!/bin/bash
# Copyright (c) 2013-2016 Alberto Otero de la Roza
# <[email protected]>, Felix Kannemann
# <[email protected]>, Erin R. Johnson <[email protected]>,
# Ross M. Dickson <[email protected]>, Hartmut Schmider
# <[email protected]>, and Axel D. Becke <[email protected]>
## modify this
chf="blyp"
c1=0.5942
c2=1.4555
basis="6-31+g*"
ecp=""
G09="g09"
POSTG="postg"
verbose=""
guesstrick=""
cat > $2.route <<EOF
%mem=2GB
%nprocs=4
#p blyp int(grid=ultrafine)
EOF
#########
if [ -f $basis ] ; then
basisfile=$basis
basis="gen"
fi
if [ -n "$ecp" ] && [ -f "$ecp" ] ; then
ecpfile=$ecp
ecp="pseudo=read"
else
ecp=""
fi
# read info from gaussian "output"
read atoms derivs charge spin < $2
# route string
sroute="units=au output=wfx"
if [ $derivs == "2" ] ; then
sroute="${sroute} freq(noraman) punch=derivatives"
do1="1"
do2="1"
elif [ $derivs == "1" ] ; then
sroute="${sroute} force punch=derivatives"
do1="1"
fi
# prepare the G09 input file ($2.in) for the single-point "Force" calculation
rm -f $2.in >& /dev/null
if [ -n "$guesstrick" ] ; then
echo "%chk=${guesstrick%.chk}.chk" >> $2.in
fi
if [ -n "$guesstrick" ] && [ -f "${guesstrick%.chk}.chk" ] ; then
guess="guess=read"
else
guess=""
fi
cat $2.route >> $2.in
cat >> $2.in <<EOF
# $sroute $basis $ecp $guess
title
$charge $spin
$(sed -n 2,$(($atoms+1))p < $2 | cut -c 1-72)
EOF
if [ "x$basis" == "xgen" ] ; then
echo "" >> $2.in
awk '/^ *$/{next}{print}' $basisfile >> $2.in
fi
if [ -n "$ecp" ] ; then
echo "" >> $2.in
awk '/^ *$/{next}{print}' $ecpfile >> $2.in
fi
cat >> $2.in <<EOF
$2.wfx
EOF
# run G09
$G09 < $2.in > $2.out
if [ -n "$do1" ] ; then
head -n $atoms fort.7 | tr D e > $2.fgauss
fi
if [ -n "$do2" ] ; then
tail -n+$((${atoms}+1)) fort.7 | tr D e > $2.qgauss
fi
rm -f fort.7 >& /dev/null
# run postG
$POSTG $c1 $c2 $2.wfx $chf > $2.outg
grep 'WARNING -- inconsistent nelec' $2.outg && exit
# energy
e=$(grep 'total energy' $2.outg | awk '{print $NF}')
printf "%20.12e%20.12e%20.12e%20.12e\n" $e 0 0 0 | tr e D > $3
# forces
if [ -n "$do1" ] ; then
grep -A $(($atoms+1)) 'dispersion forces' $2.outg | awk '{print $2, $3, $4}' | tail -n $atoms > $2.fpostg
paste $2.fgauss $2.fpostg > $2.forces
awk '{printf("%20.12e%20.12e%20.12e\n",$1-$4,$2-$5,$3-$6)}' $2.forces | tr e D >> $3
fi
# frequencies
if [ -n "$do2" ] ; then
printf "%20.12e%20.12e%20.12e\n" 0 0 0 | tr e D >> $3 # polarizability
printf "%20.12e%20.12e%20.12e\n" 0 0 0 | tr e D >> $3
for ((i=1;i<=3*${atoms};i++)) ; do
printf "%20.12e%20.12e%20.12e\n" 0 0 0 | tr e D >> $3 # dip ders
done
grep -A $((3*$atoms*(3*$atoms+1)/2+1)) 'dispersion force constant matrix' $2.outg | \
tail -n $((3*$atoms*(3*$atoms+1)/2)) | \
awk '{printf "%s ",$NF}NR%3==0{printf"\n"}' > $2.qpostg
paste $2.qgauss $2.qpostg > $2.freqs
awk '{printf("%20.12e%20.12e%20.12e\n",$1+$4,$2+$5,$3+$6)}' $2.freqs | tr e D >> $3
fi
if [ -n "$verbose" ] ; then
# output for debug
echo "#XDM# gaussian input" >> $4
cat $2.in >> $4
echo "#XDM# gaussian output" >> $4
cat $2.out >> $4
echo "#XDM# wfx file" >> $4
cat $2.wfx >> $4
echo "#XDM# postg output" >> $4
cat $2.outg >> $4
else
echo "#XDM# tail -n 50 logfile" >> $4
tail -n 50 $2.out >> $4
echo "#XDM# grep 'Delta-E' logfile" >> $4
grep 'Delta-E' $2.out >> $4
echo "#XDM# energies from postg output (hartree)" >> $4
grep 'energy' $2.outg >> $4
fi
rm -f $2.* >& /dev/null