def METHOD_NAME():
    #! A demonstration of mixed Cartesian/ZMatrix geometry specification, using variables, for
    #! the benzene-hydronium complex.  Atoms can be placed using ZMatrix coordinates, whether they belong
    #! to the same fragment or not.  Note that the Cartesian specification must come before the ZMatrix entries
    #! because the former define absolute positions, while the latter are relative.
    refENuc = 268.6171792624
    refGEOM = \
          [[   0.710500000000,    -0.794637665924,    -1.230622098778],
           [   1.421000000000,    -0.794637665924,     0.000000000000],
           [   0.710500000000,    -0.794637665924,     1.230622098778],
           [  -0.710500000000,    -0.794637665924,     1.230622098778],
           [   1.254500000000,    -0.794637665924,    -2.172857738095],
           [  -1.254500000000,    -0.794637665924,     2.172857738095],
           [  -0.710500000000,    -0.794637665924,    -1.230622098778],
           [  -1.421000000000,    -0.794637665924,     0.000000000000],
           [   2.509000000000,    -0.794637665924,     0.000000000000],
           [   1.254500000000,    -0.794637665924,     2.172857738095],
           [  -1.254500000000,    -0.794637665924,    -2.172857738095],
           [  -2.509000000000,    -0.794637665924,     0.000000000000],
           [   0.000000000000,     3.205362334076,     0.000000000000],
           [   0.494974746831,     3.555362334076,    -0.857321409974],
           [   0.494974746831,     3.555362334076,     0.857321409974],
           [  -0.989949493661,     3.555362334076,     0.000000000000]]
    dimer = qcdb.Molecule("""
    1 1
    # This part is just a normal Cartesian geometry specification for benzene
    C          0.710500000000    -0.794637665924    -1.230622098778
    C          1.421000000000    -0.794637665924     0.000000000000
    C          0.710500000000    -0.794637665924     1.230622098778
    C         -0.710500000000    -0.794637665924     1.230622098778
    H          1.254500000000    -0.794637665924    -2.172857738095
    H         -1.254500000000    -0.794637665924     2.172857738095
    C         -0.710500000000    -0.794637665924    -1.230622098778
    C         -1.421000000000    -0.794637665924     0.000000000000
    H          2.509000000000    -0.794637665924     0.000000000000
    H          1.254500000000    -0.794637665924     2.172857738095
    H         -1.254500000000    -0.794637665924    -2.172857738095
    H         -2.509000000000    -0.794637665924     0.000000000000
    # And the hydronium part is specified using a zmatrix, referencing the benzene coordinates
    X  1  CC  3  30   2  A2
    O  13 R   1  90   2  90
    H  14 OH  13 TDA  1  0
    H  14 OH  15 TDA  13 A1
    H  14 OH  15 TDA  13 -A1
    CC    = 1.421
    CH    = 1.088
    A1    = 120.0
    A2    = 180.0
    OH    = 1.05
    R     = 4.0
    units angstrom
    """)
    dimer.update_geometry()
    assert compare_values(refENuc, dimer.nuclear_repulsion_energy(), 9, "Bz-H3O+: nuclear repulsion energy")
    geom_now = qcdb.mscale(dimer.geometry(), qcdb.constants.bohr2angstroms)
    assert compare_matrices(refGEOM, geom_now, 6, "Bz-H3O+: geometry and orientation")