Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

# Copyright 2011-2017 Kwant authors. 

# 

# This file is part of Kwant. It is subject to the license terms in the file 

# LICENSE.rst found in the top-level directory of this distribution and at 

# http://kwant-project.org/license. A list of Kwant authors can be found in 

# the file AUTHORS.rst at the top-level directory of this distribution and at 

# http://kwant-project.org/authors. 

 

import keyword 

from collections import defaultdict 

 

import numpy as np 

 

import sympy 

import sympy.abc 

import sympy.physics.quantum 

from sympy.core.function import AppliedUndef 

from sympy.core.sympify import converter 

from sympy.core.core import all_classes as sympy_classes 

from sympy.physics.matrices import msigma as _msigma 

 

import warnings 

 

from .._common import reraise_warnings 

 

# TODO: remove when sympy correctly includes MutableDenseMatrix (lol). 

sympy_classes = set(sympy_classes) | {sympy.MutableDenseMatrix} 

 

momentum_operators = sympy.symbols('k_x k_y k_z', commutative=False) 

position_operators = sympy.symbols('x y z', commutative=False) 

 

pauli = [sympy.eye(2), _msigma(1), _msigma(2), _msigma(3)] 

 

extra_ns = sympy.abc._clash.copy() 

extra_ns.update({s.name: s for s in momentum_operators}) 

extra_ns.update({s.name: s for s in position_operators}) 

extra_ns.update({'kron': sympy.physics.quantum.TensorProduct, 

'eye': sympy.eye, 'identity': sympy.eye}) 

extra_ns.update({'sigma_{}'.format(c): p for c, p in zip('0xyz', pauli)}) 

 

 

# workaroud for https://github.com/sympy/sympy/issues/12060 

del extra_ns['I'] 

del extra_ns['pi'] 

 

 

################ Helpers to handle sympy 

 

def lambdify(expr, locals=None): 

"""Return a callable object for computing a continuum Hamiltonian. 

 

.. warning:: 

This function uses ``eval`` (because it calls ``sympy.sympify``), and 

thus should not be used on unsanitized input. 

 

If necessary, the given expression is sympified using 

`kwant.continuum.sympify`. It is then converted into a callable object. 

 

Parameters 

---------- 

expr : str or SymPy expression 

Expression to be converted into a callable object 

locals : dict or ``None`` (default) 

Additional definitions for `~kwant.continuum.sympify`. 

 

Examples 

-------- 

>>> f = lambdify('a + b', locals={'b': 'b + c'}) 

>>> f(1, 3, 5) 

9 

 

>>> ns = {'sigma_plus': [[0, 2], [0, 0]]} 

>>> f = lambdify('k_x**2 * sigma_plus', ns) 

>>> f(0.25) 

array([[ 0. , 0.125], 

[ 0. , 0. ]]) 

""" 

with reraise_warnings(level=4): 

expr = sympify(expr, locals) 

 

args = [s.name for s in expr.atoms(sympy.Symbol)] 

args += [str(f.func) for f in expr.atoms(AppliedUndef, sympy.Function)] 

 

return sympy.lambdify(sorted(args), expr) 

 

 

def sympify(expr, locals=None): 

"""Sympify object using special rules for Hamiltonians. 

 

If `'expr`` is already a type that SymPy understands, it will do nothing 

but return that value. Note that ``locals`` will not be used in this 

situation. 

 

Otherwise, it is sympified by ``sympy.sympify`` with a modified namespace 

such that 

 

* the position operators "x", "y" or "z" and momentum operators "k_x", 

"k_y", and "k_z" do not commute, 

* all single-letter identifiers and names of Greek letters (e.g. "pi" or 

"gamma") are treated as symbols, 

* "kron" corresponds to ``sympy.physics.quantum.TensorProduct``, and 

"identity" to ``sympy.eye``, 

* "sigma_0", "sigma_x", "sigma_y", "sigma_z" are the Pauli matrices. 

 

In addition, Python list literals are interpreted as SymPy matrices. 

 

.. warning:: 

This function uses ``eval`` (because it calls ``sympy.sympify``), and 

thus should not be used on unsanitized input. 

 

Parameters 

---------- 

expr : str or SymPy expression 

Expression to be converted to a SymPy object. 

locals : dict or ``None`` (default) 

Additional entries for the namespace under which `expr` is sympified. 

The keys must be valid Python variable names. The values may be 

strings, since they are all are sent through `continuum.sympify` 

themselves before use. (Note that this is a difference to how 

``sympy.sympify`` behaves.) 

 

.. note:: 

When a value of `locals` is already a SymPy object, it is used 

as-is, and the caller is responsible to set the commutativity of 

its symbols appropriately. This possible source of errors is 

demonstrated in the last example below. 

 

Returns 

------- 

result : SymPy object 

 

Examples 

-------- 

>>> sympify('k_x * A(x) * k_x + V(x)') 

k_x*A(x)*k_x + V(x) # as valid sympy object 

 

>>> sympify('k_x**2 + V', locals={'V': 'V_0 + V(x)'}) 

k_x**2 + V(x) + V_0 

 

>>> ns = {'sigma_plus': [[0, 2], [0, 0]]} 

>>> sympify('k_x**2 * sigma_plus', ns) 

Matrix([ 

[0, 2*k_x**2], 

[0, 0]]) 

 

>>> sympify('k_x * A(c) * k_x', locals={'c': 'x'}) 

k_x*A(x)*k_x 

>>> sympify('k_x * A(c) * k_x', locals={'c': sympy.Symbol('x')}) 

A(x)*k_x**2 

 

""" 

stored_value = None 

 

# if ``expr`` is already a ``sympy`` object we may terminate a code path 

if isinstance(expr, tuple(sympy_classes)): 

if locals: 

warnings.warn('Input expression is already SymPy object: ' 

'"locals" will not be used.', 

RuntimeWarning, 

stacklevel=2) 

 

# We assume that all present functions, like "sin", "cos", will be 

# provided by user during the final evaluation through "params". 

# Therefore we make sure they are defined as AppliedUndef, not built-in 

# sympy types. 

subs = {r: sympy.Function(str(r.func))(*r.args) 

for r in expr.atoms(sympy.Function)} 

 

return expr.subs(subs) 

 

# if ``expr`` is not a "sympy" then we proceed with sympifying process 

if locals is None: 

locals = {} 

 

for k in locals: 

if (not isinstance(k, str) 

or not k.isidentifier() or keyword.iskeyword(k)): 

raise ValueError( 

"Invalid key in 'locals': {}\nKeys must be " 

"identifiers and may not be keywords".format(repr(k))) 

 

# sympify values of locals before updating it with extra_ns 

locals = {k: sympify(v) for k, v in locals.items()} 

for k, v in extra_ns.items(): 

locals.setdefault(k, v) 

try: 

stored_value = converter.pop(list, None) 

converter[list] = lambda x: sympy.Matrix(x) 

hamiltonian = sympy.sympify(expr, locals=locals) 

# if input is for example ``[[k_x * A(x) * k_x]]`` after the first 

# sympify we are getting list of sympy objects, so we call sympify 

# second time to obtain ``sympy`` matrices. 

hamiltonian = sympy.sympify(hamiltonian) 

finally: 

195 ↛ 196line 195 didn't jump to line 196, because the condition on line 195 was never true if stored_value is not None: 

converter[list] = stored_value 

else: 

del converter[list] 

 

return hamiltonian 

 

 

def make_commutative(expr, *symbols): 

"""Make sure that specified symbols are defined as commutative. 

 

Parameters 

---------- 

expr: sympy.Expr or sympy.Matrix 

symbols: sequace of symbols 

Set of symbols that are requiered to be commutative. It doesn't matter 

of symbol is provided as commutative or not. 

 

Returns 

------- 

input expression with all specified symbols changed to commutative. 

""" 

symbols = [sympy.Symbol(s.name, commutative=False) for s in symbols] 

expr = expr.subs({s: sympy.Symbol(s.name) for s in symbols}) 

return expr 

 

 

def monomials(expr, gens=None): 

"""Parse ``expr`` into monomials in the symbols in ``gens``. 

 

Parameters 

---------- 

expr: sympy.Expr or sympy.Matrix 

Sympy expression to be parsed into monomials. 

gens: sequence of sympy.Symbol objects or strings (optional) 

Generators of monomials. If unset it will default to all 

symbols used in ``expr``. 

 

Returns 

------- 

dictionary (generator: monomial) 

 

Examples 

-------- 

>>> expr = kwant.continuum.sympify("A * (x**2 + y) + B * x + C") 

>>> monomials(expr, gens=('x', 'y')) 

{1: C, x: B, x**2: A, y: A} 

""" 

if gens is None: 

gens = expr.atoms(sympy.Symbol) 

else: 

gens = [sympify(g) for g in gens] 

 

if not isinstance(expr, sympy.MatrixBase): 

return _expression_monomials(expr, gens) 

else: 

output = defaultdict(lambda: sympy.zeros(*expr.shape)) 

for (i, j), e in np.ndenumerate(expr): 

mons = _expression_monomials(e, gens) 

for key, val in mons.items(): 

output[key][i, j] += val 

return dict(output) 

 

 

def _expression_monomials(expr, gens): 

"""Parse ``expr`` into monomials in the symbols in ``gens``. 

 

Parameters 

---------- 

expr: sympy.Expr 

Sympy expr to be parsed. 

gens: sequence of sympy.Symbol 

Generators of monomials. 

 

Returns 

------- 

dictionary (generator: monomial) 

""" 

expr = sympy.expand(expr) 

output = defaultdict(lambda: sympy.Integer(0)) 

for summand in expr.as_ordered_terms(): 

key = [] 

val = [] 

for factor in summand.as_ordered_factors(): 

symbol, exponent = factor.as_base_exp() 

if symbol in gens: 

key.append(factor) 

else: 

val.append(factor) 

output[sympy.Mul(*key)] += sympy.Mul(*val) 

 

return dict(output) 

 

 

################ general help functions 

 

def gcd(*args): 

if len(args) == 1: 

return args[0] 

 

L = list(args) 

 

while len(L) > 1: 

a = L[len(L) - 2] 

b = L[len(L) - 1] 

L = L[:len(L) - 2] 

 

while a: 

a, b = b%a, a 

 

L.append(b) 

 

return abs(b)