diff --git a/src/Adiabatic/Adiabatic.jl b/src/Adiabatic/Adiabatic.jl
new file mode 100644
index 0000000..742066f
--- /dev/null
+++ b/src/Adiabatic/Adiabatic.jl
@@ -0,0 +1,10 @@
+################################################################################
+#
+# This file simulates an adiabatic computing
+#
+################################################################################
+
+
+include("Base.jl")
+include("Hamiltonian.jl")
+include("Operator.jl")
diff --git a/src/Adiabatic/Base.jl b/src/Adiabatic/Base.jl
new file mode 100644
index 0000000..0597c13
--- /dev/null
+++ b/src/Adiabatic/Base.jl
@@ -0,0 +1,60 @@
+type AdiaComputer <: QuComput
+    ###########################################
+    #  Eternal
+    ###########################################
+    HB::AbstractMatrix      # Base Hamiltonian
+    HP::AbstractMatrix      # Problem Hamiltonian
+    maxtime::Real           # max evolution time
+    n::Int64           # number of bits
+    dt::Real                # time step
+
+    ###########################################
+    # variable for current state
+    ###########################################
+    location::Real          # time location
+    state::AbstractVector
+    eigens::AbstractVecOrMat
+
+    ###########################################
+    # measure
+    ###########################################
+
+    prob::Real
+
+    ###########################################
+    # other
+    ###########################################
+
+    nev::Int
+
+    # n is the number of bits
+    # maxtime is the max evolution time
+    function AdiaComputer{M,N}(ins::Instance{M,N},n::Int,maxtime::Real;dt=1e-2,nev=6)
+        HB = bHamilton(n)
+        HP = pHamilton(ins,n)
+        # set time location to be the beginning
+        location = 0
+        # prepare the initial state
+        state = convert(Array{Complex,1},[1/sqrt(2^n) for i=1:2^n])
+
+        # initialize states
+        eigens = eigs(HB,nev=nev,which=:SM)[1].'
+
+        # adjust nev if the number of bits is smaller than 3
+        if n<3
+            nev = 2^n
+            warn("adjusting nev to $(nev)\n")
+        end
+
+        # probility
+        prob = norm([x==findmin(diag(HP))[2]?1:0 for x=1:2^n])^2
+
+        new(HB,HP,maxtime,n,dt,location,state,eigens,prob,nev)
+    end
+end
+
+function Hamiltonian(Hs::AdiaComputer)
+    return (1-Hs.location)*Hs.HB+Hs.location*Hs.HP
+end
+
+export AdiaComputer
diff --git a/src/Adiabatic/Cooling.jl b/src/Adiabatic/Cooling.jl
new file mode 100644
index 0000000..f155401
--- /dev/null
+++ b/src/Adiabatic/Cooling.jl
@@ -0,0 +1,69 @@
+function cooler!(Hs::AdiaComputer,gamma::Real,t::Real)
+    # boundscheck(Hs,gamma,t)
+    energy = eig(full(Hamiltonian(Hs)))[1]
+    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
+    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"
+
+    Hs.state = normalize!((Hs.state-im*exp(im*gamma)*
+            trotter(-im*t*(1-Hs.location)*Hs.HB,
+                -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
+end
+
+function heater!(Hs::AdiaComputer,gamma::Real,t::Real)
+    # boundscheck(Hs,gamma,t)
+    energy = eig(full(Hamiltonian(Hs)))[1]
+    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
+    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"
+
+    Hs.state = normalize!((Hs.state+im*exp(im*gamma)*
+                trotter(-im*t*(1-Hs.location)*Hs.HB,
+                    -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
+end
+
+function daemon!(
+    Hs::AdiaComputer,
+    gamma::Real,
+    t::Real
+    )
+    dice = rand()
+    if dice <= 0.5*(1+sin(gamma))
+        #get |1> (higher energy)
+        heater!(Hs,gamma,t)
+        return 0.5*(1-sin(gamma))
+    else
+        #get |0> (lower energy)
+        cooler!(Hs,gamma,t)
+        return 0.5*(1+sin(gamma))
+    end
+end
+
+function cooling!(
+    Hs::AdiaComputer;
+    n=5)
+
+    count = 0
+    gamma,t = CoolingPara(Hs)
+
+    while count<n
+        Hs.prob *= daemon!(Hs,gamma,t)
+        count += 1
+    end
+    return Hs
+end
+
+function CoolingPara(Hs::AdiaComputer)
+    Eigen = eig(full(Hamiltonian(Hs)))[1]
+    maxEigen = maximum(Eigen)
+    minEigen = minimum(Eigen)
+
+    gamma = (maxEigen+minEigen)/(maxEigen-minEigen) * pi/2 * 0.1
+
+    if (gamma-pi/2)>0
+        t = 0.5*((gamma-pi/2)/minEigen + (gamma+pi/2)/maxEigen)
+    else
+        t = 0.5*((gamma-pi/2)/maxEigen + (gamma+pi/2)/maxEigen)
+    end
+    return gamma,t
+end
+
+export cooling!
diff --git a/src/Adiabatic/Hamiltonian.jl b/src/Adiabatic/Hamiltonian.jl
new file mode 100644
index 0000000..fe7e3db
--- /dev/null
+++ b/src/Adiabatic/Hamiltonian.jl
@@ -0,0 +1,32 @@
+# Base Hamiltonian
+function bHamilton(bitnum::Int)
+    H = sparse(0.5*(I-σ₁))
+    Iden = spdiagm([1,1])
+    for i=2:bitnum
+        H = H⊗Iden
+    end
+
+    for i=2:bitnum
+        H_BC = Iden
+        for j=2:bitnum
+            if i==j
+                H_BC = H_BC⊗sparse(0.5*(I-σ₁))
+            else
+                H_BC = H_BC⊗Iden
+            end
+        end
+        H+=H_BC
+    end
+    return H
+end
+
+# problem Hamiltonian
+function pHamilton{M,N}(ins::Instance{M,N},n::Integer)
+    sum = spzeros(2^n,2^n);
+    for clause in ins.c
+        sum += spdiagm(Int[clause(i) for i=0:2^n-1])
+    end
+    return sum
+end
+
+export bHamilton,pHamilton
diff --git a/src/Adiabatic/Operator.jl b/src/Adiabatic/Operator.jl
new file mode 100644
index 0000000..10bbf16
--- /dev/null
+++ b/src/Adiabatic/Operator.jl
@@ -0,0 +1,2 @@
+include("timeop.jl")
+include("Cooling.jl")
diff --git a/src/Adiabatic/README.md b/src/Adiabatic/README.md
new file mode 100644
index 0000000..6d9a8c4
--- /dev/null
+++ b/src/Adiabatic/README.md
@@ -0,0 +1,11 @@
+# Adiabatic
+
+This folder is for simulations on adiabatic computing.
+
+# Contents
+
+- **Adiabatic** pack adiabatic related package in this folder
+    - **Base** basic types for adiabatic computing
+    - **Hamiltonian** Hamiltonian generators for adiabatic computing
+    - **Cooling** Daemon-like Cooling accroding to [doi:10.103](http://www.nature.com/nphoton/journal/v8/n2/full/nphoton.2013.354.html)
+    - **Operator** quantum operators realted to adiabatic computing
diff --git a/src/Adiabatic/timeop.jl b/src/Adiabatic/timeop.jl
new file mode 100644
index 0000000..ac87450
--- /dev/null
+++ b/src/Adiabatic/timeop.jl
@@ -0,0 +1,21 @@
+function realtimeop!(Hs::AdiaComputer)
+    Hs.state = trotter(-im*(1-Hs.location)*Hs.HB,-im*Hs.location*Hs.HP,3)*Hs.state
+    Hs.location += Hs.dt/Hs.maxtime
+end
+
+function next_timestep!(Hs::AdiaComputer;evopercentage::Real=1/3,nev=6)
+    @assert (0<=Hs.location+evopercentage)&&( Hs.location+evopercentage-1<0.1) "evolutoin percentage out of bounds(should be in [0,1])"
+
+    const evotime = Hs.maxtime
+    const dt      = Hs.dt
+
+    # calculate eigen values
+    for i=Hs.location*evotime:dt:(Hs.location+evopercentage)*evotime
+        realtimeop!(Hs)
+        Hs.eigens = [Hs.eigens;eigs(Hamiltonian(Hs),nev=nev,which=:SM)[1].']
+    end
+
+    return Hs
+end
+
+export next_timestep!
diff --git a/src/Circuit/Base.jl b/src/Circuit/Base.jl
new file mode 100644
index 0000000..e69de29
diff --git a/src/Circuit/Circuit.jl b/src/Circuit/Circuit.jl
new file mode 100644
index 0000000..e69de29
diff --git a/src/QuCmp.jl b/src/QuCmp.jl
index 5afc76d..011860c 100644
--- a/src/QuCmp.jl
+++ b/src/QuCmp.jl
@@ -1,5 +1,10 @@
 module QuCmp
 
-# package code goes here
+abstract QuComput
+
+include("const.jl")
+include("utils/LogicExpr.jl")
+include("utils/math.jl")
+include("Adiabatic/Adiabatic.jl")
 
 end # module
diff --git a/src/State.jl b/src/State.jl
new file mode 100644
index 0000000..fad5f3e
--- /dev/null
+++ b/src/State.jl
@@ -0,0 +1,6 @@
+# abstract AbstractState
+# type State <: AbstractState
+#     s::Vector{Complex}
+# end
+#
+# function call(state::State)
diff --git a/src/const.jl b/src/const.jl
new file mode 100644
index 0000000..49bbbbf
--- /dev/null
+++ b/src/const.jl
@@ -0,0 +1,10 @@
+const σ₀=I
+const σ₁=[0 1;1 0]
+const σ₂=[0 -im;im 0]
+const σ₃=[1 0;0 -1]
+const sigmax = [0 1;1 0]
+const sigmay = [0 -im;im 0]
+const sigmaz = [1 0;0 -1]
+const hadamard = [1/sqrt(2) 1/sqrt(2);1/sqrt(2) -1/sqrt(2)]
+
+export σ₀,σ₁,σ₂,σ₃,sigmax,sigmay,sigmaz,hadamard
diff --git a/src/utils/LogicExpr.jl b/src/utils/LogicExpr.jl
new file mode 100644
index 0000000..30c4fce
--- /dev/null
+++ b/src/utils/LogicExpr.jl
@@ -0,0 +1,208 @@
+###################################################################
+#
+# This file generates SAT's Instance in classical method for tests
+# in quantum computing.
+# Ref: arxiv:quant-ph/0007071v1
+#
+###################################################################
+
+
+
+abstract AbstractBits
+abstract AbstractClause{N}
+
+immutable Bits <: AbstractBits
+    value::UInt32
+end
+
+function call(b::Bits,n::Integer)
+    @assert n>0
+    return (b.value>>(n-1))&1
+end
+
+"""
+TruthTable Example:
+
+|bits value | truth value|
+|-----------|------------|
+| 000       | 0          |
+| 001       | 1          |
+| 010       | 0          |
+| 011       | 1          |
+| 100       | 0          |
+| 101       | 1          |
+| 110       | 0          |
+| 111       | 0          |
+
+TruthTable(0b00101010)
+"""
+immutable TruthTable
+    value::UInt32
+end
+
+function call(truthtable::TruthTable,index::Integer)
+    return (truthtable.value>>index)&1
+end
+
+immutable Clause{N} <: AbstractClause{N}
+    table::TruthTable
+    ids::AbstractVector{Integer}
+
+    function Clause(table::TruthTable,ids::AbstractVector{Integer})
+        @assert length(ids)==N
+        sort!(ids)
+        new(table,ids)
+    end
+end
+
+Clause(table::TruthTable,ids::Integer...) = Clause{length(ids)}(table,[ids...])
+Clause(num::Integer,table::TruthTable,ids::Integer...) = Clause{num}(table,[ids...])
+Clause(num::Integer,table::Integer,ids::Integer...) = Clause{num}(TruthTable(table),[ids...])
+
+function call{N}(clause::Clause{N},assign::Integer)
+    return clause.table(assign)
+end
+
+################################################################################
+#
+# Clauses in exact cover problem
+#
+################################################################################
+
+type ECClause{N} <: AbstractClause{N}
+    ids::Vector{Int}
+
+    function ECClause(ids::Vector{Int})
+        @assert length(ids) == N
+        sort!(ids)
+        new(ids)
+    end
+end
+
+ECClause(ids::Vector{Int}) = ECClause{length(ids)}(ids)
+ECClause(ids::Integer...) = ECClause{length(ids)}([ids...])
+
+function call{N}(c::ECClause{N},assign::Integer)
+    res = 0
+    for i = 1:N
+        res += assign&1
+        assign = assign>>1
+    end
+    return Int(res==1)
+end
+
+function save{N}(io::IO,clause::ECClause{N})
+    write(io,clause.ids[1])
+    for id in clause.ids[2:end]
+        write(io,"\t")
+        write(io,id)
+    end
+    write(io,"\n")
+end
+
+"""
+Instance is a collection of M clauses:
+
+~~~TeX
+C_1 Λ C_2 Λ... C_M
+~~~
+
+construct an instance:
+
+- `num` is the number of bits
+- `clauses` is the collection of clause
+"""
+type Instance{M,N}
+    c::AbstractVector{AbstractClause{N}}
+end
+
+Instance{N}(num::Integer,clauses::AbstractVector{ECClause{N}}) = Instance{num,N}(clauses)
+Instance{N}(num::Integer,clause::ECClause{N},clauses::ECClause{N}...) = Instance(num,[clause,clauses...])
+Instance{N}(num::Integer,clauses::AbstractVector{Clause{N}}) = Instance{num,N}(clauses)
+Instance{N}(num::Integer,clause::Clause{N},clauses::Clause{N}...) = Instance(num,[clause,clauses...])
+
+function call{M,N}(clauses::Instance{M,N},assign::Bits)
+    res = 1
+    for clause in clauses.c
+        assignment = 0
+        digit = 0
+        for id in clause.ids
+            assignment += assign(id)<<digit
+            digit += 1
+        end
+        res &= clause(assignment)
+    end
+    return res
+end
+
+function show{M,N}(io::IO,ins::Instance{M,N})
+    for clause in ins.c[1:end-1]
+        print(io,clause.ids)
+        print(io,"Λ")
+    end
+    print(io,ins.c[end].ids)
+    print(io,"\n")
+end
+
+show{M,N}(ins::Instance{M,N}) = show(STDOUT,ins)
+
+function save{M,N}(io::IO,ins::Instance{M,N})
+    for clause in ins.c
+        save(io,clause)
+    end
+end
+
+
+function AssignNum{M,N}(instance::Instance{M,N})
+    count = 0
+    for i = 0:2^M-1
+        if instance(Bits(i)) == 1
+            count+=1
+        end
+    end
+    return count
+end
+
+import Base: push!
+
+function push!{M,N}(instance::Instance{M,N},clause::ECClause{N})
+    push!(instance.c,clause)
+end
+
+################################################################################
+#
+# generate instance with only one assignment satisfied
+#
+################################################################################
+
+function engine(n::Integer,N::Integer,maxtry=1000)
+    list = collect(1:n)
+    ids = (list|>shuffle)[1:N]
+
+    clauses = [ECClause(ids)]
+    ins = Instance(n,clauses)
+
+    for i = 1:maxtry
+        pre_ids = ids
+        ids = (list|>shuffle)[1:N]
+        if sort!(ids) != sort!(pre_ids)
+            push!(ins,ECClause(ids))
+        end
+
+        if AssignNum(ins)<2
+            return ins
+        end
+    end
+
+    return nothing
+end
+
+function generate(n::Integer,maxiter=1000;maxbits=16,N=3)
+    for i=1:maxiter
+        t = engine(n,N)
+        t === nothing || return t
+    end
+    warn("may not have an assignment\n")
+end
+
+export AbstractBits,AbstractClause,Bits,Clause,ECClause,TruthTable,save,show,generate,call,Instance
diff --git a/src/utils/math.jl b/src/utils/math.jl
new file mode 100644
index 0000000..66a78dc
--- /dev/null
+++ b/src/utils/math.jl
@@ -0,0 +1,5 @@
+include("matrix.jl")
+
+function normalize!(vec::AbstractVector)
+    return vec[:]=vec/norm(vec)
+end
diff --git a/src/utils/matrix.jl b/src/utils/matrix.jl
new file mode 100644
index 0000000..673b917
--- /dev/null
+++ b/src/utils/matrix.jl
@@ -0,0 +1,15 @@
+#trotter expansion
+function trotter(A::AbstractMatrix,B::AbstractMatrix,P::Int64)
+    return (expm(full(A/(2*P)))*expm(full(B/P))*expm(full(A/(2*P))))^P
+end
+
+function (⊗)(A::AbstractMatrix,B::AbstractMatrix)
+    @assert size(A)[1]==size(A)[2]
+    @assert size(B)[1]==size(B)[2]
+
+    return kron(A,B)
+end
+
+function (⊕)(A::AbstractMatrix,B::AbstractMatrix)
+    return full(blkdiag(sparse(A),sparse(B)))
+end