Computing Green's function for a 2D slice of a 3D crystal

Crystal: FCC lattice with harmonic nearest neighbour interactions. The domain is assumed to be cubic with periodic boundaries.

FCC Basis: 4 atom Cubic FCC basis. The lattice vectors are simply aligned along the cubic axes.

Calculation Scheme: First the 12x12 Hessian matrices are written down as a function of CUBIC lattice separation. Note: since we have only nearest neighbour interactions, we have at most 27 lattice translations whose corresponding hessian matrices have non-zero entries.

The Dynamical matrices, D(k), are then computing for all CUBIC wavevectors that can fit in the simulation box via taking the DFT of Hessian matrices, that is: D(k) = SUM{{R} H(R)*exp(i k.R)}}, where (R} is the set of all possible translations in the lattice (NOTE: R and -R must be separately counted). This DFT was not done through the FFT since the number of non-zero Hessians is much smaller than the number of lattice translations.

Each dynamical matrix was diagonalized to get the 12 polarization vectors & their corresponding eigenvalues {(P_i, L_i) : 1<= i <=12} and the Green's function in fourier space is computed using: G(k) = SUM{(i=1...12) P_i.outerProduct.P_i/L_i }

Finally the real space Green's function can be obtained for the whole crystal through an inverse Fourier Transform: G(R) = SUM{{k} G(k)*exp(i k.R)} where {k} is the set of all fourier modes that can fit in the lattice.

To obtain the Green's function for the slice, I rotated the tensor obtained and discarded all components that had a 'z' piece. Finally I was left with a 2x2 matrix which gave the 2D correlations in a slice (for a particular separation vector) from a 3D crystal.


Difficulties: Due to the fact, that G(k) is the inverse of D(k), we have problems when computing G(k) for long wavelength fourier modes (since their translational and longitudinal polarizations; eigenvalues are very small, hence these D(k) have high condition number). Specifically we can have G(k) & G(-k) differ significantly enough for G(R) to be not purely real. I solved this problem by explicitly setting: G(k)=G(-k).