-
Notifications
You must be signed in to change notification settings - Fork 9
/
constraints.py
128 lines (103 loc) · 3.3 KB
/
constraints.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
# This file is part of PyCI.
#
# PyCI is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# PyCI is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License
# along with PyCI. If not, see <http://www.gnu.org/licenses/>.
r"""PyCI RDM Constraints module."""
import numpy as np
from scipy.optimize import root
__all__ = [
"find_closest_sdp",
"calc_P",
"calc_Q",
"calc_G",
"calc_T1",
"calc_T2",
"calc_T2_prime",
]
def find_closest_sdp(dm, constraint, alpha):
r"""
Projection onto a semidefinite constraint.
Parameters
----------
dm : np.ndarray
Density matrix.
constraint : function
Positive semidefinite constraint, linear mapping.
alpha : float
Value of the correct trace.
"""
#symmetrize if necessary
constrained = constraint(dm)
L = constrained + constrained.conj().T
#find eigendecomposition
vals, vecs = np.linalg.eigh(L)
#calculate the shift, sigma0
sigma0 = calculate_shift(vals, alpha)
#calculate the closest semidefinite positive matrix with correct trace
L_closest = vals @ np.diag(vecs - sigma0) @ vecs.conj().T
# return the reconstructed density matrix
return constraint(L_closest).conj().T
def calculate_shift(eigenvalues, alpha):
r"""
Calculate the shift for density matrix for it to have the correct trace.
This step shifts the spectrum and eliminates the negative eigenvalues.
Parameters
----------
eigenvalues : np.ndarray
Eigenvalues of not shifted matrix.
alpha : float
Value of the coprrect trace.
"""
#sample code, to be confirmed
trace = lambda sigma0: np.sum(np.heaviside(eigenvalues - sigma0, 0.5)*(eigenvalues - sigma0))
constraint = lambda x: trace(x) - alpha
res = root(constraint, 0)
return res.x
def calc_P():
pass
def calc_Q():
pass
def calc_G(gamma, N, conjugate=False):
"""
Calculating G tensor
Parameters
----------
gamma: np.ndarray
1DM tensor
N: int
number of electrons in the system
conjugate: bool
conjugate or regular condition
Returns
-------
np.ndarray
Notes
-----
G is defined as:
.. math::
\mathcal{G}_1(\Gamma)_{\alpha \beta ; \gamma \delta}=\delta_{\beta \delta} \rho_{\alpha \gamma}-\Gamma_{\alpha \delta ; \gamma \beta}
\mathcal{G}^{\prime}(\Gamma)_{\alpha \beta ; \gamma \delta}=\delta_{\beta \delta} \rho_{\alpha \gamma}-\Gamma_{\alpha \delta ; \gamma \beta}-\rho_{\alpha \beta} \rho_{\gamma \delta}
"""
eye = np.eye(gamma.shape[0])
rho = 1/(N - 1) * np.einsum('abgb -> ag', gamma)
res = np.einsum('bd, ag -> abgd', eye, rho) - np.einsum('adgb -> abgd', gamma)
if not conjugate:
return res
else:
return res - np.einsum('ab, gd -> abgd', eye, eye)
def calc_T1():
pass
def calc_T2():
pass
def calc_T2_prime():
pass