Python secondquant.wicks函数代码示例

def test_fully_contracted():
    i, j, k, l = symbols("i j k l", below_fermi=True)
    a, b, c, d = symbols("a b c d", above_fermi=True)
    p, q, r, s = symbols("p q r s", cls=Dummy)

    Fock = AntiSymmetricTensor("f", (p,), (q,)) * NO(Fd(p) * F(q))
    V = (AntiSymmetricTensor("v", (p, q), (r, s)) * NO(Fd(p) * Fd(q) * F(s) * F(r))) / 4

    Fai = wicks(NO(Fd(i) * F(a)) * Fock, keep_only_fully_contracted=True, simplify_kronecker_deltas=True)
    assert Fai == AntiSymmetricTensor("f", (a,), (i,))
    Vabij = wicks(NO(Fd(i) * Fd(j) * F(b) * F(a)) * V, keep_only_fully_contracted=True, simplify_kronecker_deltas=True)
    assert Vabij == AntiSymmetricTensor("v", (a, b), (i, j))

def test_NO():
    i, j, k, l = symbols('i j k l', below_fermi=True)
    a, b, c, d = symbols('a b c d', above_fermi=True)
    p, q, r, s = symbols('p q r s', cls=Dummy)

    assert (NO(Fd(p)*F(q) + Fd(a)*F(b)) ==
       NO(Fd(p)*F(q)) + NO(Fd(a)*F(b)))
    assert (NO(Fd(i)*NO(F(j)*Fd(a))) ==
       NO(Fd(i)*F(j)*Fd(a)))
    assert NO(1) == 1
    assert NO(i) == i
    assert (NO(Fd(a)*Fd(b)*(F(c) + F(d))) ==
            NO(Fd(a)*Fd(b)*F(c)) +
            NO(Fd(a)*Fd(b)*F(d)))

    assert NO(Fd(a)*F(b))._remove_brackets() == Fd(a)*F(b)
    assert NO(F(j)*Fd(i))._remove_brackets() == F(j)*Fd(i)

    assert (NO(Fd(p)*F(q)).subs(Fd(p), Fd(a) + Fd(i)) ==
            NO(Fd(a)*F(q)) + NO(Fd(i)*F(q)))
    assert (NO(Fd(p)*F(q)).subs(F(q), F(a) + F(i)) ==
            NO(Fd(p)*F(a)) + NO(Fd(p)*F(i)))

    expr = NO(Fd(p)*F(q))._remove_brackets()
    assert wicks(expr) == NO(expr)

    assert NO(Fd(a)*F(b)) == - NO(F(b)*Fd(a))

    no = NO(Fd(a)*F(i)*F(b)*Fd(j))
    l1 = [ ind for ind in no.iter_q_creators() ]
    assert l1 == [0, 1]
    l2 = [ ind for ind in no.iter_q_annihilators() ]
    assert l2 == [3, 2]

def test_wicks():
    p,q,r,s = symbols('pqrs',above_fermi=True)

    # Testing for particles only

    str = F(p)*Fd(q)
    assert wicks(str) == NO(F(p)*Fd(q)) + KroneckerDelta(p,q)
    str = Fd(p)*F(q)
    assert wicks(str) == NO(Fd(p)*F(q))


    str = F(p)*Fd(q)*F(r)*Fd(s)
    nstr= wicks(str)
    fasit = NO(
    KroneckerDelta(p, q)*KroneckerDelta(r, s)
    + KroneckerDelta(p, q)*AnnihilateFermion(r)*CreateFermion(s)
    + KroneckerDelta(r, s)*AnnihilateFermion(p)*CreateFermion(q)
    - KroneckerDelta(p, s)*AnnihilateFermion(r)*CreateFermion(q)
    - AnnihilateFermion(p)*AnnihilateFermion(r)*CreateFermion(q)*CreateFermion(s))
    assert nstr == fasit

    assert (p*q*nstr).expand() == wicks(p*q*str)
    assert (nstr*p*q*2).expand() == wicks(str*p*q*2)


    # Testing CC equations particles and holes
    i,j,k,l = symbols('ijkl',below_fermi=True,dummy=True)
    a,b,c,d = symbols('abcd',above_fermi=True,dummy=True)
    p,q,r,s = symbols('pqrs',dummy=True)

    assert (wicks(F(a)*NO(F(i)*F(j))*Fd(b)) ==
            NO(F(a)*F(i)*F(j)*Fd(b)) +
            KroneckerDelta(a,b)*NO(F(i)*F(j)))
    assert (wicks(F(a)*NO(F(i)*F(j)*F(k))*Fd(b)) ==
            NO(F(a)*F(i)*F(j)*F(k)*Fd(b)) -
            KroneckerDelta(a,b)*NO(F(i)*F(j)*F(k)))


    expr = wicks(Fd(i)*NO(Fd(j)*F(k))*F(l))
    assert (expr ==
           -KroneckerDelta(i,k)*NO(Fd(j)*F(l)) -
            KroneckerDelta(j,l)*NO(Fd(i)*F(k)) -
            KroneckerDelta(i,k)*KroneckerDelta(j,l)+
            KroneckerDelta(i,l)*NO(Fd(j)*F(k)) +
            NO(Fd(i)*Fd(j)*F(k)*F(l)))
    expr = wicks(F(a)*NO(F(b)*Fd(c))*Fd(d))
    assert (expr ==
           -KroneckerDelta(a,c)*NO(F(b)*Fd(d)) -
            KroneckerDelta(b,d)*NO(F(a)*Fd(c)) -
            KroneckerDelta(a,c)*KroneckerDelta(b,d)+
            KroneckerDelta(a,d)*NO(F(b)*Fd(c)) +
            NO(F(a)*F(b)*Fd(c)*Fd(d)))

def test_NO():
    i, j, k, l = symbols('i j k l', below_fermi=True)
    a, b, c, d = symbols('a b c d', above_fermi=True)
    p, q, r, s = symbols('p q r s', cls=Dummy)

    assert (NO(Fd(p)*F(q) + Fd(a)*F(b)) ==
       NO(Fd(p)*F(q)) + NO(Fd(a)*F(b)))
    assert (NO(Fd(i)*NO(F(j)*Fd(a))) ==
       NO(Fd(i)*F(j)*Fd(a)))
    assert NO(1) == 1
    assert NO(i) == i
    assert (NO(Fd(a)*Fd(b)*(F(c) + F(d))) ==
            NO(Fd(a)*Fd(b)*F(c)) +
            NO(Fd(a)*Fd(b)*F(d)))

    assert NO(Fd(a)*F(b))._remove_brackets() == Fd(a)*F(b)
    assert NO(F(j)*Fd(i))._remove_brackets() == F(j)*Fd(i)

    assert (NO(Fd(p)*F(q)).subs(Fd(p), Fd(a) + Fd(i)) ==
            NO(Fd(a)*F(q)) + NO(Fd(i)*F(q)))
    assert (NO(Fd(p)*F(q)).subs(F(q), F(a) + F(i)) ==
            NO(Fd(p)*F(a)) + NO(Fd(p)*F(i)))

    expr = NO(Fd(p)*F(q))._remove_brackets()
    assert wicks(expr) == NO(expr)

    assert NO(Fd(a)*F(b)) == - NO(F(b)*Fd(a))

    no = NO(Fd(a)*F(i)*F(b)*Fd(j))
    l1 = [ ind for ind in no.iter_q_creators() ]
    assert l1 == [0, 1]
    l2 = [ ind for ind in no.iter_q_annihilators() ]
    assert l2 == [3, 2]
    no = NO(Fd(a)*Fd(i))
    assert no.has_q_creators == 1
    assert no.has_q_annihilators == -1
    assert str(no) == ':CreateFermion(a)*CreateFermion(i):'
    assert repr(no) == 'NO(CreateFermion(a)*CreateFermion(i))'
    assert latex(no) == r'\left\{a^\dagger_{a} a^\dagger_{i}\right\}'
    raises(NotImplementedError, lambda:  NO(Bd(p)*F(q)))

def main():
    print()
    print("Calculates the Coupled-Cluster energy- and amplitude equations")
    print("See 'An Introduction to Coupled Cluster Theory' by")
    print("T. Daniel Crawford and Henry F. Schaefer III")
    print("http://www.ccc.uga.edu/lec_top/cc/html/review.html")
    print()

    # setup hamiltonian
    p, q, r, s = symbols('p,q,r,s', cls=Dummy)
    f = AntiSymmetricTensor('f', (p,), (q,))
    pr = NO((Fd(p)*F(q)))
    v = AntiSymmetricTensor('v', (p, q), (r, s))
    pqsr = NO(Fd(p)*Fd(q)*F(s)*F(r))

    H = f*pr + Rational(1, 4)*v*pqsr
    print("Using the hamiltonian:", latex(H))

    print("Calculating 4 nested commutators")
    C = Commutator

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print("commutator 1...")
    comm1 = wicks(C(H, T))
    comm1 = evaluate_deltas(comm1)
    comm1 = substitute_dummies(comm1)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print("commutator 2...")
    comm2 = wicks(C(comm1, T))
    comm2 = evaluate_deltas(comm2)
    comm2 = substitute_dummies(comm2)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print("commutator 3...")
    comm3 = wicks(C(comm2, T))
    comm3 = evaluate_deltas(comm3)
    comm3 = substitute_dummies(comm3)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print("commutator 4...")
    comm4 = wicks(C(comm3, T))
    comm4 = evaluate_deltas(comm4)
    comm4 = substitute_dummies(comm4)

    print("construct Hausdoff expansion...")
    eq = H + comm1 + comm2/2 + comm3/6 + comm4/24
    eq = eq.expand()
    eq = evaluate_deltas(eq)
    eq = substitute_dummies(eq, new_indices=True,
            pretty_indices=pretty_dummies_dict)
    print("*********************")
    print()

    print("extracting CC equations from full Hbar")
    i, j, k, l = symbols('i,j,k,l', below_fermi=True)
    a, b, c, d = symbols('a,b,c,d', above_fermi=True)
    print()
    print("CC Energy:")
    print(latex(wicks(eq, simplify_dummies=True,
        keep_only_fully_contracted=True)))
    print()
    print("CC T1:")
    eqT1 = wicks(NO(Fd(i)*F(a))*eq, simplify_kronecker_deltas=True, keep_only_fully_contracted=True)
    eqT1 = substitute_dummies(eqT1)
    print(latex(eqT1))
    print()
    print("CC T2:")
    eqT2 = wicks(NO(Fd(i)*Fd(j)*F(b)*F(a))*eq, simplify_dummies=True, keep_only_fully_contracted=True, simplify_kronecker_deltas=True)
    P = PermutationOperator
    eqT2 = simplify_index_permutations(eqT2, [P(a, b), P(i, j)])
    print(latex(eqT2))

print "Setting up hamiltonian"
p,q,r,s = symbols('pqrs',dummy=True)
f = AntiSymmetricTensor('f',(p,),(q,))
pr = NO((Fd(p)*F(q)))
v = AntiSymmetricTensor('v',(p,q),(r,s))
pqsr = NO(Fd(p)*Fd(q)*F(s)*F(r))

H=f*pr +Number(1,4)*v*pqsr

print "Calculating nested commutators"
C = Commutator

T1,T2 = get_CC_operators()
T = T1+ T2
print "comm1..."
comm1 = wicks(C(H,T),simplify_dummies=True, simplify_kronecker_deltas=True)

T1,T2 = get_CC_operators()
T = T1+ T2
print "comm2..."
comm2 = wicks(C(comm1,T),simplify_dummies=True, simplify_kronecker_deltas=True)

T1,T2 = get_CC_operators()
T = T1+ T2
print "comm3..."
comm3 = wicks(C(comm2,T),simplify_dummies=True, simplify_kronecker_deltas=True)

T1,T2 = get_CC_operators()
T = T1+ T2
print "comm4..."
comm4 = wicks(C(comm3,T),simplify_dummies=True, simplify_kronecker_deltas=True)

def main():
    print
    print "Calculates the Coupled-Cluster energy- and amplitude equations"
    print "See 'An Introduction to Coupled Cluster Theory' by"
    print "T. Daniel Crawford and Henry F. Schaefer III"
    print "http://www.ccc.uga.edu/lec_top/cc/html/review.html"
    print

    # setup hamiltonian
    p, q, r, s = symbols("pqrs", dummy=True)
    f = AntiSymmetricTensor("f", (p,), (q,))
    pr = NO((Fd(p) * F(q)))
    v = AntiSymmetricTensor("v", (p, q), (r, s))
    pqsr = NO(Fd(p) * Fd(q) * F(s) * F(r))

    H = f * pr

    # Uncomment the next line to use a 2-body hamiltonian:
    # H=f*pr + Number(1,4)*v*pqsr

    print "Using the hamiltonian:", latex(H)

    print "Calculating nested commutators"
    C = Commutator

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print "comm1..."
    comm1 = wicks(C(H, T), simplify_dummies=True, simplify_kronecker_deltas=True)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print "comm2..."
    comm2 = wicks(C(comm1, T), simplify_dummies=True, simplify_kronecker_deltas=True)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print "comm3..."
    comm3 = wicks(C(comm2, T), simplify_dummies=True, simplify_kronecker_deltas=True)

    T1, T2 = get_CC_operators()
    T = T1 + T2
    print "comm4..."
    comm4 = wicks(C(comm3, T), simplify_dummies=True, simplify_kronecker_deltas=True)

    print "construct Hausdoff expansion..."
    eq = H + comm1 + comm2 / 2 + comm3 / 6 + comm4 / 24
    eq = eq.expand()
    eq = evaluate_deltas(eq)
    eq = substitute_dummies(eq, new_indices=True, reverse_order=False, pretty_indices=pretty_dummies_dict)
    print "*********************"
    print

    print "extracting CC equations from full Hbar"
    i, j, k, l = symbols("ijkl", below_fermi=True)
    a, b, c, d = symbols("abcd", above_fermi=True)
    print
    print "CC Energy:"
    print latex(wicks(eq, simplify_dummies=True, keep_only_fully_contracted=True))
    print
    print "CC T1:"
    eqT1 = wicks(NO(Fd(i) * F(a)) * eq, simplify_kronecker_deltas=True, keep_only_fully_contracted=True)
    eqT1 = substitute_dummies(eqT1, reverse_order=False)
    print latex(eqT1)
    print
    print "CC T2:"
    eqT2 = wicks(
        NO(Fd(i) * Fd(j) * F(b) * F(a)) * eq,
        simplify_dummies=True,
        keep_only_fully_contracted=True,
        simplify_kronecker_deltas=True,
    )
    P = PermutationOperator
    eqT2 = simplify_index_permutations(eqT2, [P(a, b), P(i, j)])
    print latex(eqT2)