Qi = ij*QdXi;  <span class="comment">% d(ij*Q)/dXi</span>
        dbQj = ij*QdXj + bsxfun(@times,-2*ij(:,d),Q2(d,:)); <span class="comment">% d(ij*Q)/dXj</span>
        daQi = ii*QdXi + bsxfun(@times,-2*ii(:,d),Q2(d,:)); <span class="comment">% d(ii*Q)/dXi</span>
        daQj = ii*QdXj; <span class="comment">% d(ii*Q)/dXj</span>

        dsaQi = sum(daQi.*ii,2) - 2.*aQ(:,d).*ii(:,d);
        dsaQj = sum(daQj.*ii,2);
        dsbQi = sum(dbQi.*ij,2);
        dsbQj = sum(dbQj.*ij,2) - 2.*bQ(:,d).*ij(:,d);
    