diff --git a/SeQuant/core/context.cpp b/SeQuant/core/context.cpp index c7f63c69cf..cba9b80121 100644 --- a/SeQuant/core/context.cpp +++ b/SeQuant/core/context.cpp @@ -31,6 +31,9 @@ bool operator==(const Context& ctx1, const Context& ctx2) { ctx1.canonicalization_options() == ctx2.canonicalization_options() && ctx1.braket_typesetting() == ctx2.braket_typesetting() && ctx1.braket_slot_typesetting() == ctx2.braket_slot_typesetting() && + ctx1.symmetry() == ctx2.symmetry() && + ctx1.hermiticity() == ctx2.hermiticity() && + ctx1.column_symmetry() == ctx2.column_symmetry() && *ctx1.index_space_registry() == *ctx2.index_space_registry(); } @@ -132,7 +135,10 @@ Context::Context(Options options) first_dummy_index_ordinal_(options.first_dummy_index_ordinal), canonicalization_options_(options.canonicalization_options), braket_typesetting_(options.braket_typesetting), - braket_slot_typesetting_(options.braket_slot_typesetting) {} + braket_slot_typesetting_(options.braket_slot_typesetting), + symmetry_(options.symmetry), + hermiticity_(options.hermiticity), + column_symmetry_(options.column_symmetry) {} Context Context::clone() const { Context ctx(*this); @@ -176,6 +182,12 @@ BraKetSlotTypesetting Context::braket_slot_typesetting() const { return braket_slot_typesetting_; } +Symmetry Context::symmetry() const { return symmetry_; } + +Hermiticity Context::hermiticity() const { return hermiticity_; } + +ColumnSymmetry Context::column_symmetry() const { return column_symmetry_; } + Context& Context::set(Vacuum vacuum) { vacuum_ = vacuum; return *this; @@ -229,6 +241,21 @@ Context& Context::set(BraKetSlotTypesetting bkst) { return *this; } +Context& Context::set(Symmetry symmetry) { + symmetry_ = symmetry; + return *this; +} + +Context& Context::set(Hermiticity hermiticity) { + hermiticity_ = hermiticity; + return *this; +} + +Context& Context::set(ColumnSymmetry column_symmetry) { + column_symmetry_ = column_symmetry; + return *this; +} + IndexSpace get_particle_space(const IndexSpace::QuantumNumbers& qn) { return get_default_context().index_space_registry()->particle_space(qn); } diff --git a/SeQuant/core/context.hpp b/SeQuant/core/context.hpp index 4c3ee9641d..6b1da44d4c 100644 --- a/SeQuant/core/context.hpp +++ b/SeQuant/core/context.hpp @@ -22,13 +22,11 @@ namespace sequant { /// orthonormal to their dual (bra) counterparts /// (`IndexSpaceMetric::Unit`) or not (`IndexSpaceMetric::General`); /// this affects the value of Wick contractions. -/// - `braket_symmetry`: whether the primal (ket) and dual (bra) vector space -/// _bases_ -/// are "equivalent" (homogenous to each other; `BraKetSymmetry::Symm`), -/// "conjugate to each other" (conjugate-homogenous -/// to each other; `BraKetSymmetry::Conjugate`) or are "nonequivalent" (not -/// homogeneous; `BraKetSymmetry::Nonsymm`) +/// - `symmetry`, `hermiticity`, `column_symmetry`: the symmetries given to a +/// *deserialized* tensor that does not specify them; the programmatic +/// Tensor ctors are unaffected (see Tensor::Defaults). There is no +/// `braket_symmetry` knob: a tensor's BraKetSymmetry is derived from its +/// `hermiticity` and base field. /// - `spbasis`: whether the bra/ket bases are spinor (`SPBasis::Spinor`) or /// spin-free (`SPBasis::Spinfree`). /// - `first_dummy_index_ordinal`: during its operation SeQuant will generate @@ -58,6 +56,9 @@ class Context { constexpr static auto braket_typesetting = BraKetTypesetting::ContraSub; constexpr static auto braket_slot_typesetting = BraKetSlotTypesetting::TensorPackage; + constexpr static auto symmetry = Symmetry::Nonsymm; + constexpr static auto hermiticity = Hermiticity::NonHermitian; + constexpr static auto column_symmetry = ColumnSymmetry::Nonsymm; }; /// helper for the named-parameter constructor of Context @@ -86,6 +87,14 @@ class Context { /// the BraKetSlotTypesetting object BraKetSlotTypesetting braket_slot_typesetting = Defaults::braket_slot_typesetting; + /// the default bra/ket permutational Symmetry for deserialized tensors + Symmetry symmetry = Defaults::symmetry; + /// the default Hermiticity for deserialized tensors; the braket symmetry + /// of a deserialized tensor is *derived* from this and its #base_field + Hermiticity hermiticity = Defaults::hermiticity; + /// the default ColumnSymmetry (particle-permutation symmetry) for + /// deserialized tensors + ColumnSymmetry column_symmetry = Defaults::column_symmetry; }; static Options make_default_options() { return {}; } @@ -155,6 +164,15 @@ class Context { /// \return BraKetSlotTypesetting of this context; see BraKetSlotTypesetting /// for the meaning of the possible values BraKetSlotTypesetting braket_slot_typesetting() const; + /// \return the default bra/ket permutational Symmetry for deserialized tensors + Symmetry symmetry() const; + /// \return the default Hermiticity for deserialized tensors; the braket + /// symmetry of a deserialized tensor is *derived* from this and its + /// #base_field + Hermiticity hermiticity() const; + /// \return the default ColumnSymmetry (particle-permutation symmetry) for + /// deserialized tensors + ColumnSymmetry column_symmetry() const; /// Sets the Vacuum for this context, convenient for chaining /// \param vacuum Vacuum @@ -195,6 +213,16 @@ class Context { /// \param braket_slot_typeset BraKetSlotTypesetting /// \return ref to `*this`, for chaining Context& set(BraKetSlotTypesetting braket_slot_typeset); + /// Sets the default bra/ket permutational Symmetry for deserialized tensors + /// \return ref to `*this`, for chaining + Context& set(Symmetry symmetry); + /// Sets the default Hermiticity for deserialized tensors (the braket symmetry + /// of a deserialized tensor is derived from this and its base field) + /// \return ref to `*this`, for chaining + Context& set(Hermiticity hermiticity); + /// Sets the default ColumnSymmetry for deserialized tensors + /// \return ref to `*this`, for chaining + Context& set(ColumnSymmetry column_symmetry); private: std::shared_ptr idx_space_reg_ = nullptr; @@ -207,6 +235,9 @@ class Context { BraKetTypesetting braket_typesetting_ = Defaults::braket_typesetting; BraKetSlotTypesetting braket_slot_typesetting_ = Defaults::braket_slot_typesetting; + Symmetry symmetry_ = Defaults::symmetry; + Hermiticity hermiticity_ = Defaults::hermiticity; + ColumnSymmetry column_symmetry_ = Defaults::column_symmetry; }; /// Context object equality comparison diff --git a/SeQuant/core/expressions/tensor.hpp b/SeQuant/core/expressions/tensor.hpp index bc74cb7dcd..bd59e01e21 100644 --- a/SeQuant/core/expressions/tensor.hpp +++ b/SeQuant/core/expressions/tensor.hpp @@ -54,9 +54,57 @@ DEFINE_STRONG_TYPE_FOR_RANGE_AND_RANGESIZE(ket); // strong type wrapper for objects associated with aux DEFINE_STRONG_TYPE_FOR_RANGE_AND_RANGESIZE(aux); +/// @brief named-parameter pack of a Tensor's symmetry attributes +/// +/// Passed to the Tensor constructors as a single argument (the named-parameter +/// idiom via designated initializers) instead of the order-sensitive positional +/// `Symmetry`/`BraKetSymmetry`/`ColumnSymmetry` arguments. Each field is +/// optional; an unset one falls back to the *fixed* library default (see +/// Tensor::resolve_symmetries). This makes it easy to define one symmetry pack +/// for a translation unit (e.g. `constexpr TensorSymmetries +/// particle_symmetric{.column = ColumnSymmetry::Symm};`) and reuse it at every +/// construction site. +/// +/// @note #braket and #hermiticity are two ways to specify the same underlying +/// trait: if #braket is set it is used directly (and, unless #hermiticity +/// is also set, the Hermiticity is back-filled from it); otherwise the +/// braket symmetry is *derived* from #hermiticity and the tensor's base +/// field. Prefer #hermiticity for field-agnostic physical facts. +struct TensorSymmetries { + /// bra/ket permutational symmetry + std::optional perm = std::nullopt; + /// bra<->ket exchange symmetry (an alternative to #hermiticity; see note) + std::optional braket = std::nullopt; + /// abstract (field-agnostic) adjoint symmetry + std::optional hermiticity = std::nullopt; + /// particle (column) permutation symmetry + std::optional column = std::nullopt; +}; + /// @brief a Tensor is an instance of AbstractTensor over a scalar field, i.e. /// Tensors have commutative addition and product operations class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { + public: + /// @brief the fixed library defaults used to resolve the symmetry attributes + /// a Tensor constructor leaves unspecified + /// + /// Every programmatic constructor resolves an omitted `std::optional` + /// symmetry argument (or an unset TensorSymmetries field) against these; see + /// resolve_symmetries(). They are *independent of the active + /// sequant::Context* -- unlike deserialization, whose corresponding defaults + /// come from Context. There is deliberately no BraKetSymmetry default: it is + /// *derived* from #hermiticity and the tensor's #base_field via + /// to_braket_symmetry(). + struct Defaults { + /// default bra/ket permutational Symmetry + static constexpr Symmetry symmetry = Symmetry::Nonsymm; + /// default (field-agnostic) Hermiticity; the observable BraKetSymmetry is + /// derived from this and the tensor's base field + static constexpr Hermiticity hermiticity = Hermiticity::NonHermitian; + /// default particle-exchange ColumnSymmetry + static constexpr ColumnSymmetry column_symmetry = ColumnSymmetry::Nonsymm; + }; + private: using index_container_type = container::svector; static auto make_indices(IndexList indices) { return indices; } @@ -255,56 +303,132 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { } } + /// the four symmetry attributes of a Tensor, fully resolved (no defaulting + /// left to do) + struct ResolvedSymmetries { + Symmetry symmetry; + BraKetSymmetry braket_symmetry; + Hermiticity hermiticity; + ColumnSymmetry column_symmetry; + /// whether #column_symmetry was spelled out by the caller rather than + /// defaulted; see check_symmetries() + bool column_symmetry_specified; + }; + + /// Resolves the (possibly unspecified) symmetry attributes of a tensor. + /// - #BraKetSymmetry is a *derived* property: if not given explicitly it is + /// computed from the (explicit or default) #Hermiticity and the tensor's + /// @p base_fld via to_braket_symmetry(). + /// - #Symmetry, #Hermiticity and #ColumnSymmetry fall back to the *fixed* + /// library defaults in Tensor::Defaults -- the safest, fully non-symmetric + /// / non-Hermitian choice -- when not specified. + /// @param syms the (partially specified) symmetry pack + /// @param base_fld the tensor's #base_field, used to derive the + /// #BraKetSymmetry from the #Hermiticity + /// @note Programmatic Tensor construction never consults the default + /// sequant::Context: the meaning of a ctor call is independent of + /// ambient global state (so it is predictable and lock-free). Only the + /// deserializer honors the Context-level symmetry defaults, since it is + /// the boundary that turns under-specified external input into tensors; + /// domains that need a different programmatic default (e.g. mbpt, whose + /// tensors are particle/column-symmetric) pass it explicitly. + static ResolvedSymmetries resolve_symmetries(const TensorSymmetries &syms, + Field base_fld) { + const Symmetry s_resolved = syms.perm.value_or(Defaults::symmetry); + const ColumnSymmetry ps_resolved = + syms.column.value_or(Defaults::column_symmetry); + const Hermiticity h_resolved = + syms.hermiticity.value_or(Defaults::hermiticity); + const BraKetSymmetry bks_resolved = + syms.braket.has_value() ? *syms.braket + : to_braket_symmetry(h_resolved, base_fld); + // an explicit Hermiticity is reported verbatim, to preserve traits that + // the BraKetSymmetry round-trip cannot represent (e.g. AntiHermitian) + const Hermiticity hermiticity_resolved = + syms.hermiticity.has_value() + ? *syms.hermiticity + : (syms.braket.has_value() ? to_hermiticity(*syms.braket) + : h_resolved); + return {s_resolved, bks_resolved, hermiticity_resolved, ps_resolved, + syms.column.has_value()}; + } + + // fully-resolved terminal ctor (range form) template Tensor(S &&label, const bra &bra_indices, const ket &ket_indices, const aux &aux_indices, reserved_tag, - Symmetry s = Symmetry::Nonsymm, - std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + ResolvedSymmetries rsym) : label_(toUtf16(std::forward(label))), bra_(make_indices(bra_indices)), ket_(make_indices(ket_indices)), aux_(make_indices(aux_indices)), - symmetry_(s), - braket_symmetry_(bks_opt.value_or(BraKetSymmetry::Nonsymm)), - hermiticity_(to_hermiticity(braket_symmetry_)), - column_symmetry_(ps), + symmetry_(rsym.symmetry), + braket_symmetry_(rsym.braket_symmetry), + hermiticity_(rsym.hermiticity), + column_symmetry_(rsym.column_symmetry), bra_net_rank_(ranges::count_if( bra_, [](const Index &idx) { return static_cast(idx); })), ket_net_rank_(ranges::count_if( ket_, [](const Index &idx) { return static_cast(idx); })) { validate_indices(); - check_symmetries(); + check_symmetries(rsym.column_symmetry_specified); canonicalize_slots(); } + // fully-resolved terminal ctor (move form) template Tensor(S &&label, bra &&bra_indices, ket &&ket_indices, aux &&aux_indices, reserved_tag, - Symmetry s = Symmetry::Nonsymm, - std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + ResolvedSymmetries rsym) : label_(toUtf16(std::forward(label))), bra_(std::move(bra_indices)), ket_(std::move(ket_indices)), aux_(std::move(aux_indices)), - symmetry_(s), - braket_symmetry_(bks_opt.value_or(BraKetSymmetry::Nonsymm)), - hermiticity_(to_hermiticity(braket_symmetry_)), - column_symmetry_(ps), + symmetry_(rsym.symmetry), + braket_symmetry_(rsym.braket_symmetry), + hermiticity_(rsym.hermiticity), + column_symmetry_(rsym.column_symmetry), bra_net_rank_(ranges::count_if( bra_, [](const Index &idx) { return static_cast(idx); })), ket_net_rank_(ranges::count_if( ket_, [](const Index &idx) { return static_cast(idx); })) { validate_indices(); - check_symmetries(); + check_symmetries(rsym.column_symmetry_specified); canonicalize_slots(); } + // defaulting reserved-tag ctor (range form) + template + Tensor(S &&label, const bra &bra_indices, + const ket &ket_indices, + const aux &aux_indices, reserved_tag, + const TensorSymmetries &syms = {}) + : Tensor(std::forward(label), bra_indices, ket_indices, aux_indices, + reserved_tag{}, + resolve_symmetries( + syms, sequant::base_field(make_indices(bra_indices), + make_indices(ket_indices)))) {} + + // defaulting reserved-tag ctor (move form) + template + Tensor(S &&label, bra &&bra_indices, + ket &&ket_indices, + aux &&aux_indices, reserved_tag, + const TensorSymmetries &syms = {}) + // base_field() takes its arguments by value and the indices are moved + // only in the delegated ctor's member-init list, so this reads them + // intact + : Tensor(std::forward(label), std::move(bra_indices), + std::move(ket_indices), std::move(aux_indices), reserved_tag{}, + resolve_symmetries( + syms, sequant::base_field(bra_indices, ket_indices))) {} + public: /// constructs an uninitialized Tensor /// @sa Tensor::operator bool() @@ -328,17 +452,24 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { /// to indices) /// @param ket_indices list of ket indices (or objects that can be converted /// to indices) - /// @param s the symmetry of bra or ket - /// @param bks the symmetry with respect to bra-ket exchange - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry; if unset (`std::nullopt`), + /// defaults to Tensor::Defaults::symmetry (#Symmetry::Nonsymm) + /// @param bks_opt symmetry under bra-ket exchange; if unset (`std::nullopt`) + /// it is *derived* from the default Tensor::Defaults::hermiticity + /// (#Hermiticity::NonHermitian), i.e. #BraKetSymmetry::Nonsymm + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, const bra &bra_indices, - const ket &ket_indices, Symmetry s = Symmetry::Nonsymm, + const ket &ket_indices, + std::optional s = std::nullopt, std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + std::optional ps = std::nullopt) : Tensor(std::forward(label), bra_indices, ket_indices, sequant::aux{}, - reserved_tag{}, s, bks_opt, ps) { + reserved_tag{}, + TensorSymmetries{.perm = s, .braket = bks_opt, .column = ps}) { assert_nonreserved_label(label_); } @@ -349,19 +480,26 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { /// to indices) /// @param aux_indices list of aux indices (or objects that can be /// converted to indices) - /// @param s the symmetry of bra or ket - /// @param bks the symmetry with respect to bra-ket exchange - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry; if unset (`std::nullopt`), + /// defaults to Tensor::Defaults::symmetry (#Symmetry::Nonsymm) + /// @param bks_opt symmetry under bra-ket exchange; if unset (`std::nullopt`) + /// it is *derived* from the default Tensor::Defaults::hermiticity + /// (#Hermiticity::NonHermitian), i.e. #BraKetSymmetry::Nonsymm + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, const bra &bra_indices, const ket &ket_indices, - const aux &aux_indices, Symmetry s = Symmetry::Nonsymm, + const aux &aux_indices, + std::optional s = std::nullopt, std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + std::optional ps = std::nullopt) : Tensor(std::forward(label), bra_indices, ket_indices, aux_indices, - reserved_tag{}, s, bks_opt, ps) { + reserved_tag{}, + TensorSymmetries{.perm = s, .braket = bks_opt, .column = ps}) { assert_nonreserved_label(label_); } @@ -370,18 +508,23 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { /// to indices) /// @param ket_indices list of ket indices (or objects that can be converted /// to indices) - /// @param s the symmetry of bra or ket - /// @param bks the symmetry with respect to bra-ket exchange - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry; if unset (`std::nullopt`), + /// defaults to Tensor::Defaults::symmetry (#Symmetry::Nonsymm) + /// @param bks_opt symmetry under bra-ket exchange; if unset (`std::nullopt`) + /// it is *derived* from the default Tensor::Defaults::hermiticity + /// (#Hermiticity::NonHermitian), i.e. #BraKetSymmetry::Nonsymm + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, bra &&bra_indices, ket &&ket_indices, - Symmetry s = Symmetry::Nonsymm, + std::optional s = std::nullopt, std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + std::optional ps = std::nullopt) : Tensor(std::forward(label), std::move(bra_indices), - std::move(ket_indices), sequant::aux{}, reserved_tag{}, s, - bks_opt, ps) { + std::move(ket_indices), sequant::aux{}, reserved_tag{}, + TensorSymmetries{.perm = s, .braket = bks_opt, .column = ps}) { assert_nonreserved_label(label_); } @@ -392,19 +535,24 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { /// to indices) /// @param aux_indices list of aux indices (or objects that can be /// converted to indices) - /// @param s the symmetry of bra or ket - /// @param bks the symmetry with respect to bra-ket exchange - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry; if unset (`std::nullopt`), + /// defaults to Tensor::Defaults::symmetry (#Symmetry::Nonsymm) + /// @param bks_opt symmetry under bra-ket exchange; if unset (`std::nullopt`) + /// it is *derived* from the default Tensor::Defaults::hermiticity + /// (#Hermiticity::NonHermitian), i.e. #BraKetSymmetry::Nonsymm + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, bra &&bra_indices, ket &&ket_indices, aux &&aux_indices, - Symmetry s = Symmetry::Nonsymm, + std::optional s = std::nullopt, std::optional bks_opt = std::nullopt, - ColumnSymmetry ps = ColumnSymmetry::Symm) + std::optional ps = std::nullopt) : Tensor(std::forward(label), std::move(bra_indices), std::move(ket_indices), std::move(aux_indices), reserved_tag{}, - s, bks_opt, ps) { + TensorSymmetries{.perm = s, .braket = bks_opt, .column = ps}) { assert_nonreserved_label(label_); } @@ -420,66 +568,113 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { /// @param label the tensor label /// @param bra_indices list of bra indices /// @param ket_indices list of ket indices - /// @param s the symmetry of bra or ket - /// @param h the abstract symmetry under (Hermitian) adjoint - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry + /// @param h the abstract symmetry under (Hermitian) adjoint; the observable + /// #BraKetSymmetry is derived from this and the tensor's #base_field + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, const bra &bra_indices, const ket &ket_indices, Symmetry s, Hermiticity h, - ColumnSymmetry ps = ColumnSymmetry::Symm) - // The base_field must be resolved into a BraKetSymmetry *here*, in the - // delegation, because the delegated-to ctor's body runs - // canonicalize_slots() (which keys off braket_symmetry_) -- it cannot be - // fixed up afterwards. Hence make_indices is unavoidably evaluated here - // too (the delegated ctor materializes bra_/ket_ from the same ranges - // again); the duplication is the cost of safe delegation, not an - // oversight. - : Tensor(std::forward(label), bra_indices, ket_indices, s, - to_braket_symmetry( - h, sequant::base_field(make_indices(bra_indices), - make_indices(ket_indices))), - ps) { - // Overwrite after delegation to preserve the exact trait (incl. - // AntiHermitian, which the BraKetSymmetry round-trip cannot represent). - // Stays consistent with the reserved-(anti)symmetrizer Symm->Conjugate - // demotion in the delegated ctor: that demotion only adjusts - // braket_symmetry_ (leaving hermiticity Hermitian), which h also implies. - hermiticity_ = h; + std::optional ps = std::nullopt) + // braket_symmetry_ must be derived from h before canonicalize_slots() + // runs, so let the reserved-tag ctor do it + : Tensor(std::forward(label), bra_indices, ket_indices, sequant::aux{}, + reserved_tag{}, + TensorSymmetries{.perm = s, .hermiticity = h, .column = ps}) { + assert_nonreserved_label(label_); } /// @param label the tensor label /// @param bra_indices list of bra indices /// @param ket_indices list of ket indices /// @param aux_indices list of aux indices - /// @param s the symmetry of bra or ket - /// @param h the abstract symmetry under (Hermitian) adjoint - /// @param ps the symmetry under exchange of particles + /// @param s bra/ket permutational symmetry + /// @param h the abstract symmetry under (Hermitian) adjoint; the observable + /// #BraKetSymmetry is derived from this and the tensor's #base_field + /// @param ps particle- (column-) exchange symmetry; if unset + /// (`std::nullopt`), defaults to Tensor::Defaults::column_symmetry + /// (#ColumnSymmetry::Nonsymm) template Tensor(S &&label, const bra &bra_indices, const ket &ket_indices, const aux &aux_indices, Symmetry s, Hermiticity h, - ColumnSymmetry ps = ColumnSymmetry::Symm) - // The base_field must be resolved into a BraKetSymmetry *here*, in the - // delegation, because the delegated-to ctor's body runs - // canonicalize_slots() (which keys off braket_symmetry_) -- it cannot be - // fixed up afterwards. Hence make_indices is unavoidably evaluated here - // too (the delegated ctor materializes bra_/ket_ from the same ranges - // again); the duplication is the cost of safe delegation, not an - // oversight. - : Tensor(std::forward(label), bra_indices, ket_indices, aux_indices, s, - to_braket_symmetry( - h, sequant::base_field(make_indices(bra_indices), - make_indices(ket_indices))), - ps) { - // Overwrite after delegation to preserve the exact trait (incl. - // AntiHermitian, which the BraKetSymmetry round-trip cannot represent). - // Stays consistent with the reserved-(anti)symmetrizer Symm->Conjugate - // demotion in the delegated ctor: that demotion only adjusts - // braket_symmetry_ (leaving hermiticity Hermitian), which h also implies. - hermiticity_ = h; + std::optional ps = std::nullopt) + // braket_symmetry_ must be derived from h before canonicalize_slots() + // runs, so let the reserved-tag ctor do it + : Tensor(std::forward(label), bra_indices, ket_indices, aux_indices, + reserved_tag{}, + TensorSymmetries{.perm = s, .hermiticity = h, .column = ps}) { + assert_nonreserved_label(label_); + } + + /// @} + + /// @name named-parameter symmetry constructors + /// Specify the tensor's symmetries as a single TensorSymmetries pack (the + /// named-parameter idiom) rather than positional arguments; unset fields fall + /// back to the fixed library defaults. Convenient for reusing one symmetry + /// pack (e.g. a per-TU default) across many construction sites. + /// @{ + + /// @param label the tensor label + /// @param bra_indices list of bra indices + /// @param ket_indices list of ket indices + /// @param syms the (partially specified) symmetry pack + template + Tensor(S &&label, const bra &bra_indices, + const ket &ket_indices, TensorSymmetries syms) + : Tensor(std::forward(label), bra_indices, ket_indices, sequant::aux{}, + reserved_tag{}, syms) { + assert_nonreserved_label(label_); + } + + /// @param label the tensor label + /// @param bra_indices list of bra indices + /// @param ket_indices list of ket indices + /// @param aux_indices list of aux indices + /// @param syms the (partially specified) symmetry pack + template + Tensor(S &&label, const bra &bra_indices, + const ket &ket_indices, + const aux &aux_indices, TensorSymmetries syms) + : Tensor(std::forward(label), bra_indices, ket_indices, aux_indices, + reserved_tag{}, syms) { + assert_nonreserved_label(label_); + } + + /// @param label the tensor label + /// @param bra_indices list of bra indices + /// @param ket_indices list of ket indices + /// @param syms the (partially specified) symmetry pack + template + Tensor(S &&label, bra &&bra_indices, + ket &&ket_indices, TensorSymmetries syms) + : Tensor(std::forward(label), std::move(bra_indices), + std::move(ket_indices), sequant::aux{}, reserved_tag{}, syms) { + assert_nonreserved_label(label_); + } + + /// @param label the tensor label + /// @param bra_indices list of bra indices + /// @param ket_indices list of ket indices + /// @param aux_indices list of aux indices + /// @param syms the (partially specified) symmetry pack + template + Tensor(S &&label, bra &&bra_indices, + ket &&ket_indices, + aux &&aux_indices, TensorSymmetries syms) + : Tensor(std::forward(label), std::move(bra_indices), + std::move(ket_indices), std::move(aux_indices), reserved_tag{}, + syms) { + assert_nonreserved_label(label_); } /// @} @@ -743,25 +938,37 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { std::size_t bra_net_rank_; std::size_t ket_net_rank_; - void check_symmetries() { + /// @param column_symmetry_specified whether the column symmetry was spelled + /// out by the caller rather than defaulted; a *specified* value that + /// contradicts a reserved label's defining symmetry is an error, an + /// unspecified one is simply replaced by the correct value + /// @sa make_symmetrizer(), make_antisymmetrizer() + void check_symmetries(bool column_symmetry_specified) { + if (symmetry_ == Symmetry::Symm || symmetry_ == Symmetry::Antisymm) { + // (Anti)symmetry in bra and ket indices automatically implies column + // symmetry. N.B. must run before the reserved-label check below, which + // would otherwise report the (Antisymm) antisymmetrizer as contradicting. + column_symmetry_ = ColumnSymmetry::Symm; + } + // The (anti)symmetrizer is a permutation-bookkeeping operator whose // bra<->ket orientation defines/extracts the external indices and must be // preserved; it must never be bra<->ket-symmetric (Symm), or // canonicalization would swap its bra and ket and corrupt external-index - // extraction. Over a real field a Hermitian operator becomes Symm, so - // demote Symm to Conjugate here (Conjugate is treated as no-swap by the - // canonicalizer, matching the complex-field behavior) before - // canonicalize_slots() may act on it. - if ((label_ == reserved::antisymm_label() || - label_ == reserved::symm_label()) && - braket_symmetry_ != BraKetSymmetry::Nonsymm) { - throw Exception( - "(Anti)symmetrization operators must not have braket symmetry"); - } - - if (symmetry_ == Symmetry::Symm || symmetry_ == Symmetry::Antisymm) { - // (Anti)symmetry in bra and ket indices automatically implies column - // symmetry + // extraction. + if (label_ == reserved::antisymm_label() || + label_ == reserved::symm_label()) { + if (braket_symmetry_ != BraKetSymmetry::Nonsymm) + throw Exception( + "(Anti)symmetrization operators must not have braket symmetry"); + // (Anti)symmetrization operators act on indistinguishable particles, + // hence are always column-symmetric; supplying it when unspecified makes + // them compare equal however they were built (a mismatch would silently + // prevent otherwise-equal terms from cancelling). + if (column_symmetry_specified && column_symmetry_ != ColumnSymmetry::Symm) + throw Exception( + "(Anti)symmetrization operators must be column (particle) " + "symmetric"); column_symmetry_ = ColumnSymmetry::Symm; } } @@ -941,15 +1148,100 @@ static_assert(is_tensor, using TensorPtr = std::shared_ptr; inline ExprPtr make_overlap(const Index &bra_index, const Index &ket_index) { + // overlap tensors are particle (column) symmetric by convention, and must + // compare equal to the same overlap obtained by deserialization return ex(Tensor(reserved::overlap_label(), bra{bra_index}, - ket{ket_index}, aux{}, Tensor::reserved_tag{})); + ket{ket_index}, aux{}, Tensor::reserved_tag{}, + TensorSymmetries{.perm = Symmetry::Nonsymm, + .column = ColumnSymmetry::Symm})); } inline ExprPtr make_kronecker(const Index &bra_index, const Index &ket_index) { + // as in make_overlap, a Kronecker delta is particle (column) symmetric return ex(Tensor(reserved::kronecker_label(), bra{bra_index}, - ket{ket_index}, aux{}, Tensor::reserved_tag{})); + ket{ket_index}, aux{}, Tensor::reserved_tag{}, + TensorSymmetries{.perm = Symmetry::Nonsymm, + .column = ColumnSymmetry::Symm})); +} + +/// @name (anti)symmetrization operator factories +/// The reserved (anti)symmetrization operators Ŝ/ have *defining* symmetries +/// that are not free parameters: both are braket-Nonsymm (their bra<->ket +/// orientation defines which indices are external) and column- (particle-) +/// symmetric (they act on indistinguishable particles), and  is additionally +/// Antisymm in bra and in ket. These factories spell those out, so that a +/// (anti)symmetrizer built programmatically is identical to one obtained by +/// deserialization; the Tensor ctors reject any contradicting symmetry. +/// @{ + +/// the defining symmetries of the reserved antisymmetrization operator  +inline constexpr TensorSymmetries antisymmetrizer_symmetries{ + .perm = Symmetry::Antisymm, + .braket = BraKetSymmetry::Nonsymm, + .column = ColumnSymmetry::Symm}; + +/// the defining symmetries of the reserved symmetrization operator Ŝ +inline constexpr TensorSymmetries symmetrizer_symmetries{ + .perm = Symmetry::Nonsymm, + .braket = BraKetSymmetry::Nonsymm, + .column = ColumnSymmetry::Symm}; + +/// @brief makes an antisymmetrization operator  +/// @param bra_indices the bra ("upper") external indices +/// @param ket_indices the ket ("lower") external indices +/// @return an ExprPtr to the  Tensor +template +ExprPtr make_antisymmetrizer(const bra &bra_indices, + const ket &ket_indices) { + return ex(Tensor(reserved::antisymm_label(), bra_indices, ket_indices, + antisymmetrizer_symmetries)); +} + +/// @brief makes an antisymmetrization operator  +/// @param bra_indices the bra ("upper") external indices +/// @param ket_indices the ket ("lower") external indices +/// @param aux_indices the aux indices +/// @return an ExprPtr to the  Tensor +template +ExprPtr make_antisymmetrizer(const bra &bra_indices, + const ket &ket_indices, + const aux &aux_indices) { + return ex(Tensor(reserved::antisymm_label(), bra_indices, ket_indices, + aux_indices, antisymmetrizer_symmetries)); } +/// @brief makes a symmetrization operator Ŝ +/// @param bra_indices the bra ("upper") external indices +/// @param ket_indices the ket ("lower") external indices +/// @return an ExprPtr to the Ŝ Tensor +template +ExprPtr make_symmetrizer(const bra &bra_indices, + const ket &ket_indices) { + return ex(Tensor(reserved::symm_label(), bra_indices, ket_indices, + symmetrizer_symmetries)); +} + +/// @brief makes a symmetrization operator Ŝ +/// @param bra_indices the bra ("upper") external indices +/// @param ket_indices the ket ("lower") external indices +/// @param aux_indices the aux indices +/// @return an ExprPtr to the Ŝ Tensor +template +ExprPtr make_symmetrizer(const bra &bra_indices, + const ket &ket_indices, + const aux &aux_indices) { + return ex(Tensor(reserved::symm_label(), bra_indices, ket_indices, + aux_indices, symmetrizer_symmetries)); +} + +/// @} + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_TENSOR_HPP diff --git a/SeQuant/core/io/serialization/v1/ast_conversions.hpp b/SeQuant/core/io/serialization/v1/ast_conversions.hpp index 11a6e1f7e8..e83c6a6d02 100644 --- a/SeQuant/core/io/serialization/v1/ast_conversions.hpp +++ b/SeQuant/core/io/serialization/v1/ast_conversions.hpp @@ -323,12 +323,32 @@ struct Transformer { ann(std::move(braIndices)), vac); } - // Set required symmetries for symmetrization operators + // Force the defining symmetries of the reserved (anti)symmetrization + // operators; see sequant::{anti,}symmetrizer_symmetries. + const bool is_reserved_symmetrizer = + tensor.name == reserved::antisymm_label() || + tensor.name == reserved::symm_label(); if (tensor.name == reserved::antisymm_label()) { + //  antisymmetrizes within bra and within ket, Ŝ only across the + // {bra,ket} particle columns (i.e. it is perm-Nonsymm) perm_symm = Symmetry::Antisymm; } else if (tensor.name == reserved::symm_label()) { + perm_symm = Symmetry::Nonsymm; + } + if (is_reserved_symmetrizer) { + // force it rather than passing the Context's column default through, + // which the Tensor ctor would reject as a contradicting *explicit* + // request column_symm = ColumnSymmetry::Symm; } + // likewise, force braket-Nonsymm rather than inheriting the Context's + // default Hermiticity, which could derive a non-Nonsymm braket and make a + // plain "Ŝ{...}"/"Â{...}" fail to construct; an explicit braket spec is + // left untouched so that the Tensor ctor still rejects it + if (is_reserved_symmetrizer && + (!tensor.symmetry.has_value() || + tensor.symmetry.value().braket_symm == ast::SymmetrySpec::unspecified)) + braket_symm = BraKetSymmetry::Nonsymm; // Dispatch to correct Tensor constructor (taking either BraKetSymmetry or // Hermiticity) diff --git a/SeQuant/core/io/serialization/v1/deserialize.cpp b/SeQuant/core/io/serialization/v1/deserialize.cpp index d9c07a7beb..c49c8588e6 100644 --- a/SeQuant/core/io/serialization/v1/deserialize.cpp +++ b/SeQuant/core/io/serialization/v1/deserialize.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -246,8 +247,11 @@ AST parse(const StartRule &start, std::wstring_view input, transform::DefaultSymmetries to_default_symms( const DeserializationOptions &options) { - transform::DefaultSymmetries symms{Symmetry::Nonsymm, BraKetSymmetry::Nonsymm, - ColumnSymmetry::Nonsymm}; + // unspecified symmetries default to the active Context's; the Hermiticity is + // resolved against each tensor's base_field downstream, as in the Tensor ctor + const Context &ctx = get_default_context(); + transform::DefaultSymmetries symms{ctx.symmetry(), ctx.hermiticity(), + ctx.column_symmetry()}; if (options.def_perm_symm.has_value()) { std::get<0>(symms) = options.def_perm_symm.value(); diff --git a/SeQuant/domain/mbpt/antisymmetrizer.cpp b/SeQuant/domain/mbpt/antisymmetrizer.cpp index 130df98a1d..22f06e151d 100644 --- a/SeQuant/domain/mbpt/antisymmetrizer.cpp +++ b/SeQuant/domain/mbpt/antisymmetrizer.cpp @@ -85,8 +85,10 @@ antisymm_element::antisymm_element(ExprPtr ex_) { new_kets.push_back(unique_kets_list[j].second[index_label_pos]); index_label_pos++; } - auto new_tensor = ex(label, bra(std::move(new_bras)), - ket(std::move(new_kets))); + // mbpt tensors are particle (column) symmetric + auto new_tensor = ex( + label, bra(std::move(new_bras)), ket(std::move(new_kets)), + Symmetry::Nonsymm, std::nullopt, ColumnSymmetry::Symm); new_product = new_tensor * new_product; new_product->canonicalize(); } @@ -344,9 +346,11 @@ ExprPtr max_similarity(const std::vector& original_upper, } } if (new_pairs > og_pairs) { - factor = ex(-1) * ex(factor->as().label(), - bra(std::move(current_lower)), - ket(std::move(current_upper))); + factor = ex(-1) * + ex(factor->as().label(), + bra(std::move(current_lower)), + ket(std::move(current_upper)), Symmetry::Nonsymm, + std::nullopt, ColumnSymmetry::Symm); } } else if (factor->is()) { std::vector current_upper; @@ -422,7 +426,8 @@ ExprPtr spin_sum(std::vector original_upper, new_lower.push_back(factor->as().bra()[i]); } factor = ex(L"Γ", factor->as().bra(), - factor->as().ket()); + factor->as().ket(), Symmetry::Nonsymm, + std::nullopt, ColumnSymmetry::Symm); } else if (factor->is()) { // prefactor = ex(-0.5) * // ex(factor->as().rank()) * prefactor; diff --git a/SeQuant/domain/mbpt/biorthogonalization.cpp b/SeQuant/domain/mbpt/biorthogonalization.cpp index e68850e74c..0c4b9174c8 100644 --- a/SeQuant/domain/mbpt/biorthogonalization.cpp +++ b/SeQuant/domain/mbpt/biorthogonalization.cpp @@ -651,8 +651,7 @@ ExprPtr biorthogonal_transform_pre_nnsproject_impl( auto bixs = ext_idxs | transform([](auto&& vec) { return get_bra_idx(vec); }); auto kixs = ext_idxs | transform([](auto&& vec) { return get_ket_idx(vec); }); - ExprPtr S_tensor = - ex(Tensor{reserved::symm_label(), bra(kixs), ket(bixs)}); + ExprPtr S_tensor = make_symmetrizer(bra(kixs), ket(bixs)); if (factor_out_nns_projector) { if (ext_idxs.size() > 1) { diff --git a/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp index c69c20ece5..8eb123d002 100644 --- a/SeQuant/domain/mbpt/op.cpp +++ b/SeQuant/domain/mbpt/op.cpp @@ -609,9 +609,11 @@ ExprPtr OpMaker::operator()(std::optional dep, [this, opsymm_opt, full_label, op_herm]( const auto& creidxs, const auto& annidxs, const auto& batchidxs, Symmetry opsymm) { + // mbpt operators act on indistinguishable particles, hence are + // particle (column) symmetric return ex(full_label, bra(creidxs), ket(annidxs), aux(batchidxs), opsymm_opt ? *opsymm_opt : opsymm, - op_herm); + op_herm, ColumnSymmetry::Symm); }, dep ? *dep : UseDepIdx::None, normalization); } @@ -620,8 +622,11 @@ ExprPtr OpMaker::operator()(std::optional dep, cre_spaces_, ann_spaces_, [this, opsymm_opt, full_label, op_herm]( const auto& creidxs, const auto& annidxs, Symmetry opsymm) { + // mbpt operators act on indistinguishable particles, hence are + // particle (column) symmetric return ex(full_label, bra(creidxs), ket(annidxs), - opsymm_opt ? *opsymm_opt : opsymm, op_herm); + opsymm_opt ? *opsymm_opt : opsymm, op_herm, + ColumnSymmetry::Symm); }, dep ? *dep : UseDepIdx::None, normalization); } @@ -691,8 +696,7 @@ ExprPtr F(bool use_tensor, const IndexSpace& reference_occupied) { ketidxs.push_back(m2); return ex(L"g", bra(std::move(braidxs)), ket(std::move(ketidxs)), Symmetry::Antisymm) * - ex(kronecker_label(), bra{m2}, ket{m1}, - Symmetry::Nonsymm); + make_kronecker(m2, m1); } else { // opsymm == Symmetry::Nonsymm auto braidx_J = braidxs; braidx_J.push_back(m1); @@ -703,13 +707,16 @@ ExprPtr F(bool use_tensor, const IndexSpace& reference_occupied) { auto ketidxs_K = ketidxs; using std::begin; ketidxs_K.emplace(begin(ketidxs_K), m2); + // g is an integral over indistinguishable particles, hence is + // particle (column) symmetric (which the Antisymm branch above + // gets implicitly, but Nonsymm perm does not) return (ex(L"g", bra(std::move(braidx_J)), - ket(std::move(ketidxs_J)), Symmetry::Nonsymm) - + ket(std::move(ketidxs_J)), Symmetry::Nonsymm, + std::nullopt, ColumnSymmetry::Symm) - ex(L"g", bra(std::move(braidx_K)), - ket(std::move(ketidxs_K)), - Symmetry::Nonsymm)) * - ex(kronecker_label(), bra{m2}, ket{m1}, - Symmetry::Nonsymm); + ket(std::move(ketidxs_K)), Symmetry::Nonsymm, + std::nullopt, ColumnSymmetry::Symm)) * + make_kronecker(m2, m1); } }); }; @@ -1344,9 +1351,12 @@ ExprPtr expectation_value_impl(ExprPtr expr, OpConnections connect, braidxs.size() == ketidxs.size()); // need to handle particle # violating case? const auto rank = braidxs.size(); + // an RDM over indistinguishable particles is particle (column) + // symmetric (which the Antisymm branch implies, but Nonsymm does not) return ex( rdm_label, bra(std::move(braidxs)), ket(std::move(ketidxs)), - rank > 1 && spinor ? Symmetry::Antisymm : Symmetry::Nonsymm); + rank > 1 && spinor ? Symmetry::Antisymm : Symmetry::Nonsymm, + std::nullopt, ColumnSymmetry::Symm); }; if (exptr.template is()) { diff --git a/SeQuant/domain/mbpt/rdm.cpp b/SeQuant/domain/mbpt/rdm.cpp index 1e49256f31..f3c2403e1c 100644 --- a/SeQuant/domain/mbpt/rdm.cpp +++ b/SeQuant/domain/mbpt/rdm.cpp @@ -16,6 +16,12 @@ std::wstring rdm_cumulant_label() { namespace sequant::mbpt::decompositions { +namespace { +// RDMs, cumulants and the intermediates built from them are over +// indistinguishable particles, hence are particle (column) symmetric +constexpr TensorSymmetries particle_symmetric{.column = ColumnSymmetry::Symm}; +} // namespace + ExprPtr cumu_to_density(ExprPtr ex_) { SEQUANT_ASSERT(ex_->is()); SEQUANT_ASSERT(ex_->as().rank() == 1); @@ -23,7 +29,8 @@ ExprPtr cumu_to_density(ExprPtr ex_) { auto down_0 = ex_->as().ket()[0]; auto up_0 = ex_->as().bra()[0]; - auto density = ex(rdm_label(), bra{up_0}, ket{down_0}); + auto density = + ex(rdm_label(), bra{up_0}, ket{down_0}, particle_symmetric); return density; } @@ -37,9 +44,12 @@ sequant::ExprPtr cumu2_to_density(sequant::ExprPtr ex_) { auto down_1 = ex_->as().ket()[1]; auto up_1 = ex_->as().bra()[1]; - auto density2 = ex(rdm_label(), bra{up_0, up_1}, ket{down_0, down_1}); - auto density_1 = ex(rdm_label(), bra{up_0}, ket{down_0}); - auto density_2 = ex(rdm_label(), bra{up_1}, ket{down_1}); + auto density2 = ex(rdm_label(), bra{up_0, up_1}, ket{down_0, down_1}, + particle_symmetric); + auto density_1 = + ex(rdm_label(), bra{up_0}, ket{down_0}, particle_symmetric); + auto density_2 = + ex(rdm_label(), bra{up_1}, ket{down_1}, particle_symmetric); auto d1_d2 = antisymmetrize(density_1 * density_2); return density2 + ex(-1) * d1_d2.result; @@ -57,13 +67,16 @@ ExprPtr cumu3_to_density(ExprPtr ex_) { auto down_2 = ex_->as().ket()[2]; auto up_2 = ex_->as().bra()[2]; - auto cumulant2 = - ex(rdm_cumulant_label(), bra{up_1, up_2}, ket{down_1, down_2}); - auto density_1 = ex(rdm_label(), bra{up_0}, ket{down_0}); - auto density_2 = ex(rdm_label(), bra{up_1}, ket{down_1}); - auto density_3 = ex(rdm_label(), bra{up_2}, ket{down_2}); + auto cumulant2 = ex(rdm_cumulant_label(), bra{up_1, up_2}, + ket{down_1, down_2}, particle_symmetric); + auto density_1 = + ex(rdm_label(), bra{up_0}, ket{down_0}, particle_symmetric); + auto density_2 = + ex(rdm_label(), bra{up_1}, ket{down_1}, particle_symmetric); + auto density_3 = + ex(rdm_label(), bra{up_2}, ket{down_2}, particle_symmetric); auto density3 = ex(rdm_label(), bra{up_0, up_1, up_2}, - ket{down_0, down_1, down_2}); + ket{down_0, down_1, down_2}, particle_symmetric); auto d1_d2 = antisymmetrize(density_1 * density_2 * density_3 + density_1 * cumulant2); @@ -98,7 +111,8 @@ ExprPtr one_body_sub(ExprPtr ex_) { // J. Chem. Phys. 132, 234107 (2010); auto up_0 = ex_->as().creators()[0].index(); const auto a = ex(cre({up_0}), ann({down_0})); - const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); + const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}, + particle_symmetric); auto result = a + (ex(-1) * cumu1); return (result); @@ -117,12 +131,14 @@ ExprPtr two_body_decomp( auto up_0 = ex_->as().creators()[0].index(); auto up_1 = ex_->as().creators()[1].index(); - const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); - const auto cumu2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}); + const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}, + particle_symmetric); + const auto cumu2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}, + particle_symmetric); const auto a = ex(cre{up_1}, ann{down_1}); const auto a2 = ex(cre{up_0, up_1}, ann{down_0, down_1}); - const auto double_cumu = - ex(rdm_cumulant_label(), bra{down_0, down_1}, ket{up_0, up_1}); + const auto double_cumu = ex(rdm_cumulant_label(), bra{down_0, down_1}, + ket{up_0, up_1}, particle_symmetric); auto term1 = cumu1 * a; auto term2 = cumu1 * cumu2; @@ -151,20 +167,22 @@ three_body_decomp(ExprPtr ex_, bool approx) { std::vector initial_upper{up_0, up_1, up_2}; - const auto cumulant = - ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); + const auto cumulant = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}, + particle_symmetric); const auto a = ex(cre{up_1, up_2}, ann{down_1, down_2}); auto a_cumulant = cumulant * a; - auto cumulant2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}); - auto cumulant3 = ex(rdm_cumulant_label(), bra{down_2}, ket{up_2}); + auto cumulant2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}, + particle_symmetric); + auto cumulant3 = ex(rdm_cumulant_label(), bra{down_2}, ket{up_2}, + particle_symmetric); auto cumulant_3x = cumulant * cumulant2 * cumulant3; auto a1 = ex(cre{up_0}, ann{down_0}); auto a1_cumu1_cumu2 = a1 * cumulant2 * cumulant3; - auto two_body_cumu = - ex(rdm_cumulant_label(), bra{down_1, down_2}, ket{up_1, up_2}); + auto two_body_cumu = ex(rdm_cumulant_label(), bra{down_1, down_2}, + ket{up_1, up_2}, particle_symmetric); auto a1_cumu2 = a1 * two_body_cumu; auto cumu1_cumu2 = cumulant * two_body_cumu; @@ -173,7 +191,7 @@ three_body_decomp(ExprPtr ex_, bool approx) { if (!approx) { auto cumu3 = ex(rdm_cumulant_label(), bra{down_0, down_1, down_2}, - ket{up_0, up_1, up_2}); + ket{up_0, up_1, up_2}, particle_symmetric); sum_of_terms.result = cumu3 + sum_of_terms.result; } @@ -268,11 +286,11 @@ three_body_decomposition(ExprPtr ex_, int rank, bool fast) { // make tensors which can be decomposed into the constituent pieces later // in the procedure. auto DE2 = ex(L"DE2", bra{down_0, down_1, down_2}, - ket{up_0, up_1, up_2}); + ket{up_0, up_1, up_2}, particle_symmetric); auto DDE = ex(L"DDE", bra{down_0, down_1, down_2}, - ket{up_0, up_1, up_2}); + ket{up_0, up_1, up_2}, particle_symmetric); auto D2E = ex(L"D2E", bra{down_0, down_1, down_2}, - ket{up_0, up_1, up_2}); + ket{up_0, up_1, up_2}, particle_symmetric); auto result = DE2 + D2E - ex(2) * DDE; return {result, initial_pairing}; } diff --git a/SeQuant/domain/mbpt/rules/df.cpp b/SeQuant/domain/mbpt/rules/df.cpp index e32f99dab4..b7186aaa56 100644 --- a/SeQuant/domain/mbpt/rules/df.cpp +++ b/SeQuant/domain/mbpt/rules/df.cpp @@ -31,20 +31,24 @@ ExprPtr density_fit_impl(Tensor const& tnsr, Index const& aux_idx, // when the Tensor is built. auto t1 = ex(factor_label, bra({ranges::front(tnsr.bra())}), ket({ranges::front(tnsr.ket())}), aux({aux_idx}), - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); auto t2 = ex(factor_label, bra({ranges::back(tnsr.bra())}), ket({ranges::back(tnsr.ket())}), aux({aux_idx}), - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); if (tnsr.symmetry() == Symmetry::Antisymm) { auto t3 = ex(factor_label, bra({ranges::back(tnsr.bra())}), ket({ranges::front(tnsr.ket())}), aux({aux_idx}), - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); auto t4 = ex(factor_label, bra({ranges::front(tnsr.bra())}), ket({ranges::back(tnsr.ket())}), aux({aux_idx}), - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); return t1 * t2 - t3 * t4; } diff --git a/SeQuant/domain/mbpt/rules/thc.cpp b/SeQuant/domain/mbpt/rules/thc.cpp index 7fcfbf5978..dc41140464 100644 --- a/SeQuant/domain/mbpt/rules/thc.cpp +++ b/SeQuant/domain/mbpt/rules/thc.cpp @@ -15,6 +15,12 @@ namespace sequant::mbpt { +// THC/DF factor and central tensors act on indistinguishable particles, hence +// are particle (column) symmetric +namespace { +constexpr TensorSymmetries particle_symmetric{.column = ColumnSymmetry::Symm}; +} // namespace + ExprPtr tensor_hypercontract_impl(Tensor const& tnsr, Index const& aux_idx_1, Index const& aux_idx_2, std::wstring_view factor_label, @@ -24,20 +30,21 @@ ExprPtr tensor_hypercontract_impl(Tensor const& tnsr, Index const& aux_idx_1, && tnsr.aux_rank() == 0); auto t1 = ex(factor_label, bra({ranges::front(tnsr.bra())}), ket(), - aux({aux_idx_1})); + aux({aux_idx_1}), particle_symmetric); auto t2 = ex(factor_label, bra(), ket({ranges::front(tnsr.ket())}), - aux({aux_idx_1})); + aux({aux_idx_1}), particle_symmetric); auto t3 = ex(factor_label, bra({ranges::back(tnsr.bra())}), ket(), - aux({aux_idx_2})); + aux({aux_idx_2}), particle_symmetric); auto t4 = ex(factor_label, bra(), ket({ranges::back(tnsr.ket())}), - aux({aux_idx_2})); - auto z = ex(aux_label, bra(), ket(), aux({aux_idx_1, aux_idx_2})); + aux({aux_idx_2}), particle_symmetric); + auto z = ex(aux_label, bra(), ket(), aux({aux_idx_1, aux_idx_2}), + particle_symmetric); if (tnsr.symmetry() == Symmetry::Antisymm) { auto t1a = ex(factor_label, bra({ranges::back(tnsr.bra())}), ket(), - aux({aux_idx_1})); + aux({aux_idx_1}), particle_symmetric); auto t3a = ex(factor_label, bra({ranges::front(tnsr.bra())}), ket(), - aux({aux_idx_2})); + aux({aux_idx_2}), particle_symmetric); return (t1 * t2 * z * t3 * t4) - (t1a * t2 * z * t3a * t4); } diff --git a/SeQuant/domain/mbpt/spin.cpp b/SeQuant/domain/mbpt/spin.cpp index a182e40374..6b2dff4b97 100644 --- a/SeQuant/domain/mbpt/spin.cpp +++ b/SeQuant/domain/mbpt/spin.cpp @@ -324,8 +324,8 @@ ExprPtr remove_spin(const ExprPtr& expr) { } } return ex(tensor.label(), bra(std::move(b)), ket(std::move(k)), - tensor.aux(), tensor.symmetry(), - tensor.braket_symmetry()); + tensor.aux(), tensor.symmetry(), tensor.braket_symmetry(), + tensor.column_symmetry()); }; auto remove_spin_from_product = @@ -633,7 +633,7 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { auto S = Tensor{}; if (A_is_nconserving) { S = Tensor(reserved::symm_label(), A_tensor.bra(), A_tensor.ket(), - A_tensor.aux(), Symmetry::Nonsymm); + A_tensor.aux(), symmetrizer_symmetries); } else { // A is N-nonconserving auto n = std::min(A_tensor.bra_rank(), A_tensor.ket_rank()); container::svector bra_list(A_tensor.bra().begin(), @@ -641,7 +641,8 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { container::svector ket_list(A_tensor.ket().begin(), A_tensor.ket().begin() + n); S = Tensor(reserved::symm_label(), bra(std::move(bra_list)), - ket(std::move(ket_list)), A_tensor.aux(), Symmetry::Nonsymm); + ket(std::move(ket_list)), A_tensor.aux(), + symmetrizer_symmetries); } const auto nf = rational{1, factorial(S.ket_rank())}; @@ -1231,7 +1232,9 @@ ExprPtr merge_tensors(const Tensor& O1, const Tensor& O2) { auto b = ranges::views::concat(O1.bra(), O2.bra()); auto k = ranges::views::concat(O1.ket(), O2.ket()); auto a = ranges::views::concat(O1.aux(), O2.aux()); - return ex(Tensor(O1.label(), bra(b), ket(k), aux(a), O1.symmetry())); + // preserve all of O1's symmetry attributes, not just its perm symmetry + return ex(Tensor(O1.label(), bra(b), ket(k), aux(a), O1.symmetry(), + O1.braket_symmetry(), O1.column_symmetry())); } std::vector open_shell_A_op(const Tensor& A) { diff --git a/tests/integration/srcc.cpp b/tests/integration/srcc.cpp index e1e70e83df..a0f17f0607 100644 --- a/tests/integration/srcc.cpp +++ b/tests/integration/srcc.cpp @@ -103,7 +103,9 @@ class compute_cceqvec { : AssertStrictBraKetSymmetry::Yes; auto context_resetter = sequant::set_scoped_default_context( sequant::Context({.index_space_registry_shared_ptr = so_reg, - .vacuum = Vacuum::SingleProduct}) + .vacuum = Vacuum::SingleProduct, + // mbpt works with particle-symmetric tensors + .column_symmetry = ColumnSymmetry::Symm}) .set(so_strict)); std::vector eqvec_so; switch (type) { @@ -284,7 +286,9 @@ int main(int argc, char* argv[]) { sequant::set_default_context( sequant::Context({.index_space_registry_shared_ptr = sr_reg, .vacuum = Vacuum::SingleProduct, - .spbasis = spbasis}) + .spbasis = spbasis, + // mbpt works with particle-symmetric tensors + .column_symmetry = ColumnSymmetry::Symm}) .set(strict)); TensorCanonicalizer::register_instance( std::make_shared()); diff --git a/tests/integration/stcc.cpp b/tests/integration/stcc.cpp index aa6dade25b..a477a1bf64 100644 --- a/tests/integration/stcc.cpp +++ b/tests/integration/stcc.cpp @@ -38,7 +38,9 @@ int main(int argc, char* argv[]) { .vacuum = Vacuum::SingleProduct, .canonicalization_options = CanonicalizeOptions::default_options().copy_and_set( - CanonicalizationMethod::Complete)}); + CanonicalizationMethod::Complete), + // mbpt works with particle-symmetric tensors + .column_symmetry = ColumnSymmetry::Symm}); TensorCanonicalizer::register_instance( std::make_shared()); mbpt::set_default_mbpt_context( diff --git a/tests/unit/catch2_sequant.hpp b/tests/unit/catch2_sequant.hpp index 3dc0185ef1..91573f69a2 100644 --- a/tests/unit/catch2_sequant.hpp +++ b/tests/unit/catch2_sequant.hpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -22,6 +23,15 @@ #include #include +namespace sequant { +/// Shared particle-symmetric symmetry pack for the MBPT test TUs, which must +/// spell out MBPT particle symmetry explicitly since programmatic Tensor +/// construction is Context-independent (see Tensor::Defaults). Defined `inline` +/// here rather than once per TU so that unity test builds see one definition. +inline constexpr TensorSymmetries particle_symmetric{.column = + ColumnSymmetry::Symm}; +} // namespace sequant + namespace Catch { // Make sure Catch uses proper string representation for SeQuant types diff --git a/tests/unit/data/sf_r2_direct_real_inc.hpp b/tests/unit/data/sf_r2_direct_real_inc.hpp index 1241421ac9..c36866e12a 100644 --- a/tests/unit/data/sf_r2_direct_real_inc.hpp +++ b/tests/unit/data/sf_r2_direct_real_inc.hpp @@ -13,7 +13,7 @@ namespace sequant::tests::data { inline const std::wstring& sf_r2_direct_real() { static const std::wstring s = LR"FIXTURE( -8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_3,i_4}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3,a_4;i_2,i_1}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_4,i_2}:N-N-S * t{a_3,a_4;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_1}:N-S-S * t{a_2,a_3;i_3,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_2,i_4}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;i_1,a_3}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{i_2;i_3}:N-S-S * t{a_1,a_2;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_3,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_1}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_3,i_2}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_3}:N-N-S * t{a_2,a_4;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{i_1;i_3}:N-S-S * t{a_1,a_2;i_2,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_2}:N-N-S * t{a_1,a_3;i_3,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_2,i_3}:N-N-S * t{a_3,a_4;i_4,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,i_2;i_3,i_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_3,i_1}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{a_2;a_3}:N-S-S * t{a_1,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;a_1,a_2}:N-S-S * t{a_3;i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{a_1;a_3}:N-S-S * t{a_2,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2,a_3;i_1,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;a_2,i_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_3,i_2}:N-N-S * t{a_3,a_4;i_4,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;i_2,a_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2,a_3;i_3,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_2;i_4,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_2,i_4}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{a_3;i_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_2}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_1,i_3}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_1,i_4}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2,a_3;i_4,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,i_2;a_1,a_2}:N-S-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{a_3;i_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_1,i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,i_1;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;i_1,i_2}:N-S-S * t{a_1,a_2;i_4,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_4;i_2}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3,a_4;i_1,i_2}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_4,i_2}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2,a_3;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_4}:N-N-S * t{a_2,a_4;i_2,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_4,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_3}:N-N-S * t{a_1,a_2;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_4}:N-N-S * t{a_1,a_2;i_2,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;a_1,i_3}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3,a_4;i_2,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,a_4}:N-S-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;i_1,i_2}:N-S-S * t{a_1,a_2;i_3,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_1,a_4;i_4,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_4}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_4;i_3,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_4}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;i_1,a_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,i_2;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S * t{a_2,a_4;i_4,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_1,a_4;i_3,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_2;i_2,a_3}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_3}:N-N-S * t{a_1,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_3}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_4;i_1}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;a_1,i_3}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,i_2;a_2,a_1}:N-S-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_4;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_4;i_1}:N-N-S * t{a_1,a_2;i_2,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,i_1;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_3;i_3}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_2,i_4}:N-N-S * t{a_3,a_4;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_1,i_4}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_3,i_1}:N-N-S * t{a_2,a_3;i_4,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,a_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_4,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_1,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_2,i_4}:N-N-S * t{a_2,a_3;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_1,i_3}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2,a_3;i_3,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_3}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3,a_4;i_2,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{a_3;i_3}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_4,i_3}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * f{i_3;a_3}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_3,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_4}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_4,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_1,i_3}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_2,i_3}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_2}:N-N-S * t{a_1,a_3;i_1,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;a_2,a_1}:N-S-S * t{a_3;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_4;i_3}:N-N-S * t{a_2,a_3;i_1,i_2}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_4}:N-N-S * t{a_1,a_2;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_2,i_4}:N-N-S * t{a_2,a_4;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_2;a_3,i_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S +8 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_3,i_4}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3,a_4;i_2,i_1}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_4,i_2}:N-N-S * t{a_3,a_4;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_1}:N-S-S * t{a_2,a_3;i_3,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_2,i_4}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;i_1,a_3}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * f{i_2;i_3}:N-S-S * t{a_1,a_2;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_3,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_1}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_3,i_2}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_3}:N-N-S * t{a_2,a_4;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * f{i_1;i_3}:N-S-S * t{a_1,a_2;i_2,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_2}:N-N-S * t{a_1,a_3;i_3,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_2,i_3}:N-N-S * t{a_3,a_4;i_4,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,i_2;i_3,i_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_3,i_1}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * f{a_2;a_3}:N-S-S * t{a_1,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;a_1,a_2}:N-S-S * t{a_3;i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * f{a_1;a_3}:N-S-S * t{a_2,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2,a_3;i_1,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;a_2,i_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_3,i_2}:N-N-S * t{a_3,a_4;i_4,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;i_2,a_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2,a_3;i_3,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_2;i_4,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_2,i_4}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * f{a_3;i_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_2,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_2}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_1,i_3}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_1,i_4}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S * t{a_2,a_4;i_2,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_3,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2,a_3;i_4,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,i_2;a_1,a_2}:N-S-S - 8 Ŝ{i_1,i_2;a_1,a_2} * f{a_3;i_3}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_1,i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,i_1;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;i_1,i_2}:N-S-S * t{a_1,a_2;i_4,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_4;i_2}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3,a_4;i_1,i_2}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_4,i_2}:N-N-S * t{a_2,a_3;i_3,i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2,a_3;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_2;i_3}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_2}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_4}:N-N-S * t{a_2,a_4;i_2,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_4,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_1}:N-N-S * t{a_4;i_3}:N-N-S * t{a_1,a_2;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_3;i_4}:N-N-S * t{a_1,a_2;i_2,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;a_1,i_3}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3,a_4;i_2,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,a_4}:N-S-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;i_1,i_2}:N-S-S * t{a_1,a_2;i_3,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_1,a_4;i_4,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_4}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_4;i_3,i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_4}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;i_1,a_3}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,i_2;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S * t{a_2,a_4;i_4,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_3;i_1}:N-N-S * t{a_1,a_4;i_3,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_2;i_2,a_3}:N-S-S * t{a_1,a_3;i_3,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_3}:N-N-S * t{a_1,a_3;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_3}:N-N-S * t{a_1,a_2;i_1,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_2}:N-S-S * t{a_2,a_3;i_1,i_3}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_4;i_1}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;a_1,i_3}:N-S-S * t{a_2,a_3;i_2,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,i_2;a_2,a_1}:N-S-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_4;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_3;i_3}:N-N-S * t{a_4;i_1}:N-N-S * t{a_1,a_2;i_2,i_4}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,i_1;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_4}:N-N-S * t{a_3;i_3}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_2,i_4}:N-N-S * t{a_3,a_4;i_3,i_1}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S * t{a_1,a_2;i_3,i_4}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,i_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_4;i_2}:N-N-S * t{a_2,a_3;i_1,i_4}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_3,i_1}:N-N-S * t{a_2,a_3;i_4,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,a_4}:N-S-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_1,i_3}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_1,a_3;i_4,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_2}:N-N-S * t{a_2,a_4;i_1,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_2,i_4}:N-N-S * t{a_2,a_3;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_2;i_4}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_4;i_1,i_3}:N-N-S * t{a_2,a_3;i_2,i_4}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;i_2,a_3}:N-S-S * t{a_2,a_3;i_3,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_4,i_3}:N-N-S * t{a_2,a_4;i_1,i_2}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3,a_4;i_2,i_1}:N-N-S - 8 Ŝ{i_1,i_2;a_1,a_2} * f{a_3;i_3}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_2;i_4,i_3}:N-N-S * t{a_3,a_4;i_1,i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,a_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * f{i_3;a_3}:N-S-S * t{a_3;i_2}:N-N-S * t{a_1,a_2;i_3,i_1}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_2;i_3}:N-N-S * t{a_3;i_4}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S * t{a_2,a_4;i_4,i_2}:N-N-S + 8 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_1,i_3}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_3;i_3}:N-N-S * t{a_1,a_4;i_1,i_2}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;i_3,i_4}:N-S-S * t{a_2;i_4}:N-N-S * t{a_1,a_3;i_2,i_3}:N-N-S + 16 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_2}:N-S-S * t{a_4;i_2}:N-N-S * t{a_1,a_3;i_1,i_3}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;a_2,a_1}:N-S-S * t{a_3;i_2}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;i_3,a_1}:N-S-S * t{a_4;i_3}:N-N-S * t{a_2,a_3;i_1,i_2}:N-N-S - 16 Ŝ{i_1,i_2;a_1,a_2} * g{i_2,a_3;i_3,i_4}:N-S-S * t{a_3;i_4}:N-N-S * t{a_1,a_2;i_1,i_3}:N-N-S + 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,i_4;a_3,a_4}:N-S-S * t{a_1,a_3;i_2,i_4}:N-N-S * t{a_2,a_4;i_1,i_3}:N-N-S - 2 Ŝ{i_1,i_2;a_1,a_2} * g{a_3,a_4;a_1,a_2}:N-S-S * t{a_3;i_2}:N-N-S * t{a_4;i_1}:N-N-S - 4 Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_2;a_3,i_2}:N-S-S * t{a_1;i_3}:N-N-S * t{a_3;i_1}:N-N-S )FIXTURE"; return s; } diff --git a/tests/unit/test_biorthogonalization.cpp b/tests/unit/test_biorthogonalization.cpp index 8d5320d5a9..0cd04f12c8 100644 --- a/tests/unit/test_biorthogonalization.cpp +++ b/tests/unit/test_biorthogonalization.cpp @@ -27,20 +27,19 @@ TEST_CASE("biorthogonalization", "[Biorthogonalization]") { L"1/2 t{a1;i1}", L"1/6 (2 g{a1,a2;i1,i2} + g{a2,a1;i1,i2})", // cmp. Wang & Knizia (2018), DOI: arXiv:1805.00565 L"1/120 (" - "-7 t{i_1,i_2,i_3;a_2,a_3,a_1}:N-C-S " - "- 7 t{i_1,i_2,i_3;a_3,a_1,a_2}:N-C-S " - "- t{i_1,i_2,i_3;a_2,a_1,a_3}:N-C-S " - "- t{i_1,i_2,i_3;a_3,a_2,a_1}:N-C-S " - "- t{i_1,i_2,i_3;a_1,a_3,a_2}:N-C-S " - "+ 17t{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S)"}; + "-7 t{i_1,i_2,i_3;a_2,a_3,a_1}:N-N-S " + "- 7 t{i_1,i_2,i_3;a_3,a_1,a_2}:N-N-S " + "- t{i_1,i_2,i_3;a_2,a_1,a_3}:N-N-S " + "- t{i_1,i_2,i_3;a_3,a_2,a_1}:N-N-S " + "- t{i_1,i_2,i_3;a_1,a_3,a_2}:N-N-S " + "+ 17t{i_1,i_2,i_3;a_1,a_2,a_3}:N-N-S)"}; REQUIRE(inputs.size() == expected_outputs.size()); for (std::size_t i = 0; i < inputs.size(); ++i) { CAPTURE(i); - ExprPtr input_expr = - deserialize(inputs.at(i), {.def_col_symm = ColumnSymmetry::Symm}); + ExprPtr input_expr = deserialize(inputs.at(i)); auto externals = external_indices(input_expr); @@ -81,8 +80,7 @@ TEST_CASE("biorthogonalization", "[Biorthogonalization]") { container::svector expressions; container::svector expected; for (std::size_t k = 0; k < inputs.at(i).size(); ++k) { - ResultExpr parsed = deserialize( - inputs.at(i).at(k), {.def_col_symm = ColumnSymmetry::Symm}); + ResultExpr parsed = deserialize(inputs.at(i).at(k)); expressions.push_back(parsed); expected.push_back( @@ -110,8 +108,7 @@ TEST_CASE("biorthogonalization", "[Biorthogonalization]") { container::svector expressions; for (const std::wstring &str : current_inputs) { - expressions.push_back(deserialize( - str, {.def_col_symm = ColumnSymmetry::Symm})); + expressions.push_back(deserialize(str)); } REQUIRE_THROWS_WITH( diff --git a/tests/unit/test_canonicalize.cpp b/tests/unit/test_canonicalize.cpp index 5b32c3e2ba..be93d48d62 100644 --- a/tests/unit/test_canonicalize.cpp +++ b/tests/unit/test_canonicalize.cpp @@ -28,6 +28,9 @@ #include #include +// the `particle_symmetric` symmetry pack (column = Symm) is defined in +// catch2_sequant.hpp and shared across the MBPT test TUs + TEST_CASE("canonicalization", "[algorithms]") { using namespace sequant; @@ -42,25 +45,25 @@ TEST_CASE("canonicalization", "[algorithms]") { SECTION("Tensors") { { auto op = ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("g{p1,p2;p3,p4}")); } { auto op = ex(L"g", bra{L"p_2", L"p_1"}, ket{L"p_3", L"p_4"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("g{p1,p2;p4,p3}")); } { auto op = ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_4", L"p_3"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("g{p1,p2;p4,p3}")); } { auto op = ex(L"g", bra{L"p_2", L"p_1"}, ket{L"p_4", L"p_3"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("g{p1,p2;p3,p4}")); } @@ -80,13 +83,13 @@ TEST_CASE("canonicalization", "[algorithms]") { // aux indices { auto op = ex(L"B", bra{L"p_1"}, ket{L"p_2"}, aux{L"p_3"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("B{p1;p2;p3}")); } { auto op = ex(L"B", bra{L"p_1", L"p_2"}, ket{L"p_4", L"p_3"}, - aux{L"p_5"}, Symmetry::Nonsymm); + aux{L"p_5"}, particle_symmetric); canonicalize(op); REQUIRE_THAT(op, SimplifiesTo("B{p1,p2;p4,p3;p5}")); } @@ -113,11 +116,11 @@ TEST_CASE("canonicalization", "[algorithms]") { { auto input = ex(reserved::symm_label(), bra{L"a_1", L"a_2"}, - ket{L"i_1", L"i_2"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"a_5"}, ket{L"i_5"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"i_5"}, ket{L"a_1"}, Symmetry::Nonsymm) * + ket{L"i_1", L"i_2"}, particle_symmetric) * + ex(L"f", bra{L"a_5"}, ket{L"i_5"}, particle_symmetric) * + ex(L"t", bra{L"i_5"}, ket{L"a_1"}, particle_symmetric) * ex(L"t", bra{L"i_1", L"i_2"}, ket{L"a_5", L"a_2"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(input); REQUIRE_THAT( input, @@ -126,11 +129,11 @@ TEST_CASE("canonicalization", "[algorithms]") { { auto input = ex(reserved::symm_label(), bra{L"a_1", L"a_2"}, - ket{L"i_1", L"i_2"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"a_5"}, ket{L"i_5"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"i_1"}, ket{L"a_5"}, Symmetry::Nonsymm) * + ket{L"i_1", L"i_2"}, particle_symmetric) * + ex(L"f", bra{L"a_5"}, ket{L"i_5"}, particle_symmetric) * + ex(L"t", bra{L"i_1"}, ket{L"a_5"}, particle_symmetric) * ex(L"t", bra{L"i_5", L"i_2"}, ket{L"a_1", L"a_2"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(input); REQUIRE_THAT( input, @@ -176,13 +179,13 @@ TEST_CASE("canonicalization", "[algorithms]") { q2->adjoint(); auto input = ex(reserved::symm_label(), bra{L"a_1", L"a_2"}, - ket{L"i_1", L"i_2"}, Symmetry::Nonsymm) * - q2 * ex(L"f", bra{L"a_5"}, ket{L"i_5"}, Symmetry::Nonsymm) * + ket{L"i_1", L"i_2"}, particle_symmetric) * + q2 * ex(L"f", bra{L"a_5"}, ket{L"i_5"}, particle_symmetric) * ex(L"p") * - ex(L"t", bra{L"i_1"}, ket{L"a_5"}, Symmetry::Nonsymm) * + ex(L"t", bra{L"i_1"}, ket{L"a_5"}, particle_symmetric) * ex(L"q1") * ex(L"t", bra{L"i_5", L"i_2"}, ket{L"a_1", L"a_2"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(input); REQUIRE_THAT(input, SimplifiesTo("p q1 q2^* Ŝ{a_1,a_2;i_1,i_2} f{a_3;i_3} " @@ -190,22 +193,23 @@ TEST_CASE("canonicalization", "[algorithms]") { } { // Product containing adjoint of a Tensor auto f2 = ex(L"f", bra{L"a_1", L"a_2"}, ket{L"i_5", L"i_2"}, - Symmetry::Nonsymm, BraKetSymmetry::Nonsymm); + Symmetry::Nonsymm, BraKetSymmetry::Nonsymm, + ColumnSymmetry::Symm); f2->adjoint(); auto input1 = ex(reserved::symm_label(), bra{L"a_1", L"a_2"}, - ket{L"i_1", L"i_2"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"a_5"}, ket{L"i_5"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"i_1"}, ket{L"a_5"}, Symmetry::Nonsymm) * f2; + ket{L"i_1", L"i_2"}, particle_symmetric) * + ex(L"f", bra{L"a_5"}, ket{L"i_5"}, particle_symmetric) * + ex(L"t", bra{L"i_1"}, ket{L"a_5"}, particle_symmetric) * f2; canonicalize(input1); REQUIRE_THAT(input1, SimplifiesTo("Ŝ{a_1,a_2;i_1,i_2} f{a_3;i_3} " "f⁺{i_1,i_3;a_1,a_2}:N-N-S t{i_2;a_3}")); auto input2 = ex(reserved::symm_label(), bra{L"a_1", L"a_2"}, - ket{L"i_1", L"i_2"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"a_5"}, ket{L"i_5"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"i_1"}, ket{L"a_5"}, Symmetry::Nonsymm) * f2 * + ket{L"i_1", L"i_2"}, particle_symmetric) * + ex(L"f", bra{L"a_5"}, ket{L"i_5"}, particle_symmetric) * + ex(L"t", bra{L"i_1"}, ket{L"a_5"}, particle_symmetric) * f2 * ex(L"w") * ex(rational{1, 2}); canonicalize(input2); REQUIRE_THAT(input2, @@ -217,11 +221,11 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(rational{1, 2}) * ex(L"B", bra{L"p_2"}, ket{L"p_4"}, aux{L"p_5"}, - Symmetry::Nonsymm) * + particle_symmetric) * ex(L"B", bra{L"p_1"}, ket{L"p_3"}, aux{L"p_5"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm); + particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric) * + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric); canonicalize(input); // because bra and ket are in same space dummy renaming flips the bra and // ket even though the tensors are not bra-ket symmetric @@ -238,9 +242,11 @@ TEST_CASE("canonicalization", "[algorithms]") { // whether bra/ket swap should occur for each tensor auto input = ex(rational{1, 2}) * ex(L"B", bra{L"p_2"}, ket{L"p_1"}, aux{L"p_5"}, - Symmetry::Nonsymm, BraKetSymmetry::Symm) * + Symmetry::Nonsymm, BraKetSymmetry::Symm, + ColumnSymmetry::Symm) * ex(L"B", bra{L"p_1"}, ket{L"p_2"}, aux{L"p_5"}, - Symmetry::Nonsymm, BraKetSymmetry::Symm); + Symmetry::Nonsymm, BraKetSymmetry::Symm, + ColumnSymmetry::Symm); REQUIRE_THAT(input, EquivalentTo("1/2 B{p1;p2;p5}:N-S B{p1;p2;p5}:N-S")); } // SF R2 ±pair extracted from the real-field CCSD doubles. Under @@ -268,9 +274,9 @@ TEST_CASE("canonicalization", "[algorithms]") { .spbasis = SPBasis::Spinfree}) .set(AssertStrictBraKetSymmetry::No)); auto input = deserialize( - L"8 * Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_1}:N-S-S " + L"8 * Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_1}:N-S-S " L"* t{a_2,a_3;i_2,i_3}:N-N-S " - L"+ 8 * Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;a_1,i_3}:N-S-S " + L"+ 8 * Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;a_1,i_3}:N-S-S " L"* t{a_2,a_3;i_2,i_3}:N-N-S"); simplify(input); const std::size_t n = @@ -291,10 +297,10 @@ TEST_CASE("canonicalization", "[algorithms]") { ctx_min.set(AssertStrictBraKetSymmetry::No); auto resetter = set_scoped_default_context(ctx_min); auto exA = deserialize( - L"8 * Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3,a_1;a_3,i_1}:N-S-S " + L"8 * Ŝ{i_1,i_2;a_1,a_2} * g{i_3,a_1;a_3,i_1}:N-S-S " L"* t{a_2,a_3;i_2,i_3}:N-N-S"); auto exB = deserialize( - L"8 * Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_1,a_3;a_1,i_3}:N-S-S " + L"8 * Ŝ{i_1,i_2;a_1,a_2} * g{i_1,a_3;a_1,i_3}:N-S-S " L"* t{a_2,a_3;i_2,i_3}:N-N-S"); REQUIRE(exA); REQUIRE(exB); @@ -476,14 +482,14 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(rational{1, 2}) * ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm) + + particle_symmetric) * + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric) + ex(rational{1, 2}) * ex(L"g", bra{L"p_2", L"p_1"}, ket{L"p_4", L"p_3"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm); + particle_symmetric) * + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric); simplify(input); canonicalize(input); REQUIRE_THAT(input, @@ -496,13 +502,13 @@ TEST_CASE("canonicalization", "[algorithms]") { ex(rational{1, 2}) * ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"}, Symmetry::Symm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm) + + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric) + ex(rational{1, 2}) * ex(L"g", bra{L"p_2", L"p_1"}, ket{L"p_4", L"p_3"}, Symmetry::Symm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm); + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric); canonicalize(input); REQUIRE_THAT(input, EquivalentTo("g{p2,p3;p1,p4}:S t{p1;p2} t{p4;p3}")); } @@ -513,13 +519,13 @@ TEST_CASE("canonicalization", "[algorithms]") { ex(rational{1, 2}) * ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"}, Symmetry::Antisymm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm) + + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric) + ex(rational{1, 2}) * ex(L"g", bra{L"p_2", L"p_1"}, ket{L"p_4", L"p_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"p_3"}, ket{L"p_1"}, Symmetry::Nonsymm) * - ex(L"t", bra{L"p_4"}, ket{L"p_2"}, Symmetry::Nonsymm); + ex(L"t", bra{L"p_3"}, ket{L"p_1"}, particle_symmetric) * + ex(L"t", bra{L"p_4"}, ket{L"p_2"}, particle_symmetric); canonicalize(input); REQUIRE_THAT(input, EquivalentTo("g{p2,p3;p1,p4}:A t{p1;p2} t{p4;p3}")); } @@ -530,13 +536,13 @@ TEST_CASE("canonicalization", "[algorithms]") { ex(rational{4, 3}) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"a_3", L"i_1"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_3"}, Symmetry::Nonsymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_3"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_4", L"i_2"}, Symmetry::Antisymm) - ex(rational{1, 3}) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_1", L"a_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}, Symmetry::Nonsymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_3", L"i_2"}, Symmetry::Antisymm); canonicalize(input); @@ -550,16 +556,16 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(rational{4, 3}) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"a_3", L"i_1"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_3"}, Symmetry::Nonsymm) * + particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_3"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_4", L"i_2"}, - Symmetry::Nonsymm) - + particle_symmetric) - ex(rational{1, 3}) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_1", L"a_3"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}, Symmetry::Nonsymm) * + particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_3", L"i_2"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(input); REQUIRE(input->size() == 1); @@ -573,16 +579,16 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(-4) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_3", L"i_2", L"i_4"}, Symmetry::Nonsymm) + + ket{L"i_3", L"i_2", L"i_4"}, particle_symmetric) + ex(-4) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_2", L"i_4", L"i_3"}, Symmetry::Nonsymm); + ket{L"i_2", L"i_4", L"i_3"}, particle_symmetric); canonicalize(input); REQUIRE_THAT( input, @@ -594,17 +600,17 @@ TEST_CASE("canonicalization", "[algorithms]") { auto term1 = ex(-4) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_3", L"i_2", L"i_4"}, Symmetry::Nonsymm); + ket{L"i_3", L"i_2", L"i_4"}, particle_symmetric); auto term2 = ex(-4) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_2", L"i_4", L"i_3"}, Symmetry::Nonsymm); + ket{L"i_2", L"i_4", L"i_3"}, particle_symmetric); canonicalize(term1); canonicalize(term2); REQUIRE_THAT(term1, @@ -625,16 +631,16 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(2) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_3", L"i_4", L"i_2"}, Symmetry::Nonsymm) + + ket{L"i_3", L"i_4", L"i_2"}, particle_symmetric) + ex(2) * ex(reserved::symm_label(), bra{L"i_1", L"i_2", L"i_3"}, - ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}) * + ket{L"a_1", L"a_2", L"a_3"}, particle_symmetric) * + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, - ket{L"i_2", L"i_3", L"i_4"}, Symmetry::Nonsymm); + ket{L"i_2", L"i_3", L"i_4"}, particle_symmetric); canonicalize(input); REQUIRE_THAT( input, EquivalentTo( @@ -647,20 +653,20 @@ TEST_CASE("canonicalization", "[algorithms]") { auto input = ex(rational{4, 3}) * ex(L"B", bra{L"i_3"}, ket{L"a_3"}, aux{L"p_5"}, - Symmetry::Nonsymm) * + particle_symmetric) * ex(L"B", bra{L"i_4"}, ket{L"i_1"}, aux{L"p_5"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_3"}, Symmetry::Nonsymm) * + particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_3"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_4", L"i_2"}, - Symmetry::Nonsymm) - + particle_symmetric) - ex(rational{1, 3}) * ex(L"B", bra{L"i_3"}, ket{L"i_1"}, aux{L"p_5"}, - Symmetry::Nonsymm) * + particle_symmetric) * ex(L"B", bra{L"i_4"}, ket{L"a_3"}, aux{L"p_5"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}, Symmetry::Nonsymm) * + particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric) * ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_3", L"i_2"}, - Symmetry::Nonsymm); + particle_symmetric); canonicalize(input); simplify(input); diff --git a/tests/unit/test_eval_expr.cpp b/tests/unit/test_eval_expr.cpp index 4d8339404f..d0002e8677 100644 --- a/tests/unit/test_eval_expr.cpp +++ b/tests/unit/test_eval_expr.cpp @@ -192,9 +192,12 @@ TEST_CASE("eval_expr", "[EvalExpr]") { res = deserialize(L"Amplitude{i1;a1} = t{a1;i1}"); root_expr = binarize(res)->expr(); REQUIRE(root_expr.is()); - REQUIRE(root_expr.as() == Tensor(L"Amplitude", - bra(IndexList{L"i_1"}), - ket(IndexList{L"a_1"}))); + // the deserialized ResultExpr's Amplitude picks up the Context's column + // symmetry (Symm), so the programmatic reference must request it too -- + // programmatic ctors are Context-independent (see Tensor::Defaults) + REQUIRE(root_expr.as() == + Tensor(L"Amplitude", bra(IndexList{L"i_1"}), ket(IndexList{L"a_1"}), + TensorSymmetries{.column = ColumnSymmetry::Symm})); } SECTION("Adjoint op") { diff --git a/tests/unit/test_main.cpp b/tests/unit/test_main.cpp index 53d43a326e..be040855ce 100644 --- a/tests/unit/test_main.cpp +++ b/tests/unit/test_main.cpp @@ -48,7 +48,12 @@ int main(int argc, char* argv[]) { .braket_typesetting = BraKetTypesetting::ContraSub, // to_latex() reference outputs predominantly assume the original // (naive) convention - .braket_slot_typesetting = BraKetSlotTypesetting::Naive}); + .braket_slot_typesetting = BraKetSlotTypesetting::Naive, + // mbpt works with particle-symmetric tensors, so default new/parsed + // tensors to particle (column) symmetry; bra-ket and permutational + // symmetries stay at the safe Nonsymm/NonHermitian library defaults and + // are specified explicitly where needed (e.g. Hermitian integrals) + .column_symmetry = ColumnSymmetry::Symm}); TensorCanonicalizer::set_cardinal_tensor_labels( sequant::mbpt::cardinal_tensor_labels()); // uncomment to enable verbose output ... diff --git a/tests/unit/test_mbpt.cpp b/tests/unit/test_mbpt.cpp index 0e3e8f4e10..37c1e9ef84 100644 --- a/tests/unit/test_mbpt.cpp +++ b/tests/unit/test_mbpt.cpp @@ -1039,8 +1039,8 @@ SECTION("SRSF") { auto scalar = expr->as().scalar(); REQUIRE(scalar == 1); REQUIRE_THAT(expr, - EquivalentTo("Ŝ{a_3,a_4;i_3,i_4}:N-C-S * " - "f{a_1,a_2;i_1,i_2}:N-C-S * ã{i_3,i_4;a_3,a_4}")); + EquivalentTo("Ŝ{a_3,a_4;i_3,i_4} * " + "f{a_1,a_2;i_1,i_2}:N-N-S * ã{i_3,i_4;a_3,a_4}")); } } // SECTION("SRSF") @@ -1135,9 +1135,9 @@ SECTION("rules") { const std::vector expected = { L"t{a1,a2;i1,i2} t{a3;i3}", L"t{a1,a2;i1,i2} g{a3;i3}", - L"t{a1,a2;i1,i2} B{i1;a1;x_1} B{i2;a2;x_1}", - L"t{a1,a2;i1,i2} (B{i1;a1;x_1} B{i2;a2;x_1} " - "- B{i2;a1;x_1} B{i1;a2;x_1})", + L"t{a1,a2;i1,i2} B{i1;a1;x_1}:N-C-S B{i2;a2;x_1}:N-C-S", + L"t{a1,a2;i1,i2} (B{i1;a1;x_1}:N-C-S B{i2;a2;x_1}:N-C-S " + "- B{i2;a1;x_1}:N-C-S B{i1;a2;x_1}:N-C-S)", }; REQUIRE(inputs.size() == expected.size()); @@ -1363,7 +1363,7 @@ SECTION("avoided-connections") { // only one term with no A-{f,g} connection REQUIRE(res2.is()); const std::wstring expected2 = - L"-1 Â{i_1;a_1}:A-C-S t{a_1,a_2;i_2,i_1}:A-C-S f{i_2;a_2}:A-C-S"; + L"-1 Â{i_1;a_1} t{a_1,a_2;i_2,i_1}:A-C-S f{i_2;a_2}:A-C-S"; REQUIRE_THAT(sequant::simplify(res2), EquivalentTo(expected2)); // same test as above but from Operator level and labels for connectivity @@ -1383,7 +1383,7 @@ SECTION("avoided-connections") { REQUIRE(res4_full.size() == 4); REQUIRE(res4.is()); // only single term survives const std::wstring expected4 = - L"Â{i_1;a_2}:A-C-S Â{a_1;i_2}:A-C-S g{i_3,i_2;a_3,a_1}:A-C-S " + L"Â{i_1;a_2} Â{a_1;i_2} g{i_3,i_2;a_3,a_1}:A-C-S " L"t{a_3,a_2;i_3,i_1}:A-C-S"; REQUIRE_THAT(simplify(res4), EquivalentTo(expected4)); } diff --git a/tests/unit/test_parse.cpp b/tests/unit/test_parse.cpp index f4256ad14f..2cf135abec 100644 --- a/tests/unit/test_parse.cpp +++ b/tests/unit/test_parse.cpp @@ -89,6 +89,30 @@ TEST_CASE("serialization", "[serialization]") { ctx.set(mbpt::make_sr_spaces()); auto ctx_resetter = set_scoped_default_context(ctx); + SECTION("default symmetries follow the Context") { + // the deserializer sources unspecified tensor symmetries from the active + // Context (unlike the programmatic Tensor ctor, which uses fixed + // defaults -- see test_tensor "programmatic ctor defaults are + // Context-independent"). Flip the Context column-symmetry default and + // confirm parsed tensors follow it. + auto ctx_symm = get_default_context(); + ctx_symm.set(ColumnSymmetry::Symm); + { + auto resetter = set_scoped_default_context(ctx_symm); + REQUIRE(deserialize(L"t{i1,i2;a1,a2}") + ->as() + .column_symmetry() == ColumnSymmetry::Symm); + } + auto ctx_nonsymm = get_default_context(); + ctx_nonsymm.set(ColumnSymmetry::Nonsymm); + { + auto resetter = set_scoped_default_context(ctx_nonsymm); + REQUIRE(deserialize(L"t{i1,i2;a1,a2}") + ->as() + .column_symmetry() == ColumnSymmetry::Nonsymm); + } + } + SECTION("Scalar tensor") { auto expr = deserialize(L"t{}"); REQUIRE(expr->is()); @@ -500,6 +524,60 @@ TEST_CASE("serialization", "[serialization]") { io::serialization::SerializationError, serializationErrorMatches(9, 1, "Invalid symmetry specifier")); } + + SECTION("(anti)symmetrization operators reject braket symmetry") { + // a reserved (anti)symmetrizer must be braket-Nonsymm; an explicit + // non-Nonsymm braket spec is rejected by the Tensor ctor + REQUIRE_THROWS(deserialize(L"Ŝ{i1,i2;a1,a2}:N-C-S")); + REQUIRE_THROWS(deserialize(L"Â{i1,i2;a1,a2}:A-C-S")); + // the default (braket-Nonsymm) form is accepted + REQUIRE_NOTHROW(deserialize(L"Ŝ{i1,i2;a1,a2}")); + REQUIRE_NOTHROW(deserialize(L"Â{i1,i2;a1,a2}")); + // ... and stays accepted even when the Context defaults to Hermitian: + // reserved (anti)symmetrizers must not inherit the Context braket + // default (which would otherwise derive a non-Nonsymm braket and throw) + { + auto ctx = get_default_context(); + ctx.set(Hermiticity::Hermitian); + auto resetter = set_scoped_default_context(ctx); + REQUIRE_NOTHROW(deserialize(L"Ŝ{i1,i2;a1,a2}")); + REQUIRE_NOTHROW(deserialize(L"Â{i1,i2;a1,a2}")); + } + } + + SECTION("(anti)symmetrization operators have fixed perm symmetry") { + // the defining bra/ket permutational symmetry of a reserved + // (anti)symmetrizer is not a free parameter either, so it must not be + // inherited from the Context; a parsed operator must match the one + // built by make_(anti)symmetrizer() + auto ctx = get_default_context(); + ctx.set(Symmetry::Antisymm); + auto resetter = set_scoped_default_context(ctx); + const auto S = deserialize(L"Ŝ{i1,i2;a1,a2}"); + REQUIRE(S->as().symmetry() == Symmetry::Nonsymm); + REQUIRE(S == make_symmetrizer(bra{Index(L"i_1"), Index(L"i_2")}, + ket{Index(L"a_1"), Index(L"a_2")})); + const auto A = deserialize(L"Â{i1,i2;a1,a2}"); + REQUIRE(A->as().symmetry() == Symmetry::Antisymm); + REQUIRE(A == make_antisymmetrizer(bra{Index(L"i_1"), Index(L"i_2")}, + ket{Index(L"a_1"), Index(L"a_2")})); + } + + SECTION("(anti)symmetrization operators are always column symmetric") { + // ... including under a Context whose column default is the library + // default Nonsymm: the deserializer must force Symm rather than pass + // the (contradicting) Context default to the Tensor ctor + auto ctx = get_default_context(); + ctx.set(ColumnSymmetry::Nonsymm); + auto resetter = set_scoped_default_context(ctx); + for (const auto& input : + {L"Ŝ{i1,i2;a1,a2}", L"Â{i1,i2;a1,a2}", L"Ŝ{i1,i2;a1,a2}:N-N-N", + L"Â{i1,i2;a1,a2}:A-N-N"}) { + ExprPtr expr; + REQUIRE_NOTHROW(expr = deserialize(input)); + REQUIRE(expr->as().column_symmetry() == ColumnSymmetry::Symm); + } + } } } diff --git a/tests/unit/test_spin.cpp b/tests/unit/test_spin.cpp index 3660245d1e..7d35a1ee05 100644 --- a/tests/unit/test_spin.cpp +++ b/tests/unit/test_spin.cpp @@ -33,6 +33,9 @@ #include #include +// the particle-symmetric default symmetry pack `particle_symmetric` (column = +// Symm) is defined in catch2_sequant.hpp and shared across the MBPT test TUs + TEST_CASE("spin", "[spin]") { using namespace sequant; using namespace sequant::mbpt; @@ -61,8 +64,8 @@ TEST_CASE("spin", "[spin]") { Index i1(L"i_1"); Index a1(L"a_1", {i1}); - const auto expr = - ex(L"t", bra{i1}, ket{a1}) * ex(L"F", bra{a1}, ket{i1}); + const auto expr = ex(L"t", bra{i1}, ket{a1}, particle_symmetric) * + ex(L"F", bra{a1}, ket{i1}, particle_symmetric); REQUIRE_NOTHROW(spintrace(expr)); { // assume spin-free spaces auto expr_st = spintrace(expr); @@ -195,7 +198,7 @@ TEST_CASE("spin", "[spin]") { auto p3 = Index(L"p↑_3"); auto p4 = Index(L"p↓_4"); - auto input = ex(L"t", bra{p1, p2}, ket{p3, p4}); + auto input = ex(L"t", bra{p1, p2}, ket{p3, p4}, particle_symmetric); REQUIRE(can_expand(input->as()) == true); REQUIRE(ms_conserving_columns(input->as()) == true); @@ -206,7 +209,7 @@ TEST_CASE("spin", "[spin]") { for (auto& i : result->as().const_braket_indices()) REQUIRE(i.space().base_key() == L"p"); - input = ex(L"t", bra{p1, p3}, ket{p2, p4}); + input = ex(L"t", bra{p1, p3}, ket{p2, p4}, particle_symmetric); REQUIRE_THAT(swap_spin(input), EquivalentTo("t{p↓1,p↓3;p↑2,p↑4}")); REQUIRE(can_expand(input->as()) == false); REQUIRE(ms_conserving_columns(input->as()) == false); @@ -214,7 +217,7 @@ TEST_CASE("spin", "[spin]") { SECTION("Tensor: expand_antisymm") { // 1-body - auto input = ex(L"t", bra{L"a_1"}, ket{L"i_1"}); + auto input = ex(L"t", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric); auto result = expand_antisymm(input->as()); REQUIRE(input->as() == result->as()); REQUIRE(!result->is()); @@ -293,17 +296,18 @@ TEST_CASE("spin", "[spin]") { REQUIRE_THAT( result, EquivalentTo( - "1/4 g{p↑_1,p↑_2;p↑_3,p↑_4}:N-C-S - 1/4 " - "g{p↑_2,p↑_1;p↑_3,p↑_4}:N-C-S + 1/4 g{p↑_1,p↓_2;p↑_3,p↓_4}:N-C-S " - "+ 1/4 g{p↑_2,p↓_1;p↑_4,p↓_3}:N-C-S + 1/4 " - "g{p↓_1,p↓_2;p↓_3,p↓_4}:N-C-S - 1/4 " - "g{p↓_2,p↓_1;p↓_3,p↓_4}:N-C-S")); + "1/4 g{p↑_1,p↑_2;p↑_3,p↑_4}:N-N-S - 1/4 " + "g{p↑_2,p↑_1;p↑_3,p↑_4}:N-N-S + 1/4 g{p↑_1,p↓_2;p↑_3,p↓_4}:N-N-S " + "+ 1/4 g{p↑_2,p↓_1;p↑_4,p↓_3}:N-N-S + 1/4 " + "g{p↓_1,p↓_2;p↓_3,p↓_4}:N-N-S - 1/4 " + "g{p↓_2,p↓_1;p↓_3,p↓_4}:N-N-S")); } } SECTION("Product") { - const auto expr = ex(L"f", bra{L"i_1"}, ket{L"a_1"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_1"}); + const auto expr = + ex(L"f", bra{L"i_1"}, ket{L"a_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric); auto result = spintrace(expr); canonicalize(result); REQUIRE_THAT(result, EquivalentTo("2 f{i1;a1} t{a1;i1}")); @@ -312,11 +316,12 @@ TEST_CASE("spin", "[spin]") { SECTION("Scaled Product") { { // 1/2 * g * t1 * t1 - const auto expr = ex(rational{1, 2}) * - ex(L"g", bra{L"i_1", L"i_2"}, - ket{L"a_1", L"a_2"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_1"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}); + const auto expr = + ex(rational{1, 2}) * + ex(L"g", bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric); auto result = spintrace(expr); canonicalize(result); REQUIRE_THAT(result, @@ -349,13 +354,15 @@ TEST_CASE("spin", "[spin]") { SECTION("Sum") { // f * t1 + 1/2 * g * t1 * t1 + 1/4 * g * t2 - const auto ex1 = ex(L"f", bra{L"i_1"}, ket{L"a_1"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_1"}); - const auto ex2 = ex(rational{1, 2}) * - ex(L"g", bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, - Symmetry::Antisymm) * - ex(L"t", bra{L"a_1"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}); + const auto ex1 = + ex(L"f", bra{L"i_1"}, ket{L"a_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric); + const auto ex2 = + ex(rational{1, 2}) * + ex(L"g", bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric); const auto ex3 = ex(rational{1, 4}) * ex(L"g", bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, Symmetry::Antisymm) * @@ -421,8 +428,8 @@ TEST_CASE("spin", "[spin]") { Symmetry::Antisymm) * ex(L"g", bra{L"a_1", L"a_2"}, ket{L"a_3", L"a_4"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_4"}, ket{L"i_2"}); + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_4"}, ket{L"i_2"}, particle_symmetric); result = expand_A_op(input); REQUIRE(result->is()); REQUIRE(result->size() == 4); @@ -437,10 +444,10 @@ TEST_CASE("spin", "[spin]") { Symmetry::Antisymm) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"a_3", L"a_4"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_4"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_3"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}); + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_4"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_3"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric); result = expand_A_op(input); REQUIRE_THAT( result, @@ -560,11 +567,12 @@ SECTION("Expand Symmetrizer") { ex(symm_label(), bra{L"i_1", L"i_2", L"i_3"}, ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Nonsymm) * ex(L"g", bra{L"i_4", L"i_5"}, ket{L"a_4", L"a_5"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"a_3"}, ket{L"i_4"}) * - ex(L"t", bra{L"a_5"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_4"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_1", L"a_2"}, ket{L"i_5", L"i_3"}); + particle_symmetric) * + ex(L"t", bra{L"a_3"}, ket{L"i_4"}, particle_symmetric) * + ex(L"t", bra{L"a_5"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_4"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_1", L"a_2"}, ket{L"i_5", L"i_3"}, + particle_symmetric); auto result = S_maps(input); REQUIRE(result->is()); REQUIRE(result->size() == 6); @@ -589,23 +597,23 @@ SECTION("partial expansion + S_maps = full expansion") { Symmetry::Antisymm); auto result = symmetrize_expr(input); REQUIRE_THAT( - result, SimplifiesTo( - "1/2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * t{a_1,a_2;i_1,i_2}:A-C-S " - "- 1/2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * t{a_1,a_2;i_2,i_1}:A-C-S")); + result, + SimplifiesTo("1/2 Ŝ{i_1,i_2;a_1,a_2} * t{a_1,a_2;i_1,i_2}:A-N-S " + "- 1/2 Ŝ{i_1,i_2;a_1,a_2} * t{a_1,a_2;i_2,i_1}:A-N-S")); result = S_maps(result); REQUIRE_THAT( result, SimplifiesTo( - "1/4 t{a_1,a_2;i_1,i_2}:A-C-S + 1/4 t{a_2,a_1;i_2,i_1}:A-C-S " - "- 1/4 t{a_1,a_2;i_2,i_1}:A-C-S - 1/4 t{a_2,a_1;i_1,i_2}:A-C-S")); + "1/4 t{a_1,a_2;i_1,i_2}:A-N-S + 1/4 t{a_2,a_1;i_2,i_1}:A-N-S " + "- 1/4 t{a_1,a_2;i_2,i_1}:A-N-S - 1/4 t{a_2,a_1;i_1,i_2}:A-N-S")); result = expand_A_op(input); REQUIRE_THAT( result, SimplifiesTo( - "1/4 t{a_1,a_2;i_1,i_2}:A-C-S - 1/4 t{a_1,a_2;i_2,i_1}:A-C-S " - "- 1/4 t{a_2,a_1;i_1,i_2}:A-C-S + 1/4 t{a_2,a_1;i_2,i_1}:A-C-S")); + "1/4 t{a_1,a_2;i_1,i_2}:A-N-S - 1/4 t{a_1,a_2;i_2,i_1}:A-N-S " + "- 1/4 t{a_2,a_1;i_1,i_2}:A-N-S + 1/4 t{a_2,a_1;i_2,i_1}:A-N-S")); } SECTION("partial spintracing + S_maps = full spintracing") { @@ -617,30 +625,31 @@ SECTION("partial spintracing + S_maps = full spintracing") { input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}); REQUIRE_THAT( result, - EquivalentTo("-2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * t{a_1,a_2;i_2,i_1}:N-C-S " - "+ 4 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * t{a_1,a_2;i_1,i_2}:N-C-S")); + EquivalentTo("-2 Ŝ{i_1,i_2;a_1,a_2} * t{a_1,a_2;i_2,i_1}:N-N-S " + "+ 4 Ŝ{i_1,i_2;a_1,a_2} * t{a_1,a_2;i_1,i_2}:N-N-S")); result = S_maps(result); REQUIRE_THAT( result, EquivalentTo( - "-1 t{a_1,a_2;i_2,i_1}:N-C-S - 1 t{a_2,a_1;i_1,i_2}:N-C-S " - "+ 2 t{a_1,a_2;i_1,i_2}:N-C-S + 2 t{a_2,a_1;i_2,i_1}:N-C-S")); + "-1 t{a_1,a_2;i_2,i_1}:N-N-S - 1 t{a_2,a_1;i_1,i_2}:N-N-S " + "+ 2 t{a_1,a_2;i_1,i_2}:N-N-S + 2 t{a_2,a_1;i_2,i_1}:N-N-S")); result = closed_shell_spintrace( input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}, true); REQUIRE_THAT( result, EquivalentTo( - "-1 t{a_1,a_2;i_2,i_1}:N-C-S - 1 t{a_2,a_1;i_1,i_2}:N-C-S " - "+ 2 t{a_1,a_2;i_1,i_2}:N-C-S + 2 t{a_2,a_1;i_2,i_1}:N-C-S")); + "-1 t{a_1,a_2;i_2,i_1}:N-N-S - 1 t{a_2,a_1;i_1,i_2}:N-N-S " + "+ 2 t{a_1,a_2;i_1,i_2}:N-N-S + 2 t{a_2,a_1;i_2,i_1}:N-N-S")); } SECTION("Symmetrize expression") { { // g * t1 + g * t1 - auto input = ex(L"g", bra{L"a_1", L"a_2"}, ket{L"i_1", L"a_3"}, - Symmetry::Symm) * - ex(L"t", bra{L"a_3"}, ket{L"i_2"}) + - ex(L"g", bra{L"a_2", L"a_1"}, ket{L"i_2", L"a_3"}, - Symmetry::Symm) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}); + auto input = + ex(L"g", bra{L"a_1", L"a_2"}, ket{L"i_1", L"a_3"}, + Symmetry::Symm) * + ex(L"t", bra{L"a_3"}, ket{L"i_2"}, particle_symmetric) + + ex(L"g", bra{L"a_2", L"a_1"}, ket{L"i_2", L"a_3"}, + Symmetry::Symm) * + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric); auto const ext_idxs = external_indices(input); auto bixs = ext_idxs | ranges::views::transform( @@ -648,7 +657,7 @@ SECTION("Symmetrize expression") { auto kixs = ext_idxs | ranges::views::transform( [](auto&& vec) { return get_ket_idx(vec); }); auto result = ex(Tensor{symm_label(), bra(std::move(kixs)), - ket(std::move(bixs))}) * + ket(std::move(bixs)), particle_symmetric}) * input; REQUIRE_THAT(result, EquivalentTo("2 Ŝ{i1,i2;a1,a2} g{a1,a2;i2,a3}:S t{a3;i1}")); @@ -664,22 +673,23 @@ SECTION("Symmetrize expression") { auto kixs = ext_idxs | ranges::views::transform( [](auto&& vec) { return get_ket_idx(vec); }); auto result = ex(Tensor{symm_label(), bra(std::move(kixs)), - ket(std::move(bixs))}) * + ket(std::move(bixs)), particle_symmetric}) * input; REQUIRE_THAT(result, EquivalentTo("Ŝ{i1,i2;a1,a2} * g{a1,a2;i1,i2}:S")); } { // g * t1 * t1 * t1 + g * t1 * t1 * t1 - auto input = ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_1", L"a_3"}, - Symmetry::Symm) * - ex(L"t", bra{L"a_1"}, ket{L"i_3"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}) * - ex(L"t", bra{L"a_3"}, ket{L"i_2"}) + - ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_2", L"a_3"}, - Symmetry::Symm) * - ex(L"t", bra{L"a_2"}, ket{L"i_3"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_4"}) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}); + auto input = + ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_1", L"a_3"}, + Symmetry::Symm) * + ex(L"t", bra{L"a_1"}, ket{L"i_3"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric) * + ex(L"t", bra{L"a_3"}, ket{L"i_2"}, particle_symmetric) + + ex(L"g", bra{L"i_3", L"i_4"}, ket{L"i_2", L"a_3"}, + Symmetry::Symm) * + ex(L"t", bra{L"a_2"}, ket{L"i_3"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_4"}, particle_symmetric) * + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric); auto const ext_idxs = external_indices(input); auto bixs = ext_idxs | ranges::views::transform( @@ -688,7 +698,7 @@ SECTION("Symmetrize expression") { [](auto&& vec) { return get_ket_idx(vec); }); auto result = ex(Tensor{reserved::symm_label(), bra(std::move(kixs)), - ket(std::move(bixs))}) * + ket(std::move(bixs)), particle_symmetric}) * input; REQUIRE_THAT( result, @@ -701,15 +711,17 @@ SECTION("Symmetrize expression") { ex(2) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"a_3", L"a_4"}, Symmetry::Symm) * - ex(L"t", bra{L"a_3"}, ket{L"i_3"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_4"}) * - ex(L"t", bra{L"a_1", L"a_4"}, ket{L"i_1", L"i_2"}) + + ex(L"t", bra{L"a_3"}, ket{L"i_3"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_4"}, particle_symmetric) * + ex(L"t", bra{L"a_1", L"a_4"}, ket{L"i_1", L"i_2"}, + particle_symmetric) + ex(2) * ex(L"g", bra{L"i_3", L"i_4"}, ket{L"a_3", L"a_4"}, Symmetry::Symm) * - ex(L"t", bra{L"a_3"}, ket{L"i_3"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_4"}) * - ex(L"t", bra{L"a_2", L"a_4"}, ket{L"i_2", L"i_1"}); + ex(L"t", bra{L"a_3"}, ket{L"i_3"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_4"}, particle_symmetric) * + ex(L"t", bra{L"a_2", L"a_4"}, ket{L"i_2", L"i_1"}, + particle_symmetric); auto const ext_idxs = external_indices(input); auto bixs = ext_idxs | ranges::views::transform( @@ -718,7 +730,7 @@ SECTION("Symmetrize expression") { [](auto&& vec) { return get_ket_idx(vec); }); auto result = ex(Tensor{reserved::symm_label(), bra(std::move(kixs)), - ket(std::move(bixs))}) * + ket(std::move(bixs)), particle_symmetric}) * input; simplify(result); REQUIRE(result->is() == false); @@ -739,7 +751,7 @@ SECTION("Swap bra kets") { // Tensor { auto input = ex(L"g", bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, - Symmetry::Nonsymm); + particle_symmetric); auto result = swap_bra_ket(input); REQUIRE_THAT(result, EquivalentTo("g{a1,a2;i1,i2}")); } @@ -747,18 +759,19 @@ SECTION("Swap bra kets") { // Product { auto input = ex(L"g", bra{L"a_5", L"a_6"}, ket{L"i_5", L"i_6"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"i_2"}, ket{L"a_6"}); + particle_symmetric) * + ex(L"t", bra{L"i_2"}, ket{L"a_6"}, particle_symmetric); auto result = swap_bra_ket(input); REQUIRE_THAT(result, EquivalentTo("g{i5,i6;a5,a6} t{a6;i2}")); } // Sum { - auto input = ex(L"f", bra{L"i_1"}, ket{L"i_5"}) + - ex(L"g", bra{L"a_5", L"a_6"}, ket{L"i_5", L"i_6"}, - Symmetry::Nonsymm) * - ex(L"t", bra{L"i_2"}, ket{L"a_6"}); + auto input = + ex(L"f", bra{L"i_1"}, ket{L"i_5"}, particle_symmetric) + + ex(L"g", bra{L"a_5", L"a_6"}, ket{L"i_5", L"i_6"}, + particle_symmetric) * + ex(L"t", bra{L"i_2"}, ket{L"a_6"}, particle_symmetric); auto result = swap_bra_ket(input); // TODO: This should be EquivalentTo but canonicalization currently doesn't // permit expressions that break co/contra variance of index contractions. @@ -801,9 +814,9 @@ SECTION("Closed-shell spintrace CCD") { closed_shell_CC_spintrace_v1(pno_ccd_energy_so_as_sum); REQUIRE_THAT(pno_ccd_energy_sf, EquivalentTo("2 g{a1,a2;i1,i2}:N-C " - "t{i1,i2;a1,a2}:N-C - " + "t{i1,i2;a1,a2}:N-N - " "g{a1,a2;i1,i2}:N-C " - "t{i1,i2;a2,a1}:N-C")); + "t{i1,i2;a2,a1}:N-N")); } { // CSV (aka PNO) for more compact equations const auto pno_ccd_energy_so = deserialize( @@ -817,9 +830,9 @@ SECTION("Closed-shell spintrace CCD") { closed_shell_CC_spintrace_v2(pno_ccd_energy_so_as_sum); REQUIRE_THAT(pno_ccd_energy_sf, EquivalentTo("2 g{a1,a2;i1,i2}:N-C " - "t{i1,i2;a1,a2}:N-C - " + "t{i1,i2;a1,a2}:N-N - " "g{a1,a2;i1,i2}:N-C " - "t{i1,i2;a2,a1}:N-C")); + "t{i1,i2;a2,a1}:N-N")); } } } @@ -828,8 +841,9 @@ SECTION("Closed-shell spintrace CCSD") { // These terms from CCSD R1 equations { // A * f - const auto input = ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"f", bra{L"a_1"}, ket{L"i_1"}); + const auto input = + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"f", bra{L"a_1"}, ket{L"i_1"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -840,10 +854,11 @@ SECTION("Closed-shell spintrace CCSD") { { // - A * f * t1 - const auto input = ex(-1) * - ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"f", bra{L"i_2"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_2"}); + const auto input = + ex(-1) * + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"f", bra{L"i_2"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_2"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -852,9 +867,10 @@ SECTION("Closed-shell spintrace CCSD") { { // A * f * t1 - const auto input = ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"f", bra{L"a_1"}, ket{L"a_2"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_1"}); + const auto input = + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"f", bra{L"a_1"}, ket{L"a_2"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_1"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -911,11 +927,12 @@ SECTION("Closed-shell spintrace CCSD") { { // A * g * t1 * t1 - const auto input = ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"g", bra{L"i_2", L"a_1"}, - ket{L"a_2", L"a_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}); + const auto input = + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"g", bra{L"i_2", L"a_1"}, ket{L"a_2", L"a_3"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -925,11 +942,12 @@ SECTION("Closed-shell spintrace CCSD") { { // A * g * t2 * t2 - const auto input = ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"g", bra{L"i_2", L"i_3"}, - ket{L"i_1", L"a_2"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_3"}); + const auto input = + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"g", bra{L"i_2", L"i_3"}, ket{L"i_1", L"a_2"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_3"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -939,11 +957,12 @@ SECTION("Closed-shell spintrace CCSD") { { // A * f * t1 * t1 - const auto input = ex(-1) * - ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"f", bra{L"i_2"}, ket{L"a_2"}) * - ex(L"t", bra{L"a_2"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_2"}); + const auto input = + ex(-1) * + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"f", bra{L"i_2"}, ket{L"a_2"}, particle_symmetric) * + ex(L"t", bra{L"a_2"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_2"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -952,13 +971,14 @@ SECTION("Closed-shell spintrace CCSD") { { // -1/2 * A * g * t1 * t2 - const auto input = ex(rational{-1, 2}) * - ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"g", bra{L"i_2", L"i_3"}, - ket{L"a_2", L"a_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_1"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_2", L"a_3"}, - ket{L"i_1", L"i_3"}, Symmetry::Antisymm); + const auto input = + ex(rational{-1, 2}) * + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"g", bra{L"i_2", L"i_3"}, ket{L"a_2", L"a_3"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_1"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_2", L"a_3"}, ket{L"i_1", L"i_3"}, + Symmetry::Antisymm); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -969,13 +989,14 @@ SECTION("Closed-shell spintrace CCSD") { { // -1/2 * A * g * t1 * t2 - const auto input = ex(rational{-1, 2}) * - ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"g", bra{L"i_2", L"i_3"}, - ket{L"a_2", L"a_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_1", L"a_3"}, - ket{L"i_2", L"i_3"}, Symmetry::Antisymm); + const auto input = + ex(rational{-1, 2}) * + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"g", bra{L"i_2", L"i_3"}, ket{L"a_2", L"a_3"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_2", L"i_3"}, + Symmetry::Antisymm); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -986,13 +1007,14 @@ SECTION("Closed-shell spintrace CCSD") { { // A * g * t1 * t2 - const auto input = ex(1) * - ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * - ex(L"g", bra{L"i_2", L"i_3"}, - ket{L"a_2", L"a_3"}, Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_1", L"a_3"}, - ket{L"i_1", L"i_3"}, Symmetry::Antisymm); + const auto input = + ex(1) * + ex(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) * + ex(L"g", bra{L"i_2", L"i_3"}, ket{L"a_2", L"a_3"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_1", L"a_3"}, ket{L"i_1", L"i_3"}, + Symmetry::Antisymm); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -1005,12 +1027,13 @@ SECTION("Closed-shell spintrace CCSD") { { // - A * g * t1 * t1 * t1 - auto input = ex(-1) * - ex(L"g", bra{L"i_2", L"i_3"}, ket{L"a_2", L"a_3"}, - Symmetry::Antisymm) * - ex(L"t", bra{L"a_2"}, ket{L"i_2"}) * - ex(L"t", bra{L"a_3"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_1"}, ket{L"i_3"}); + auto input = + ex(-1) * + ex(L"g", bra{L"i_2", L"i_3"}, ket{L"a_2", L"a_3"}, + Symmetry::Antisymm) * + ex(L"t", bra{L"a_2"}, ket{L"i_2"}, particle_symmetric) * + ex(L"t", bra{L"a_3"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1"}, ket{L"i_3"}, particle_symmetric); auto result = ex(rational{1, 2}) * spintrace(input, IdxGroupList{{L"i_1", L"a_1"}}); expand(result); @@ -1072,7 +1095,7 @@ SECTION("Closed-shell spintrace CCSDT terms") { ket{L"a_1", L"a_2", L"a_3"}, Symmetry::Antisymm) * ex(L"t", bra{L"a_1", L"a_2", L"a_3"}, ket{L"i_2", L"i_3", L"i_4"}, Symmetry::Antisymm) * - ex(L"f", bra{L"i_4"}, ket{L"i_1"}); + ex(L"f", bra{L"i_4"}, ket{L"i_1"}, particle_symmetric); auto result = expand_A_op(input); REQUIRE(result->size() == 36); @@ -1083,14 +1106,14 @@ SECTION("Closed-shell spintrace CCSDT terms") { REQUIRE_THAT( result, - EquivalentTo("24 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * f{i_4;i_3}:N-C-S * " - "t{a_1,a_2,a_3;i_1,i_2,i_4}:N-C-S - 12" - " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * f{i_4;i_3}:N-C-S * " - "t{a_1,a_2,a_3;i_2,i_1,i_4}:N-C-S + 12" - " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * f{i_4;i_1}:N-C-S * " - "t{a_1,a_2,a_3;i_2,i_3,i_4}:N-C-S - 24" - " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * f{i_4;i_2}:N-C-S * " - "t{a_1,a_2,a_3;i_1,i_3,i_4}:N-C-S ")); + EquivalentTo("24 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * f{i_4;i_3}:N-N-S * " + "t{a_1,a_2,a_3;i_1,i_2,i_4}:N-N-S - 12" + " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * f{i_4;i_3}:N-N-S * " + "t{a_1,a_2,a_3;i_2,i_1,i_4}:N-N-S + 12" + " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * f{i_4;i_1}:N-N-S * " + "t{a_1,a_2,a_3;i_2,i_3,i_4}:N-N-S - 24" + " Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * f{i_4;i_2}:N-N-S * " + "t{a_1,a_2,a_3;i_1,i_3,i_4}:N-N-S ")); } SECTION( @@ -1116,25 +1139,25 @@ SECTION("Closed-shell spintrace CCSDT terms") { REQUIRE_THAT( result_1, EquivalentTo( - " 8 g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-C-S + " + " 8 g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-N-S + " "2" - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4" - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_3,i_1,i_2}:N-C-S - 4" - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_1,i_3,i_2}:N-C-S - 4 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_2,i_1,i_3}:N-C-S - 4 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_2,i_1,i_3}:N-C-S + 2 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_3,i_1,i_2}:N-C-S - 4" - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_3,i_2,i_1}:N-C-S + 2 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_3,i_2}:N-C-S - 4 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_2,i_3}:N-C-S + 8 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_1,i_2,i_3}:N-C-S + 8 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_2,i_1,i_3}:N-C-S + 2 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_3,i_2,i_1}:N-C-S - 4 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_3,i_2,i_1}:N-C-S + 2 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-C-S + 2" - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-C-S")); + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4" + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_3,i_1,i_2}:N-N-S - 4" + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_1,i_3,i_2}:N-N-S - 4 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_2,i_1,i_3}:N-N-S - 4 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_2,i_1,i_3}:N-N-S + 2 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_3,i_1,i_2}:N-N-S - 4" + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_3,i_2,i_1}:N-N-S + 2 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_1,i_3,i_2}:N-N-S - 4 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_1,i_2,i_3}:N-N-S + 8 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_1,i_2,i_3}:N-N-S + 8 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_2,i_1,i_3}:N-N-S + 2 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_3,i_2,i_1}:N-N-S - 4 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_3,i_2,i_1}:N-N-S + 2 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-N-S + 2" + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-N-S")); // the new efficient method, does spintracing with partial expansion, then // expanding by S_map (this method is used in @@ -1150,24 +1173,24 @@ SECTION("Closed-shell spintrace CCSDT terms") { REQUIRE_THAT( result_2, EquivalentTo( - "8 g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-C-S + 2" - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_3,i_1,i_2}:N-C-S - 4 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_1,i_3,i_2}:N-C-S - 4 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_2,i_1,i_3}:N-C-S - 4" - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_2,i_1,i_3}:N-C-S+ 2 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_3,i_1,i_2}:N-C-S - 4 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_3,i_2,i_1}:N-C-S + 2 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4" - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_3,i_2}:N-C-S - 4" - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_2,i_3}:N-C-S + 8 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_1,i_2,i_3}:N-C-S + 8 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_2,i_1,i_3}:N-C-S + 2 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_3,i_2,i_1}:N-C-S - 4 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_2,i_3,i_1}:N-C-S - 4 " - "g{a_2,a_3;a_4,a_5}:N-C-S * t{a_1,a_4,a_5;i_3,i_2,i_1}:N-C-S + 2 " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-C-S + 2 " - "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-C-S")); + "8 g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-N-S + 2" + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_3,i_1,i_2}:N-N-S - 4 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_1,i_3,i_2}:N-N-S - 4 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_2,i_1,i_3}:N-N-S - 4" + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_2,i_1,i_3}:N-N-S+ 2 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_3,i_1,i_2}:N-N-S - 4 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_3,i_2,i_1}:N-N-S + 2 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4" + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_1,i_3,i_2}:N-N-S - 4" + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_1,i_2,i_3}:N-N-S + 8 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_1,i_2,i_3}:N-N-S + 8 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_2,i_1,i_3}:N-N-S + 2 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_3,i_2,i_1}:N-N-S - 4 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_2,i_3,i_1}:N-N-S - 4 " + "g{a_2,a_3;a_4,a_5}:N-N-S * t{a_1,a_4,a_5;i_3,i_2,i_1}:N-N-S + 2 " + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-N-S + 2 " + "g{a_1,a_3;a_4,a_5}:N-N-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-N-S")); } SECTION("the most expensive terms in CCSDT in v2") { // results in 1 term @@ -1181,7 +1204,7 @@ SECTION("Closed-shell spintrace CCSDT terms") { result, EquivalentTo( L"3 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-C-S")); + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_3,i_1,i_2}:N-N-S")); } SECTION("most expensive CCSDT term in v1") { // results in 4 terms @@ -1196,13 +1219,13 @@ SECTION("Closed-shell spintrace CCSDT terms") { result, EquivalentTo( L"-6/5 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * " - "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-C-S + " - "3 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * g{a_1,a_2;a_4,a_5}:N-C-S *" - " t{a_3,a_4,a_5;i_3,i_1,i_2}:N-C-S -" - " 3/5 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * g{a_1,a_2;a_4,a_5}:N-C-S *" - " t{a_3,a_4,a_5;i_3,i_2,i_1}:N-C-S -" - " 6/5 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3}:N-C-S * g{a_1,a_2;a_4,a_5}:N-C-S *" - " t{a_3,a_4,a_5;i_2,i_1,i_3}:N-C-S")); + "g{a_1,a_2;a_4,a_5}:N-N-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-N-S + " + "3 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * g{a_1,a_2;a_4,a_5}:N-N-S *" + " t{a_3,a_4,a_5;i_3,i_1,i_2}:N-N-S -" + " 3/5 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * g{a_1,a_2;a_4,a_5}:N-N-S *" + " t{a_3,a_4,a_5;i_3,i_2,i_1}:N-N-S -" + " 6/5 Ŝ{i_1,i_2,i_3;a_1,a_2,a_3} * g{a_1,a_2;a_4,a_5}:N-N-S *" + " t{a_3,a_4,a_5;i_2,i_1,i_3}:N-N-S")); } SECTION("f * t3") { @@ -1478,8 +1501,9 @@ SECTION("Open-shell spin-tracing") { // Tensor canonicalize { - auto t3 = ex(Tensor(L"t", bra{a3A, a2B, a2A}, ket{i1A, i2B, i3A})); - auto f = ex(Tensor(L"f", bra{a1A}, ket{a2A})); + auto t3 = ex(Tensor(L"t", bra{a3A, a2B, a2A}, ket{i1A, i2B, i3A}, + particle_symmetric)); + auto f = ex(Tensor(L"f", bra{a1A}, ket{a2A}, particle_symmetric)); auto ft3 = f * t3; ft3->canonicalize(); REQUIRE_THAT(ft3, EquivalentTo("f{a↑1;a↑2} t{a↑3,a↑2,a↓2;i↑1,i↑3,i↓2}")); @@ -1500,10 +1524,11 @@ SECTION("Open-shell spin-tracing") { // f_oo * t2 { - auto input = ex(rational{1, 2}) * - ex(L"f", bra{L"i_3"}, ket{L"i_1"}) * - ex(L"t", bra{L"a_1", L"a_2"}, ket{L"i_2", L"i_3"}, - Symmetry::Antisymm); + auto input = + ex(rational{1, 2}) * + ex(L"f", bra{L"i_3"}, ket{L"i_1"}, particle_symmetric) * + ex(L"t", bra{L"a_1", L"a_2"}, ket{L"i_2", L"i_3"}, + Symmetry::Antisymm); auto result = open_shell_spintrace( input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}); @@ -1538,10 +1563,11 @@ SECTION("Open-shell spin-tracing") { // f * t3 { - auto input = ex(rational{1, 12}) * - ex(L"f", bra{L"a_1"}, ket{L"a_4"}) * - ex(L"t", bra{L"a_2", L"a_3", L"a_4"}, - ket{L"i_1", L"i_2", L"i_3"}, Symmetry::Antisymm); + auto input = + ex(rational{1, 12}) * + ex(L"f", bra{L"a_1"}, ket{L"a_4"}, particle_symmetric) * + ex(L"t", bra{L"a_2", L"a_3", L"a_4"}, + ket{L"i_1", L"i_2", L"i_3"}, Symmetry::Antisymm); auto result = open_shell_spintrace( input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}}); @@ -1566,26 +1592,26 @@ SECTION("Open-shell spin-tracing") { Symmetry::Antisymm); auto g = Tensor(L"g", bra{i3A, i4A}, ket{i1A, i2A}, Symmetry::Antisymm); - auto t3 = - Tensor(L"t", bra{a1A, a2A, a3B}, ket{i3A, i4A, i3B}, Symmetry::Nonsymm); + auto t3 = Tensor(L"t", bra{a1A, a2A, a3B}, ket{i3A, i4A, i3B}, + particle_symmetric); auto input = ex(A2_aab) * ex(g) * ex(t3); auto result = expand_A_op(input); result->visit(reset_idx_tags); REQUIRE_THAT(result, - EquivalentTo("-1 g{i↑_3,i↑_4;i↑_1,i↑_2}:A-C-S * " - "t{a↑_1,a↑_2,a↓_3;i↑_4,i↑_3,i↓_3}:N-C-S")); + EquivalentTo("-1 g{i↑_3,i↑_4;i↑_1,i↑_2}:A-N-S * " + "t{a↑_1,a↑_2,a↓_3;i↑_4,i↑_3,i↓_3}:N-N-S")); g = Tensor(L"g", bra{i4A, i5A}, ket{i1A, i2A}, Symmetry::Antisymm); - t3 = - Tensor(L"t", bra{a1A, a2A, a3B}, ket{i4A, i5A, i3B}, Symmetry::Nonsymm); + t3 = Tensor(L"t", bra{a1A, a2A, a3B}, ket{i4A, i5A, i3B}, + particle_symmetric); input = ex(A2_aab) * ex(g) * ex(t3); result = expand_A_op(input); result->visit(reset_idx_tags); REQUIRE_THAT(result, - EquivalentTo("-1 g{i↑_3,i↑_4;i↑_1,i↑_2}:A-C-S * " - "t{a↑_1,a↑_2,a↓_3;i↑_4,i↑_3,i↓_3}:N-C-S")); + EquivalentTo("-1 g{i↑_3,i↑_4;i↑_1,i↑_2}:A-N-S * " + "t{a↑_1,a↑_2,a↓_3;i↑_4,i↑_3,i↓_3}:N-N-S")); } // CCSDT R3 10 aaa, bbb diff --git a/tests/unit/test_tensor.cpp b/tests/unit/test_tensor.cpp index aa9ecb6198..8637062564 100644 --- a/tests/unit/test_tensor.cpp +++ b/tests/unit/test_tensor.cpp @@ -25,6 +25,9 @@ #include #include +// the `particle_symmetric` symmetry pack (column = Symm) is defined in +// catch2_sequant.hpp and shared across the MBPT test TUs + TEST_CASE("tensor", "[elements]") { using namespace sequant; @@ -59,8 +62,12 @@ TEST_CASE("tensor", "[elements]") { REQUIRE(ranges::distance(t2.const_indices().begin(), t2.const_indices().end()) == 2); REQUIRE(t2.symmetry() == Symmetry::Nonsymm); - REQUIRE(t2.braket_symmetry() == BraKetSymmetry::Conjugate); - REQUIRE(t2.column_symmetry() == ColumnSymmetry::Symm); + // programmatic ctor defaults are fixed and Context-independent: braket + // symmetry derives from the default NonHermitian -> Nonsymm, and column + // symmetry defaults to the safe Nonsymm (the Context's Symm default applies + // only to deserialized tensors) + REQUIRE(t2.braket_symmetry() == BraKetSymmetry::Nonsymm); + REQUIRE(t2.column_symmetry() == ColumnSymmetry::Nonsymm); REQUIRE(t2.label() == L"F"); // bra/kets of different rank @@ -82,8 +89,8 @@ TEST_CASE("tensor", "[elements]") { REQUIRE(ranges::distance(t3.const_indices().begin(), t3.const_indices().end()) == 2); REQUIRE(t3.symmetry() == Symmetry::Nonsymm); - REQUIRE(t3.braket_symmetry() == BraKetSymmetry::Conjugate); - REQUIRE(t3.column_symmetry() == ColumnSymmetry::Symm); + REQUIRE(t3.braket_symmetry() == BraKetSymmetry::Nonsymm); + REQUIRE(t3.column_symmetry() == ColumnSymmetry::Nonsymm); REQUIRE(t3.label() == L"N"); REQUIRE_NOTHROW(Tensor(L"g", bra{Index{L"i_1"}, Index{L"i_2"}}, @@ -110,14 +117,14 @@ TEST_CASE("tensor", "[elements]") { SECTION("null indices") { // null indices ok in asymmetric bra or ket REQUIRE_NOTHROW(Tensor(L"N", bra{L"i_2", L"", L"i_3"}, - ket{L"", L"i_1", L""}, aux{}, Symmetry::Nonsymm)); + ket{L"", L"i_1", L""}, aux{}, particle_symmetric)); REQUIRE_NOTHROW(Tensor(L"N", bra{L"", L"i_1", L""}, ket{L"i_2", L"", L"i_3"}, aux{}, - Symmetry::Nonsymm)); + particle_symmetric)); REQUIRE_NOTHROW( - Tensor(L"N", bra{L""}, ket{L"i_1"}, aux{}, Symmetry::Nonsymm)); - Tensor t(L"N", bra{L""}, ket{L"i_1"}, aux{}, Symmetry::Nonsymm); + Tensor(L"N", bra{L""}, ket{L"i_1"}, aux{}, particle_symmetric)); + Tensor t(L"N", bra{L""}, ket{L"i_1"}, aux{}, particle_symmetric); REQUIRE(t.bra_rank() == 0); REQUIRE(t.ket_rank() == 1); REQUIRE(t.bra_net_rank() == 0); @@ -127,29 +134,29 @@ TEST_CASE("tensor", "[elements]") { // in fact slots of asymmetric particle-symmetric tensors are kept in // canonical order - REQUIRE(Tensor(L"N", bra{L""}, ket{L"i_1"}, aux{}, Symmetry::Nonsymm) == - Tensor(L"N", bra{}, ket{L"i_1"}, aux{}, Symmetry::Nonsymm)); + REQUIRE(Tensor(L"N", bra{L""}, ket{L"i_1"}, aux{}, particle_symmetric) == + Tensor(L"N", bra{}, ket{L"i_1"}, aux{}, particle_symmetric)); REQUIRE(Tensor(L"N", bra{L"i_2", L"", L"i_3"}, ket{L"", L"i_1", L""}, - aux{}, Symmetry::Nonsymm) == + aux{}, particle_symmetric) == Tensor(L"N", bra{L"i_2", L"i_3"}, ket{L"", L"", L"i_1"}, aux{}, - Symmetry::Nonsymm)); + particle_symmetric)); REQUIRE(Tensor(L"N", bra{L"", L"i_1", L""}, ket{L"i_2", L"", L"i_3"}, - aux{}, Symmetry::Nonsymm) == + aux{}, particle_symmetric) == Tensor(L"N", bra{L"i_1", L"", L""}, ket{L"", L"i_2", L"i_3"}, - aux{}, Symmetry::Nonsymm)); + aux{}, particle_symmetric)); REQUIRE( Tensor(L"N", bra{L"", L"i_1", L"", L"i_4"}, - ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, Symmetry::Nonsymm) == + ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, particle_symmetric) == Tensor(L"N", bra{L"i_4", L"i_1", L"", L""}, - ket{L"i_5", L"", L"i_2", L"i_3"}, aux{}, Symmetry::Nonsymm)); + ket{L"i_5", L"", L"i_2", L"i_3"}, aux{}, particle_symmetric)); // in fact unnecessary null indices are dropped in canonicalization REQUIRE( Tensor(L"N", bra{L"", L"i_1", L"", L"i_4"}, - ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, Symmetry::Nonsymm) == + ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, particle_symmetric) == Tensor(L"N", bra{L"i_4", L"i_1"}, ket{L"i_5", L"", L"i_2", L"i_3"}, - aux{}, Symmetry::Nonsymm)); + aux{}, particle_symmetric)); Tensor t5(L"N", bra{L"", L"i_1", L"", L"i_4"}, - ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, Symmetry::Nonsymm); + ket{L"i_2", L"", L"i_3", L"i_5"}, aux{}, particle_symmetric); REQUIRE(t5.bra_rank() == 2); REQUIRE(t5.bra_net_rank() == 2); REQUIRE(t5.bra()[0] == L"i_4"); @@ -344,7 +351,7 @@ TEST_CASE("tensor", "[elements]") { auto t2 = Tensor(L"F", bra{L"i_1"}, ket{L"i_2"}, aux{L"i_3"}); auto t3 = Tensor(L"F", bra{L"i_1"}, ket{L"i_2"}, aux{L"i_3", L"i_4"}); auto t4 = Tensor(L"F", bra{Index(L"i_1", {L"i_5", L"i_6"}), Index{}}, - ket{L"", L"i_2"}, aux{L"i_3", L"i_4"}, Symmetry::Nonsymm); + ket{L"", L"i_2"}, aux{L"i_3", L"i_4"}, particle_symmetric); auto h1 = ex(L"F", bra{L"i_1"}, ket{L"i_2"}) * ex(cre({L"i_1"}), ann({L"i_2"})); @@ -432,7 +439,9 @@ TEST_CASE("tensor", "[elements]") { SECTION("adjoint") { auto f1 = Tensor(L"F", bra{L"i_1", L"i_2"}, ket{L"i_3", L"i_4"}); REQUIRE_NOTHROW(f1.adjoint()); - REQUIRE(to_latex(f1) == L"{F^{{i_1}{i_2}}_{{i_3}{i_4}}}"); + // F is now non-Hermitian by default (braket Nonsymm), so its adjoint is + // marked with the conjugation superscript + REQUIRE(to_latex(f1) == L"{F⁺^{{i_1}{i_2}}_{{i_3}{i_4}}}"); auto t1 = Tensor(L"t", bra{L"a_1"}, ket{L"i_1"}, Symmetry::Nonsymm, BraKetSymmetry::Nonsymm); @@ -445,7 +454,7 @@ TEST_CASE("tensor", "[elements]") { ex(cre{L"i_1"}, ann{L"i_2"}); h1 = adjoint(h1); REQUIRE(to_latex(h1) == - L"{{\\tilde{a}^{{i_2}}_{{i_1}}}{F^{{i_1}}_{{i_2}}}}"); + L"{{\\tilde{a}^{{i_2}}_{{i_1}}}{F⁺^{{i_1}}_{{i_2}}}}"); h1 = adjoint(h1); REQUIRE(to_latex(h1) == L"{{F^{{i_2}}_{{i_1}}}{\\tilde{a}^{{i_1}}_{{i_2}}}}"); @@ -551,9 +560,11 @@ TEST_CASE("tensor_hermiticity", "[elements]") { }; auto make_diff = [&idx](Field f) { auto a = ex(L"g", bra{idx(L"i_1", f)}, ket{idx(L"i_2", f)}, - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); auto b = ex(L"g", bra{idx(L"i_2", f)}, ket{idx(L"i_1", f)}, - Symmetry::Nonsymm, Hermiticity::Hermitian); + Symmetry::Nonsymm, Hermiticity::Hermitian, + ColumnSymmetry::Symm); ExprPtr diff = a - b; simplify(diff); return diff; @@ -564,4 +575,95 @@ TEST_CASE("tensor_hermiticity", "[elements]") { // complex spaces -> Conjugate -> the two orientations remain distinct REQUIRE_FALSE(make_diff(Field::Complex) == ex(0)); } + + SECTION("programmatic ctor defaults are Context-independent") { + // The Context's default symmetries govern *deserialization* only (see + // test_parse). A programmatic Tensor ctor always resolves unspecified + // attributes against fixed library defaults -- fully non-symmetric / + // non-Hermitian -- regardless of the active Context. Set the Context + // defaults to non-default values and confirm they are ignored here. + auto ctx = get_default_context(); + ctx.set(Symmetry::Antisymm) + .set(Hermiticity::Hermitian) + .set(ColumnSymmetry::Symm); + auto resetter = set_scoped_default_context(ctx); + + auto g = Tensor(L"g", bra{L"i_1", L"i_2"}, ket{L"i_3", L"i_4"}); + REQUIRE(g.symmetry() == Symmetry::Nonsymm); + REQUIRE(g.hermiticity() == Hermiticity::NonHermitian); + REQUIRE(g.column_symmetry() == ColumnSymmetry::Nonsymm); + REQUIRE(g.braket_symmetry() == + to_braket_symmetry(Hermiticity::NonHermitian, g.base_field())); + + // explicitly-specified symmetries override the fixed defaults + auto t = Tensor(L"t", bra{L"a_1"}, ket{L"i_1"}, Symmetry::Nonsymm, + BraKetSymmetry::Nonsymm, ColumnSymmetry::Symm); + REQUIRE(t.braket_symmetry() == BraKetSymmetry::Nonsymm); + REQUIRE(t.hermiticity() == Hermiticity::NonHermitian); + REQUIRE(t.column_symmetry() == ColumnSymmetry::Symm); + } +} + +TEST_CASE("(anti)symmetrizer factories", "[elements]") { + using namespace sequant; + + SECTION("defining symmetries") { + const auto S = make_symmetrizer(bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}); + REQUIRE(S->as().label() == reserved::symm_label()); + REQUIRE(S->as().symmetry() == Symmetry::Nonsymm); + REQUIRE(S->as().braket_symmetry() == BraKetSymmetry::Nonsymm); + REQUIRE(S->as().column_symmetry() == ColumnSymmetry::Symm); + + const auto A = + make_antisymmetrizer(bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}); + REQUIRE(A->as().label() == reserved::antisymm_label()); + REQUIRE(A->as().symmetry() == Symmetry::Antisymm); + REQUIRE(A->as().braket_symmetry() == BraKetSymmetry::Nonsymm); + REQUIRE(A->as().column_symmetry() == ColumnSymmetry::Symm); + + // the factories agree with what the deserializer produces + REQUIRE(*S == *deserialize(L"Ŝ{i_1,i_2;a_1,a_2}")); + REQUIRE(*A == *deserialize(L"Â{i_1,i_2;a_1,a_2}")); + + // ... and with the aux-index overloads + REQUIRE_NOTHROW(make_symmetrizer(bra{L"i_1"}, ket{L"a_1"}, aux{L"x_1"})); + REQUIRE_NOTHROW( + make_antisymmetrizer(bra{L"i_1"}, ket{L"a_1"}, aux{L"x_1"})); + } + + SECTION("ctor rejects contradicting symmetries") { + // braket symmetry, like column symmetry below, is a defining property of a + // reserved (anti)symmetrizer, not a free parameter + REQUIRE_THROWS_AS( + Tensor(reserved::symm_label(), bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, + TensorSymmetries{.braket = BraKetSymmetry::Symm}), + Exception); + REQUIRE_THROWS_AS( + Tensor(reserved::symm_label(), bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}, + TensorSymmetries{.column = ColumnSymmetry::Nonsymm}), + Exception); + // ... but for  a Nonsymm column request is *promoted* rather than + // rejected: Antisymm bra/ket permutational symmetry implies column + // symmetry for any tensor, reserved or not, and that promotion runs first + REQUIRE(Tensor(reserved::antisymm_label(), bra{L"i_1", L"i_2"}, + ket{L"a_1", L"a_2"}, + TensorSymmetries{.perm = Symmetry::Antisymm, + .column = ColumnSymmetry::Nonsymm}) + .column_symmetry() == ColumnSymmetry::Symm); + REQUIRE_THROWS_AS( + Tensor(reserved::antisymm_label(), bra{L"i_1", L"i_2"}, + ket{L"a_1", L"a_2"}, + TensorSymmetries{.perm = Symmetry::Nonsymm, + .column = ColumnSymmetry::Nonsymm}), + Exception); + + // an *unspecified* column symmetry is silently supplied, so that the + // library's Context-independent column default (Nonsymm) does not produce + // a symmetrizer that differs from the deserialized one + const auto S = Tensor(reserved::symm_label(), bra{L"i_1", L"i_2"}, + ket{L"a_1", L"a_2"}); + REQUIRE(S.column_symmetry() == ColumnSymmetry::Symm); + REQUIRE(S == make_symmetrizer(bra{L"i_1", L"i_2"}, ket{L"a_1", L"a_2"}) + ->as()); + } } diff --git a/tests/unit/test_tensor_network.cpp b/tests/unit/test_tensor_network.cpp index f4be9861eb..97fb9aafe2 100644 --- a/tests/unit/test_tensor_network.cpp +++ b/tests/unit/test_tensor_network.cpp @@ -283,10 +283,10 @@ TEMPLATE_TEST_CASE("tensor_network_shared", "[elements]", TensorNetworkV1, std::vector>> tests{ {L"G{;;a1,a2,a3,a4} T{;;i3,i2,a3,a4}", - v3 ? idxvec_t{L"i_3", L"i_2", L"a_2", L"a_1"} + v3 ? idxvec_t{L"i_2", L"i_3", L"a_2", L"a_1"} : idxvec_t{L"i_2", L"i_3", L"a_1", L"a_2"}}, {L"G{;;a1,a2,a3,a4} T{;;i2,i3,a3,a4}", - v3 ? idxvec_t{L"i_2", L"i_3", L"a_2", L"a_1"} + v3 ? idxvec_t{L"i_3", L"i_2", L"a_2", L"a_1"} : idxvec_t{L"i_3", L"i_2", L"a_1", L"a_2"}}, }; @@ -1182,12 +1182,11 @@ TEST_CASE("tensor_network_v2", "[elements][valgrind_skip]") { } SECTION("special") { - auto factors = - deserialize( - L"Ŝ{i_1;a_1}:N-C-S g{i_2,a_1;a_2,i_1}:N-C-S " - L"t{a_2;i_2}:N-C-S") - ->as() - .factors(); + auto factors = deserialize( + L"Ŝ{i_1;a_1} g{i_2,a_1;a_2,i_1}:N-C-S " + L"t{a_2;i_2}:N-C-S") + ->as() + .factors(); TensorNetworkV2 tn(factors); diff --git a/tests/unit/test_wick.cpp b/tests/unit/test_wick.cpp index ac71bca4cd..ee45ff5784 100644 --- a/tests/unit/test_wick.cpp +++ b/tests/unit/test_wick.cpp @@ -34,6 +34,10 @@ #include namespace sequant { + +// the particle-symmetric default symmetry pack `particle_symmetric` (column = +// Symm) is defined in catch2_sequant.hpp and shared across the MBPT test TUs + struct WickAccessor {}; template <> @@ -1015,9 +1019,9 @@ TEST_CASE("wick", "[algorithms][wick][valgrind_skip]") { // multiply tensor factors and expand auto wick_result_2 = ex(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"}, - Symmetry::Nonsymm) * + particle_symmetric) * ex(L"t", bra{L"a_4", L"a_5"}, ket{L"i_4", L"i_5"}, - Symmetry::Nonsymm) * + particle_symmetric) * wick_result; expand(wick_result_2); REQUIRE(wick_result_2->size() == 4); // still 4 terms @@ -1305,23 +1309,25 @@ TEST_CASE("wick", "[algorithms][wick][valgrind_skip]") { // sequant::wprintf(to_latex_align(Ld_H2N_L), L" = \n", // to_latex_align(result2, 0, 2), L"\n"); REQUIRE(result2.as().size() == 5); - // Sum-term ordering shifted vs master after the removal of - // Context::braket_symmetry: the canonical sort key for the two - // -1/2 ã·v̄·w terms swapped. The expression is mathematically - // unchanged. + // Sum-term ordering depends on the tensors' braket and column + // symmetries; the canonical sort key for these ã·v̄·w terms reflects the + // fixed (Context-independent) Tensor ctor defaults -- notably w, which + // carries no bra/ket and so is column-Nonsymm by default -- and shifted + // vs the previous (Context-Symm) ordering. The expression is + // mathematically unchanged (same five terms, reordered). REQUIRE( result2.to_latex() == L"{ \\bigl( - " - L"{{{\\frac{1}{4}}}{\\tilde{a}^{{p_3}{p_4}{p_5}}_{{p_1}{p_2}{p_5}" - L"}}{\\bar{v}^{{p_1}{p_2}}_{{p_3}{p_4}}}{w^{}_{}[{p_5}]}} - " - L"{{{\\frac{1}{2}}}{\\tilde{a}^{{p_2}{p_3}}_{{e_1}{p_1}}}{\\bar{" - L"v}^{{e_1}{p_1}}_{{p_2}{p_3}}}{w^{}_{}[{e_1}]}} + " - L"{{\\tilde{a}^{{p_2}}_{{p_1}}}{\\bar{v}^{{e_1}{p_1}}_{{e_1}{p_2}" - L"}}{w^{}_{}[{e_1}]}} - " - L"{{{\\frac{1}{2}}}{\\tilde{a}^{{e_1}{p_3}}_{{p_1}{p_2}}}{\\bar{" - L"v}^{{p_1}{p_2}}_{{e_1}{p_3}}}{w^{}_{}[{e_1}]}} + " - L"{{{\\frac{1}{4}}}{\\tilde{a}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{" - L"v}^{{p_1}{p_2}}_{{p_3}{p_4}}}{w^{}_{}[{e_1}]}}\\bigr) }"); + L"{{{\\frac{1}{4}}}{\\tilde{a}^{{p_3}{p_4}{p_5}}_{{p_1}{p_2}{p_5}}}" + L"{\\bar{v}^{{p_1}{p_2}}_{{p_3}{p_4}}}{w^{}_{}[{p_5}]}} - " + L"{{{\\frac{1}{2}}}{\\tilde{a}^{{e_1}{p_3}}_{{p_1}{p_2}}}{\\bar{v}^" + L"{{p_1}{p_2}}_{{e_1}{p_3}}}{w^{}_{}[{e_1}]}} + " + L"{{{\\frac{1}{4}}}{\\tilde{a}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{v}^" + L"{{p_1}{p_2}}_{{p_3}{p_4}}}{w^{}_{}[{e_1}]}} + " + L"{{\\tilde{a}^{{p_2}}_{{p_1}}}{\\bar{v}^{{e_1}{p_1}}_{{e_1}{p_2}}}" + L"{w^{}_{}[{e_1}]}} - " + L"{{{\\frac{1}{2}}}{\\tilde{a}^{{p_2}{p_3}}_{{e_1}{p_1}}}{\\bar{v}^" + L"{{e_1}{p_1}}_{{p_2}{p_3}}}{w^{}_{}[{e_1}]}}\\bigr) }"); } // simplified example with "diagonal" operator from the paper, inspired by @@ -1340,12 +1346,15 @@ TEST_CASE("wick", "[algorithms][wick][valgrind_skip]") { REQUIRE(result.as().size() == 5); // clang-format off + // h is aux-only (no bra/ket particle columns) -> column-Nonsymm, like + // the programmatic h in `input`; pin it so the parsed reference matches + // the fixed-default programmatic ctor (overlaps/deltas stay Symm) auto expected = deserialize( - "- h{;;p_3} ã{p_1,p_3;p_2,p_3}" - "+ h{;;p_3} δ{p_1;a_1} δ{a_2;p_2} ã{p_3;p_3} s{a_1;a_2} " - "- h{;;a_1} δ{a_2;p_2} ã{p_1;a_1} s{a_1;a_2} " - "- h{;;a_2} δ{p_1;a_1} s{a_1;a_2} ã{a_2;p_2} " - "+ h{;;a_3} δ{p_1;a_1} δ{a_2;p_2} s{a_1;a_3} s{a_3;a_2} "); + "- h{;;p_3}:N-N-N ã{p_1,p_3;p_2,p_3}" + "+ h{;;p_3}:N-N-N δ{p_1;a_1} δ{a_2;p_2} ã{p_3;p_3} s{a_1;a_2} " + "- h{;;a_1}:N-N-N δ{a_2;p_2} ã{p_1;a_1} s{a_1;a_2} " + "- h{;;a_2}:N-N-N δ{p_1;a_1} s{a_1;a_2} ã{a_2;p_2} " + "+ h{;;a_3}:N-N-N δ{p_1;a_1} δ{a_2;p_2} s{a_1;a_3} s{a_3;a_2} "); // clang-format on simplify(expected);