!{\src2tex{textfont=tt}}
!!****f* ABINIT/symdm9
!!
!! NAME
!! symdm9
!!
!! FUNCTION
!! Use the set of special k points calculated by the Monkhorst &
!! Pack Technique.
!! Check if all the informations for the k points are present in
!! the DDB to determine their dynamical matrices.
!! Generate the dynamical matrices of the set of k points which
!! samples homogeneously the entire Brillouin zone.
!!
!! COPYRIGHT
!! Copyright (C) 1999-2012 ABINIT group (JCC,XG)
!! This file is distributed under the terms of the
!! GNU General Public Licence, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! INPUTS
!! blkflg(nsize,nblok)= flag of existence for each element of the DDB
!! blknrm(1,nblok)=norm of qpt providing normalization
!! blkqpt(1<ii<9,nblok)=q vector of a phonon mode (ii=1,2,3)
!! blktyp(nblok)=1 or 2 depending on non-stationary or stationary block
!!  3 for third order derivatives
!! blkval(2,3*mpert*3*mpert,nblok)= all the dynamical matrices
!! gprim(3,3)=dimensionless primitive translations in reciprocal space
!! indsym = mapping of atoms under symops
!! mpert =maximum number of ipert
!! natom=number of atoms in unit cell
!! nblok=number of blocks in the DDB
!! nqpt=number of special q points
!! nsym=number of space group symmetries
!! rfmeth = 1 if non-stationary block
!!  2 if stationary block
!!  3 if third order derivatives
!! rprim(3,3)=dimensionless primitive translations in real space
!! spqpt(3,nqpt)=set of special q points generated by the Monkhorst
!!  & Pack Method
!! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
!! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
!!
!! OUTPUT
!! dynmat(2,3,natom,3,natom,nqpt)=dynamical matrices relative to the q
!!  points of the B.Z. sampling
!!
!! TODO
!! A full description of the inputs should be included
!!
!! PARENTS
!!      mkifc9
!!
!! CHILDREN
!!      leave_new,wrtout
!!
!! SOURCE

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

#include "abi_common.h"


subroutine symdm9(blkflg,blknrm,blkqpt,blktyp,blkval,&
& dynmat,gprim,indsym,mpert,natom,nblok,nqpt,nsym,rfmeth,&
& rprim,spqpt,symrec,symrel)

 use m_profiling

 use defs_basis

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'symdm9'
 use interfaces_14_hidewrite
 use interfaces_16_hideleave
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: mpert,natom,nblok,nqpt,nsym,rfmeth
!arrays
 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),blktyp(nblok)
 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
 real(dp),intent(in) :: blknrm(3,nblok),blkqpt(9,nblok)
 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert,nblok),gprim(3,3),rprim(3,3)
 real(dp),intent(in) :: spqpt(3,nqpt)
 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt)

!Local variables -------------------------
!tol sets tolerance for equality of q points between those of
!the DDB and those of the sampling grid
!scalars
 integer :: ia,ib,iblok,idir1,idir2,ii,ipert1,ipert2,iqpt,isym,jj,kk,ll
 integer :: mu,nu,q1,q2
 real(dp),parameter :: tol=2.d-8
 real(dp) :: arg1,arg2,im,re,sumi,sumr
 character(len=500) :: message
!arrays
 integer,allocatable :: qtest(:,:)
 real(dp) :: qq(3),qsym(6),ss(3,3)
 real(dp),allocatable :: ddd(:,:,:,:,:)

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

 ABI_ALLOCATE(ddd,(2,3,natom,3,natom))
!Check if the blkqpt points and their symmetrics are sufficient
!in the DDB to retrieve all the q points of the B.Z. sampling

!Initialization of a test variable
 ABI_ALLOCATE(qtest,(nqpt,3))
 do iqpt=1,nqpt
   qtest(iqpt,1)=0
 end do

!Q points coming from the DDB
!write(std_out,*)' Nbr. of Blocks -> ',nblok

 do iblok=1,nblok
   if (blktyp(iblok)==rfmeth) then
     qq(1)=blkqpt(1,iblok)/blknrm(1,iblok)
     qq(2)=blkqpt(2,iblok)/blknrm(1,iblok)
     qq(3)=blkqpt(3,iblok)/blknrm(1,iblok)
!    Calculation of the symmetric points (including Time Reversal)
     do isym=1,nsym
       qsym(1)=qq(1)*symrec(1,1,isym)+qq(2)*symrec(1,2,isym)+&
&       qq(3)*symrec(1,3,isym)
       qsym(2)=qq(1)*symrec(2,1,isym)+qq(2)*symrec(2,2,isym)+&
&       qq(3)*symrec(2,3,isym)
       qsym(3)=qq(1)*symrec(3,1,isym)+qq(2)*symrec(3,2,isym)+&
&       qq(3)*symrec(3,3,isym)
!      Dont forget the Time Reversal symmetry
       qsym(4)=-qq(1)*symrec(1,1,isym)-qq(2)*symrec(1,2,isym)&
&       -qq(3)*symrec(1,3,isym)
       qsym(5)=-qq(1)*symrec(2,1,isym)-qq(2)*symrec(2,2,isym)&
&       -qq(3)*symrec(2,3,isym)
       qsym(6)=-qq(1)*symrec(3,1,isym)-qq(2)*symrec(3,2,isym)&
&       -qq(3)*symrec(3,3,isym)
!      Comparison between the q points and their symmetric points
!      and the set of q points which samples the entire Brillouin zone
       do iqpt=1,nqpt
         if (mod(abs(spqpt(1,iqpt)-qsym(1))+tol,1._dp)<2*tol)then
           if (mod(abs(spqpt(2,iqpt)-qsym(2))+tol,1._dp)<2*tol)then
             if (mod(abs(spqpt(3,iqpt)-qsym(3))+tol,1._dp)<2*tol)then

!              write(std_out,*)' q point from the DDB ! '
!              write(std_out,*)' block -> ',iblok
!              write(std_out,*)' sym.  -> ',isym
!              write(std_out,*)' No Time Reversal '
!              write(std_out,*)'(',qsym(1),',',qsym(2),',',qsym(3),')'
!              write(std_out,*)' '

               qtest(iqpt,1)=iblok
               qtest(iqpt,2)=isym
               qtest(iqpt,3)=0
             end if
           end if
         end if
         if (mod(abs(spqpt(1,iqpt)-qsym(4))+tol,1._dp)<2*tol)then
           if (mod(abs(spqpt(2,iqpt)-qsym(5))+tol,1._dp)<2*tol)then
             if (mod(abs(spqpt(3,iqpt)-qsym(6))+tol,1._dp)<2*tol)then

!              write(std_out,*)' q point from the DDB ! '
!              write(std_out,*)' block -> ',iblok
!              write(std_out,*)' sym.  -> ',isym
!              write(std_out,*)' Time Reversal '
!              write(std_out,*)'(',qsym(4),',',qsym(5),',',qsym(6),')'
!              write(std_out,*)' '

               qtest(iqpt,1)=iblok
               qtest(iqpt,2)=isym
               qtest(iqpt,3)=1
             end if
           end if
         end if

!        End of the loop on the q points of the sampling
       end do

!      End of the loop on the symmetries
     end do

   end if

!  End of the loop on the q points of the DDB
 end do

!Check if all the informations relatives to the q points sampling
!are found in the DDB; if not => stop message
 do iqpt=1,nqpt
   if (qtest(iqpt,1)==0) then
     write(message, '(a,a,a)' )&
&     ' symdm9 : the bloks found in the DDB are characterized',ch10,&
&     '  by the following wavevectors :'
     call wrtout(std_out,message,'COLL')
     do iblok=1,nblok
       write(message, '(a,4d20.12)')&
&       ' ',blkqpt(1,iblok),blkqpt(2,iblok),&
&       blkqpt(3,iblok),blknrm(1,iblok)
       call wrtout(std_out,message,'COLL')
     end do
     write(message, '(a,a,a,a,a,i6,a,a,a,3es16.6,a,a,a,a)' )&
&     ' symdm9 : ERROR -',ch10,&
&     '  Informations are missing in the DDB.',ch10,&
&     '  The dynamical matrix number',iqpt,' cannot be built,',ch10,&
&     '  since no blok with wavevector',spqpt(1:3,iqpt),ch10,&
&     '  has been found.',ch10,&
&     '  Action : add the required blok in the DDB, or modify your input file.'
     call wrtout(std_out,message,'COLL')
     call leave_new('COLL')
   end if
 end do

!Generation of the dynamical matrices relative to the q points
!of the set which samples the entire Brillouin zone
 do iqpt=1,nqpt
   q1=qtest(iqpt,1)
   q2=qtest(iqpt,2)
!  Check if the symmetry accompagnied with time reversal : q <- -q
   if (qtest(iqpt,3)==0) then
     do ii=1,3
       qq(ii)=blkqpt(ii,q1)/blknrm(1,q1)
       do jj=1,3
         ss(ii,jj)=0.0_dp
         do kk=1,3
           do ll=1,3
             ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2)
           end do
         end do
       end do
     end do

!    DEBUG
!    write(std_out,*)ss(1,1)
!    write(std_out,*)ss(2,1)
!    write(std_out,*)ss(3,1)
!    write(std_out,*)ss(1,2)
!    write(std_out,*)ss(2,2)
!    write(std_out,*)ss(3,2)
!    write(std_out,*)ss(1,3)
!    write(std_out,*)ss(2,3)
!    write(std_out,*)ss(3,3)
!    ENDDEBUG

   else

     do ii=1,3
       qq(ii)=-blkqpt(ii,q1)/blknrm(1,q1)
       do jj=1,3
         ss(ii,jj)=0.0_dp
         do kk=1,3
           do ll=1,3
             ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2)
           end do
         end do
       end do
     end do

   end if
!  
!  Check whether all the information is contained
!  in the DDB
   do ipert2=1,natom
     do idir2=1,3
       do ipert1=1,natom
         do idir1=1,3
           if(blkflg(idir1,ipert1,idir2,ipert2,q1)/=1)then
             write(message, '(a,a,a,a,a,i6,a,a,a,4i4,a,a,a,a)' )&
&             ' symdm9 : ERROR -',ch10,&
&             '  Informations are missing in the DDB.',ch10,&
&             '  In blok',q1,' the following element is missing :',ch10,&
&             '  idir1,ipert1,idir2,ipert2=',idir1,ipert1,idir2,ipert2,ch10,&
&             '  Action : add the required information in the DDB,',ch10,&
&             '  or modify your input file.'
             call wrtout(std_out,message,'COLL')
             call leave_new('COLL')
           end if
         end do
       end do
     end do
   end do

!  Read the dynamical matrices in the DDB
   do ipert2=1,natom
     do idir2=1,3
       do ipert1=1,natom
         do idir1=1,3
           ddd(:,idir1,ipert1,idir2,ipert2)=&
&           blkval(:,idir1+3*(ipert1-1+mpert*(idir2-1+3*(ipert2-1))),q1)
         end do
       end do
     end do
   end do

!  Calculation of the dynamical matrix of a symmetrical q point
   do ia=1,natom
     do ib=1,natom
!      write(std_out,*)'atom-> ',ia,indsym(4,q2,ia)
!      write(std_out,*)'atom-> ',ib,indsym(4,q2,ib)
       arg1=two_pi*(qq(1)*indsym(1,q2,ia)+qq(2)*indsym(2,q2,ia)&
&       +qq(3)*indsym(3,q2,ia))
       arg2=two_pi*(qq(1)*indsym(1,q2,ib)+qq(2)*indsym(2,q2,ib)&
&       +qq(3)*indsym(3,q2,ib))
       re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
       im=cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)
!      write(std_out,*)'re : ',re
!      write(std_out,*)'im : ',im
       do mu=1,3
         do nu=1,3
!          write(std_out,*)' '
           sumr=0.0_dp
           sumi=0.0_dp
           do ii=1,3
             do jj=1,3
!              If there is Time Reversal : D.M. <- Complex Conjugate D.M.
               if (qtest(iqpt,3)==0) then
                 sumr=sumr+ss(mu,ii)*ss(nu,jj)*&
&                 ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
                 sumi=sumi+ss(mu,ii)*ss(nu,jj)*&
&                 ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
               else
                 sumr=sumr+ss(mu,ii)*ss(nu,jj)*&
&                 ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
                 sumi=sumi-ss(mu,ii)*ss(nu,jj)*&
&                 ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
               end if
             end do
           end do

!          Dynmat -> Dynamical Matrix for the q point of the sampling
!          write(std_out,*)' Sumr -> ',mu,nu,sumr
!          write(std_out,*)' Sumi -> ',mu,nu,sumi
           dynmat(1,mu,ia,nu,ib,iqpt)=re*sumr-im*sumi
           dynmat(2,mu,ia,nu,ib,iqpt)=re*sumi+im*sumr

!          DEBUG
!          if((ia==2 .or. ia==3) .and. ib==1)then
!          write(std_out,'(5i3,2es16.8)' )&
!          &       mu,ia,nu,ib,iqpt,dynmat(1:2,mu,ia,nu,ib,iqpt)
!          end if
!          ENDDEBUG

!          End loop on the coordinates
         end do
       end do

!      End loop on the ia atoms
     end do

!    End loop on the ib atoms
   end do

!  End loop on the q points of the sampling
 end do

 ABI_DEALLOCATE(ddd)
 ABI_DEALLOCATE(qtest)

end subroutine symdm9
!!***
