def METHOD_NAME():
    atoms = load_atoms("structures/carbon_chain.xyz")
    morse1 = MorsePotentialPerAtom(r0=3, epsilon=2, rho0=6)
    morse2 = MorsePotentialPerAtom(r0=3, epsilon=4, rho0=6)
    mcfm_pot = create_mcfm_potential(atoms,
                                     classical_calculator=morse1,
                                     qm_calculator=morse2,
                                     special_atoms_list=None,
                                     double_bonded_atoms_list=None)
    dynamics = get_dynamics(atoms, mcfm_pot, T=100)
    # ------ Add constraints
    fixed = 29
    moving = 1
    direction = np.array([0, 0, 1])
    velocity = 5e-4 * units.Ang
    c_fixed = ase.constraints.FixAtoms([fixed])
    c_relax = ase.constraints.FixAtoms([moving])
    c_moving = LinearConstraint(moving, direction, velocity)
    atoms.set_constraint([c_fixed, c_moving])
    # ------ Adjust flagging energies
    mcfm_pot.qm_cluster.flagging_module.ema_parameter = 0.003
    mcfm_pot.qm_cluster.flagging_module.qm_flag_potential_energies =\
        np.ones((len(atoms), 2), dtype=float)
    mcfm_pot.qm_cluster.flagging_module.qm_flag_potential_energies[:, 0] *= -5
    mcfm_pot.qm_cluster.flagging_module.qm_flag_potential_energies[:, 1] *= -6
    atoms.arrays["qm_flag_potential_energies[in_out]"] =\
        mcfm_pot.qm_cluster.flagging_module.qm_flag_potential_energies
    # ------ Run dynamics
    dynamics.run(14000)
    atoms.set_constraint([c_fixed, c_relax])
    dynamics.run(3000)