QCLS: Cluster Expansion Based Solver for Open Quantum System
-
C++ Build environment (Simply download VSBuildTools and install C++ build tools if you don't know how to configure it.)
-
Python 3.7 and later version
-
Scipy 1.4.1 and later version
-
numpy 1.16.5 and later version
Quick install
pip install -r requirements.txt
Run the following code under the path ./cluster_py
python ./setup.py install
pip install QCLSolver
conda install -c yesunhuang qclsolver
from QCLSolver.solver import Solve
from QCLSolver.data import Data
from QCLSolver.tool import*
Using QCLS
is quite similar to using the mesolve
in QuTiP, a quite simple procedure can be followed:
-
Create lists for representing the Hamiltonian, collapse operators and the initial state. Determine the tracking operators and the required cluster-expansion order.
(We are using the upper class English letters to represent creation operators and lower class ones to represent the annihilation operators. The system will automatically detect how many different letters u have inputted and assume the same number of modes in the system. )
-
Derive the equations for the systems using the construct function of
Data
class and store the instance. -
Feed the
Data
instance into theSolve
method and evolve the system. -
Process the and visualize the result.
The following code shows how to use QCLS
to simulate the Janes-Cumming model. (See also the same title notebooks from Quantum mechanics lectures with QuTiP)
The Hamiltonian can be expressed as
#parameters
wc = 1.0 * 2 * pi # cavity frequency
wa = 1.0 * 2 * pi # atom frequency
g = 0.05 * 2 * pi # coupling strength
kappa = 0.005 # cavity dissipation rate
gamma = 0.05
#Create the time list
tlist = np.linspace(0,25,101) #numpy has been imported as np
#Build up the Hamiltonian operators
Hamilton = [['Aa',wc],['Bb',wa],['Ab',g],['aB',g]]
#Build up the Collapse operators
Co_ps = [['a',kappa ],['b',gamma ]]
#Define the inital states (Only folk states are supported now)
psi0 = [0,1]
#Add the tracking operators (The result will be returned at the same order of the operators in T_o)
T_o = ['Aa','Bb']
#Derive the equations under the expanson order of 2
data = Data(Hamilton ,Co_ps , T_o , 2)
#Evolve the system (We are using the ODE solver `solve_vip` from scipy, thus the parameters and the output forms are almost identity )
sol = Solve(data , psi0, (0,25), t_eval = tlist )
#Visualize the result via matplotlib
fig, axes = plt.subplots(1, 1, figsize=(10,6))
axes.plot(tlist, np.real(sol.y[0]),color='red',linestyle='--',label="Cavity(QCLS)")
axes.plot(tlist, np.real(sol.y[1]),color='blue',linestyle='--',label="Atom excited state(QCLS)")
axes.legend(loc=0)
axes.set_xlabel('Time')
axes.set_ylabel('Occupation probability')
plt.show()
-
Prototype
class Data: def __init__(self, H, Co_ps, trackList, maxOPLen)
-
Parameters
- H: The list of Hamiltonian, the form should be of [['Operators','Coefficient']], Coefficient can be a complex number or a function
- Co_ps: Collapse operators, the form is the same as Hamiltonian.
- trackOp: The list of the operators whose mean value you are going to track.
- maxOPLen: The order of cluster expansion
-
Prototype
def UpdateInitialState(self, initialState)
-
Parameters
- initialState: A list of initial folk states. The order of mode is following the alphabet order of the letters u used to represent the operators.
-
Return
A list of new current mean value of tracking operators calculated based on the initial state.
-
Prototype
Calculate(self)
-
Parameters
None
-
Return
A list of the slope values using to evolve the mean value of operators
Warning: this function is for advanced use and might still have bugs.
-
Prototype
SetCoefHOList(self, coefHOList, args=None, t=0, ow=True)
-
Parameters
- coefHoList: New list of Hamiltonian' s coefficient.
Other default parameters is for helper use. Plz do not modify them.
-
Important
Run
UpdateCoef
withForceUpdate=true
after running this method.
Warning: this function is for advanced use and might still have bugs.
-
Prototype
SetCoefCOList(self, coefCOList, args=None, t=0, ow=True)
-
Parameters
- coefCOList: New list of collapse operators' coefficient.
Other default parameters is for helper use. Plz do not modify them.
-
Important
Run
UpdateCoef
withForceUpdate=true
after running this method.
Warning: this function is for advanced use and might still have bugs.
-
Prototype
def SetCurrentValue(self, curValueList)
-
Parameters
- curValueList: New list of the mean value of tracking operators.
-
Important
Please be noted that after the deriving process, the
T_o
list not only contains the initial tracking operators but all the operators needed to calculate the evolution.
Warning: this function is for advanced use and might still have bugs.
-
Prototype
def GetCurrentValue(self)
-
Important
Please be noted that after the deriving process, the
T_o
list not only contains the initial tracking operators but all the operators needed to calculate the evolution.
Warning: this function is for advanced use and might still have bugs.
-
Prototype
def UpdateCoef(self, t, args=None, ForceUpdate=False)
-
Parameters
1.t: current time point.
2.ForceUpdate: Force an advanced update to the coefficient list.
This function is simply connecting QCLS
to outer solver.
-
Prototype
def Solve(ddata, Initial_State, t_span, user_args=None, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)
-
Parameters
1.ddata: The
Data
instance. 2.Initial_State: A list for representing the initial state.
3.other paramters: the same as the parameters requested by
solve_ivp
from SciPy -
Return:
The same as the parameters requested by
solve_ivp
from SciPy. The result of y will be returned at the same order of the operators in T_o. -
Advanced
- If u wanna use another solver, u need to implement a convert function following the examples as follows.
def ConvertToSolver(t,y,data): data.SetCurrentValue(y.tolist()) dy=np.asarray(data.Calculate()) return dy
- Or if you are building your own solver.
def NaiveSolver(data,pace,t_range,initialState): data.UpdateInitialState(initialState) y0=data.GetCurrentValue() n=int((t_range[1]-t_range[0])//pace) tlist=np.linspace(t_range[0],t_range[1],n) y=np.zeros([len(y0),n],dtype=complex) for i in range(0,len(y0)): y[i,0]=y0[i] for i in range(1,n): k1=np.asarray(data.Calculate()) y_s=y[...,i-1]+k1*pace*0.5 data.SetCurrentValue(y_s.tolist()) k2=np.asarray(data.Calculate()) y_s=y[...,i-1]+k2*pace*0.5 data.SetCurrentValue(y_s.tolist()) k3=np.asarray(data.Calculate()) y_s=y[...,i-1]+k3*pace data.SetCurrentValue(y_s.tolist()) k4=np.asarray(data.Calculate()) dy=(k1+2*k2+2*k3+k4)/6 y[...,i]=y[...,i-1]+dy*pace data.SetCurrentValue(y[...,i].tolist()) return (y,tlist)
This function is used to investigate cluster expansion.
-
Prototype
def cluster_expansion(op: str, max_op_len: int = 5) -> Tuple[List[int], Dict[str, int], List[Tuple[complex, Tuple[int]]]]:
-
Parameters
1.op: The string representation of operator.
2.max_op_len:order of cluster expansion(please be noted that the operator inputted will be expanded for at least one time.)
-
Return:
The information of the cluster expansion form of inputted operator.
-
Check out the source codes and the jupyter notebook files.
-
See also the published paper "Classical-to-quantum transition in multimode nonlinear systems with strong photon-photon coupling" on Physical Review A.
-
If you are interesting in the implement, check up the ImplementDetail pdfs as well as the old public archive. (Warm prompt: those files are so messy.)