@@ -934,6 +934,96 @@ def contract_ss(self, fcivec, norb, nelec):
         with self.make_space(norb, nelec, compress_links=False) as gas:
             return _GasSpinPlan(gas).contract(fcivec)
 
+    def transform_ci_within_gas(self, fcivec, norb, nelec, u):
+        """Return CI coefficients for ``mo_new = mo_old @ u`` within each GAS.
+
+        ``u`` must be real, orthogonal and block diagonal in the GAS partition.
+        Adjacent Givens rotations act on local string axes of each legal GAS
+        block; no full-CAS embedding or dense determinant rotation is formed.
+        Reflections and orbital permutations are included. The input is not
+        modified. This explicit method is separate from the unrestricted
+        ``transform_ci_for_orbital_rotation`` hook probed by native Newton.
+        """
+
+        gas_orbs, _, _ = self._space_spec(norb, nelec)
+        u = _as_c_double(u, (int(norb), int(norb)))
+        if not numpy.all(numpy.isfinite(u)) or not numpy.allclose(
+                u.T @ u, numpy.eye(norb), atol=1e-10, rtol=0):
+            raise ValueError("GAS orbital rotation must be finite and orthogonal")
+        allowed = numpy.zeros(u.shape, dtype=bool)
+        factors = []
+        offset = 0
+        for size in gas_orbs:
+            sl = slice(offset, offset + size)
+            allowed[sl, sl] = True
+            work = u[sl, sl].copy()
+            rotations = []
+            # L_m ... L_1 U = D, so the passive CI transform applies the
+            # exterior representations of L_1, ..., L_m, then diagonal D.
+            # Adjacent orbital pairs have no intervening fermionic sign.
+            for col in range(size - 1):
+                for q in range(size - 1, col, -1):
+                    p = q - 1
+                    a, b = work[p, col], work[q, col]
+                    if b == 0:
+                        continue
+                    norm = numpy.hypot(a, b)
+                    c, s = a / norm, b / norm
+                    row_p, row_q = work[p].copy(), work[q].copy()
+                    work[p] = c * row_p + s * row_q
+                    work[q] = -s * row_p + c * row_q
+                    rotations.append((p, q, c, s))
+            factors.append((rotations, numpy.where(work.diagonal() < 0, -1., 1.)))
+            offset += size
+        if numpy.any(numpy.abs(u[~allowed]) > 1e-12):
+            raise ValueError("orbital rotation must not mix different GAS subspaces")
+
+        source = _as_c_double(fcivec)
+        with self.make_space(norb, nelec, compress_links=False) as gas:
+            if source.size != gas.ndet:
+                raise ValueError("CI vector size does not match GAS determinant count")
+            result = source.reshape(-1).copy()
+            ngas = len(gas_orbs)
+            local = {}
+            for block in gas.block_descriptors():
+                occupations = [int(gas._gas.sector_occ[sid * ngas + g])
+                               for sid in (block['sa'], block['sb'])
+                               for g in range(ngas)]
+                shape = tuple(cistring.num_strings(gas_orbs[g % ngas], n)
+                              for g, n in enumerate(occupations))
+                begin = block['offset']
+                size = block['na'] * block['nb']
+                tensor = result[begin:begin + size].reshape(shape)
+                # C GAS sector strides have the last subspace varying fastest:
+                # (alpha GAS axes..., beta GAS axes...).
+                for axis, n in enumerate(occupations):
+                    g = axis % ngas
+                    key = (g, n)
+                    if key not in local:
+                        strings = cistring.make_strings(range(gas_orbs[g]), n)
+                        rotations, phases = factors[g]
+                        pairs = {}
+                        for p, q, _, _ in rotations:
+                            if (p, q) not in pairs:
+                                left = numpy.flatnonzero(
+                                    ((strings >> p) & 1) & ~((strings >> q) & 1))
+                                right = numpy.searchsorted(
+                                    strings, strings[left] ^ ((1 << p) | (1 << q)))
+                                pairs[p, q] = (left, right)
+                        signs = numpy.ones(strings.size)
+                        for p in numpy.flatnonzero(phases < 0):
+                            signs *= 1 - 2 * ((strings >> p) & 1)
+                        local[key] = pairs, signs
+                    pairs, signs = local[key]
+                    view = numpy.moveaxis(tensor, axis, 0)
+                    for p, q, c, s in factors[g][0]:
+                        left, right = pairs[p, q]
+                        old_left, old_right = view[left].copy(), view[right].copy()
+                        view[left] = c * old_left + s * old_right
+                        view[right] = -s * old_left + c * old_right
+                    view *= signs.reshape((-1,) + (1,) * (view.ndim - 1))
+        return result.reshape(source.shape)
+
     def transform_ci_for_orbital_rotation(self, fcivec, norb, nelec, u):
         raise NotImplementedError(
             "arbitrary active-orbital rotations do not preserve a GAS space")
 
@@ -27,10 +27,11 @@
 ``gas_restr`` and ``gas_restr_type`` are normalized by the same GAS helper
 routines.  This module supports ordinary state averaging, including
 zero-weight roots, GAS-safe canonicalization of inactive/external orbitals,
+optional within-GAS pseudo-natural orbitals with synchronized CI rotation,
 density fitting through PySCF's ``mcscf.df`` machinery, energy-only scanners
 and the GASCI-native spin-penalty Hamiltonian.  It does not implement
 state-average-mix, state-specific excited-state wrappers,
-active-space natural-orbital rotations, analytic gradients/NACs or the
+unrestricted active-space natural-orbital rotations, analytic gradients/NACs or the
 legacy two-step CASSCF driver.
 """
 
@@ -756,6 +757,10 @@ class GASSCF(newton_casscf.CASSCF):
         returns occupations as one array per subspace. These methods leave
         the computational orbitals and CI vectors unchanged.
 
+        To rotate computational orbitals and CI together, explicitly call
+        ``mc.canonicalize_(gas_pseudo_natorb=True)``. The variant without the
+        trailing underscore returns the transformed results without writeback.
+
         Analysis orbitals can be exported with PySCF's Molden writer::
 
             from pyscf.tools import molden
@@ -1255,30 +1260,112 @@ def casci(self, mo_coeff=None, ci0=None, eris=None,
 
     def canonicalize(self, mo_coeff=None, ci=None, eris=None, sort=False,
                      gas_natorb=False, gasdm1=None, verbose=None,
-                     cas_natorb=None, **kwargs):
-        """Canonicalize core/external orbitals without rotating GAS subspaces.
-
-        Native CASSCF calls this method with its ``natorb`` flag as the fifth
-        positional argument.  In GASSCF, any automatic active-space natural
-        orbital rotation would generally mix GAS subspaces and invalidate the
-        restricted determinant space, so active-space natural orbitals are
-        explicitly guarded.  The active GAS block is otherwise kept fixed;
-        only inactive and external orbitals are canonicalized by the inherited
-        PySCF machinery.
+                     cas_natorb=None, *, gas_pseudo_natorb=False, **kwargs):
+        """Return ``(mo_coeff, ci, mo_energy)`` without writing to this object.
+
+        By default only core/external orbitals are canonicalized. Opt in with
+        ``gas_pseudo_natorb=True`` to diagonalize the density separately within
+        each GAS subspace and transform every CI root into the new basis,
+        including zero-weight roots. SA uses the weighted density unless
+        ``gasdm1`` is supplied in the input active-orbital basis.
+
+        Active occupations are ordered from largest to smallest within each
+        unfrozen symmetry block of a GAS subspace. Frozen orbitals and orbital
+        symmetry/``extrasym`` labels are respected. ``sort`` retains its native
+        meaning for core/external orbital energies. ``mo_energy`` contains
+        Fock diagonal elements, not pseudo-natural occupations.
+
+        These are pseudo-natural orbitals: density between GAS subspaces need
+        not vanish. Rotations across GAS subspaces remain unsupported, as do
+        ``gas_natorb=True`` and ``cas_natorb=True``. Use ``canonicalize_`` to
+        write the returned orbitals, CI and orbital energies to this object.
         """
 
         if gas_natorb or cas_natorb:
             _unsupported("GAS natural-orbital rotation")
-        return gasci.GASCI.canonicalize(
-            self, mo_coeff, ci, eris, sort=sort, gas_natorb=False,
+        if not gas_pseudo_natorb:
+            return gasci.GASCI.canonicalize(
+                self, mo_coeff, ci, eris, sort=sort, gas_natorb=False,
+                gasdm1=gasdm1, verbose=verbose, **kwargs)
+
+        mo_coeff = self.mo_coeff if mo_coeff is None else mo_coeff
+        ci = self.ci if ci is None else ci
+        if mo_coeff is None or ci is None:
+            raise ValueError("pseudo-natural canonicalization requires orbitals and CI")
+        if gasdm1 is None:
+            gasdm1 = self.make_gasdm1(
+                ci=ci, state=None if self._has_state_weights() else 0)
+        gasdm1 = numpy.asarray(gasdm1)
+        if (gasdm1.shape != (self.ncas, self.ncas)
+                or numpy.iscomplexobj(gasdm1)
+                or not numpy.all(numpy.isfinite(gasdm1))):
+            raise ValueError("gasdm1 must be a finite real (ncas, ncas) matrix")
+
+        nmo = mo_coeff.shape[1]
+        frozen = numpy.zeros(nmo, dtype=bool)
+        if isinstance(self.frozen, (int, numpy.integer)):
+            frozen[:self.frozen] = True
+        elif self.frozen is not None:
+            frozen[self.frozen] = True
+        orbsym = getattr(mo_coeff, 'orbsym', numpy.zeros(nmo, dtype=int))
+        extrasym = getattr(self, 'extrasym', None)
+        if extrasym is None:
+            extrasym = numpy.zeros(nmo, dtype=int)
+        rotation = numpy.eye(self.ncas)
+        offset = 0
+        for size in self.gas_orbs:
+            groups = {}
+            for p in range(offset, offset + size):
+                i = self.ncore + p
+                if not frozen[i]:
+                    groups.setdefault((orbsym[i], extrasym[i]), []).append(p)
+            for indices in groups.values():
+                block = numpy.ix_(indices, indices)
+                _, vectors = self._natural_eigensystem(gasdm1[block], sort=True)
+                rotation[block] = vectors
+            offset += size
+
+        # Transform CI explicitly, without exposing a general rotation hook to
+        # Newton. In particular, never solve GASCI again in the rotated basis.
+        def transform(vector):
+            return self.fcisolver.transform_ci_within_gas(
+                vector, self.ncas, self.nelecas, rotation)
+
+        if isinstance(ci, (list, tuple)):
+            ci_new = type(ci)(transform(vector) for vector in ci)
+        elif numpy.ndim(ci) == 2 and self.fcisolver.nroots > 1:
+            ci_new = numpy.asarray([transform(vector) for vector in ci])
+        else:
+            ci_new = transform(ci)
+        # The native routine may reorder orbital labels when sort=True;
+        # own the metadata as well as the MO array to keep this call nonmutating.
+        mo_input = mo_coeff.copy()
+        if getattr(mo_coeff, 'orbsym', None) is not None:
+            mo_input = lib.tag_array(mo_input, orbsym=numpy.array(orbsym, copy=True))
+        # Both calls use the original basis, so a supplied ERIS remains valid.
+        fock = self.get_fock(mo_coeff, ci, eris, gasdm1, verbose)
+        mo_new, _, mo_energy = gasci.GASCI.canonicalize(
+            self, mo_input, ci, eris, sort=sort, gas_natorb=False,
             gasdm1=gasdm1, verbose=verbose, **kwargs)
+        active = slice(self.ncore, self.ncore + self.ncas)
+        mo_new[:, active] = mo_coeff[:, active] @ rotation
+        mo_energy[active] = numpy.einsum(
+            'pi,pi->i', mo_new[:, active], fock @ mo_new[:, active])
+        return mo_new, ci_new, mo_energy
 
     def canonicalize_(self, mo_coeff=None, ci=None, eris=None, sort=False,
                       gas_natorb=False, gasdm1=None, verbose=None,
-                      cas_natorb=None, **kwargs):
+                      cas_natorb=None, *, gas_pseudo_natorb=False, **kwargs):
+        """Write the canonicalized MO, CI and orbital energies to this object.
+
+        Returns the same three values as :meth:`canonicalize`; the method
+        without the trailing underscore leaves these result attributes intact.
+        """
+
         mo_coeff, ci, mo_energy = self.canonicalize(
             mo_coeff, ci, eris, sort=sort, gas_natorb=gas_natorb,
-            gasdm1=gasdm1, verbose=verbose, cas_natorb=cas_natorb, **kwargs)
+            gasdm1=gasdm1, verbose=verbose, cas_natorb=cas_natorb,
+            gas_pseudo_natorb=gas_pseudo_natorb, **kwargs)
         self.mo_coeff = mo_coeff
         self.ci = ci
         self.mo_energy = mo_energy
 
@@ -275,6 +275,72 @@ def test_gas_fci_vector_converters(self):
             numpy.testing.assert_array_equal(projected_twice, projected)
             self.assertEqual(numpy.count_nonzero(projected), gas.ndet)
 
+    def test_ci_rotation_within_gas_matches_full_fci(self):
+        rng = numpy.random.default_rng(211)
+        norb, nelec = 5, (3, 2)
+        gas_orbs = (3, 2)
+        bounds = ((2, 4), (5, 5))
+        _, blocks = addons_gas.normalize_gas_spec(
+            gas_orbs, nelec, bounds, 'cumulative-occ')
+        supergroups = numpy.unique(blocks[:, :2] + blocks[:, 2:], axis=0)
+        cases = [
+            (gas_orbs, nelec, bounds, 'cumulative-occ'),
+            (gas_orbs, nelec, supergroups, 'supergroup'),
+            (gas_orbs, nelec, blocks[::2], 'spin-supergroup'),
+            ((5,), nelec, ((3, 2),), 'spin-supergroup'),
+            (gas_orbs, (0, 0), ((0, 0, 0, 0),), 'spin-supergroup'),
+            (gas_orbs, (5, 5), ((3, 2, 3, 2),), 'spin-supergroup'),
+        ]
+        for sizes, electrons, restriction, kind in cases:
+            solver = fci_gas.FCISolver(
+                gas_orbs=sizes, gas_restr=restriction, gas_restr_type=kind)
+            u = numpy.zeros((norb, norb))
+            offset = 0
+            for size in sizes:
+                q = numpy.linalg.qr(rng.normal(size=(size, size)))[0]
+                # Explicitly include determinant -1, not just proper rotations.
+                if numpy.linalg.det(q) > 0:
+                    q[:, 0] *= -1
+                u[offset:offset+size, offset:offset+size] = q
+                offset += size
+            permutation = numpy.eye(norb)
+            permutation[:, [0, 1]] = permutation[:, [1, 0]]
+            with solver.make_space(norb, electrons) as gas:
+                ci = rng.normal(size=gas.ndet)
+                ci /= numpy.linalg.norm(ci)
+                before = ci.copy()
+                full = fci_gas.gas2fci(ci, gas)
+                for rotation in (u, permutation):
+                    with self.subTest(kind=kind, nelec=electrons, sizes=sizes):
+                        expected = fci_addons.transform_ci_for_orbital_rotation(
+                            full, norb, electrons, rotation)
+                        with mock.patch.object(fci_gas, 'gas2fci',
+                                               side_effect=AssertionError('no CAS embedding')):
+                            actual = solver.transform_ci_within_gas(
+                                ci, norb, electrons, rotation)
+                        numpy.testing.assert_allclose(
+                            fci_gas.gas2fci(actual, gas), expected, atol=2e-12, rtol=0)
+                        recovered = solver.transform_ci_within_gas(
+                            actual, norb, electrons, rotation.T)
+                        numpy.testing.assert_allclose(recovered, ci, atol=2e-12, rtol=0)
+                        self.assertAlmostEqual(numpy.linalg.norm(actual), 1., places=12)
+                        numpy.testing.assert_array_equal(ci, before)
+
+    def test_ci_rotation_within_gas_rejects_invalid_transforms(self):
+        solver = fci_gas.FCISolver(gas_orbs=(2, 2),
+                                  gas_restr=((1, 1, 1, 1),))
+        ci = numpy.ones(16) / 4
+        with self.assertRaisesRegex(ValueError, 'orthogonal'):
+            solver.transform_ci_within_gas(ci, 4, (2, 2), numpy.eye(4) * 2)
+        mixed = numpy.eye(4)
+        mixed[:, [1, 2]] = mixed[:, [2, 1]]
+        with self.assertRaisesRegex(ValueError, 'mix different GAS'):
+            solver.transform_ci_within_gas(ci, 4, (2, 2), mixed)
+        with self.assertRaisesRegex(TypeError, 'real-valued'):
+            solver.transform_ci_within_gas(ci, 4, (2, 2), numpy.eye(4, dtype=complex))
+        with self.assertRaisesRegex(ValueError, 'CI vector size'):
+            solver.transform_ci_within_gas(ci[:-1], 4, (2, 2), numpy.eye(4))
+
     def test_restricted_hamiltonian_projection(self):
         cases = (
             {