[anaddb] FYI: plot eigenvalues of IFC matrix

MULTIBINIT, aTDEP, ANADDB, optics, cut3d, mrg* postprocessors

Moderators: MMNSchmitt, gonze

Locked
t-nissie
Posts: 12
Joined: Sat Jul 10, 2010 4:21 am
Contact:

[anaddb] FYI: plot eigenvalues of IFC matrix

Post by t-nissie » Fri Aug 19, 2011 6:09 pm

Hi,

Do you want to plot a dispersion relation of eigenvalues of
inter-atomic force constant (IFC) matrix? It may be more
useful for discussing distortions of crystals than phonon dispersion.
For examples, Fig. 3 in [K. Leung, E. Cockayne, and A. F. Wright, Phys. Rev. B 65,
214111 (2002) http://link.aps.org/doi/10.1103/PhysRevB.65.214111 ] and
Fig. 5 (B) in http://loto.sourceforge.net/feram/parameters/ .

Try this patch for abinit-6.8.1/src/72_response/phfrq3.F90,
then _B2EPS.freq file will contain the eigenvalues, not phonon frequencies.

Code: Select all

--- phfrq3.F90~ 2011-06-26 09:32:46.000000000 +0900
+++ phfrq3.F90  2011-07-06 17:56:52.621968140 +0900
@@ -229,7 +229,7 @@
 !Include the masses in the dynamical matrix
  do ipert1=1,natom
    do ipert2=1,natom
-     fac=1.0_dp/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
+     fac=1.0_dp   !/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
      do idir1=1,3
        do idir2=1,3
          i1=idir1+(ipert1-1)*3
@@ -239,9 +239,9 @@
          displ(2*index  )=displ(2*index  )*fac*nearidentity(idir1,idir2)
 !        This is to break slightly the translation invariance, and make
 !        the automatic tests more portable
-         if(ipert1==ipert2 .and. idir1==idir2)then
-           displ(2*index-1)=displ(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp
-         end if
+         !if(ipert1==ipert2 .and. idir1==idir2)then
+         !  displ(2*index-1)=displ(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp
+         !end if
        end do
      end do
    end do
@@ -288,15 +288,16 @@
 
 !Get the phonon frequencies (negative by convention, if
 !the eigenvalue of the dynamical matrix is negative)
- do imode=1,3*natom
-   if(eigval(imode)>=1.0d-16)then
-     phfrq(imode)=sqrt(eigval(imode))
-   else if(eigval(imode)>=-1.0d-16)then
-     phfrq(imode)=zero
-   else
-     phfrq(imode)=-sqrt(-eigval(imode))
-   end if
- end do
+!  do imode=1,3*natom
+!    if(eigval(imode)>=1.0d-16)then
+!      phfrq(imode)=sqrt(eigval(imode))
+!    else if(eigval(imode)>=-1.0d-16)then
+!      phfrq(imode)=zero
+!    else
+!      phfrq(imode)=-sqrt(-eigval(imode))
+!    end if
+!  end do
+phfrq(:)=eigval(:)
 
 !Fix the phase of the eigenvectors
  allocate(dum(2,0))


If you want to use GNUPLOT, Python, Ruby, etc,
this patch for abinit-6.8.1/src/77_ddb/sortph.F90 is useful.

Code: Select all

--- sortph.F90~ 2011-06-26 10:13:13.000000000 +0900
+++ sortph.F90  2011-07-06 10:56:28.000000000 +0900
@@ -134,13 +134,13 @@
 
 
 !Write frequencies in a file
- write(fmt_phfrq,'(a,i3,a)') '(', 3*natom, 'd14.6)'
+ write(fmt_phfrq,'(a,i3,a)') '(', 3*natom, 'e14.6)'
  write(ufreq,fmt_phfrq) (phfrqNew(j),j=1,3*natom)
 
 !write displacements in a file
  do imode=1,3*natom
    do iatom=1,natom
-     write(udispl,'(d14.6)') &
+     write(udispl,'(e14.6)') &
      sqrt(       displNew(3*(iatom-1)+1,imode)*   &
 &     conjg(displNew(3*(iatom-1)+1,imode)) + &
 &     displNew(3*(iatom-1)+2,imode)*   &


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

Locked