[anaddb] 77_ddb/sortph.F90 fails to connect branches

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] 77_ddb/sortph.F90 fails to connect branches

Post by t-nissie » Sat Jul 10, 2010 5:00 am

Hi all,

77_ddb/sortph.F90 sometimes fails to connect branches correctly.
Here are my quick-hacked 77_ddb/sortph.F90 and a
small patch for 72_response/phfrq3.F90.
Understand what they are, then use them.
But, you do NOT have to understand GROUP THEORY !!!

This technique may be also used to plot electronic band structures.

Ciao, ciao,
Takeshi Nishimatsu
http://loto.sourceforge.net/feram/
Fast MD program for perovskite-type ferroelectrics

Code: Select all

!{\src2tex{textfont=tt}}
!!****f* ABINIT/sortph
!! NAME
!! sortph
!!
!! FUNCTION
!! Sort the energies in order to have fine phonon
!! dispersion curves
!! It is best not to include the gamma point in the list
!!
!! COPYRIGHT
!! Copyright (C) 2002-2010 ABINIT group (FDortu,MVeithen)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! MODIFIED
!! Takeshi Nishimatsu
!!
!! INPUTS
!!  displ(2,3*natom,3*natom)= contain
!!   the displacements of atoms in cartesian coordinates.
!!   The first index means either the real or the imaginary part,
!!   The second index runs on the direction and the atoms displaced
!!   The third index runs on the modes.
!!  filnam=name of output files
!!   hacmm1,hartev,harthz,xkb= different conversion factors
!!  natom= number of atom
!!  phfrq(3*natom)= phonon frequencies in Hartree
!!  qphon(3)=phonon wavevector
!!  udispl=unit number for output of phonon eigendisplacements
!!  ufreq=unit number for output of phonon frequencies
!!
!! OUTPUT
!!  (only writing ?)
!!
!! NOTES
!! Called by one processor only
!!
!! PARENTS
!!      mkphbs,thm9
!!
!! CHILDREN
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine sortph(displ,filnam, natom,phfrq,qphon,udispl,ufreq)

use defs_basis
use defs_datatypes
use defs_abitypes

implicit none

!Arguments -----------------------------------
!scalars
integer,intent(in) :: natom,udispl,ufreq
character(len=fnlen),intent(in) :: filnam
!arrays
complex(dp),intent(in) :: displ(3*natom,3*natom)
real(dp),intent(in)    :: phfrq(3*natom),qphon(3)

!Local variables-------------------------------
!scalars
integer :: iatom,imode,j,i(1)
character(len=fnlen) :: file_displ,file_freq
!arrays
logical     ::               mask(3*natom)
real(dp)    ::           phfrqNew(3*natom)
complex(dp) ::           displNew(3*natom,3*natom)
complex(dp) ::    transpose_displ(3*natom,3*natom)
real(dp)    ::     abs_similarity(3*natom,3*natom)  !|<displNew|displLast>|
complex(dp),allocatable,save :: displLast(:,:)


! *********************************************************************

write(6, '(a)' )' sortph : enter '

if(.not.allocated(displLast)) then
   file_freq  = trim(filnam)//".freq" !---------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_freq)
   open(ufreq,FILE=trim(file_freq),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_freq),' opened '
   file_displ = trim(filnam)//".displ" !--------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_displ)
   open(udispl,FILE=trim(file_displ),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_displ),' opened '
   allocate(displLast(3*natom,3*natom)) !---------------------------------------
   phfrqNew(:) = phfrq(:)
   displNew(:,:) = displ(:,:)
else
   !Avoid gfortran 4.2.1 bug, with which you CANNOT conjg(transpose(displ))
   transpose_displ = transpose(displ)
   abs_similarity = abs(matmul(conjg(transpose_displ),displLast))
   mask(:) = .true.
   do j = 1, 3*natom
      i(:) = maxloc( abs_similarity(:,j), mask(:) )
      mask(i(1)) = .false.
      phfrqNew(j) = phfrq(i(1))
      displNew(:,j) = displ(:,i(1))
   end do
endif


!  Write frequencies in a file
!  CE FORMAT DOIT ETRE REECRIT CAR CAS PARTICULIER 3*natom=15
write(ufreq,'(15(5d14.6))')  (phfrqNew(j),j=1,3*natom)

!  write displacements in a file
do imode=1,3*natom
   do iatom=1,natom
      write(udispl,'(d14.6)') &
           sqrt(       displNew(3*(iatom-1)+1,imode)*   &
           &     conjg(displNew(3*(iatom-1)+1,imode)) + &
           &           displNew(3*(iatom-1)+2,imode)*   &
           &     conjg(displNew(3*(iatom-1)+2,imode)) + &
           &           displNew(3*(iatom-1)+3,imode)*   &
           &     conjg(displNew(3*(iatom-1)+3,imode)) )
   end do
end do

displLast(:,:) = displNew(:,:)
write(6, '(a)' )' sortph : exit '

end subroutine sortph
!!***


Code: Select all

--- 72_response/phfrq3.F90~     2010-05-16 05:54:52.000000000 +0900
+++ 72_response/phfrq3.F90      2010-07-10 11:07:29.421562936 +0900
@@ -331,20 +331,20 @@
      end do
    end do
  end do
-
-!Get the phonon displacements
- do imode=1,3*natom
-   do idir1=1,3
-     do ipert1=1,natom
-       i1=idir1+(ipert1-1)*3
-       index=i1+3*natom*(imode-1)
-       displ(2*index-1)=eigvec(2*index-1)&
-&       /  sqrt(amu(typat(ipert1))*amu_emass)
-       displ(2*index  )=eigvec(2*index  )&
-&       /  sqrt(amu(typat(ipert1))*amu_emass)
-     end do
-   end do
- end do
+displ=eigvec
+! !Get the phonon displacements
+!  do imode=1,3*natom
+!    do idir1=1,3
+!      do ipert1=1,natom
+!        i1=idir1+(ipert1-1)*3
+!        index=i1+3*natom*(imode-1)
+!        displ(2*index-1)=eigvec(2*index-1)&
+! &       /  sqrt(amu(typat(ipert1))*amu_emass)
+!        displ(2*index  )=eigvec(2*index  )&
+! &       /  sqrt(amu(typat(ipert1))*amu_emass)
+!      end do
+!    end do
+!  end do
 
 end subroutine phfrq3
 !!***

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

Re: [anaddb] 77_ddb/sortph.F90 fails to connect branches

Post by t-nissie » Sat Jul 10, 2010 5:42 am

Simplified:

Code: Select all

!{\src2tex{textfont=tt}}
!!****f* ABINIT/sortph
!! NAME
!! sortph
!!
!! FUNCTION
!! Sort the energies in order to have fine phonon
!! dispersion curves
!! It is best not to include the gamma point in the list
!!
!! COPYRIGHT
!! Copyright (C) 2002-2010 ABINIT group (FDortu,MVeithen)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! MODIFIED
!! Takeshi Nishimatsu
!!
!! INPUTS
!!  displ(2,3*natom,3*natom)= contain
!!   the displacements of atoms in cartesian coordinates.
!!   The first index means either the real or the imaginary part,
!!   The second index runs on the direction and the atoms displaced
!!   The third index runs on the modes.
!!  filnam=name of output files
!!   hacmm1,hartev,harthz,xkb= different conversion factors
!!  natom= number of atom
!!  phfrq(3*natom)= phonon frequencies in Hartree
!!  qphon(3)=phonon wavevector
!!  udispl=unit number for output of phonon eigendisplacements
!!  ufreq=unit number for output of phonon frequencies
!!
!! OUTPUT
!!  (only writing ?)
!!
!! NOTES
!! Called by one processor only
!!
!! PARENTS
!!      mkphbs,thm9
!!
!! CHILDREN
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine sortph(displ,filnam, natom,phfrq,qphon,udispl,ufreq)

use defs_basis
use defs_datatypes
use defs_abitypes

implicit none

!Arguments -----------------------------------
!scalars
integer,intent(in) :: natom,udispl,ufreq
character(len=fnlen),intent(in) :: filnam
!arrays
complex(dp),intent(in) :: displ(3*natom,3*natom)
real(dp),intent(in)    :: phfrq(3*natom),qphon(3)

!Local variables-------------------------------
!scalars
integer :: iatom,imode,j,i(1)
character(len=fnlen) :: file_displ,file_freq
!arrays
logical     ::               mask(3*natom)
real(dp)    ::         phfrqLocal(3*natom)
complex(dp) ::    transpose_displ(3*natom,3*natom)
real(dp)    ::     abs_similarity(3*natom,3*natom)  ! |<displ|displSaved>|
complex(dp),allocatable,save :: displSaved(:,:)


! *********************************************************************

write(6, '(a)' )' sortph : enter '

if(allocated(displSaved)) then
   transpose_displ = transpose(displ)   !Avoid gfortran 4.2.1 bug, with which you CANNOT conjg(transpose(displ))
   abs_similarity = abs(matmul(conjg(transpose_displ),displSaved))   ! |<displ|displSaved>|
   mask(:) = .true.
   do j = 1, 3*natom
      i(:) = maxloc( abs_similarity(:,j), mask(:) )
      mask(i(1)) = .false.
      phfrqLocal(j) = phfrq(i(1))
      displSaved(:,j) = displ(:,i(1))
   end do
else
   allocate(displSaved(3*natom,3*natom)) !-------------------------------------------------
   phfrqLocal(:) = phfrq(:)
   displSaved(:,:) = displ(:,:)
   file_freq  = trim(filnam)//".freq" !---------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_freq)
   open(ufreq,FILE=trim(file_freq),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_freq),' opened '
   file_displ = trim(filnam)//".displ" !--------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_displ)
   open(udispl,FILE=trim(file_displ),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_displ),' opened '
endif

!  Write frequencies in a file
!  CE FORMAT DOIT ETRE REECRIT CAR CAS PARTICULIER 3*natom=15
write(ufreq,'(15(5d14.6))')  (phfrqLocal(j),j=1,3*natom)

!  write displacements in a file
do imode=1,3*natom
   do iatom=1,natom
      write(udispl,'(d14.6)') &
           sqrt(       displSaved(3*(iatom-1)+1,imode)*   &
           &     conjg(displSaved(3*(iatom-1)+1,imode)) + &
           &           displSaved(3*(iatom-1)+2,imode)*   &
           &     conjg(displSaved(3*(iatom-1)+2,imode)) + &
           &           displSaved(3*(iatom-1)+3,imode)*   &
           &     conjg(displSaved(3*(iatom-1)+3,imode)) )
   end do
end do

write(6, '(a)' )' sortph : exit '
end subroutine sortph
!!***
Takeshi Nishimatsu
love && peace && free_software
http://loto.sourceforge.net/feram/
Fast MD program for perovskite-type ferroelectrics

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

Re: [anaddb] 77_ddb/sortph.F90 fails to connect branches

Post by t-nissie » Sat Aug 07, 2010 12:31 pm

Here is the final(?) patch for files in 6.2.1:

* Keep 72_response/phfrq3.F90 as in 6.2.1.
* 77_ddb/sortph.F90: file attached.
* 77_ddb/interfaces_77_ddb.F90: patch attached.
* 77_ddb/mkphbs.F90: patch attached.
* 77_ddb/thm9.F90: patch attached.
* I put http://loto.sourceforge.net/loto/PbTiO3 ... ddb.tar.gz , an example for this patch.
Find and view PbTiO3-new-anaddb/PbTiO3-phonon-band.pdf in the tar ball.
MD5 (PbTiO3-new-anaddb.tar.gz) = d829556780f7e943ea7750691c6f8530

Note that, even with this patch, sortph.F90 fails to connect
branches from the Gamma point when the system has LO-TO
splitting. In that case, you should move data column in your
B2EPS.freq file by your hands.

77_ddb/sortph.F90:

Code: Select all

!{\src2tex{textfont=tt}}
!!****f* ABINIT/sortph
!! NAME
!! sortph
!!
!! FUNCTION
!! Sort the energies in order to have fine phonon
!! dispersion curves
!! It is best not to include the gamma point in the list
!!
!! COPYRIGHT
!! Copyright (C) 2002-2010 ABINIT group (FDortu,MVeithen)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! MODIFIED
!! Takeshi Nishimatsu
!!
!! INPUTS
!!  eigvec(2*3*natom*3*natom)= contain
!!  the eigenvectors of the dynamical matrix.
!!  displ(2,3*natom,3*natom)= contain
!!   the displacements of atoms in cartesian coordinates.
!!   The first index means either the real or the imaginary part,
!!   The second index runs on the direction and the atoms displaced
!!   The third index runs on the modes.
!!  filnam=name of output files
!!   hacmm1,hartev,harthz,xkb= different conversion factors
!!  natom= number of atom
!!  phfrq(3*natom)= phonon frequencies in Hartree
!!  qphon(3)=phonon wavevector
!!  udispl=unit number for output of phonon eigendisplacements
!!  ufreq=unit number for output of phonon frequencies
!!
!! OUTPUT
!!  (only writing ?)
!!
!! NOTES
!! Called by one processor only
!!
!! PARENTS
!!      mkphbs,thm9
!!
!! CHILDREN
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine sortph(eigvec,displ,filnam, natom,phfrq,qphon,udispl,ufreq)

use defs_basis
use defs_datatypes
use defs_abitypes

implicit none

!Arguments -----------------------------------
!scalars
integer,intent(in) :: natom,udispl,ufreq
character(len=fnlen),intent(in) :: filnam
!arrays
complex(dp),intent(in) :: eigvec(3*natom,3*natom)
complex(dp),intent(in) :: displ(3*natom,3*natom)
real(dp),intent(in)    :: phfrq(3*natom),qphon(3)

!Local variables-------------------------------
!scalars
integer :: iatom,imode,j,i(1)
character(len=fnlen) :: file_displ,file_freq
!arrays
logical     ::               mask(3*natom)
real(dp)    ::           phfrqNew(3*natom)
complex(dp) ::           displNew(3*natom,3*natom)
complex(dp) ::          eigvecNew(3*natom,3*natom)
complex(dp) ::   transpose_eigvec(3*natom,3*natom)
real(dp)    ::     abs_similarity(3*natom,3*natom)  !|<displNew|displLast>|
complex(dp),allocatable,save :: eigvecLast(:,:)


! *********************************************************************

write(6, '(a)' )' sortph : enter '

if(.not.allocated(eigvecLast)) then
   file_freq  = trim(filnam)//".freq" !---------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_freq)
   open(ufreq,FILE=trim(file_freq),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_freq),' opened '
   file_displ = trim(filnam)//".displ" !--------------------------------------------------
   write(6, '(a,a)' )' sortph : opening file ',trim(file_displ)
   open(udispl,FILE=trim(file_displ),STATUS='replace',ACCESS='sequential',ACTION='write')
   write(6, '(a,a,a)' )' sortph : file ',trim(file_displ),' opened '
   allocate(eigvecLast(3*natom,3*natom)) !---------------------------------------
    phfrqNew(:)   =  phfrq(:)
    displNew(:,:) =  displ(:,:)
   eigvecNew(:,:) = eigvec(:,:)
else
   !Avoid gfortran 4.2.1 bug, with which you CANNOT conjg(transpose(displ))
   transpose_eigvec = transpose(eigvec)
   abs_similarity = abs(matmul(conjg(transpose_eigvec),eigvecLast))
   mask(:) = .true.
   do j = 1, 3*natom
      i(:) = maxloc( abs_similarity(:,j), mask(:) )
      mask(i(1)) = .false.
      phfrqNew(j)   =    phfrq(i(1))
      displNew(:,j) =  displ(:,i(1))
     eigvecNew(:,j) = eigvec(:,i(1))
   end do
endif


!  Write frequencies in a file
!  CE FORMAT DOIT ETRE REECRIT CAR CAS PARTICULIER 3*natom=15
write(ufreq,'(15(5d14.6))')  (phfrqNew(j),j=1,3*natom)

!  write displacements in a file
do imode=1,3*natom
   do iatom=1,natom
      write(udispl,'(d14.6)') &
           sqrt(       displNew(3*(iatom-1)+1,imode)*   &
           &     conjg(displNew(3*(iatom-1)+1,imode)) + &
           &           displNew(3*(iatom-1)+2,imode)*   &
           &     conjg(displNew(3*(iatom-1)+2,imode)) + &
           &           displNew(3*(iatom-1)+3,imode)*   &
           &     conjg(displNew(3*(iatom-1)+3,imode)) )
   end do
end do

eigvecLast(:,:) = eigvecNew(:,:)
write(6, '(a)' )' sortph : exit '

end subroutine sortph
!!***


Code: Select all

--- interfaces_77_ddb.F90~      2010-08-07 18:14:44.000000000 +0900
+++ interfaces_77_ddb.F90       2010-08-07 18:20:11.000000000 +0900
@@ -1855,7 +1855,7 @@
 end interface
 
 interface
- subroutine sortph(displ,filnam,& 
+ subroutine sortph(eigvec,displ,filnam,& 
   &  natom,phfrq,qphon,udispl,ufreq)
   use defs_basis
   implicit none
@@ -1863,6 +1863,7 @@
   integer,intent(in) :: udispl
   integer,intent(in) :: ufreq
   character(len=fnlen),intent(in) :: filnam
+  real(dp),intent(in) :: eigvec(2,3*natom,3*natom)
   real(dp),intent(in) :: displ(2,3*natom,3*natom)
   real(dp),intent(in) :: phfrq(3*natom)
   real(dp),intent(in) :: qphon(3)


Code: Select all

--- mkphbs.F90~ 2010-08-07 18:03:05.000000000 +0900
+++ mkphbs.F90  2010-08-07 18:05:42.000000000 +0900
@@ -286,7 +286,7 @@
 !  (visualization of phonon band structures)
    if (anaddb_dtset%eivec == 4) then
      tmpfilename = trim(outfile_radix)//"_B2EPS"
-     call sortph(displ,tmpfilename,&
+     call sortph(eigvec,displ,tmpfilename,&
 &     natom,phfrq,qphon,udispl,ufreq)
    end if
 


Code: Select all

--- thm9.F90~   2010-08-07 18:13:17.000000000 +0900
+++ thm9.F90    2010-08-07 18:16:04.000000000 +0900
@@ -308,7 +308,7 @@
 
      if (present(themflag)) then
        if (themflag ==2)then
-         call sortph(displ, trim(outfilename_radix), natom,phfrq,qphon,udispl,ufreq)
+         call sortph(eigvec, displ, trim(outfilename_radix), natom,phfrq,qphon,udispl,ufreq)
        end if
      end if
 
Takeshi Nishimatsu
love && peace && free_software
http://loto.sourceforge.net/feram/
Fast MD program for perovskite-type ferroelectrics

mverstra
Posts: 655
Joined: Wed Aug 19, 2009 12:01 pm

Re: [anaddb] 77_ddb/sortph.F90 fails to connect branches

Post by mverstra » Sat Sep 04, 2010 11:28 am

Hi Takeshi - this is great. I have converted your sortph to the abinit coding style and will try it out.

thanks a lot for your help

Matthieu
Matthieu Verstraete
University of Liege, Belgium

harada57
Posts: 4
Joined: Wed Aug 09, 2017 6:24 am

Re: [anaddb] 77_ddb/sortph.F90 fails to connect branches

Post by harada57 » Wed Aug 09, 2017 10:29 am

I am new to abinit and trying to find the lifetimes using GWA. Did you find answer to your question about where the abinit outputs the imaginary part of self energy ?


วิธีเล่นบาคาร่า

Locked