Abstract A variational Monte Carlo calculation is performed to investigate the origin of large isotope effects in ferroelectric or antiferroelectric phase transitions of hydrogen bonded materials. Non-adiabatic coupling between an electron and a nucleus is especially considered in the variational wavefunction. A characteristic feature of the hydrogen bond is that, in a certain range of oxygen-oxygen distance, two stable states for a hydrogen (deuteron) nucleus, i.e., a localized state and a tunneling state, is energetically in subtle balance with each other. It is demonstrated that the balance is strongly affected by isotope substitution. This can easily lead to the large isotope effect. However, the difference of the charge state between H and D reported by an x-ray experiment cannot be explained even by including the non-adiabatic effect.