Bond_style Gaussian question in simulation (n=1,compared with harmonic)

If you check out this post, the documentation for bond_style gaussian is missing a factor of 2 in the exponent (i.e., it should be \exp(-2(r-r_i)^2/w_i)). Though, I’m not sure if that accounts for your huge difference.