The ABINIT Users Mailing List ( CLOSED )

["alexandra.carvalho@epfl.ch" <alexandra.carvalho@epfl.ch>]

Dear Takeshi,

Thank you very much for answering my question and
also for providing the code which clarified the units
used in the DDB file.

Kind regards,

Alexandra

Quoting 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
>
>
>
>