-
Notifications
You must be signed in to change notification settings - Fork 0
/
turbo2nw.py
218 lines (158 loc) · 5.18 KB
/
turbo2nw.py
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
#!/usr/bin/python
import re, sys, math, logging, random
sys.dont_write_bytecode = True
from qcldm.turbomole_format.control_format import ControlFormat
from qcldm.turbomole_format.turbo_basis import TurboBasis
from qcldm.gauss_functions.gauss_function import GaussFunction, GaussFunctionNormed, GaussFunctionContracted
from qcldm.util.log_colorizer import init_log
from qcldm.util.units import Units
from qcldm.util.elements import ELEMENTS
from qcldm.structures.atom_vector import AtomKeys
from qcldm.atom.shells import Shells
TEMPLATE = '''start cluster
title "cluster"
geometry "full-cluster" units au noautoz
symmetry c1
%%COORD%%
end
basis spherical
%%BASIS%%
end
ecp
%%ECP%%
end
so
%%SOECP%%
end
set geometry "full-cluster"
charge %%CHARGE%%
scf
vectors input hcore
vectors output ./cluster.movecs
uhf
maxiter 100
end
task scf energy
dft
XC pbe0
CONVERGENCE fast
vectors input ./cluster.movecs
vectors output ./cluster.movecs
iterations 200
grid lebedev 350 17 becke;tolerances accCoul 20 tol_rho 20
end
task dft gradient
'''
def dummy_basis():
gc = GaussFunctionContracted()
gc.fs.append([1, GaussFunctionNormed(100, 0)])
return TurboBasis('', [gc])
def dummy_ecp(name):
return TurboBasis.TurboEcp('', ELEMENTS[name].number, TurboBasis.EcpPart([[0, GaussFunction(1, 2)]]), [TurboBasis.EcpPart([[0.00000001, GaussFunction(1, 2)]])], [])
def build_basis(c, specs, placeholder):
basis_str = ''
ecp_str = ''
so_str = ''
# print specs
elbasmap = {}
elecpmap = {}
for spec in specs:
# print(spec)
el = None
bas = None
if spec[0] == 'none' and spec[1] == None:
# print(spec)
# continue
el = placeholder
bas = dummy_basis()
elecpmap[el] = dummy_ecp(placeholder)
elif spec[0] == 'none':
el = spec[1].split()[0]
bas = dummy_basis()
else:
el = spec[0].split()[0]
bas = c.bases[spec[0]]
el = el[0].upper() + el[1:].lower()
assert el not in list(elbasmap.keys()) or bas.name == elbasmap[el].name, 'multiple bases for %s: %s AND %s' % (el, bas.name, elbasmap[el].name)
elbasmap[el] = bas
if spec[1] != None:
assert el not in list(elecpmap.keys()), 'multiple ecps for %s' % el
elecpmap[el] = c.ecps[spec[1]]
for el in sorted(elbasmap.keys()):
bas = elbasmap[el]
for f in bas.functions:
l = Shells.SHELLS[f.fs[0][1].l]
basis_str += el + ' ' * (6 - len(el)) + l + '\n'
for c, g in f.fs:
basis_str += "%18.10f %15.8f" % (g.a, c) + '\n'
for el in sorted(elecpmap.keys()):
ecp = elecpmap[el]
ecp_str += el + ' ' * (2 - len(el)) + ' nelec ' + str(ecp.ncore) + '\n'
ecp_str += el + ' ' * (2 - len(el)) + ' ul\n'
for c, g in ecp.local.functions:
ecp_str += ' %2d %15.8f %15.8f\n' % (g.l, g.a, c)
for l, part in enumerate(ecp.semilocal):
ecp_str += el + ' ' * (2 - len(el)) + ' ' + Shells.SHELLS[l] + '\n'
for c, g in part.functions:
ecp_str += ' %2d %15.8f %15.8f\n' % (g.l, g.a, c)
for el in sorted(elecpmap.keys()):
ecp = elecpmap[el]
if not ecp.spinorbit:
continue
for l, part in enumerate(ecp.spinorbit):
so_str += el + ' ' * (2 - len(el)) + ' ' + Shells.SHELLS[l] + '\n'
for c, g in part.functions:
so_str += ' %2d %15.8f %15.8f\n' % (g.l, g.a, c)
return basis_str, ecp_str, so_str
init_log(sys.argv)
c = ControlFormat.from_path('.')
#for a in c.cell.atoms:
# print a
placeholder = 'He'
with open('test.nw', 'w') as f:
coord_block = ''
specs = set()
charge = 0.
embcount = {}
bqbases = []
smap = c.species_map()
for i, a in enumerate(c.cell.atoms):
embedded = AtomKeys.ESTIMATED_CHARGE in list(a.data().keys())
name = a.name()
if name.lower() == 'q':
# name = 'bq'
name = placeholder
a.data()[AtomKeys.ESTIMATED_CHARGE] += ELEMENTS[name].number
if embedded:
if name not in list(embcount.keys()):
embcount[name] = 0
embcount[name] = embcount[name] + 1
if name != 'bq':
name += str(embcount[name])
real_basisname = smap[i + 1][0]
if name == 'bq' and real_basisname != 'none':
bas_to_copy = c.bases[real_basisname]
if real_basisname not in bqbases:
bqbases.append(real_basisname)
bqbasisname = 'bq%d %s' % (bqbases.index(real_basisname) + 1, real_basisname.split()[1])
bqbasis = TurboBasis(bqbasisname, bas_to_copy.functions)
c.bases[bqbasisname] = bqbasis
smap[i + 1] = bqbasisname, smap[i + 1][1]
name = 'bq%d' % (bqbases.index(real_basisname) + 1)
coord_block += "%3s %15.10f %15.10f %15.10f" % (name, a.position()[0] / Units.BOHR, a.position()[1] / Units.BOHR, a.position()[2] / Units.BOHR)
if embedded:
ncore = c.ecps[smap[i + 1][1]].ncore if smap[i + 1][1] else 0
emcharge = a.data()[AtomKeys.ESTIMATED_CHARGE] + ncore
real_basisname = smap[i + 1][0]
if 'bq' in name and real_basisname != 'none':
if emcharge == 0:
emcharge = 1e-8
charge += emcharge
coord_block += ' charge %13.8f' % emcharge
coord_block += "\n"
specs.add(smap[i + 1])
charge = charge - round(charge)
print(charge)
chargestr = "%13.8f" % charge
basis_block, ecp_block, so_block = build_basis(c, specs, placeholder)
f.write(TEMPLATE.replace('%%COORD%%', coord_block[:-1]).replace('%%BASIS%%', basis_block[:-1]).replace('%%ECP%%', ecp_block[:-1]).replace('%%SOECP%%', so_block[:-1]).replace('%%CHARGE%%', chargestr))