Please make sure to send follow-up questions to the list.
In general, you can’t model a potential with extra terms without changing the code.
In this specific case, however, since all you’re adding is a constant term, it’s not clear why it’s even necessary. It doesn’t affect either the energy or the dynamics. Couldn’t you just add in the extra –N^2 term during postprocessing?