Rotations

From time to time the wavefunctions need to be rotated by a band $\times$ band matrix $M$, e.g.

$$\psi^\prime_{pb} = \sum_k \psi_{pk}M_{kb}$$ (4.3)

Clearly this transformation, in general, mixes data between all of the nodes. Rather than storing an entire wavefunction we transform just the parts needed for one node at a time.

! construct the data for node `node'
do node=1,bnd_nodes
  transform local data to give contribution to `node'
  reduce over bnd-group
end do

This is simplest if we encapsulate all such operations into a single subroutine, overloaded for the various wavefunction types. In fact such a routine already exists, wave_rotate.

There are also occasions when we wish to transform a wavefunction by a matrix that is not square, so that the resulting number of bands is not the same as the initial number of bands. This occurs when, for example, we orthogonalise a slice to an existing wavefunction; we copy the wavefunction to a temporary slice, transform the data by multiplying by the overlap matrix, which changes the number of bands to the same as the slice, and then subtract it from slice. This is not a `rotation', and goes against the specification for wave_rotate. For this reason we coded a new, more general subroutine wave_transform which allows any linear transformation of the bands.

The orthogonalisation subroutines can now be coded compactly as:

! orthogonalise slice to all bands of wvfn
call wave_dot_all(wvfn,slice,overlap)

call wave_copy(wvfn,tmp_slice)
call wave_transform(tmp_slice,overlap)

! Subtract overlapping components
call wave_add(tmp_slice,slice,c1=-cmplx_1)