[anaddb] FYI: plot eigenvalues of IFC matrix
Posted: 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.
If you want to use GNUPLOT, Python, Ruby, etc,
this patch for abinit-6.8.1/src/77_ddb/sortph.F90 is useful.
Enjoy,
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,