A general purpose direct BEM code has been developed for three-dimensional crack problems in piezoelectric structures. Special 3D non-continuous crack tip elements and several techniques for determining crack tip parameters were implemented. To calculate the electromechanical energy release rate for a virtual crack extension, the θ-method is employed, which was originally suggested by BONNET for linear elastic materials. The paper presents the generalization and numerical realization of the θ-method to 3D piezoelectric cracks. The great advantage of the θ-method is the direct computation of energy release rate, whereas the way via K-factors and the IRWIN matrix is more complicated. The efficiency and accuracy of the technique are shown for various example problems by comparing with analytical solutions.