The ABINIT Users Mailing List ( CLOSED )

[NISHIMATSU Takeshi <t_nissie@yahoo.co.jp>]

alexandra,

> Also, could you please tell me if the dynamical matrix in the DDB file  
> uses the default database of isotope masses, or is given for unitary  
> masses.

I guess the matrix in DDB file is the *reduced* force constant matrix.
It is not weighted by 1/sqrt(M_iM_j).
Atomic mass is given by src/11util/atmdata.F90.

With attached code, you can read (and diagonalize) the force constant matrix.

===================================================================
! diagonalize15x15.F -*-f90-*-
! Time-stamp: <2008-12-24 13:02:29 takeshi>
! Author: Takeshi NISHIMATSU
! Usage: tail -225 perovskite_o_DS5_DDB | SOMEWHERE/diagonalize15x15
!!
program diagonalize15x15
  implicit none
  integer, parameter :: N_ATOM = 5
  integer, parameter :: N_FREEDOM = N_ATOM*3
  integer, parameter :: LWORK = N_FREEDOM*20   ! for LAPACK zheev
  integer I, J, alpha, beta, io, info
  complex*16 Force_Matrix_Hermetian(N_FREEDOM,N_FREEDOM)
  real*8     Force_Matrix_real_symmetric(N_FREEDOM,N_FREEDOM)
  ! The force matrix (or the second derivative of total energy) becomes
  ! real symmetric, if its k-point has a (very) good symmetry.
  ! The force matrix recorded in the end of each DDB file of ABINIT has a
  ! unit of Hartree/(lattice constant)^2. NOT Hartree/Bohr^2.
  complex*16 work(LWORK)                       ! for LAPACK zheev
  real*8    rwork(LWORK)                       ! for LAPACK zheev
  real*8    eigenvalues(N_FREEDOM)
  real*8    tmp_r, tmp_i
  character(len=50) :: fmt_eigen
  logical   all_real

  all_real = .true.
  do
     read(5, *, IOSTAT=io)    alpha, I, beta, J, tmp_r, tmp_i
     !write(6,'(4i4,2f10.4)') alpha, I, beta, J, tmp_r, tmp_i
     if (io.ne.0) exit
     if (I<=N_ATOM .and. J<=N_ATOM) then
        write(6,'(4i4,2i5,2f15.9)') alpha, I, beta, J, 3*(I-1)+alpha, 
3*(J-1)+beta , tmp_r, tmp_i
        Force_Matrix_real_symmetric(3*(I-1)+alpha, 3*(J-1)+beta) = tmp_r
        Force_Matrix_Hermetian(     3*(I-1)+alpha, 3*(J-1)+beta) = 
cmplx(tmp_r,tmp_i)
        if ( abs(tmp_i) > 1.0d-4 ) all_real = .false.
     end if
  end do

  if (all_real) then
     !LAPACK dsyev(JOBZ, UPLO, N,         A,                        LDA,      
 W,           WORK, 
LWORK, INFO)
     call    dsyev('V',  'L',  N_FREEDOM, Force_Matrix_real_symmetric, 
N_FREEDOM, eigenvalues, rwork,
LWORK, info)
     write(fmt_eigen,'(a,i3,a)') '(i2,f11.6,', N_ATOM, '(f10.3,f7.3,f7.3))'
     do j=1,N_FREEDOM
        write(6,fmt_eigen) j, eigenvalues(j), 
(Force_Matrix_real_symmetric(i,j), i=1,N_FREEDOM)
     end do
  else
     !LAPACK zheev(JOBZ,UPLO,N,        A,                  LXL,      W,       
  
WORK,LWORK,RWORK,INFO)
     call    zheev('V', 'L',
N_FREEDOM,Force_Matrix_Hermetian,N_FREEDOM,eigenvalues,work,LWORK,rwork,info)
     write(fmt_eigen,'(a,i3,a)') "(i2,f11.6,", N_ATOM, 
"((f8.2,f5.2)(f6.2,f5.2)(f6.2,f5.2)))"
     do j=1,N_FREEDOM
        write(6,fmt_eigen) j, eigenvalues(j), (Force_Matrix_Hermetian(i,j), 
i=1,N_FREEDOM)
     end do
  end if
end program diagonalize15x15
===================================================================

Good luck,
-- 
 Takeshi Nishimatsu
 love && peace && free_software
 http://loto.sourceforge.net/feram/   Fast MD program for perovskite-type 
ferroelectrics