[gmx-users] Phonon calculations in periodic crystals

[Anne Kelley]

I am trying to use GROMACS to calculate the phonons (normal modes) of a bulk crystal, CdSe.  I have found a simple force field, Coulomb + Lennard-Jones, in the literature (Rabani, J. Chem. Phys. 116, 258, 2002) which the author showed reproduced the phonon dispersion curves and other mechanical properties of bulk wurtzite CdSe quite well.  A number of other workers have used this force field in molecular dynamics simulations.  But when I use Rabani's force field with GROMACS I get phonon frequencies that are much too high, up to about 2.2 times the experimental ones.

I am doing all of my calculations with the double precision version of GROMACS.  I have made a .top file for CdSe using Rabani's Lennard-Jones parameters and ionic charges, and a .gro file containing an integer number of unit cells with the known lattice constants.  I first do an energy minimization until the maximum forces are around 1.e-4, and get the right crystal structure and lattice constants.  I am using periodic boundary conditions with PME.  I then use the "nm" integrator (with the -t option to read in the more precise .trr structure file) to calculate the Hessian, and then the g_nmeig_d program to diagonalize the Hessian and get the normal modes.  This all seems to work fine, but I don't get the literature values for the frequencies (calculated maximum about 450 cm-1, literature and experimental about 215 cm-1).  I have checked that when I enter the correct masses and known harmonic force constant for the H2 molecule, I get back the right vibrational frequency.  I have tried changing the size of the system (5, 7, or 9 unit cells in each direction) and it has almost no effect on the frequencies.  I have tried things like changing the Coulomb and Lennard-Jones cutoffs, and even tried regular Ewald rather than PME (which took a very long time), but these had no significant effect on my results.  I also tried calculating the phonon spectrum for a different material, AgBr, using a Coulomb + Buckingham potential from the literature (J. Phys. Chem. 99, 14344, 1995).  This gave me a better result, but still the distribution of frequencies is not correct and the maximum phonon frequency is about 15% higher than what the authors got with the same force field.

Are you aware of any issues with GROMACS in doing normal mode calculations on periodic systems?  Can you suggest any likely things I'm doing wrong?

Anne Kelley

Anne Myers Kelley
Professor of Chemistry, School of Natural Sciences
Secretary-Treasurer, APS Division of Laser Science
University of California, Merced
5200 North Lake Road, Merced, CA 95343
Tel. 209-228-4345
amkelley at ucmerced.edu
http://faculty.ucmerced.edu/amkelley/