I manage to find how position and ADP constraints are handled  
by constraint objects featuring methods to return the independent  
variables and the gradient with respect to them. Then cctbx/cctbx/
examples/structure_factor_derivatives_4.py shows how those features  
can be used. Would Hydrogen geometry constraint be best handled by  
mimicking those special positions constraint?

I am not sure. The handling of hydrogens seems more complicated than
dealing with the constraints due to symmetry. I just had a quick look
at the SHELX manual. I see many special cases.

Aplenty indeed. Right. Let me precise my question then. The two classes dealing with position and ADP symmetry constraints feature the same methods:
- n_independent_params
- independent_params
- all_params
- independent_gradients
- independent_curvatures
- gradient_sum_matrix

I was trying to figure out whether this interface would be sufficient and/or necessary for the Hydrogen geometry constraints. I am not quite sure I understand what is gradient_sum_matrix exactly in fact. 

If I had a chance to work on this, I'd start with a Python script
emulating some of the cases handled by SHELX, refining a few small
structures. After getting a few concrete examples to work, I'd try to
generalize the procedures for handling the special cases. I'd switch to
C++ only after having fully conquered the math and the bookkeeping in
Python, and being sure about the final interfaces.

It is the philosophy of the cctbx, isn't it? You have constructed the cctbx so as to use it as Python library, relegating C++ to some sort of assembly language for efficiency. The lack of abstraction in the C++ code of the cctbx (hardly any inheritance, no advanced genericity) would require wrapping a lot of the cctbx classes behind proxies inheriting from abstract types. A typical example for me are those classes dealing with position constraints and ADP constraints. They are unrelated in C++, cctbx::sgtbx::site_constraints and cctbx::sgtbx::tensor_rank_2::constraints, although they have both have the same member functions listed above. Of course, from Python, it does not matter, thanks to duck typing: if two objects answer the same method calls, then they are by all practical means of the same type.

To implement the core refinement, I'd copy and modify
cctbx/cctbx/xray/minimization.py. See also
phenix/phenix/substructure/hyss/minimization.py, which you can find
in the CCI Apps bundle (http://phenix-online.org/download/cci_apps/).

I have already started to study the former. Thanks for pointing out the latter. 

Hope this helps. As you can tell, I am not an expert in constrained
hydrogen refinement, but I'm very interested in seeing such algorithms
in the cctbx (we'd probably want to use them in phenix.refine). Let me
know if you feel I can do something to help out.

Our group would be more than happy to contribute them to the cctbx since we absolutely need them for our project. 
Cheerio,

Luc