I don't think that's currently possible. The trace goes GeometryOptimizer.kernel -> geometric.optimize.run_optimizer -> CASCI_GradScanner.__call__ -> CASCI_Scanner.__call__ and the second step is, as far as PySCF is concerned, a black box: geomeTRIC is third-party software. I don't think you can pass kwargs all the way down. This would have to be implemented as a behavior change of CASCI_Scanner and/or CASCI_GradScanner, and would probably require adding new attributes to the class(es). Or you could just monkeypatch it, I guess (c.f.):

myopt = mycas.Gradients().optimizer(solver='geomeTRIC')
def mycall (self, mol_or_geom, mo_coeff=None, ci0=None):
    if isinstance(mol_or_geom, gto.MoleBase):
        mol = mol_or_geom
    else:
        mol = self.mol.set_geom_(mol_or_geom, inplace=False)

    self.reset (mol)

    if mo_coeff is None:
        mf_scanner = self._scf
        mf_scanner(mol)
        mo_coeff = mf_scanner.mo_coeff
    mo_coeff = mycas.sort_mo (cas_orb, mo_coeff=mo_coeff)
    if ci0 is None:
        ci0 = self.ci
    self.mol = mol
    e_tot = self.kernel(mo_coeff, ci0)[0]
    return e_tot
myopt.method.base.__call__ = mycall
mol_eq = myopt.kernel ()

Zooming out a bit, I wonder what your scientific goal is here, because I think this has the same fundamental problem as state-specific excited-state MC-SCF geometry optimization: orbitals will sometimes change order as the geometry changes. (1) Since you don't yet know what the final geometry is, how do you know which are the orbitals you want? Alternatively, one could compare the HF orbitals at step n to the stored CASCI mo_coeff (i.e., from step n-1), and select the ncas orbitals of the former with the highest overlap of the latter. That's probably more robust, but it's still not perfectly robust, since nothing can guarantee that the orbitals at the two geometries are qualitatively comparable so there might be a bunch of overlaps like 0.499999, 0.5000001, etc (and we aren't allowed to make linear combinations since this is CASCI! (2)). This all also increases the danger that the CPHF step of the CASCI gradient evaluation is unsolvable (since there are orbital-energy-difference denominators in those expressions), and you would have to be very careful never to sort any HF-occupied orbital above any HF-virtual orbital because then the expressions in grad.casci.grad_elec become incorrect.

I'm sure I'm overthinking the issue as it pertains to your specific case. I just wanted to offer some insight as to why I think PySCF hasn't implemented this behavior yet. Developers would need to know what general types of behaviors are most useful to the largest number of users.

1. Of course, this problem can also affect CASCI @ HF frontier orbitals, but it's probably less severe if all the active orbitals are clumped together around the frontier.
2. Some people do things like "CASCI @ MP2 natural orbitals" or "CASCI @ Boys localized orbitals". But anything of this type would require a totally different gradient implementation.

I don't think that's currently possible. The trace goes GeometryOptimizer.kernel -> geometric.optimize.run_optimizer -> CASCI_GradScanner.__call__ -> CASCI_Scanner.__call__ and the second step is, as far as PySCF is concerned, a black box: geomeTRIC is third-party software. I don't think you can pass kwargs all the way down. This would have to be implemented as a behavior change of CASCI_Scanner and/or CASCI_GradScanner, and would probably require adding new attributes to the class(es). Or you could just monkeypatch it, I guess (c.f.):

myopt = mycas.Gradients().optimizer(solver='geomeTRIC')
def mycall (self, mol_or_geom, mo_coeff=None, ci0=None):
    if isinstance(mol_or_geom, gto.MoleBase):
        mol = mol_or_geom
    else:
        mol = self.mol.set_geom_(mol_or_geom, inplace=False)

    self.reset (mol)

    if mo_coeff is None:
        mf_scanner = self._scf
        mf_scanner(mol)
        mo_coeff = mf_scanner.mo_coeff
    mo_coeff = mycas.sort_mo (cas_orb, mo_coeff=mo_coeff)
    if ci0 is None:
        ci0 = self.ci
    self.mol = mol
    e_tot = self.kernel(mo_coeff, ci0)[0]
    return e_tot
myopt.method.base.__call__ = mycall
mol_eq = myopt.kernel ()

Zooming out a bit, I wonder what your scientific goal is here, because I think this has the same fundamental problem as state-specific excited-state MC-SCF geometry optimization: orbitals will sometimes change order as the geometry changes. (1) Since you don't yet know what the final geometry is, how do you know which are the orbitals you want? Alternatively, one could compare the HF orbitals at step n to the stored CASCI mo_coeff (i.e., from step n-1), and select the ncas orbitals of the former with the highest overlap of the latter. That's probably more robust, but it's still not perfectly robust, since nothing can guarantee that the orbitals at the two geometries are qualitatively comparable so there might be a bunch of overlaps like 0.499999, 0.5000001, etc (and we aren't allowed to make linear combinations since this is CASCI! (2)). This all also increases the danger that the CPHF step of the CASCI gradient evaluation is unsolvable (since there are orbital-energy-difference denominators in those expressions), and you would have to be very careful never to sort any HF-occupied orbital above any HF-virtual orbital because then the expressions in grad.casci.grad_elec become incorrect.

I'm sure I'm overthinking the issue as it pertains to your specific case. I just wanted to offer some insight as to why I think PySCF hasn't implemented this behavior yet. Developers would need to know what general types of behaviors are most useful to the largest number of users.

1. Of course, this problem can also affect CASCI @ HF frontier orbitals, but it's probably less severe if all the active orbitals are clumped together around the frontier.
2. Some people do things like "CASCI @ MP2 natural orbitals" or "CASCI @ Boys localized orbitals". But anything of this type would require a totally different gradient implementation.

Thank you for your serious reply. I think you've thought about this issue very comprehensively and given me some inspiration. Perhaps the above method can only be used if the initial geometrie is close to the desired geometrie and the orbitals order is the same. I'll reconsider it based on what I'm trying to solve in conjunction with your input. Thank you very much for your reply!