!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_dynmat
!! NAME
!!  m_dynmat
!!
!! FUNCTION
!!  This module provides low-level tools to operate on the dynamical matrix
!!
!! COPYRIGHT
!!  Copyright (C) 2014 ABINIT group (XG, JCC, MJV, NH, RC, MVeithen, MM)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! TODO
!!  Use more explicative names for the procedures!
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

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

#include "abi_common.h"

module m_dynmat
    
 use defs_basis
 use m_profiling
 use m_errors
 use m_linalg_interfaces

 use m_numeric_tools,   only : wrap2_pmhalf
 use m_ewald,           only : ewald9

 implicit none

 private

 public :: asria_calc           ! Calculate the correction for the Acoustic sum rule on 
                                !   the InterAtomic Forces or on the dynamical matrix directly
 public :: asria_corr           ! Imposition of the Acoustic sum rule on the InterAtomic Forces
                                !   or on the dynamical matrix directly from the previously calculated d2asr
 public :: asrprs               ! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry
 public :: cart29               ! Transform a second-derivative matrix from reduced  coordinates to cartesian coordinates, and also
                                !   1) add the ionic part of the effective charges,
                                !   2) normalize the electronic dielectric tensor, and add the vacuum polarisation
 public :: cart39               ! Transform a vector from reduced coordinates to cartesian coordinates,
                                !   taking into account the perturbation from which it was derived,
                                !   and also check the existence of the new values.
 public :: chkph3               ! Check the completeness of the dynamical matrix
 public :: chneu9               ! Imposition of the Acoustic sum rule on the Effective charges
 public :: d2sym3               ! Build (nearly) all the other matrix elements that can be build using symmetries.
 public :: q0dy3_apply          ! Takes care of the inclusion of the ewald q=0 term in the dynamical matrix
 public :: q0dy3_calc           ! Calculate the q=0 correction term to the dynamical matrix
 public :: symdyma              ! Symmetrize the dynamical matrices
 public :: wings3               ! Suppress the wings of the cartesian 2DTE for which the diagonal element is not known
 public :: asrif9               ! Imposes the Acoustic Sum Rule to Interatomic Forces
 public :: make_bigbox          ! Generates a Big Box of R points for the Fourier Transforms the dynamical matrix
 public :: canat9               ! From reduced to canonical coordinates
 public :: canct9               ! Convert from canonical coordinates to cartesian coordinates
 public :: chkrp9               ! Check if the rprim used for the definition of the unit cell (in the
                                ! inputs) are consistent with the rprim used in the routine generating  the Big Box 
 public :: dist9                ! Compute the distance between atoms in the big box
 public :: ftifc_q2r            ! Fourier transform of the dynamical matrices to obtain interatomic forces (real space).
 public :: ftifc_r2q            ! Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space).
 public :: ifclo9               ! Convert from cartesian coordinates to local coordinates
 public :: wght9                ! Generates a weight to each R points of the Big Box and for each pair of atoms
 public :: d3sym                ! Given a set of calculated elements of the 3DTE matrix,
                                ! build (nearly) all the other matrix elements that can be build using symmetries.
 public :: symdm9               ! Generate the dynamical matrix in the full BZ from the irreducible q-points.
 public :: axial9               ! Generates the local coordinates system from the  knowledge of the first vector (longitudinal) and
                                !   the ifc matrix in cartesian coordinates
 public :: dymfz9               ! Multiply the dynamical matrix by a phase shift to account for normalized canonical coordinates.
 public :: nanal9               ! Subtract/Add the non-analytical part from one dynamical matrix with number iqpt.
 public :: gtdyn9               ! Generates a dynamical matrix from interatomic force constants and long-range electrostatic interactions.


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

contains 
!!***

!!****f* m_dynmat/asria_calc
!! NAME
!! asria_calc
!!
!! FUNCTION
!! Calculate the correction for the Acoustic sum rule on the InterAtomic Forces
!! or on the dynamical matrix directly
!!
!! INPUTS
!! asr=(0 => no ASR, 1 or 2=> the diagonal element is modified to give the ASR,
!!      5 => impose hermitian solution using lapack call)
!! d2cart=matrix of second derivatives of total energy, in cartesian coordinates
!! mpert =maximum number of ipert
!! natom=number of atom
!!
!! OUTPUT
!! d2asr=matrix used to store the correction needed to fulfill
!! the acoustic sum rule.
!!
!! PARENTS
!!      anaddb,gath3
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine asria_calc(asr,d2asr,d2cart,mpert,natom)


!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 'asria_calc'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: asr,mpert,natom
!arrays
 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
 real(dp),intent(out) :: d2asr(2,3,natom,3,natom)

!Local variables-------------------------------
!scalars
 integer :: idir1,idir2,ii,ipert1,ipert2
 integer :: constrank, imatelem, iconst, nconst, nd2_packed, info
 character(len=500) :: message
!arrays
 integer, allocatable :: packingindex(:,:,:,:)
 real(dp), allocatable :: constraints(:,:,:)
 real(dp), allocatable :: d2cart_packed(:,:)
 real(dp), allocatable :: singvals(:)
 real(dp), allocatable :: constr_rhs(:,:)
 real(dp), allocatable :: work(:,:),rwork(:)

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

 d2asr = zero

 if (asr==0) return

 write(message, '(a)' )&
& ' asria_calc : calculation of the correction to the ASR for the interatomic forces.'
 call wrtout(std_out,message,'COLL')
 do ipert1=1,natom
   do idir1=1,3
     do idir2=1,3

!      Compute d2asr
       do ipert2=1,natom
         d2asr(:,idir1,ipert1,idir2,ipert1)=&
&         d2asr(:,idir1,ipert1,idir2,ipert1)+&
&         d2cart(:,idir1,ipert1,idir2,ipert2)
       end do
     end do
   end do
 end do

!holistic method: overwrite d2asr with hermitian solution
 if (asr == 5) then
   nconst = 9*natom
   nd2_packed = 3*natom*(3*natom+1)/2
   ABI_ALLOCATE(constraints,(2,nconst, nd2_packed))
   ABI_ALLOCATE(d2cart_packed,(2,nd2_packed))
   ABI_ALLOCATE(constr_rhs,(2,nd2_packed))
   ABI_ALLOCATE(singvals,(nconst))
   ABI_ALLOCATE(work,(2,3*nd2_packed))
   ABI_ALLOCATE(rwork,(5*nd2_packed))
   ABI_ALLOCATE(packingindex,(3,natom,3,natom))
   ii=1
   packingindex=-1
   do ipert2=1,natom
     do idir2=1,3
       do ipert1=1,ipert2-1
         do idir1=1,3
           packingindex(idir1,ipert1,idir2,ipert2) = ii
           ii = ii+1
         end do
       end do
       do idir1=1,idir2
         packingindex(idir1,ipert2,idir2,ipert2) = ii
         ii = ii+1
       end do
     end do
   end do
!  setup constraint matrix
   constraints = zero
   do ipert1=1,natom
     do idir1=1,3
       do idir2=1,3
         iconst = idir2+3*(idir1-1 + 3*(ipert1-1))
!        set all atom forces, this component
         do ipert2=1,natom
           imatelem = packingindex(idir1,ipert1,idir2,ipert2)
           if (imatelem == -1) then
             imatelem = packingindex(idir2,ipert2,idir1,ipert1)
           end if
           constraints(1,iconst,imatelem) = one
         end do
       end do
     end do
   end do

   d2cart_packed = -999.0d0
   do ipert2=1,natom
     do idir2=1,3
       do ipert1=1,natom
         do idir1=1,3
           imatelem = packingindex(idir1,ipert1,idir2,ipert2)
           if (imatelem == -1) cycle
           d2cart_packed(:,imatelem) = d2cart(:,idir1,ipert1,idir2,ipert2)
         end do
       end do
     end do
   end do
   constr_rhs = zero
   constr_rhs(1,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(1,:))
   constr_rhs(2,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(2,:))

!  lwork = 3*nd2_packed
   call zgelss (nconst,nd2_packed,1,constraints,nconst,constr_rhs,nd2_packed,&
&   singvals,-one,constrank,work,3*nd2_packed,rwork,info)

   write(message,"(a,i0)")'zgelss returned info: ', info
   ABI_CHECK(info /=0, message)

!  unpack 
   do ipert2=1,natom
     do idir2=1,3
       do ipert1=1,natom
         do idir1=1,3
           imatelem = packingindex(idir1,ipert1,idir2,ipert2)
           if (imatelem == -1) then
             imatelem = packingindex(idir2,ipert2,idir1,ipert1)
!            NOTE: should complex conjugate the correction below.
           end if
           d2asr(:,idir1,ipert1,idir2,ipert2) = constr_rhs(:,imatelem)
         end do
       end do
     end do
   end do

   ABI_DEALLOCATE(constraints)
   ABI_DEALLOCATE(d2cart_packed)
   ABI_DEALLOCATE(singvals)
   ABI_DEALLOCATE(constr_rhs)
   ABI_DEALLOCATE(work)
   ABI_DEALLOCATE(rwork)
   ABI_DEALLOCATE(packingindex)
 end if

end subroutine asria_calc
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/asria_corr
!! NAME
!! asria_corr
!!
!! FUNCTION
!! Imposition of the Acoustic sum rule on the InterAtomic Forces
!!  or on the dynamical matrix directly from the previously calculated d2asr
!!
!! INPUTS
!! asr=(0 => no ASR, 1 or 2=> the diagonal element is modified to give the ASR,
!!      5 => impose hermitian solution using lapack call)
!! d2asr=matrix used to store the correction needed to fulfill
!! the acoustic sum rule.
!! mpert =maximum number of ipert
!! natom=number of atom
!!
!! OUTPUT
!! Input/Output:
!! d2cart=matrix of second derivatives of total energy, in cartesian coordinates
!!
!! SIDE EFFECTS
!!
!! PARENTS
!!      anaddb,elast9,gath3,instr9,mkphbs,thmeig
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine asria_corr(asr,d2asr,d2cart,mpert,natom)


!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 'asria_corr'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: asr,mpert,natom
!arrays
 real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)

!Local variables-------------------------------
!scalars
 integer :: idir1,idir2,ipert1,ipert2

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

 if (asr==0) return
 !call wrtout(std_out,' asria_corr: imposition of the ASR for the interatomic forces.','COLL')

! Remove d2asr
 do ipert2=1,natom
   do idir2=1,3
     do ipert1=1,natom
       do idir1=1,3
         d2cart(:,idir1,ipert1,idir2,ipert2)= d2cart(:,idir1,ipert1,idir2,ipert2) - d2asr(:,idir1,ipert1,idir2,ipert2)
       end do
     end do
   end do
 end do

end subroutine asria_corr
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/asrprs
!! NAME
!! asrprs
!!
!! FUNCTION
!! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry
!!
!! INPUTS
!!  asr=(3 => 1D systems, all elements are modified to give ASR and 
!!            rotational symmetry)
!!      (4 => 0D systems, all elements are modified to give ASR and 
!!            rotational symmetry)
!!  asrflg=(1 => the correction to enforce asr is computed from
!!           d2cart, but NOT applied;
!!          2 => one uses the previously determined correction)
!!  minvers=previously calculated inverted coefficient matrix
!!  mpert =maximum number of ipert
!!  natom=number of atom
!!  rotinv=(1,2,3 => for linear systems along x,y,z
!!          4 => non-linear molecule 
!!  xcart=cartesian coordinates of the ions
!!
!! OUTPUT
!!  (see side effects)
!!
!! SIDE EFFECTS
!! Input/Output:
!! d2cart=matrix of second derivatives of total energy, in cartesian coordinates
!! minvers=inverse of the supermatrix for future application of the corrections
!!
!! PARENTS
!!      anaddb,mkphbs
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine asrprs(asr,asrflag,rotinv,uinvers,vtinvers,singular,d2cart,mpert,natom,xcart)


!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 'asrprs'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: asr,asrflag,mpert,natom,rotinv
!arrays
 real(dp),intent(in) :: xcart(3,natom)
 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
 real(dp),intent(inout) :: singular(1:3*natom*(3*natom-1)/2)
 real(dp),intent(inout) :: uinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)
 real(dp),intent(inout) :: vtinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)

!Local variables-------------------------------
!scalars
 integer :: column,idir1,idir2,ii,info,ipert1,ipert2,jj,n3,row,superdim
 real(dp) :: rcond,test
! real(dp) :: tau ! tau is present but commented out in this routine
 character(len=500) :: message
!arrays
 integer :: check(3,natom,3)
 real(dp) :: tmp(natom,3,3),weightf(1:natom,1:natom)
 real(dp),allocatable :: d2cartold(:,:,:,:,:),d2vecc(:),d2veccnew(:),d2vecr(:)
 real(dp),allocatable :: d2vecrnew(:),superm(:,:),umatrix(:,:),vtmatrix(:)
 real(dp),allocatable :: work(:)

! *********************************************************************
 
 if(asr/=3 .and. asr/=4)then
   write(message,'(3a,i0)')&
&   'The argument asr should be 3 or 4,',ch10,&
&   'however, asr = ',asr
   MSG_BUG(message)
 end if
 
 if (asr==3.or.asr==4)then
   write(message, '(a,a)' ) ch10, &
&   'asrprs: imposition of the ASR for the interatomic forces and rotational invariance'
   call wrtout(std_out,message,'COLL')
 end if
 
 write(message,'(a,i6)')' asrflag is', asrflag 
 call wrtout(std_out,message,'COLL')
 call wrtout(ab_out,message,'COLL')

!variables for the dimensions of the matrices

!n1=3*natom*(3*natom-1)/2
!n2=9*natom
 n3=3*natom

 superdim=9*natom*(natom-1)/2+n3

 ABI_ALLOCATE(d2vecr,(1:superdim))
 ABI_ALLOCATE(d2vecc,(1:superdim))
 d2vecr=0d0
 d2vecc=0d0

!should be changed set to delta function for debugging
 weightf=1d0
!tau=1d-10
 do ii=1, natom
!  do jj=1, ii-1
!  weightf(ii,jj)= & 
!  &     ((xcart(1,ii)-xcart(1,jj))**2+(xcart(2,ii)-xcart(2,jj))**2+(xcart(3,ii)-xcart(3,jj))**2)**tau
!  enddo
   weightf(ii,ii)=0d0
 end do

 ABI_ALLOCATE(d2cartold,(2,3,mpert,3,mpert))

 d2cartold=d2cart

!setup vector with uncorrected derivatives

 do ipert1=1, natom
   do ipert2=1, ipert1-1
     do idir1=1,3
       do idir2=1,3
         row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
         if(abs(d2cart(1,idir1,ipert1,idir2,ipert2))<1d-6)then
           d2cart(1,idir1,ipert1,idir2,ipert2)=0d0
         else
           d2vecr(row)=4*weightf(ipert1,ipert2)*d2cart(1,idir1,ipert1,idir2,ipert2)
         end if
         if(abs(d2cart(2,idir1,ipert1,idir2,ipert2))<1d-6) then
           d2cart(2,idir1,ipert1,idir2,ipert2)=0d0
         else
           d2vecc(row)=4*weightf(ipert1,ipert2)*d2cart(2,idir1,ipert1,idir2,ipert2)
         end if
       end do
     end do
   end do
 end do

 if(asrflag==1) then !calculate the pseudo-inverse of the supermatrix
   ABI_ALLOCATE(superm,(1:superdim,1:superdim))
   
   superm=0d0

!  Setting up the supermatrix containing G, A, D

   do ipert1=1, natom
     do idir1=1, 3
!      Setting up G 
       idir2=mod(idir1,3)+1
       row=3*(ipert1-1)+idir1
       do ipert2=1, ipert1-1
         column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
         superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
         column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir2
         superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
       end do
       do ipert2=ipert1+1, natom
         column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir1-1)+rotinv
         superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
         column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir2-1)+rotinv
         superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
       end do
     end do
     do idir1=1, 3
!      Setting up D
       idir2=mod(idir1,3)+1
       ii=mod(idir1+1,3)+1
       do ipert2=1, ipert1-1
         row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
         column=9*natom*(natom-1)/2+3*(ipert1-1)+idir1
         superm(column,row)=superm(column,row)+xcart(idir2,ipert2)-xcart(idir2,ipert1)
         column=9*natom*(natom-1)/2+3*(ipert1-1)+ii
         superm(column,row)=superm(column,row)+xcart(ii,ipert1)-xcart(ii,ipert2)
         row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+rotinv
         column=9*natom*(natom-1)/2+3*(ipert2-1)+idir1
         superm(column,row)=superm(column,row)+xcart(idir2,ipert1)-xcart(idir2,ipert2)
         column=9*natom*(natom-1)/2+3*(ipert2-1)+ii
         superm(column,row)=superm(column,row)+xcart(ii,ipert2)-xcart(ii,ipert1)
       end do 
!      Setting up A
       do idir2=1, 3
         do ipert2=1, ipert1-1
           column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
           row=n3+column
           superm(column,row)=4*weightf(ipert1,ipert2)
         end do
       end do
     end do
   end do 

!  calculate the pseudo-inverse of the supermatrix

   ABI_ALLOCATE(work,(1:6*superdim))
   ABI_ALLOCATE(vtmatrix,(1:superdim))
   ABI_ALLOCATE(umatrix,(1:superdim,1:superdim))

!  singular value decomposition of superm

   call dgesvd('A','O',superdim,superdim,superm,superdim,singular,umatrix,superdim, &
&   vtmatrix, 1, work,6*superdim,info)

   ABI_DEALLOCATE(vtmatrix)
   ABI_DEALLOCATE(work)

   write(message, '(a,es16.8,es16.8)' )&
&   ' Largest and smallest values from svd', singular(1), singular(superdim) 
   call wrtout(std_out,message,'COLL')
   call wrtout(ab_out,message,'COLL')

!  Invert U and V**T, orthogonal matrices
   
   do ii=1, superdim
     do jj=1, superdim
       uinvers(ii,jj)=umatrix(jj,ii)
       vtinvers(ii,jj)=superm(jj,ii)
     end do
   end do
   
   ABI_DEALLOCATE(umatrix)
   ABI_DEALLOCATE(superm)
   
   write(message,'(a,a)')' asrprs: done with asrflag 1', ch10 
   call wrtout(std_out,message,'COLL')

 end if !asrflag=1
 
 if(asrflag==2) then 
   
   ABI_ALLOCATE(d2vecrnew,(1:superdim))
   ABI_ALLOCATE(d2veccnew,(1:superdim))

!  Calculate V**T**-1 Sigma**-1 U**-1 *rhs

   do ii=1, superdim
     d2vecrnew(ii)=0d0
     d2veccnew(ii)=0d0
     do jj=1, superdim
       d2vecrnew(ii)=d2vecrnew(ii)+uinvers(ii,jj)*d2vecr(jj)
       d2veccnew(ii)=d2veccnew(ii)+uinvers(ii,jj)*d2vecc(jj)
     end do
   end do

   rcond=1d-10*singular(1)
   do ii=1, superdim
     if(singular(ii)>rcond) then
       d2vecrnew(ii)=d2vecrnew(ii)/singular(ii)
       d2veccnew(ii)=d2veccnew(ii)/singular(ii)
     else
       d2vecrnew(ii)=0d0
       d2veccnew(ii)=0d0
     end if
   end do 

   do ii=1, superdim
     d2vecr(ii)=0d0
     d2vecc(ii)=0d0
     do jj=1, superdim
       d2vecr(ii)=d2vecr(ii)+vtinvers(ii,jj)*d2vecrnew(jj)
       d2vecc(ii)=d2vecc(ii)+vtinvers(ii,jj)*d2veccnew(jj)
     end do
   end do

!  Store vector back into the matrix of 2nd order derivates
   
   do ipert1=1, natom
     do ipert2=1, ipert1-1
       do idir1=1,3
         do idir2=1,3
           row=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
           d2cart(1,idir1,ipert1,idir2,ipert2)=d2vecr(row)
           d2cart(2,idir1,ipert1,idir2,ipert2)=d2vecc(row)
           d2cart(1,idir2,ipert2,idir1,ipert1)=d2vecr(row)
           d2cart(2,idir2,ipert2,idir1,ipert1)=d2vecc(row)
         end do
       end do
     end do
   end do
   
   ABI_DEALLOCATE(d2vecrnew)
   ABI_DEALLOCATE(d2veccnew)
   
   check=0

   do ipert1=1, natom
     do idir1=1, 3
       do idir2=1, 3
         d2cart(1,idir1,ipert1,idir2,ipert1)=0d0
         d2cart(2,idir1,ipert1,idir2,ipert1)=0d0
         tmp(ipert1,idir1,idir2)=0d0
         do ipert2=1, natom
           if(ipert2/=ipert1) then
             tmp(ipert1,idir1,idir2)=tmp(ipert1,idir1,idir2) & 
&             -d2cart(1,idir1,ipert1,idir2,ipert2) & 
&             -d2cart(1,idir2,ipert2,idir1,ipert1)
           end if
         end do
       end do
     end do
   end do 
   
   do ipert1=1, natom
     do idir1=1, 3
       do idir2=1, 3
         d2cart(1,idir1,ipert1,idir2,ipert1)=tmp(ipert1,idir1,idir2)/2
         d2cart(1,idir2,ipert1,idir1,ipert1)=d2cart(1,idir1,ipert1,idir2,ipert1)
       end do
     end do
   end do
   
   write(std_out,*) 'this should all be zero'
   
   do ipert1=1, natom
     do idir1=1, 3
       do idir2=1, 3
         test=0d0
         do ipert2=1, natom
           test=test+d2cart(1,idir1,ipert1,idir2,ipert2)+d2cart(1,idir2,ipert2,idir1,ipert1)
         end do
         write(std_out,'(i3,i3,i3,es11.3)') idir1,ipert1,idir2,test
         
         write(message, '(i3,i3,i3,es11.3)' ) idir1,ipert1,idir2,test 
         call wrtout(ab_out,message,'COLL')

       end do
     end do
   end do

   write(std_out,*) 'these as well'
   do ipert2=1, natom
     do idir1=1, 3
       do idir2=1, 3
         test=0d0
         do ipert1=1, natom
           test=test+d2cart(1,idir1,ipert1,idir2,ipert2)
         end do
         write(std_out,'(i3,i3,i3,i3,es11.3)') idir1,ipert1,idir2,ipert2,test
       end do
     end do
   end do

   write(message,'(a,a)')' asrprs: done with asrflag 2', ch10 
   call wrtout(std_out,message,'COLL')
   
 end if !ends asrflag=2

 ABI_DEALLOCATE(d2vecr)
 ABI_DEALLOCATE(d2vecc)
 ABI_DEALLOCATE(d2cartold)

end subroutine asrprs
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/cart29
!! NAME
!! cart29
!!
!!
!! FUNCTION
!! Transform a second-derivative matrix from reduced
!! coordinates to cartesian coordinates, and also
!! 1) add the ionic part of the effective charges,
!! 2) normalize the electronic dielectric tensor, and
!!    add the vacuum polarisation
!!
!! INPUTS
!!  blkflg(3,mpert,3,mpert,nblok)=
!!   ( 1 if the element of the dynamical matrix has been calculated ;
!!     0 otherwise )
!!  blkval(2,3,mpert,3,mpert,nblok)=DDB values
!!  gprimd(3,3)=basis vector in the reciprocal space
!!  iblok=number of the blok that will be transformed
!!  mpert =maximum number of ipert
!!  natom=number of atom
!!  nblok=number of blocks (dimension of blkflg and blkval)
!!  ntypat=number of atom types
!!  rprimd(3,3)=basis vector in the real space
!!  typat(natom)=integer label of each type of atom (1,2,...)
!!  ucvol=unit cell volume
!!  zion(ntypat)=charge corresponding to the atom type
!!
!! OUTPUT
!!  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
!!  2DTE matrix has been calculated correctly ; 0 otherwise )
!!  d2cart(2,3,mpert,3,mpert)=
!!    dynamical matrix, effective charges, dielectric tensor,....
!!    all in cartesian coordinates
!!
!! PARENTS
!!      gath3,m_ddb_blk
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine cart29(blkflg,blkval,carflg,d2cart,&
& gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)


!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 'cart29'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: iblok,mpert,natom,nblok,ntypat
 real(dp),intent(in) :: ucvol
!arrays
 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),typat(natom)
 integer,intent(out) :: carflg(3,mpert,3,mpert)
 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3)
 real(dp),intent(in) :: zion(ntypat)
 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: idir1,idir2,ii,ipert1,ipert2
!arrays
 integer :: flg1(3),flg2(3)
 real(dp) :: vec1(3),vec2(3)

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

!First, copy the data blok in place.
 d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok)

!Cartesian coordinates transformation (in two steps)
!First step
 do ipert1=1,mpert
   do ipert2=1,mpert
     do ii=1,2
       do idir1=1,3
         do idir2=1,3
           vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2)
!          Note here blkflg
           flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok)
         end do
         call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2)
         do idir2=1,3
           d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2)
!          And here carflg
           carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2)
         end do
       end do
     end do
   end do
 end do

!Second step
 do ipert1=1,mpert
   do ipert2=1,mpert
     do ii=1,2
       do idir2=1,3
         do idir1=1,3
           vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2)
!          Note here carflg
           flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2)
         end do
         call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2)
         do idir1=1,3
           d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1)
!          And here carflg again
           carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1)
         end do
       end do
     end do
   end do
 end do

!For the dielectric tensor, takes into account the volume
!of the unit cell, and add the unit matrix (polarization of the vacuum)
 do idir1=1,3
   do idir2=1,3
     do ii=1,2
       d2cart(ii,idir1,natom+2,idir2,natom+2)=&
&       -four_pi/ucvol*d2cart(ii,idir1,natom+2,idir2,natom+2)
     end do
   end do
 end do

 do idir1=1,3
   d2cart(1,idir1,natom+2,idir1,natom+2)=&
&   1.0_dp+d2cart(1,idir1,natom+2,idir1,natom+2)
 end do

!Add the ionic charges to delta z to get the effective charges
 do ipert1=1,natom
   do idir1=1,3
     d2cart(1,idir1,ipert1,idir1,natom+2)=&
&     zion(typat(ipert1))+d2cart(1,idir1,ipert1,idir1,natom+2)
   end do
 end do
 do ipert2=1,natom
   do idir2=1,3
     d2cart(1,idir2,natom+2,idir2,ipert2)=&
&     zion(typat(ipert2))+d2cart(1,idir2,natom+2,idir2,ipert2)
   end do
 end do

!For the piezoelectric tensor, takes into account the volume of the unit cell
 do ipert2=natom+3,natom+4
   do idir1=1,3
     do idir2=1,3
       do ii=1,2
         d2cart(ii,idir1,natom+2,idir2,ipert2)=&
&         (1.0_dp/ucvol)*d2cart(ii,idir1,natom+2,idir2,ipert2)
       end do
     end do
   end do
 end do

end subroutine cart29
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/cart39
!! NAME
!! cart39
!!
!!
!! FUNCTION
!! Transform a vector from reduced coordinates to cartesian coordinates,
!! taking into account the perturbation from which it was derived,
!! and also check the existence of the new values.
!!
!! INPUTS
!!  flg1(3)=tell if information of each component of vec1 is valid
!!  gprimd(3,3)=basis vector in the reciprocal space
!!  ipert=number of the perturbation
!!  natom=number of atom
!!  rprimd(3,3)=basis vector in the real space
!!  vec1(3)=input vector, in reduced coordinates
!!
!! OUTPUT
!!  flg2(3)=tell if information of each component of vec2 is valid
!!  vec2(3)=output vector, in cartesian coordinates
!!
!! PARENTS
!!      gath3,m_ddb_blk,m_dynmat
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine cart39(flg1,flg2,gprimd,ipert,natom,rprimd,vec1,vec2)


!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 'cart39'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: ipert,natom
!arrays
 integer,intent(in) :: flg1(3)
 integer,intent(out) :: flg2(3)
 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3),vec1(3)
 real(dp),intent(out) :: vec2(3)

!Local variables -------------------------
!scalars
 integer :: idir,ii

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

!Treat phonon-type perturbation
 if(ipert>=1.and.ipert<=natom)then

   do idir=1,3
     vec2(idir)=0.0_dp
     flg2(idir)=1
     do ii=1,3
       if(abs(gprimd(idir,ii))>1.0d-10)then
         if(flg1(ii)==1)then
           vec2(idir)=vec2(idir)+gprimd(idir,ii)*vec1(ii)
         else
           flg2(idir)=0
         end if
       end if
     end do
     if(flg2(idir)==0)vec2(idir)=0.0_dp
   end do

!  Treat electric field perturbation
 else if(ipert==natom+2) then
!  OCL SCALAR
   do idir=1,3
     vec2(idir)=0.0_dp
     flg2(idir)=1
!    OCL SCALAR
     do ii=1,3
       if(abs(rprimd(idir,ii))>1.0d-10)then
         if(flg1(ii)==1)then
           vec2(idir)=vec2(idir)+rprimd(idir,ii)*vec1(ii)/two_pi
         else
           flg2(idir)=0
         end if
       end if
     end do
     if(flg2(idir)==0)vec2(idir)=0.0_dp
   end do

!  Treat other perturbations
 else
   do idir=1,3
     vec2(idir)=vec1(idir)
     flg2(idir)=flg1(idir)
   end do
 end if

end subroutine cart39
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/chkph3
!! NAME
!! chkph3
!!
!! FUNCTION
!! Check the completeness of the dynamical matrix and eventually send a warning
!!
!! INPUTS
!!  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
!!  2DTE matrix has been calculated correctly ; 0 otherwise )
!!  idir = direction of the eventual electric field
!!  mpert =maximum number of ipert
!!  natom=number of atoms in unit cell
!!
!! OUTPUT
!!  eventually send a warning message
!!
!! PARENTS
!!      respfn
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine chkph3(carflg,idir,mpert,natom)


!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 'chkph3'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: idir,mpert,natom
!arrays
 integer,intent(in) :: carflg(3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: idir1,idir2,ipert1,ipert2,send
 character(len=500) :: message

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

 send=0

!Check the elements of the analytical part of the dynamical matrix
 do ipert2=1,natom
   do idir2=1,3
     do ipert1=1,natom
       do idir1=1,3
         if(carflg(idir1,ipert1,idir2,ipert2)==0)then
           send=1
         end if
       end do
     end do
   end do
 end do

!If some electric field is present
 if(idir/=0)then

!  Check the dielectric constant
   if(carflg(idir,natom+2,idir,natom+2)==0)then
     send=1
   end if

!  Check the effective charges
   do ipert1=1,natom
     do idir1=1,3
       if(carflg(idir1,ipert1,idir,natom+2)==0)then
         send=1
       end if
     end do
   end do

 end if

!If needed, send the message
 if(send==1)then
   write(message, '(a,a,a,a)' )&
&   ' chkph3 : WARNING -',ch10,&
&   '  The dynamical matrix was incomplete :',&
&   ' phonon frequencies may be wrong ...'
   call wrtout(std_out,message,'COLL')
   call wrtout(ab_out,message,'COLL')
 end if

end subroutine chkph3
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/chneu9
!! NAME
!! chneu9
!!
!! FUNCTION
!! Imposition of the Acoustic sum rule on the Effective charges
!! and suppress the imaginary part of the dynamical matrix
!!
!! INPUTS
!!  chneut=(0 => no ASR, 1 => equal repartition,2 => weighted repartition )
!!  mpert =maximum number of ipert
!!  natom=number of atom
!!  ntypat=number of types of atoms in unit cell
!!  selectz=selection of some parts of the effective charge tensor attached to one atom.
!!    (0=> no selection, 1=> trace only, 2=> symmetric part only)
!!  typat(natom)=type of the atom
!!  zion(ntypat)=atomic charge for every type of atom
!!
!! SIDE EFFECTS
!!  Input/Output
!!  d2cart=matrix of second derivatives of total energy, in cartesian
!!       coordinates
!!
!! PARENTS
!!      anaddb,gath3
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion)


!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 'chneu9'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: chneut,mpert,natom,ntypat,selectz
!arrays
 integer,intent(in) :: typat(natom)
 real(dp),intent(in) :: zion(ntypat)
 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: idir1,idir2,ii,ipert1,ipert2
 character(len=500) :: message
!arrays
 real(dp) :: sum(2)
 real(dp),allocatable :: wghtat(:)

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

 ABI_ALLOCATE(wghtat,(natom))

!In case of acoustic sum rule imposition, compute the weights on each atom.
 if (chneut==1)then

!  The weight is the same for all atom
   do ipert1=1,natom
     wghtat(ipert1)=1./natom
   end do

 else if (chneut==2) then

!  The weight is proportional to the diagonal electronic screening charge of the atom
   sum(1)=0.0_dp
   do ipert1=1,natom
     wghtat(ipert1)=0.0_dp
     do idir1=1,3
       wghtat(ipert1)=wghtat(ipert1)+&
&       d2cart(1,idir1,ipert1,idir1,natom+2)+&
&       d2cart(1,idir1,natom+2,idir1,ipert1)-2*zion(typat(ipert1))
     end do
     sum(1)=sum(1)+wghtat(ipert1)
   end do

!  Normalize the weights to unity
   do ipert1=1,natom
     wghtat(ipert1)=wghtat(ipert1)/sum(1)
   end do
 end if

!Calculation of the violation of the charge neutrality
!and imposition of the charge neutrality condition
 if (chneut/=0)then
   write(message, '(a,a,a,a,a,a,a)' )&
&   ' The violation of the charge neutrality conditions',ch10,&
&   ' by the effective charges is as follows :',ch10,&
&   '    atom        electric field',ch10,&
&   ' displacement     direction   '
   call wrtout(ab_out,message,'COLL')
   do idir1=1,3
     do idir2=1,3
       do ii=1,2
         sum(ii)=0.0_dp
         do ipert1=1,natom
           sum(ii)=sum(ii)+d2cart(ii,idir1,ipert1,idir2,natom+2)
         end do
         do ipert1=1,natom
           d2cart(ii,idir1,ipert1,idir2,natom+2)=&
&           d2cart(ii,idir1,ipert1,idir2,natom+2)-sum(ii)*wghtat(ipert1)
         end do
       end do
       write(message, '(i8,i16,2f16.6)' ) idir1,idir2,sum(1),sum(2)
       call wrtout(ab_out,message,'COLL')
     end do
   end do
   write(message, '(a)' )' '
   call wrtout(ab_out,message,'COLL')

!  The same for the symmetrical part
   do idir1=1,3
     do idir2=1,3
       do ii=1,2
         sum(ii)=0.0_dp
         do ipert2=1,natom
           sum(ii)=sum(ii)+d2cart(ii,idir1,natom+2,idir2,ipert2)
         end do
         do ipert2=1,natom
           d2cart(ii,idir1,natom+2,idir2,ipert2)=&
&           d2cart(ii,idir1,natom+2,idir2,ipert2)-sum(ii)*wghtat(ipert2)
         end do
       end do
     end do
   end do
 end if

!Selection of the trace of the effective charge tensor attached to each atom
 if(selectz==1)then
   do ipert1=1,natom
     do ii=1,2
       sum(ii)=0.0_dp
       do idir1=1,3
         sum(ii)=sum(ii)+d2cart(ii,idir1,ipert1,idir1,natom+2)
       end do
       do idir1=1,3
         do idir2=1,3
           d2cart(ii,idir1,ipert1,idir2,natom+2)=0.0_dp
         end do
       end do
       do idir1=1,3
         d2cart(ii,idir1,ipert1,idir1,natom+2)=sum(ii)/3.0_dp
       end do
     end do
   end do
!  Do the same for the symmetrical part of d2cart
   do ipert2=1,natom
     do ii=1,2
       sum(ii)=0.0_dp
       do idir1=1,3
         sum(ii)=sum(ii)+d2cart(ii,idir1,natom+2,idir1,ipert2)
       end do
       do idir1=1,3
         do idir2=1,3
           d2cart(ii,idir1,natom+2,idir2,ipert2)=0.0_dp
         end do
       end do
       do idir1=1,3
         d2cart(ii,idir1,natom+2,idir1,ipert2)=sum(ii)/3.0_dp
       end do
     end do
   end do
 end if

!Selection of the symmetric part of the effective charge tensor attached to each atom
 if(selectz==2)then
   do ipert1=1,natom
     do ii=1,2
       do idir1=1,3
         do idir2=1,3
           sum(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
&           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
           d2cart(ii,idir1,ipert1,idir2,natom+2)=sum(ii)
           d2cart(ii,idir2,ipert1,idir1,natom+2)=sum(ii)
         end do
       end do
     end do
   end do
!  Do the same for the symmetrical part of d2cart
   do ipert1=1,natom
     do ii=1,2
       do idir1=1,3
         do idir2=1,3
           sum(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
&           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
           d2cart(ii,idir1,ipert1,idir2,natom+2)=sum(ii)
           d2cart(ii,idir2,ipert1,idir1,natom+2)=sum(ii)
         end do
       end do
     end do
   end do
 end if

!Write the effective charge tensor
 write(message, '(a,a,a,a,a,a,a)' )&
& ' Effective charge tensors after ',ch10,&
& ' imposition of the charge neutrality,',ch10,&
& ' and eventual restriction to some part :',ch10,&
& '   atom    displacement  '
 call wrtout(ab_out,message,'COLL')

 do ipert1=1,natom
   do idir1=1,3
     write(message, '(2i10,3es16.6)' )ipert1,idir1,&
&     (d2cart(1,idir1,ipert1,idir2,natom+2),idir2=1,3)
     call wrtout(ab_out,message,'COLL')
   end do
 end do

!Zero the imaginary part of the dynamical matrix
 write(message, '(a)' )&
& ' Now, the imaginary part of the dynamical matrix is zeroed '
 call wrtout(ab_out,message,'COLL')
 call wrtout(std_out,message,'COLL')

 do ipert1=1,natom
   do ipert2=1,natom
     do idir1=1,3
       do idir2=1,3
         d2cart(2,idir1,ipert1,idir2,ipert2)=0.0_dp
       end do
     end do
   end do
 end do

 ABI_DEALLOCATE(wghtat)

end subroutine chneu9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/d2sym3
!! NAME
!! d2sym3
!!
!! FUNCTION
!! Given a set of calculated elements of the 2DTE matrix d2,
!! build (nearly) all the other matrix elements that can be build using symmetries.
!!
!! 1. Perform first some completion by symmetrisation (exchange) 
!!    over the two defining perturbations
!! 2. For each element, uses every symmetry, and build the element, in case
!!    EITHER all the needed elements are available,
!!    OR the only missing is itself
!!    OR the perturbation is the electric field, in a diamond
!!    symmetry (the last case was coded rather dirty)
!!
!! INPUTS
!!  indsym(4,nsym,natom)=indirect indexing array : for each
!!   isym,iatom, fourth element is label of atom into which iatom is sent by
!!   INVERSE of symmetry operation isym; first three elements are the primitive
!!   translations which must be subtracted after the transformation to get back
!!   to the original unit cell.
!!  mpert =maximum number of ipert
!!  natom= number of atoms
!!  nsym=number of space group symmetries
!!  qpt(3)=wavevector of the perturbation
!!  symq(4,2,nsym)= (integer) three first numbers define the G vector ;
!!   fourth number is zero if the q-vector is not preserved,
!!              is 1 otherwise
!!   second index is one without time-reversal symmetry,
!!                two with time-reversal symmetry
!!  symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
!!  symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
!!  timrev=1 if the time-reversal symmetry preserves the wavevector,
!!     modulo a reciprocal lattice vector, timrev=0 otherwise
!!
!! SIDE EFFECTS
!!  Input/Output
!!  d2(2,3,mpert,3,mpert)= matrix of the 2DTE
!!  blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
!!     matrix has been calculated ; 0 otherwise)
!!
!! NOTES
!!   The complete search would be to have the possibility
!!   of a set of missing elements. See notes of July 2, 1994,
!!   in the blue notebook 'computer codes'
!!   The partial solution adopted here takes into
!!   account some mirror symmetries
!!   as well as the tetrahedral symmetry of the diamond lattice
!!   On 010331, replaced the loops up to mpert by loops up to
!!   natom+2, because of a crash bug under Windows. However,
!!   the problem lies likely in the use of the indsym array.
!!
!! PARENTS
!!      completeperts,m_ddb_blk,respfn
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine d2sym3(blkflg,d2,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev)


!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 'd2sym3'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: mpert,natom,nsym,timrev
!arrays
 integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym),symrec(3,3,nsym)
 integer,intent(in) :: symrel(3,3,nsym)
 integer,intent(inout) :: blkflg(3,mpert,3,mpert)
 real(dp),intent(in) :: qpt(3)
 real(dp),intent(inout) :: d2(2,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: exch12,found,idir1,idir2,idisy1,idisy2,ipert1,ipert2
 integer :: ipesy1,ipesy2,isign,isym,ithree,itirev,noccur,quit,quit1
! integer :: ii,jj ! these appear below in commented out DEBUG sections
 real(dp) :: arg1,arg2,im,norm,re,sumi,sumr,xi,xr
!arrays
 integer :: sym1(3,3),sym2(3,3)

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

!Here look after exchange of 1 and 2 axis,
!for electric field in diamond symmetry
 exch12=0
 if ( qpt(1)**2+qpt(2)**2+qpt(3)**2 < tol16 ) then

   do isym=1,nsym
     exch12=1
     if(symrel(1,1,isym)/=0)exch12=0
     if(symrel(1,2,isym)/=1)exch12=0
     if(symrel(1,3,isym)/=0)exch12=0
     if(symrel(2,1,isym)/=1)exch12=0
     if(symrel(2,2,isym)/=0)exch12=0
     if(symrel(2,3,isym)/=0)exch12=0
     if(symrel(3,1,isym)/=0)exch12=0
     if(symrel(3,2,isym)/=0)exch12=0
     if(symrel(3,3,isym)/=1)exch12=0
!    if(exch12==1) write(std_out,*)' d2sym3 : found exchange 1 2 =',isym
     if(exch12==1)exit
   end do

 end if

!Exchange of perturbations

!Consider two cases : either time-reversal symmetry
!conserves the wavevector, or not
 if(timrev==0)then

!  do ipert1=1,mpert  See notes
   do ipert1=1,min(natom+2,mpert)
     do idir1=1,3
!      write(std_out,*) 'd2sym3: Im', d2(2,idir1,ipert1,idir1,ipert1)

!      Since the matrix is hermitian, the diagonal elements are real
       d2(2,idir1,ipert1,idir1,ipert1)=0.0_dp
!      
!      do ipert2=1,mpert See notes
       do ipert2=1,min(natom+2,mpert)
         do idir2=1,3

!          If an element exists
           if(blkflg(idir1,ipert1,idir2,ipert2)==1)then

!            Either complete the symmetric missing element
             if(blkflg(idir2,ipert2,idir1,ipert1)==0)then

               d2(1,idir2,ipert2,idir1,ipert1)= d2(1,idir1,ipert1,idir2,ipert2)
               d2(2,idir2,ipert2,idir1,ipert1)=-d2(2,idir1,ipert1,idir2,ipert2)

               blkflg(idir2,ipert2,idir1,ipert1)=1

!              Or symmetrize (the matrix is hermitian) in case both exists
!              (Note : this opportunity has been disabled for more
!              obvious search for bugs in the code )
!              else
!              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
!              sumi=d2(1,idir2,ipert2,idir1,ipert1)-d2(1,idir1,ipert1,idir2,ipert2)
!              d2(1,idir2,ipert2,idir1,ipert1)=sumr/2.0_dp
!              d2(1,idir1,ipert1,idir2,ipert2)=sumr/2.0_dp
!              d2(2,idir2,ipert2,idir1,ipert1)=sumi/2.0_dp
!              d2(2,idir1,ipert1,idir2,ipert2)=-sumi/2.0_dp
             end if
           end if

         end do
       end do

     end do
   end do

!  Here, case with time-reversal symmetry
 else

!  do ipert1=1,mpert See notes
   do ipert1=1,min(natom+2,mpert)
     do idir1=1,3
!      do ipert2=1,mpert See notes
       do ipert2=1,min(natom+2,mpert)
         do idir2=1,3
           d2(2,idir1,ipert1,idir2,ipert2)=0.0_dp

!          If an element exists
           if(blkflg(idir1,ipert1,idir2,ipert2)==1)then

!            Either complete the symmetric missing element
             if(blkflg(idir2,ipert2,idir1,ipert1)==0)then

               d2(1,idir2,ipert2,idir1,ipert1)=d2(1,idir1,ipert1,idir2,ipert2)
               blkflg(idir2,ipert2,idir1,ipert1)=1

!              Or symmetrize (the matrix is hermitian) in case both exists
!              (Note : this opportunity has been disabled for more
!              obvious search for bugs in the code )
!              else
!              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
!              d2(1,idir2,ipert2,idir1,ipert1)=sumr/2.0_dp
!              d2(1,idir1,ipert1,idir2,ipert2)=sumr/2.0_dp
             end if

           end if
         end do
       end do
     end do
   end do

 end if
!Big Big Loop : symmetrize three times, because
!of some cases in which one element is not yet available
!at the first pass, and even at the second one !
 do ithree=1,3

!  Big loop on all elements
!  do ipert1=1,mpert See notes
   do ipert1=1,min(natom+2,mpert)
     do idir1=1,3
!      do ipert2=1,mpert See notes
       do ipert2=1,min(natom+2,mpert)
         do idir2=1,3

!          Will get element (idir1,ipert1,idir2,ipert2)
!          so this element should not yet be present ...
           if(blkflg(idir1,ipert1,idir2,ipert2)/=1)then

             d2(1,idir1,ipert1,idir2,ipert2)=0.0_dp
             d2(2,idir1,ipert1,idir2,ipert2)=0.0_dp

!            Loop on all symmetries, including time-reversal
             quit1=0
             do isym=1,nsym
               do itirev=1,2
                 isign=3-2*itirev

                 if(symq(4,itirev,isym)/=0)then
                   found=1

!                  Here select the symmetric of ipert1
                   if(ipert1<=natom)then
                     ipesy1=indsym(4,isym,ipert1)
                     sym1(:,:)=symrec(:,:,isym)
                   else if(ipert1==(natom+2) .and. qpt(1)**2+qpt(2)**2+qpt(3)**2 < tol16)then
                     ipesy1=ipert1
                     sym1(:,:)=symrel(:,:,isym)
                   else
                     found=0
                   end if

!                  Here select the symmetric of ipert2
                   if(ipert2<=natom)then
                     ipesy2=indsym(4,isym,ipert2)
                     sym2(:,:)=symrec(:,:,isym)
                   else if(ipert2==(natom+2) .and. qpt(1)**2+qpt(2)**2+qpt(3)**2 < tol16)then
                     ipesy2=ipert2
                     sym2(:,:)=symrel(:,:,isym)
                   else
                     found=0
                   end if

!                  Now that a symmetric perturbation has been obtained,
!                  including the expression of the symmetry matrix, see
!                  if the symmetric values are available
                   if( found==1 ) then

                     sumr=0.0_dp
                     sumi=0.0_dp
                     noccur=0
                     quit=0
                     do idisy1=1,3
                       do idisy2=1,3
                         if(sym1(idir1,idisy1)/=0 .and. sym2(idir2,idisy2)/=0 )then
                           if(blkflg(idisy1,ipesy1,idisy2,ipesy2)==1)then
                             sumr=sumr+sym1(idir1,idisy1)*sym2(idir2,idisy2)*&
&                             d2(1,idisy1,ipesy1,idisy2,ipesy2)
                             sumi=sumi+sym1(idir1,idisy1)*sym2(idir2,idisy2)*&
&                             d2(2,idisy1,ipesy1,idisy2,ipesy2)

!                            Here, in case the symmetric of the element
!                            is the element, or the symmetric with
!                            respect to permutation of perturbations
!                            (some more conditions on the time-reversal
!                            symmetry must be fulfilled although)
                           else if(  idisy1==idir1 .and. ipesy1==ipert1&
&                             .and. idisy2==idir2 .and. ipesy2==ipert2&
&                             .and.(isign==1 .or. timrev==1 &
&                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
&                             then
                             noccur=noccur+sym1(idir1,idisy1)*sym2(idir2,idisy2)
                           else if(  idisy1==idir2 .and. ipesy1==ipert2&
&                             .and. idisy2==idir1 .and. ipesy2==ipert1&
&                             .and.(isign==-1 .or. timrev==1&
&                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
&                             then
                             noccur=noccur+sym1(idir1,idisy1)*sym2(idir2,idisy2)

!                            Here, electric field case
                           else if( exch12==1 .and. &
&                             ipert1==natom+2 .and. ipert2==natom+2&
&                             .and.(( idisy1+idir1 ==3 &
&                             .and. idisy2==3 .and. idir2==3)&
&                             .or. ( idisy1+idir2 ==3&
&                             .and. idisy2==3 .and. idir1==3)&
&                             .or. ( idisy2+idir2 ==3&
&                             .and. idisy1==3 .and. idir1==3)&
&                             .or. ( idisy2+idir1 ==3&
&                             .and. idisy1==3 .and. idir2==3)))&
&                             then
                             noccur=noccur+sym1(idir1,idisy1)*sym2(idir2,idisy2)

                           else
!                            Not found
                             found=0
                             quit=1
                             exit
                           end if

                         end if
                       end do
                       if(quit==1)exit
                     end do
                   end if

!                  Now, if still found, put the correct value into array d2
                   if(found==1)then

!                    In case of phonons, need to take into account the
!                    time-reversal symmetry, and the shift back to the unit cell
!                    
!                    XG990712 : I am not sure this must be kept for electric field ...
!                    1) Consider time-reversal symmetry
                     sumi=isign*sumi

                     if(ipert1<=natom .and. ipert2<=natom)then
!                      2) Shift the atoms back to the unit cell.
                       arg1=two_pi*( qpt(1)*indsym(1,isym,ipert1)&
&                       +qpt(2)*indsym(2,isym,ipert1)&
&                       +qpt(3)*indsym(3,isym,ipert1) )
                       arg2=two_pi*( qpt(1)*indsym(1,isym,ipert2)&
&                       +qpt(2)*indsym(2,isym,ipert2)&
&                       +qpt(3)*indsym(3,isym,ipert2) )
                       re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
!                      XG010117 Must use isign
                       im=isign*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))
                     else
                       re=1.0_dp
                       im=0.0_dp
                     end if

!                    Final check, could still fail if the
!                    element was its own symmetric
                     if (abs(1.0_dp-re*noccur) < 1.0d-6 .and. abs(im*noccur)  < 1.0d-6) then
                       found=0

!                      DEBUG
!                      write(std_out,*)' element is its own symmetric ...'
!                      ENDDEBUG
                     end if

                   end if

                   if(found==1)then

!                    DEBUG
!                    write(std_out,*)' all found !  isym, isign= ',isym,isign
!                    write(std_out,'(9i4)' )((sym1(ii,jj),ii=1,3),jj=1,3)
!                    write(std_out,'(9i4)' )((sym2(ii,jj),ii=1,3),jj=1,3)
!                    write(std_out,*)sumr,sumi
!                    ENDDEBUG

                     if(noccur==0)then
                       d2(1,idir1,ipert1,idir2,ipert2)=re*sumr-im*sumi
                       d2(2,idir1,ipert1,idir2,ipert2)=re*sumi+im*sumr
                     else
!                      See page July 2, 1994 in computer codes notebook
                       xr=re*sumr-im*sumi
                       xi=re*sumi+im*sumr
                       norm=1.0_dp+noccur**2-2.0_dp*re*noccur
                       xr=xr/norm
                       xi=xi/norm
                       d2(1,idir1,ipert1,idir2,ipert2)=&
&                       (1.0_dp-re*noccur)*xr-im*noccur*xi
                       d2(2,idir1,ipert1,idir2,ipert2)=&
&                       (1.0_dp-re*noccur)*xi+im*noccur*xr
                     end if

!                    The element has been constructed !
                     blkflg(idir1,ipert1,idir2,ipert2)=1

                     quit1=1
                     exit ! Exit loop on symmetry operations
                   end if

!                  End loop on all symmetries + time-reversal
                 end if
               end do
               if(quit1==1)exit
             end do

           end if
         end do ! End big loop on all elements
       end do
     end do
   end do

 end do !  End Big Big Loop

!DEBUG
!write(std_out,*)' d2sym3 : exit d2sym3 '
!write(std_out,*)' Write blkflg '
!do ipert1=9,12
!do idir1=1,3
!do ipert2=9,12
!do idir2=1,3
!write(std_out,*)'idir1,ipert1,idir2,ipert2,blkflg=',idir1,ipert1,idir2,ipert2,blkflg(idir1,ipert1,idir2,ipert2)
!enddo
!enddo
!enddo
!enddo
!write(std_out,*)' d2sym3 : timrev=',timrev
!ENDDEBUG

end subroutine d2sym3
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/q0dy3_apply
!! NAME
!! q0dy3_apply
!!
!! FUNCTION
!! Takes care of the inclusion of the ewald q=0 term in the dynamical
!! matrix - corrects the dyew matrix provided as input
!!
!! INPUTS
!!  dyewq0(3,3,natom) = part needed to correct
!!    the dynamical matrix for atom self-interaction.
!!  natom= number of atom in the unit cell
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!  dyew(2,3,natom,3,natom)= dynamical matrix corrected on output
!!
!! NOTES
!! Should be used just after each call to ewald3, for both
!! q==0 and the real wavelength.
!!
!! The q0dy3_apply should be used in conjunction with the subroutine
!! ewald3 (or ewald9):
!! First, the call of ewald3 with q==0 should be done,
!!   then the call to q0dy3_calc will produce
!!   the dyewq0 matrix from the (q=0) dyew matrix
!! Second, the call of ewald3 with the real q (either =0 or diff 0)
!!   should be done, then the call to q0dy3_apply
!!   will produce the correct dynamical matrix dyew starting from
!!   the previously calculated dyewq0 and the bare(non-corrected)
!!   dyew matrix
!!
!! PARENTS
!!      m_dynmat,m_ifc,respfn
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine q0dy3_apply(natom,dyewq0,dyew)


!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 'q0dy3_apply'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom
!arrays
 real(dp),intent(in) :: dyewq0(3,3,natom)
 real(dp),intent(inout) :: dyew(2,3,natom,3,natom)

!Local variables -------------------------
!scalars
 integer :: ia,mu,nu

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

 do mu=1,3
   do nu=1,3
     do ia=1,natom
       dyew(1,mu,ia,nu,ia)=dyew(1,mu,ia,nu,ia)-dyewq0(mu,nu,ia)
     end do
   end do
 end do

end subroutine q0dy3_apply
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/q0dy3_calc
!! NAME
!! q0dy3_calc
!!
!! FUNCTION
!! Calculate the q=0 correction term to the dynamical matrix
!!
!! INPUTS
!!  dyew(2,3,natom,3,natom)= dynamical matrix
!!    input, non-corrected, for q=0 if option=1 or 2
!!  natom= number of atom in the unit cell
!!  option= either 1 or 2:
!!     1: use dyew to calculate dyewq0 symmetrical form
!!     2: use dyew to calculate dyewq0 symmetrical form
!!
!! OUTPUT
!!  dyewq0(3,3,natom) = part needed to correct
!!    the dynamical matrix for atom self-interaction.
!!
!! NOTES
!! Should be used just after each call to ewald3, for both
!! q==0 and the real wavelength.
!!
!! If option=1 or 2, q0dy3_calc uses an Ewald dynamical matrix at q=0,
!! called dyew, to produce a contracted form called dyewq0 :
!! either:
!!   in an unsymmetrical form (if option=1), or
!!   in a symmetrical form (if option=2).
!!
!! The q0dy3_calc should be used in conjunction with the subroutine
!! ewald3 (or ewald9).
!! First, the call of ewald3 with q==0 should be done ,
!!   then the call to q0dy3_calc will produce
!!   the dyewq0 matrix from the (q=0) dyew matrix
!! Second, the call of ewald3 with the real q (either =0 or diff 0)
!!   should be done, then the call to q0dy3_apply
!!   will produce the correct dynamical matrix dyew starting from
!!   the previously calculated dyewq0 and the bare(non-corrected)
!!   dyew matrix
!!
!! PARENTS
!!      hybrid9,m_ifc,respfn
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine q0dy3_calc(natom,dyewq0,dyew,option)


!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 'q0dy3_calc'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom,option
!arrays
 real(dp),intent(in) :: dyew(2,3,natom,3,natom)
 real(dp),intent(out) :: dyewq0(3,3,natom)

!Local variables -------------------------
!scalars
 integer :: ia,ib,mu,nu
 character(len=500) :: message

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

 if(option==1.or.option==2)then
   do mu=1,3
     do nu=1,3
       do ia=1,natom
         dyewq0(mu,nu,ia)=0.0_dp
         do ib=1,natom
           dyewq0(mu,nu,ia)=dyewq0(mu,nu,ia)+dyew(1,mu,ia,nu,ib)
         end do
       end do
     end do
   end do
 else
   write (message, '(3a)')&
&   'option should be 1 or 2.',ch10,&
&   'action: correct calling routine'
   MSG_BUG(message)
 end if

 if(option==2)then
   do ia=1,natom
     do mu=1,3
       do nu=mu,3
         dyewq0(mu,nu,ia)=(dyewq0(mu,nu,ia)+dyewq0(nu,mu,ia))/2
         dyewq0(nu,mu,ia)=dyewq0(mu,nu,ia)
       end do
     end do
   end do
 end if

end subroutine q0dy3_calc
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/symdyma
!!
!! NAME
!! symdyma
!!
!! FUNCTION
!! Symmetrize the dynamical matrices
!!
!! INPUTS
!! indsym(4,msym*natom)=indirect indexing array : for each
!!   isym,iatom, fourth element is label of atom into which iatom is sent by
!!   INVERSE of symmetry operation isym; first three elements are the primitive
!!   translations which must be subtracted after the transformation to get back
!!   to the original unit cell.
!! natom=number of atoms in unit cell
!! nsym=number of space group symmetries
!! qptn(3)=normalized phonon wavevector
!! rprimd(3,3)=dimensional primitive translations (bohr)
!! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
!!
!! SIDE EFFECTS
!! Input/Output
!! dmati(2*3*natom*3*natom)=dynamical matrices relative to the q
!!  points of the B.Z. sampling
!!
!! NOTES
!! the procedure of the symmetrization of the dynamical matrix follows the
!! equations in: Hendrikse et al., Computer Phys. Comm. 86, 297 (1995)
!!
!! TODO
!! A full description of the equations should be included
!!
!! PARENTS
!!      phfrq3,relaxpol
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine symdyma(dmati,indsym,natom,nsym,qptn,rprimd,symrel)


!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 'symdyma'
 use interfaces_32_util
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom,nsym
!arrays
 integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
 real(dp),intent(in) :: qptn(3),rprimd(3,3)
 real(dp),intent(inout) :: dmati(2*3*natom*3*natom)

!Local variables -------------------------
!scalars
 integer :: i1,i2,iat,idir,ii,index,isign,isym,itirev,jat,jdir,jj,kk,ll
 integer :: niat,njat,timrev
 real(dp) :: arg1,arg2,dmint,im,re,sumi,sumr
!arrays
 integer :: indij(natom,natom),symq(4,2,nsym),symrec(3,3,nsym)
 real(dp) :: TqR(3,3),TqS_(3,3),dynmat(2,3,natom,3,natom)
 real(dp) :: dynmatint(2*nsym,2,3,natom,3,natom),gprimd(3,3)
 real(dp) :: symcart(3,3,nsym)

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

!0) initializations
 call matr3inv(rprimd,gprimd)
 do isym=1,nsym
   call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
 end do

 TqR=zero
 TqS_=zero
 dynmat=zero

!Note: dynmat is used as work space here
 i1=0
 do iat=1,natom
   do idir=1,3
     i1=i1+1
     i2=0
     do jat=1,natom
       do jdir=1,3
         i2=i2+1
         index=i1+3*natom*(i2-1)
         dynmat(1,idir,iat,jdir,jat)=dmati(2*index-1)
         dynmat(2,idir,iat,jdir,jat)=dmati(2*index  )
       end do
     end do
   end do
 end do

!Transform symrel to cartesian coordinates (RC coding)
!do isym=1,nsym
!symcart(:,:,isym)=matmul(rprimd,matmul(dble(symrel(:,:,isym)),gprimd))
!end do

!Coding from symdm9
 do isym=1,nsym
   do jj=1,3
     symcart(:,jj,isym)=zero
     do kk=1,3
       do ll=1,3
         symcart(:,jj,isym)=symcart(:,jj,isym)+rprimd(:,kk)*gprimd(jj,ll)*symrel(kk,ll,isym)
       end do
     end do
   end do
 end do


!Get the symq of the CURRENT Q POINT
!mjv: set prtvol=0 for production runs.
 call symq3(nsym,qptn,symq,symrec,timrev,prtvol=0)

 indij(:,:)=0
 dynmatint=zero

 do isym=1,nsym                                 ! loop over all the symmetries
!  write(std_out,*) 'current symmetry',isym
   do itirev=1,2                                 ! loop over the time-reversal symmetry
     isign=3-2*itirev                             ! to take into accont the time-reversal
!    write(std_out,*) 'timereversal',isign
     if (symq(4,itirev,isym)==1) then             ! isym belongs to the wave vector point group
!      write(std_out,*) 'isym belongs to the wave vector point group'
       do iat=1,natom                              ! loop over the atoms
         do jat=1,natom                            ! loop over the atoms
           niat=indsym(4,isym,iat)                   ! niat={R|t}iat
           njat=indsym(4,isym,jat)                   ! njat={R|t}jat
           indij(niat,njat)=indij(niat,njat)+1
!          write(std_out,'(a,5i5)') 'current status:',iat,jat,niat,njat,indij(niat,njat)
!          phase calculation, arg1 and arg2 because of two-atom derivative
           arg1=two_pi*( qptn(1)*indsym(1,isym,iat)+&
&           qptn(2)*indsym(2,isym,iat)+&
&           qptn(3)*indsym(3,isym,iat) )
           arg2=two_pi*( qptn(1)*indsym(1,isym,jat)+&
&           qptn(2)*indsym(2,isym,jat)+&
&           qptn(3)*indsym(3,isym,jat) )
           re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
           im=isign*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))

           do idir=1,3                               ! loop over displacements
             do jdir=1,3                              ! loop over displacements
!              we pick the (iat,jat) (3x3) block of the dyn.mat.
               sumr=zero
               sumi=zero
               do ii=1,3
                 do jj=1,3
                   sumr=sumr+symcart(idir,ii,isym)*dynmat(1,ii,niat,jj,njat)*symcart(jdir,jj,isym)
                   sumi=sumi+symcart(idir,ii,isym)*dynmat(2,ii,niat,jj,njat)*symcart(jdir,jj,isym)
                 end do
               end do
               sumi=isign*sumi

               dynmatint(nsym*(itirev-1)+isym,1,idir,iat,jdir,jat)=re*sumr-im*sumi
               dynmatint(nsym*(itirev-1)+isym,2,idir,iat,jdir,jat)=re*sumi+im*sumr

             end do
           end do
         end do
       end do ! end treatment of the (iat,jat) (3x3) block of dynmat
     end if ! symmetry check
   end do ! time-reversal
 end do ! symmetries


!4) make the average, get the final symmetric dynamical matrix
 do iat=1,natom
   do jat=1,natom
     do idir=1,3
       do jdir=1,3
         dmint=zero
         do isym=1,2*nsym
           dmint=dmint+dynmatint(isym,1,idir,iat,jdir,jat)
         end do
         dynmat(1,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
         dmint=zero
         do isym=1,2*nsym
           dmint=dmint+dynmatint(isym,2,idir,iat,jdir,jat)
         end do
         dynmat(2,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
       end do
     end do
   end do
 end do

 i1=0
 do iat=1,natom
   do idir=1,3
     i1=i1+1
     i2=0
     do jat=1,natom
       do jdir=1,3
         i2=i2+1
         index=i1+3*natom*(i2-1)
         dmati(2*index-1)=dynmat(1,idir,iat,jdir,jat)
         dmati(2*index  )=dynmat(2,idir,iat,jdir,jat)
       end do
     end do
   end do
 end do

end subroutine symdyma
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/wings3
!! NAME
!! wings3
!!
!! FUNCTION
!!  Suppress the wings of the cartesian 2DTE for which
!!  the diagonal element is not known
!!
!! INPUTS
!!  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
!!  2DTE matrix has been calculated correctly ; 0 otherwise )
!!  d2cart(2,3,mpert,3,mpert)=
!!   dynamical matrix, effective charges, dielectric tensor,....
!!   all in cartesian coordinates
!!  mpert =maximum number of ipert
!!
!! OUTPUT
!!  d2cart(2,3,mpert,3,mpert) without the wings
!!
!! PARENTS
!!      respfn
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine wings3(carflg,d2cart,mpert)


!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 'wings3'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: mpert
!arrays
 integer,intent(inout) :: carflg(3,mpert,3,mpert)
 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: idir,idir1,ipert,ipert1

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

 do ipert=1,mpert
   do idir=1,3
     if(carflg(idir,ipert,idir,ipert)==0)then
       do ipert1=1,mpert
         do idir1=1,3
           carflg(idir,ipert,idir1,ipert1)=0
           carflg(idir1,ipert1,idir,ipert)=0
           d2cart(1,idir,ipert,idir1,ipert1)=0.0_dp
           d2cart(2,idir,ipert,idir1,ipert1)=0.0_dp
           d2cart(1,idir1,ipert1,idir,ipert)=0.0_dp
           d2cart(2,idir1,ipert1,idir,ipert)=0.0_dp
         end do
       end do
     end if
   end do
 end do

end subroutine wings3
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/asrif9
!!
!! NAME
!! asrif9
!!
!! FUNCTION
!! Imposes the Acoustic Sum Rule to Interatomic Forces
!!
!! INPUTS
!! asr= Option for the imposition of the ASR
!!  0 => no ASR,
!!  1 => modify "asymmetrically" the diagonal element
!!  2 => modify "symmetrically" the diagonal element
!! natom= Number of atoms in the unit cell
!! nrpt= Number of R points in the Big Box
!! rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
!!  These coordinates are normalized (=> * acell(3)!!)
!! wghatm(natom,natom,nrpt)= Weight associated to the couple of atoms and the R vector
!! atmfrc(2,3,natom,3,natom,nrpt)= Interatomic Forces
!!
!! OUTPUT
!! atmfrc(2,3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces
!!
!! TODO
!! List of ouput should be included.
!!
!! PARENTS
!!      hybrid9,m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm)


!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 'asrif9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: asr,natom,nrpt
!arrays
 real(dp),intent(in) :: rpt(3,nrpt),wghatm(natom,natom,nrpt)
 real(dp),intent(inout) :: atmfrc(2,3,natom,3,natom,nrpt)

!Local variables -------------------------
!scalars
 integer :: found,ia,ib,irpt,izero,mu,nu
 real(dp) :: sum

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

 DBG_ENTER("COLL")

 if(asr==1.or.asr==2)then

   found=0
!  Search for the R vector which is equal to ( 0 , 0 , 0 )
!  This vector leaves the atom a on itself !
   do irpt=1,nrpt
     if (all(abs(rpt(:,irpt))<=1.0d-10)) then
       found=1
       izero=irpt
     end if
     if (found==1) exit
   end do

   if(found==0)then
     MSG_BUG('Not able to find the vector R=(0,0,0).')
   end if

   do mu=1,3
     do nu=1,3
       do ia=1,natom
         sum=0.0_dp
         do ib=1,natom

!          Get the sum of interatomic forces acting on the atom ia,
!          either in a symmetrical manner, or an unsymmetrical one.
           if(asr==1)then
             do irpt=1,nrpt
               sum=sum+wghatm(ia,ib,irpt)*atmfrc(1,mu,ia,nu,ib,irpt)
             end do
           else if(asr==2)then
             do irpt=1,nrpt
               sum=sum+&
&               (wghatm(ia,ib,irpt)*atmfrc(1,mu,ia,nu,ib,irpt)+&
&               wghatm(ia,ib,irpt)*atmfrc(1,nu,ia,mu,ib,irpt))/2
             end do
           end if
         end do

!        Correct the self-interaction in order to fulfill the ASR
         atmfrc(1,mu,ia,nu,ia,izero)=atmfrc(1,mu,ia,nu,ia,izero)-sum
         if(asr==2)then
           atmfrc(1,nu,ia,mu,ia,izero)=atmfrc(1,mu,ia,nu,ia,izero)
         end if

       end do
     end do
   end do
 end if

 DBG_EXIT("COLL")

end subroutine asrif9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/make_bigbox
!! NAME
!! make_bigbox
!!
!! FUNCTION
!! Helper functions that faciliates the generation  of a Big Box containing 
!! all the R points in the cartesian real space needed to Fourier Transform
!! the dynamical matrix into its corresponding interatomic force.
!! See bigbx9 for the algorithm.
!!
!! INPUTS
!! brav= Bravais Lattice (1=S.C.;2=F.C.C.;3=BCC;4=Hex.)
!! ngqpt(3)= Numbers used to generate the q points to sample the
!!   Brillouin zone using an homogeneous grid
!! nqshft= number of q-points in the repeated cell for the Brillouin zone sampling
!!  When nqshft is not 1, but 2 or 4 (only other allowed values), 
!!  the limits for the big box have to be extended by a factor of 2.
!! rprim(3,3)= Normalized coordinates in real space  !!! IS THIS CORRECT?
!!
!! OUTPUT
!! nprt= Number of R points in the Big Box
!! rpt(3,nrpt)= Canonical coordinates of the R points in the unit cell
!!  These coordinates are normalized (=> * acell(3)!!)
!!  The array is allocated here wit the proper dimension. Client code is responsible
!!  for the deallocation.
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine make_bigbox(brav,ngqpt,nqshft,rprim,nrpt,rpt)


!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 'make_bigbox'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: brav,nqshft
 integer,intent(out) :: nrpt
!arrays
 integer,intent(in) :: ngqpt(3)
 real(dp),intent(in) :: rprim(3,3)
 real(dp),allocatable,intent(out) :: rpt(:,:)

!Local variables -------------------------
!scalars
 integer :: choice,mrpt
!arrays
 real(dp) :: dummy_rpt(3,1)

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

 DBG_ENTER("COLL")

 ! Compute the number of points (cells) in real space
 choice=0
 call bigbx9(brav,choice,1,ngqpt,nqshft,mrpt,rprim,dummy_rpt)

 ! Now we can allocate and calculate the points and the weights.
 nrpt = mrpt
 ABI_ALLOCATE(rpt,(3,nrpt))

 choice=1
 call bigbx9(brav,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)

 DBG_EXIT("COLL")

end subroutine make_bigbox
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/bigbx9
!! NAME
!! bigbx9
!!
!! FUNCTION
!! Generation of a Big Box containing all the R points in the
!! cartesian real space needed to Fourier Transforms the dynamical
!! matrix into its corresponding interatomic force.
!!
!! INPUTS
!! brav= Bravais Lattice (1=S.C.;2=F.C.C.;3=BCC;4=Hex.)
!! choice= if 0, simply count nrpt ; if 1, checks that the input mrpt
!!   is the same as nrpt, and generate rpt(3,mrpt)
!! mrpt=dimension of rpt
!! ngqpt(3)= Numbers used to generate the q points to sample the
!!  Brillouin zone using an homogeneous grid
!! nqshft= number of q-points in the repeated cell for the Brillouin zone sampling
!!  When nqshft is not 1, but 2 or 4 (only other allowed values), 
!!  the limits for the big box have to be extended by a factor of 2.
!! rprim(3,3)= Normalized coordinates in real space  !!! IS THIS CORRECT?
!!
!! OUTPUT
!! nprt= Number of R points in the Big Box
!! rpt(3,mrpt)= Canonical coordinates of the R points in the unit cell
!!  These coordinates are normalized (=> * acell(3)!!)
!!  (output only if choice=1)
!!
!! PARENTS
!!      m_dynmat
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine bigbx9(brav,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)


!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 'bigbx9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: brav,choice,mrpt,nqshft
 integer,intent(out) :: nrpt
!arrays
 integer,intent(in) :: ngqpt(3)
 real(dp),intent(in) :: rprim(3,3)
 real(dp),intent(out) :: rpt(3,mrpt)

!Local variables -------------------------
!In some cases, the atoms coordinates are not packed in the
! [0,1]^3 cube. Then, the parameter "buffer" might be increased,
!to search relevant pairs of atoms in bigger boxes than usual.
!scalars
 integer,parameter :: buffer=1
 integer :: irpt,lim1,lim2,lim3,lqshft,r1,r2,r3
 character(len=500) :: msg

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

 lqshft=1
 if(nqshft/=1)lqshft=2

!Simple Cubic Lattice
 if (brav==1) then
   lim1=((ngqpt(1)/2)+1)*lqshft+buffer
   lim2=((ngqpt(2)/2)+1)*lqshft+buffer
   lim3=((ngqpt(3)/2)+1)*lqshft+buffer
   nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
   if(choice/=0)then
     if (nrpt/=mrpt) then
       write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
       MSG_BUG(msg)
     end if
     irpt=0
     do r1=-lim1,lim1
       do r2=-lim2,lim2
         do r3=-lim3,lim3
           irpt=irpt+1
           rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
           rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
           rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
         end do
       end do
     end do
   end if

!  Face Centered Cubic Lattice
 else if (brav==2) then
   lim1=((ngqpt(1)+3)/4)*lqshft+buffer
   lim2=((ngqpt(2)+3)/4)*lqshft+buffer
   lim3=((ngqpt(3)+3)/4)*lqshft+buffer
   nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*4
   if(choice/=0)then
     if (nrpt/=mrpt) then
       write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
       MSG_BUG(msg)
     end if
     irpt=0
     do r1=-lim1,lim1
       do r2=-lim2,lim2
         do r3=-lim3,lim3
           irpt=irpt+4
           rpt(1,irpt-3)=r1
           rpt(2,irpt-3)=r2
           rpt(3,irpt-3)=r3
           rpt(1,irpt-2)=r1
           rpt(2,irpt-2)=r2+0.5
           rpt(3,irpt-2)=r3+0.5
           rpt(1,irpt-1)=r1+0.5
           rpt(2,irpt-1)=r2
           rpt(3,irpt-1)=r3+0.5
           rpt(1,irpt)=r1+0.5
           rpt(2,irpt)=r2+0.5
           rpt(3,irpt)=r3
         end do
       end do
     end do
   end if

!  Body Centered Cubic Lattice
 else if (brav==3) then
   lim1=((ngqpt(1)+3)/4)*lqshft+buffer
   lim2=((ngqpt(2)+3)/4)*lqshft+buffer
   lim3=((ngqpt(3)+3)/4)*lqshft+buffer
   nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*2
   if(choice/=0)then
     if(nrpt/=mrpt) then
       write(msg,'(2(a,i0))')' nrpt= ',nrpt,' is not equal to mrpt= ',mrpt
       MSG_BUG(msg)
     end if
     irpt=0
     do r1=-lim1,lim1
       do r2=-lim2,lim2
         do r3=-lim3,lim3
           irpt=irpt+2
           rpt(1,irpt-1)=r1
           rpt(2,irpt-1)=r2
           rpt(3,irpt-1)=r3
           rpt(1,irpt)=r1+0.5
           rpt(2,irpt)=r2+0.5
           rpt(3,irpt)=r3+0.5
         end do
       end do
     end do
   end if

!  Hexagonal Lattice
 else if (brav==4) then
   lim1=(ngqpt(1)+1)*lqshft+buffer
   lim2=(ngqpt(2)+1)*lqshft+buffer
   lim3=((ngqpt(3)/2)+1)*lqshft+buffer
   nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
   if(choice/=0)then
     if(nrpt/=mrpt)then
       write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt=',mrpt
       MSG_BUG(msg)
     end if
     irpt=0
     do r1=-lim1,lim1
       do r2=-lim2,lim2
         do r3=-lim3,lim3
           irpt=irpt+1
           rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
           rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
           rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
         end do
       end do
     end do
   end if

 else
   write(msg,'(a,i0,a)')' The value of brav= ',brav,' is not allowed (should be 1, 2 or 4).'
   MSG_BUG(msg)
 end if

end subroutine bigbx9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/canat9
!! NAME
!! canat9
!!
!! FUNCTION
!! Transforms an atom whose coordinates (xred*rprim) would not be
!! in the chosen unit cell used to generate the interatomic forces
!! to its correspondent (rcan) in canonical coordinates.
!!
!! INPUTS
!! brav= Bravais Lattice (1=S.C.;2=F.C.C.;3=BCC;4=Hex.)
!! natom= Number of atoms in the unit cell
!! rprim(3,3)= Normalized coordinates  of primitive vectors
!!
!! OUTPUT
!! rcan(3,natom)  = Atomic position in canonical coordinates
!! trans(3,natom) = Atomic translations : xred = rcan + trans
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE


subroutine canat9(brav,natom,rcan,rprim,trans,xred)


!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 'canat9'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: brav,natom
!arrays
 real(dp),intent(in) :: rprim(3,3),xred(3,natom)
 real(dp),intent(out) :: rcan(3,natom),trans(3,natom)

!Local variables -------------------------
!scalars
 integer :: found,iatom,ii
 character(len=500) :: message
!arrays
 real(dp) :: dontno(3,4),rec(3),rok(3),shift(3),tt(3)

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

 DBG_ENTER("COLL")


!Normalization of the cartesian atomic coordinates
!If not normalized : rcan(i) <- rcan(i) * acell(i)

 do iatom=1,natom
   rcan(:,iatom)=xred(1,iatom)*rprim(:,1)+&
&   xred(2,iatom)*rprim(:,2)+&
&   xred(3,iatom)*rprim(:,3)
 end do

!Study of the different cases for the Bravais lattice :

!Simple Cubic Lattice
 if (brav==1) then

   do iatom=1,natom

!    Canon will produces these coordinate transformations
!    (Note : here we still use reduced coordinates )
     call wrap2_pmhalf(xred(1,iatom),rok(1),shift(1))
     call wrap2_pmhalf(xred(2,iatom),rok(2),shift(2))
     call wrap2_pmhalf(xred(3,iatom),rok(3),shift(3))

!    New coordinates : rcan
     rcan(:,iatom)=rok(1)*rprim(:,1)+&
&     rok(2)*rprim(:,2)+&
&     rok(3)*rprim(:,3)
!    Translations between New and Old coordinates
     tt(:)=xred(1,iatom)*rprim(:,1)+&
&     xred(2,iatom)*rprim(:,2)+&
&     xred(3,iatom)*rprim(:,3)
     trans(:,iatom)=tt(:)-rcan(:,iatom)
   end do

 else if (brav==2) then
!  Face Centered Lattice
!  Special possible translations in the F.C.C. case
   dontno(:,:)=0.0_dp
   dontno(2,2)=0.5_dp
   dontno(3,2)=0.5_dp
   dontno(1,3)=0.5_dp
   dontno(3,3)=0.5_dp
   dontno(1,4)=0.5_dp
   dontno(2,4)=0.5_dp
   do iatom=1,natom
     found=0
     do ii=1,4
       if (found==1) exit
!      Canon will produces these coordinate transformations
       call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
       call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
       call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
!      In the F.C.C., ABS[ Ri ] + ABS[ Rj ] < or = 1/2
!      The equal signs has been treated using a tolerance parameter
!      not to have twice the same point in the unit cell !
       rok(1)=rok(1)-1.0d-10
       rok(2)=rok(2)-2.0d-10
       rok(3)=rok(3)-5.0d-10
       if (abs(rok(1))+abs(rok(2))<=0.5_dp) then
         if (abs(rok(1))+abs(rok(3))<=0.5_dp) then
           if (abs(rok(2))+abs(rok(3))<=0.5_dp) then
             tt(:)=rcan(:,iatom)
!            New coordinates : rcan
             rcan(1,iatom)=rok(1)+1.0d-10
             rcan(2,iatom)=rok(2)+2.0d-10
             rcan(3,iatom)=rok(3)+5.0d-10
!            Translations between New and Old coordinates
             trans(:,iatom)=tt(:)-rcan(:,iatom)
             found=1
           end if
         end if
       end if
     end do
   end do

 else if (brav==3) then
!  Body Centered Cubic Lattice
!  Special possible translations in the B.C.C. case

   dontno(:,1)=0.0_dp
   dontno(:,2)=0.5_dp
   do iatom=1,natom
     found=0
     do ii=1,2
       if (found==1) exit
!      Canon will produce these coordinate transformations
       call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
       call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
       call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
!      In the F.C.C., ABS[ Ri ] < or = 1/2
!      and    ABS[ R1 ] + ABS[ R2 ] + ABS[ R3 ] < or = 3/4
!      The equal signs has been treated using a tolerance parameter
!      not to have twice the same point in the unit cell !
       rok(1)=rok(1)-1.0d-10
       rok(2)=rok(2)-2.0d-10
       rok(3)=rok(3)-5.0d-10
       if(abs(rok(1))+abs(rok(2))+abs(rok(3))<=0.75_dp)then
         if ( abs(rok(1))<=0.5_dp .and.&
&         abs(rok(2))<=0.5_dp .and.&
&         abs(rok(3))<=0.5_dp       ) then
           tt(:)=rcan(:,iatom)
!          New coordinates : rcan
           rcan(1,iatom)=rok(1)+1.0d-10
           rcan(2,iatom)=rok(2)+2.0d-10
           rcan(3,iatom)=rok(3)+5.0d-10
!          Translations between New and Old coordinates
           trans(:,iatom)=tt(:)-rcan(:,iatom)
           found=1
         end if
       end if
     end do
   end do

 else if (brav==4) then
!  Hexagonal Lattice
!  In this case, it is easier first to work in reduced coordinates space !
   do iatom=1,natom
!    Passage from the reduced space to the "lozenge" cell
     rec(1)=xred(1,iatom)-0.5_dp
     rec(2)=xred(2,iatom)-0.5_dp
     rec(3)=xred(3,iatom)
!    Canon will produces these coordinate transformations
     call wrap2_pmhalf(rec(1),rok(1),shift(1))
     call wrap2_pmhalf(rec(2),rok(2),shift(2))
     call wrap2_pmhalf(rec(3),rok(3),shift(3))
     rec(1)=rok(1)+0.5_dp
     rec(2)=rok(2)+0.5_dp
     rec(3)=rok(3)
!    Passage in Cartesian Normalized Coordinates
     rcan(:,iatom)=rec(1)*rprim(:,1)+&
&     rec(2)*rprim(:,2)+&
&     rec(3)*rprim(:,3)
!    Use of a tolerance parameter not to have twice the same point
!    in the unit cell !
     rcan(1,iatom)=rcan(1,iatom)-1.0d-10
     rcan(2,iatom)=rcan(2,iatom)-2.0d-10
!    Passage to the honey-com hexagonal unit cell !
     if (rcan(1,iatom)>0.5_dp) then
       rcan(1,iatom)=rcan(1,iatom)-1.0_dp
     end if
     if (rcan(1,iatom)>0.0_dp.and.rcan(1,iatom)+sqrt(3.0_dp)*rcan&
&     (2,iatom)>1.0_dp) then
       rcan(1,iatom)=rcan(1,iatom)-0.5_dp
       rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
     end if
     if (rcan(1,iatom)<=0.0_dp.and.sqrt(3.0_dp)*rcan(2,iatom)-rcan&
&     (1,iatom)>1.0_dp) then
       rcan(1,iatom)=rcan(1,iatom)+0.5_dp
       rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
     end if
!    Translations between New and Old coordinates
     tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)&
&     +xred(3,iatom)*rprim(:,3)
     trans(:,iatom)=tt(:)-rcan(:,iatom)
   end do

!  End of the possible cases for brav : 1, 2, 4.
 else

   write(message, '(a,i0,a,a,a)' )&
&   'The required value of brav=',brav,' is not available.',ch10,&
&   'It should be 1,2 or 4 .'
   MSG_BUG(message)
 end if

 call wrtout(std_out,' Canonical Atomic Coordinates ',"COLL")

 do iatom=1,natom
   write(message, '(a,i5,3es18.8)' )' atom',iatom,rcan(1,iatom),rcan(2,iatom),rcan(3,iatom)
   call wrtout(std_out,message,'COLL')
 end do

 DBG_EXIT("COLL")

end subroutine canat9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/canct9
!!
!! NAME
!! canct9
!!
!! FUNCTION
!! Convert from canonical coordinates to cartesian coordinates
!! a vector defined by its index=ib+natom*(irpt-1)
!!
!! INPUTS
!! acell(3)=length scales by which rprim is to be multiplied
!! gprim(3,3)=dimensionless primitive translations in reciprocal space
!! index= index of the atom
!! natom=number of atoms in unit cell
!! nrpt= Number of R points in the Big Box
!! rcan(3,natom)=canonical coordinates of atoms
!! rprim(3,3)=dimensionless primitive translations in real space
!! rpt(3,nrpt)=canonical coordinates of the points in the BigBox.
!!
!! OUTPUT
!! ib=number of the atom in the unit cell
!! irpt= number of the unit cell to which belong the atom
!! rcart(3)=cartesian coordinate of the atom indexed by index.
!!
!! PARENTS
!!      hybrid9,m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine canct9(acell,gprim,ib,index,irpt,natom,nrpt,rcan,rcart,rprim,rpt)


!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 'canct9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: index,natom,nrpt
 integer,intent(out) :: ib,irpt
!arrays
 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3)
 real(dp),intent(in) :: rpt(3,nrpt)
 real(dp),intent(out) :: rcart(3)

!Local variables -------------------------
!scalars
 integer :: jj
!arrays
 real(dp) :: xred(3)

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

 irpt=(index-1)/natom+1
 ib=index-natom*(irpt-1)

!Transform the canonical coordinates to reduced coord.
 do jj=1,3
   xred(jj)=gprim(1,jj)*(rpt(1,irpt)+rcan(1,ib))&
&   +gprim(2,jj)*(rpt(2,irpt)+rcan(2,ib))&
&   +gprim(3,jj)*(rpt(3,irpt)+rcan(3,ib))
 end do

!Then to cartesian coordinates (here the position of
!the atom b)

 do jj=1,3
   rcart(jj)=xred(1)*acell(1)*rprim(jj,1)+&
&   xred(2)*acell(2)*rprim(jj,2)+&
&   xred(3)*acell(3)*rprim(jj,3)
 end do

end subroutine canct9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/chkrp9
!! NAME
!! chkrp9
!!
!! FUNCTION
!! Check if the rprim used for the definition of the unit cell (in the
!! inputs) are consistent with the rprim used in the routine generating
!! the Big Box needed to generate the interatomic forces.
!!
!! INPUTS
!! brav=bravais lattice (1=simple lattice,2=face centered lattice,
!!  3=centered lattice,4=hexagonal lattice)
!! rprimd(3,3)=dimensional primitive translations for real space (bohr)
!!
!! OUTPUT
!!  (only checking)
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine chkrp9(brav,rprim)


!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 'chkrp9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: brav
!arrays
 real(dp),intent(in) :: rprim(3,3)

!Local variables -------------------------
!scalars
 integer :: ii,jj
 character(len=500) :: message

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

 if (brav==1) then
!  Simple Cubic Lattice No condition in this case !
   continue

 else if (brav==2) then
!  Face Centered Lattice
   do ii=1,3
     do jj=1,3
       if (  ( ii==jj .and. abs(rprim(ii,jj))>tol10) .or.&
&       ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
         write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
&         'The input variable rprim does not correspond to the',ch10,&
&         'fixed rprim to be used with brav=2 and ifcflag=1 :',ch10,&
&         '   0  1/2  1/2',ch10,&
&         '  1/2  0   1/2',ch10,&
&         '  1/2 1/2   0 ',ch10,&
&         'Action: rebuild your DDB by using the latter rprim.'
         MSG_ERROR(message)
       end if
     end do
   end do

 else if (brav==3) then
!  Body Centered Cubic Lattice
   do ii=1,3
     do jj=1,3
       if (  ( ii==jj .and. abs(rprim(ii,jj)+.5_dp)>tol10) .or.&
&       ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
         write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
&         'The input variable rprim does not correspond to the',ch10,&
&         'fixed rprim to be used with brav=3 and ifcflag=1 :',ch10,&
&         '  -1/2  1/2  1/2',ch10,&
&         '   1/2 -1/2  1/2',ch10,&
&         '   1/2  1/2 -1/2',ch10,&
&         'Action: rebuild your DDB by using the latter rprim.'
         MSG_ERROR(message)
       end if
     end do
   end do

 else if (brav==4) then
!  Hexagonal Lattice
   if (abs(rprim(1,1)-1.0_dp)>tol10 .or. &
&   abs(rprim(3,3)-1.0_dp)>tol10 .or. &
&   abs(rprim(2,1)      )>tol10 .or. &
&   abs(rprim(3,1)      )>tol10 .or. &
&   abs(rprim(1,3)      )>tol10 .or. &
&   abs(rprim(2,3)      )>tol10 .or. &
&   abs(rprim(3,2)      )>tol10 .or. &
&   abs(rprim(1,2)+0.5_dp)>tol10 .or. &
&   abs(rprim(2,2)-0.5_dp*sqrt(3.0_dp))>tol10 ) then
     write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
&     'The input variable rprim does not correspond to the',ch10,&
&     'fixed rprim to be used with brav=4 and ifcflag=1 :',ch10,&
&     '   1      0      0',ch10,&
&     '  -1/2 sqrt[3]/2 0',ch10,&
&     '   0      0      1',ch10,&
&     'Action: rebuild your DDB by using the latter rprim.'
     MSG_ERROR(message)
   end if

 else

   write(message, '(a,i4,a,a,a,a,a)' )&
&   'The value of brav=',brav,' is not allowed.',ch10,&
&   'Only  1,2,3 or 4 are allowed.',ch10,&
&   'Action: change the value of brav in your input file.'
   MSG_ERROR(message)
 end if

end subroutine chkrp9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/dist9
!! NAME
!! dist9
!!
!! FUNCTION
!! Compute the distance between atoms
!!
!! INPUTS
!! acell(3)=length scales by which rprim is to be multiplied
!! dist(natom,natom,nrpt)=distances between atoms
!! gprim(3,3)=dimensionless primitive translations in reciprocal space
!! natom=number of atoms in unit cell
!! nrpt= Number of R points in the Big Box
!! rcan(3,natom)=canonical coordinates of atoms
!! rprim(3,3)=dimensionless primitive translations in real space
!! rpt(3,nrpt)=cartesian coordinates of the points in the BigBox.
!!
!! OUTPUT
!! dist(natom,natom,nrpt)=distances between atoms
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt)


!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 'dist9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
! real(dp),intent(in) :: ucvol
!scalars
 integer,intent(in) :: natom,nrpt
!arrays
 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3)
 real(dp),intent(in) :: rpt(3,nrpt)
 real(dp),intent(out) :: dist(natom,natom,nrpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,ii,irpt
!arrays
 real(dp) :: ra(3),rb(3),rdiff(3),red(3),rptcar(3),xred(3)

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

!BIG loop on all generic atoms
 do ia=1,natom

!  First transform canonical coordinates to reduced coordinates
   do ii=1,3
     xred(ii)=gprim(1,ii)*rcan(1,ia)+gprim(2,ii)*rcan(2,ia)&
&     +gprim(3,ii)*rcan(3,ia)
   end do

!  Then to cartesian coordinates
   ra(:)=xred(1)*acell(1)*rprim(:,1)+xred(2)*acell(2)*rprim(:,2)+&
&   xred(3)*acell(3)*rprim(:,3)

   do ib=1,natom
     do ii=1,3
       xred(ii)=gprim(1,ii)*rcan(1,ib)+gprim(2,ii)*rcan(2,ib)&
&       +gprim(3,ii)*rcan(3,ib)
     end do
     do ii=1,3
       rb(ii)=xred(1)*acell(1)*rprim(ii,1)+&
&       xred(2)*acell(2)*rprim(ii,2)+&
&       xred(3)*acell(3)*rprim(ii,3)
     end do

     do irpt=1,nrpt
!      First transform it to reduced coordinates
       do ii=1,3
         red(ii)=gprim(1,ii)*rpt(1,irpt)+gprim(2,ii)*rpt(2,irpt)&
&         +gprim(3,ii)*rpt(3,irpt)
       end do
!      Then to cartesian coordinates
       do ii=1,3
         rptcar(ii)=red(1)*acell(1)*rprim(ii,1)+&
&         red(2)*acell(2)*rprim(ii,2)+&
&         red(3)*acell(3)*rprim(ii,3)
       end do
       do ii=1,3
         rdiff(ii)=-rptcar(ii)+ra(ii)-rb(ii)
       end do
       dist(ia,ib,irpt)=(rdiff(1)**2+rdiff(2)**2+rdiff(3)**2)**0.5

     end do

   end do
 end do

end subroutine dist9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/ftifc_q2r
!!
!! NAME
!! ftifc_q2r
!!
!! FUNCTION
!!   Generates the Fourier transform of the dynamical matrices
!!   to obtain interatomic forces (real space).
!!
!! INPUTS
!! dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
!! gprim(3,3)= Normalized coordinates in reciprocal space
!! natom= Number of atoms in the unit cell
!! nqpt= Number of q points in the Brillouin zone
!! nrpt= Number of R points in the Big Box
!! rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
!!           These coordinates are normalized (=> * acell(3)!!)
!! spqpt(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
!!
!! OUTPUT
!! atmfrc(2,3,natom,3,natom,nrpt)= Interatomic Forces in real space !!
!!  We used the imaginary part just for debugging !
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine ftifc_q2r(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt)


!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 'ftifc_q2r'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom,nqpt,nrpt
!arrays
 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
 real(dp),intent(out) :: atmfrc(2,3,natom,3,natom,nrpt)
 real(dp),intent(in) :: dynmat(2,3,natom,3,natom,nqpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,iqpt,irpt,mu,nu
 real(dp) :: im,kr,re
!arrays
 real(dp) :: kk(3)

! *********************************************************************
 
 DBG_ENTER("COLL")

!Interatomic Forces from Dynamical Matrices
 atmfrc(:,:,:,:,:,:)=0.0_dp
 do irpt=1,nrpt
   do iqpt=1,nqpt

!    Calculation of the k coordinates in Normalized Reciprocal coordinates
     kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
     kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
     kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)

!    Product of k and r
     kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)

!    Get the phase factor
     re=cos(two_pi*kr)
     im=sin(two_pi*kr)

!    Now, big inner loops on atoms and directions
!    The indices are ordered to give better speed
     do ib=1,natom
       do nu=1,3
         do ia=1,natom
           do mu=1,3
!            Real and imaginary part of the interatomic forces
             atmfrc(1,mu,ia,nu,ib,irpt)=atmfrc(1,mu,ia,nu,ib,irpt)&
&             +re*dynmat(1,mu,ia,nu,ib,iqpt)&
&             +im*dynmat(2,mu,ia,nu,ib,iqpt)
!            The imaginary part should be equal to zero !!!!!!
!            atmfrc(2,mu,ia,nu,ib,irpt)=atmfrc(2,mu,ia,nu,ib,irpt)
!            &          +re*dynmat(2,mu,ia,nu,ib,iqpt)
!            &          -im*dynmat(1,mu,ia,nu,ib,iqpt)
           end do
         end do
       end do
     end do

   end do
 end do

!The sum has to be weighted by a normalization factor of 1/nqpt
 atmfrc(:,:,:,:,:,:)=atmfrc(:,:,:,:,:,:)/nqpt

 DBG_EXIT("COLL")

end subroutine ftifc_q2r
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/ftifc_r2q
!!
!! NAME
!! ftifc_r2q
!!
!! FUNCTION
!! (r->q):
!!   Generates the Fourier transform of the interatomic forces
!!   to obtain dynamical matrices (reciprocal space).
!!
!! INPUTS
!! atmfrc(2,3,natom,3,natom,nrpt)= Interatomic Forces in real space !!
!!  We used the imaginary part just for debugging !
!! gprim(3,3)= Normalized coordinates in reciprocal space
!! natom= Number of atoms in the unit cell
!! nqpt= Number of q points in the Brillouin zone
!! nrpt= Number of R points in the Big Box
!! rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
!!           These coordinates are normalized (=> * acell(3)!!)
!! spqpt(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
!! wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
!!
!! OUTPUT
!! dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
!!
!! PARENTS
!!      m_dynmat
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine ftifc_r2q(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt,wghatm)


!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 'ftifc_r2q'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom,nqpt,nrpt
!arrays
 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
 real(dp),intent(in) :: wghatm(natom,natom,nrpt)
 real(dp),intent(in) :: atmfrc(2,3,natom,3,natom,nrpt)
 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,iqpt,irpt,mu,nu
 real(dp) :: facti,factr,im,kr,re
!arrays
 real(dp) :: kk(3)

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

 dynmat(:,:,:,:,:,:)=0.0_dp
 do iqpt=1,nqpt
   do irpt=1,nrpt

!    Calculation of the k coordinates in Normalized Reciprocal coordinates
     kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)* gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
     kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)* gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
     kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)* gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)

!    Product of k and r
     kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)

!    Get phase factor
     re=cos(two_pi*kr)
     im=sin(two_pi*kr)

!    Inner loop on atoms and directions
     do ib=1,natom
       do ia=1,natom
         if(abs(wghatm(ia,ib,irpt))>1.0d-10)then
           factr=re*wghatm(ia,ib,irpt)
           facti=im*wghatm(ia,ib,irpt)
           do nu=1,3
             do mu=1,3
!              Real and imaginary part of the dynamical matrices
               dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)&
&               +factr*atmfrc(1,mu,ia,nu,ib,irpt)
!              Atmfrc should be real
!              &       -im*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt)
               dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)&
&               +facti*atmfrc(1,mu,ia,nu,ib,irpt)
!              Atmfrc should be real
!              &        +re*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt)
             end do
           end do
         end if
       end do
     end do
   end do
 end do

end subroutine ftifc_r2q
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/ifclo9
!! NAME
!! ifclo9
!!
!! FUNCTION
!! Convert from cartesian coordinates to local coordinates
!! the 3*3 interatomic force constant matrix
!!
!! INPUTS
!! ifccar(3,3)= matrix of interatomic force constants in cartesian
!!  coordinates
!! vect1(3)= cartesian coordinates of the first local vector
!! vect2(3)= cartesian coordinates of the second local vector
!! vect3(3)= cartesian coordinates of the third local vector
!!
!! OUTPUT
!! ifcloc(3,3)= matrix of interatomic force constants in local coordinates
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine ifclo9(ifccar,ifcloc,vect1,vect2,vect3)


!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 'ifclo9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!arrays
 real(dp),intent(in) :: ifccar(3,3),vect1(3),vect2(3),vect3(3)
 real(dp),intent(out) :: ifcloc(3,3)

!Local variables -------------------------
!scalars
 integer :: ii,jj
!arrays
 real(dp) :: work(3,3)

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

 do jj=1,3
   do ii=1,3
     work(jj,ii)=0.0_dp
   end do
   do ii=1,3
     work(jj,1)=work(jj,1)+ifccar(jj,ii)*vect1(ii)
     work(jj,2)=work(jj,2)+ifccar(jj,ii)*vect2(ii)
     work(jj,3)=work(jj,3)+ifccar(jj,ii)*vect3(ii)
   end do
 end do

 do jj=1,3
   do ii=1,3
     ifcloc(ii,jj)=0.0_dp
   end do
   do ii=1,3
     ifcloc(1,jj)=ifcloc(1,jj)+vect1(ii)*work(ii,jj)
     ifcloc(2,jj)=ifcloc(2,jj)+vect2(ii)*work(ii,jj)
     ifcloc(3,jj)=ifcloc(3,jj)+vect3(ii)*work(ii,jj)
   end do
 end do

end subroutine ifclo9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/wght9
!! NAME
!! wght9
!!
!! FUNCTION
!! Generates a weight to each R points of the Big Box and for each pair of atoms
!! For each R points included in the space generates by moving
!! the unit cell around each atom; the weight will be one.
!! Border conditions are provided.
!! The R points outside the chosen space will have a 0 weight.
!!
!! INPUTS
!! brav = Bravais lattice (1=S.C.;2=F.C.C.;4=Hex.)
!! gprim(3,3)= Normalized coordinates in reciprocal space
!! natom= Number of atoms in the unit cell
!! ngqpt(6)= Numbers used to sample the Brillouin zone
!! nqpt= Number of q points used in the homogeneous grid
!!  sampling the Brillouin zone
!! nqshft=number of shift vectors in the repeated cell
!! nrpt=Number of R points in the Big Box
!! qshft(3,nqshft)=vectors that will be used to determine
!!  the shifts from (0. 0. 0.)
!! rcan(3,natom)=Atomic position in canonical coordinates
!! rpt(3,nprt)=Canonical coordinates of the R points in the unit cell
!!  These coordinates are normalized (=> * acell(3))
!!
!! OUTPUT
!! wghatm(natom,natom,nrpt)= Weight associated to the couple of atoms and the R vector
!!  The vector r(atom2)-r(atom1)+rpt should be inside the moving box
!! ngqpt(6)= can be modified
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine wght9(brav,gprim,natom,ngqpt,nqpt,nqshft,nrpt,qshft,rcan,rpt,wghatm)


!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 'wght9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: brav,natom,nqpt,nqshft,nrpt
!arrays
 integer,intent(inout) :: ngqpt(9)
 real(dp),intent(in) :: gprim(3,3),qshft(3,4),rcan(3,natom),rpt(3,nrpt)
 real(dp),intent(out) :: wghatm(natom,natom,nrpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,ii,iqshft,irpt,jqshft,nbordh,tok
 real(dp) :: factor,sum
 character(len=500) :: message
!arrays
 integer :: nbord(9)
 real(dp) :: rdiff(9),red(3,3)

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

 DBG_ENTER("COLL")


!First analyze the vectors qshft
 if(nqshft/=1)then

   if(brav==4)then
     write(message,'(a,a,a,i5,a,a,a)' )&
&     'For the time being, only nqshft=1',ch10,&
&     'is allowed with brav=4, while it is nqshft=',nqshft,'.',ch10,&
&     'Action : in the input file, correct either brav or nqshft.'
     MSG_ERROR(message)
   end if

   if(nqshft==2)then
!    Make sure that the q vectors form a BCC lattice
     do ii=1,3
       if(abs(abs(qshft(ii,1)-qshft(ii,2))-.5_dp)>1.d-10)then
         write(message, '(a,a,a,a,a,a,a)' )&
&         'The test of the q1shft vectors shows that they',ch10,&
&         'do not generate a body-centered lattice, which',ch10,&
&         'is mandatory for nqshft=2.',ch10,&
&         'Action: change the q1shft vectors in your input file.'
         MSG_ERROR(message)
       end if
     end do
   else if(nqshft==4)then
!    Make sure that the q vectors form a FCC lattice
     do iqshft=1,3
       do jqshft=iqshft+1,4
         tok=0
         do ii=1,3
!          Test on the presence of a +-0.5 difference
           if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-.5_dp) <1.d-10) tok=tok+1
!          Test on the presence of a 0 or +-1.0 difference
           if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-1._dp) <1.d-10  .or.&
&             abs(qshft(ii,iqshft)-qshft(ii,jqshft)) < 1.d-10) tok=tok+4
         end do
!        Test 1 should be satisfied twice, and test 2 once
         if(tok/=6)then
           write(message, '(7a)' )&
&           'The test of the q1shft vectors shows that they',ch10,&
&           'do not generate a face-centered lattice, which',ch10,&
&           'is mandatory for nqshft=4.',ch10,&
&           'Action: change the q1shft vectors in your input file.'
           MSG_ERROR(message)
         end if
       end do
     end do
   else
     write(message, '(a,i4,3a)' )&
&     'nqshft must be 1, 2 or 4. It is nqshft=',nqshft,'.',ch10,&
&     'Action: change nqshft in your input file.'
     MSG_ERROR(message)
   end if

 end if

 factor=0.5_dp
 if(brav==2 .or. brav==3)then
   factor=0.25_dp
 end if
 if(nqshft/=1)factor=factor*2

!DEBUG
!write(std_out,*)'factor,ngqpt',factor,ngqpt(1:3)
!ENDDEBUG


!Begin the big loop on ia and ib
 do ia=1,natom
   do ib=1,natom

!    Simple Lattice
     if (brav==1) then
!      In this case, it is better to work in reduced coordinates
!      As rcan is in canonical coordinates, => multiplication by gprim
       do ii=1,3
         red(1,ii)=  rcan(1,ia)*gprim(1,ii) +rcan(2,ia)*gprim(2,ii) +rcan(3,ia)*gprim(3,ii)
         red(2,ii)=  rcan(1,ib)*gprim(1,ii) +rcan(2,ib)*gprim(2,ii) +rcan(3,ib)*gprim(3,ii)
       end do
     end if

     do irpt=1,nrpt

!      Initialization of the weights to 1.0
       wghatm(ia,ib,irpt)=1.0_dp

!      Compute the difference vector

!      Simple Cubic Lattice
       if (brav==1) then
!        Change of rpt to reduced coordinates
         do ii=1,3
           red(3,ii)=  rpt(1,irpt)*gprim(1,ii) +rpt(2,irpt)*gprim(2,ii) +rpt(3,irpt)*gprim(3,ii)
           rdiff(ii)=red(2,ii)-red(1,ii)+red(3,ii)
         end do
!        Other lattices
       else
         do ii=1,3
           rdiff(ii)=rcan(ii,ib)-rcan(ii,ia)+rpt(ii,irpt)
         end do
!        DEBUG
!        if(ia==1 .and. ib==10)&
!        &     write(std_out,'(a,i5,a,3es16.6)' )' irpt=',irpt,'  rdiff=',rdiff(1:3)
!        ENDDEBUG
       end if

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

!      Assignement of weights

       if(nqshft==1 .and. brav/=4)then

         do ii=1,3
!          If the rpt vector is greater than
!          the allowed space => weight = 0.0
           if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
             wghatm(ia,ib,irpt)=0.0_dp
!            If the point is in a boundary position => weight/2
           else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
!            If the point is in a boundary position => weight/2
             wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
           end if
         end do

!        DEBUG
!        if(ia==1 .and. ib==10)&
!        &     write(std_out,*)' wghatm(ia,ib,irpt)=',wghatm(ia,ib,irpt)
!        ENDDEBUG

       else if(brav==4)then
!        Hexagonal

!        Examination of the X and Y boundaries in order to form an hexagon
!        First generate the relevant boundaries
         rdiff(4)=0.5_dp*( rdiff(1)+sqrt(3.0_dp)*rdiff(2) )
         ngqpt(4)=ngqpt(1)
         rdiff(5)=0.5_dp*( rdiff(1)-sqrt(3.0_dp)*rdiff(2) )
         ngqpt(5)=ngqpt(1)

!        Test the four inequalities
         do ii=1,5
           if(ii/=2)then

             nbord(ii)=0
!            If the rpt vector is greater than
!            the allowed space => weight = 0.0
             if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
               wghatm(ia,ib,irpt)=0.0_dp
             else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
!              If the point is in a boundary position increment nbord(ii)
               nbord(ii)=1
             end if

           end if
         end do

!        Computation of weights
         nbordh=nbord(1)+nbord(4)+nbord(5)
         if (nbordh==1) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
         else if (nbordh==2) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
         else if (nbordh/=0) then
           message = 'There is a problem of borders and weights (hex).'
           MSG_BUG(message)
         end if
         if (nbord(3)==1)then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
         end if

!        BCC packing of k-points
       else if(nqshft==2 .and. brav/=4)then

!        First, generate the relevant boundaries
         rdiff(4)= rdiff(1)+rdiff(2)
         rdiff(5)= rdiff(1)-rdiff(2)
         rdiff(6)= rdiff(1)+rdiff(3)
         rdiff(7)= rdiff(1)-rdiff(3)
         rdiff(8)= rdiff(3)+rdiff(2)
         rdiff(9)= rdiff(3)-rdiff(2)
         if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
           write(message, '(a,a,a,3i6,a,a,a,a)' )&
&           'In the BCC case, the three ngqpt numbers ',ch10,&
&           '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
&           'should be equal.',ch10,&
&           'Action : use identical ngqpt(1:3) in your input file.'
           MSG_ERROR(message)
         end if
         do ii=4,9
           ngqpt(ii)=ngqpt(1)
         end do

!        Test the relevant inequalities
         nbord(1)=0
         do ii=4,9
!          If the rpt vector is greater than the allowed space => weight = 0.0
           if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
             wghatm(ia,ib,irpt)=0.0_dp
           else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
!            If the point is in a boundary position increment nbord(1)
             nbord(1)=nbord(1)+1
           end if
         end do

!        Computation of weights
         if (nbord(1)==1) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
         else if (nbord(1)==2) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
         else if (nbord(1)==3) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
         else if (nbord(1)==4) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/6
         else if (nbord(1)/=0) then
           message = ' There is a problem of borders and weights (BCC).'
           MSG_ERROR(message)
         end if

!        FCC packing of k-points
       else if(nqshft==4 .and. brav/=4)then

!        First, generate the relevant boundaries
         rdiff(4)= (rdiff(1)+rdiff(2)+rdiff(3))*2._dp/3._dp
         rdiff(5)= (rdiff(1)-rdiff(2)+rdiff(3))*2._dp/3._dp
         rdiff(6)= (rdiff(1)+rdiff(2)-rdiff(3))*2._dp/3._dp
         rdiff(7)= (rdiff(1)-rdiff(2)-rdiff(3))*2._dp/3._dp
         if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
           write(message, '(a,a,a,3i6,a,a,a,a)' )&
&           'In the FCC case, the three ngqpt numbers ',ch10,&
&           '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
&           'should be equal.',ch10,&
&           'Action : use identical ngqpt(1:3) in your input file.'
           MSG_ERROR(message)
         end if
         do ii=4,7
           ngqpt(ii)=ngqpt(1)
         end do

!        Test the relevant inequalities
         nbord(1)=0
         do ii=1,7

!          If the rpt vector is greater than the allowed space => weight = 0.0
           if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
             wghatm(ia,ib,irpt)=0.0_dp
!            If the point is in a boundary position increment nbord(1)
           else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
!            DEBUG
!            write(std_out,*)'ia,ib,irpt,ii,nbord(1)'
!            write(std_out,*)ia,ib,irpt,ii,nbord(1)
!            ENDDEBUG
             nbord(1)=nbord(1)+1
           end if
         end do

!        Computation of weights
         if (nbord(1)==1) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
         else if (nbord(1)==2) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
         else if (nbord(1)==3) then
           wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
         else if (nbord(1)/=0 .and. wghatm(ia,ib,irpt)>1.d-10) then
!          Interestingly nbord(1)==4 happens for some points outside of the volume
           message = ' There is a problem of borders and weights (FCC).'
           MSG_BUG(message)
         end if

       else
         write(message, '(3a,i0,a)' )&
&         'One should not arrive here ... ',ch10,&
&         'The value nqshft ',nqshft,' is not available'
         MSG_BUG(message)
       end if
     end do ! Assignement of weights is done

   end do ! End of the double loop on ia and ib
 end do

! Check the results
 do ia=1,natom
   do ib=1,natom
     sum=zero
     do irpt=1,nrpt
!      Check if the sum of the weights is equal to the number of q points
       sum=sum+wghatm(ia,ib,irpt)
!      DEBUG
!      write(std_out,'(a,3i5)' )' atom1, atom2, irpt ; rpt ; wghatm ',ia,ib,irpt
!      write(std_out,'(3es16.6,es18.6)' )&
!      &    rpt(1,irpt),rpt(2,irpt),rpt(3,irpt),wghatm(ia,ib,irpt)
!      ENDDEBUG
     end do
     if (abs(sum-nqpt)>1.0d-10) then
       write(message, '(a,a,a,2i4,a,a,es14.4,a,a,i4,a,a,a,a,a,a)' )&
&       'The sum of the weight is not equal to nqpt.',ch10,&
&       'atoms :',ia,ib,ch10,&
&       'The sum of the weights is : ',sum,ch10,&
&       'The number of q points is : ',nqpt,ch10,&
&       'You might increase "buffer" in bigbx9.f, and recompile.',ch10,&
&       'Actually, this can also happen when ngqpt is 0 0 0,',ch10,&
&       'if brav/=1, in which case you should change brav to 1.'
       MSG_BUG(message)
     end if
   end do
 end do

 DBG_EXIT("COLL")

end subroutine wght9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/d3sym
!! NAME
!! d3sym
!!
!! FUNCTION
!! Given a set of calculated elements of the 3DTE matrix,
!! build (nearly) all the other matrix elements that can be build using symmetries.
!!
!! INPUTS
!!  indsym(4,nsym,natom)=indirect indexing array : for each
!!   isym,iatom, fourth element is label of atom into which iatom is sent by
!!   INVERSE of symmetry operation isym; first three elements are the primitive
!!   translations which must be subtracted after the transformation to get back
!!   to the original unit cell.
!!  mpert =maximum number of ipert
!!  natom= number of atoms
!!  nsym=number of space group symmetries
!!  symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
!!  symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
!!
!! SIDE EFFECTS
!!  Input/Output
!!  blkflg(3,mpert,3,mpert,3,mpert)= matrix that indicates if an
!!   element of d3 is available (1 if available, 0 otherwise)
!!  d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
!!
!! PARENTS
!!      m_ddb_blk,nonlinear
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine d3sym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel)


!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 'd3sym'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: mpert,natom,nsym
!arrays
 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
 real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3
 integer :: ipesy1,ipesy2,ipesy3,isym,ithree
 real(dp) :: sumi,sumr
!arrays
 integer :: sym1(3,3),sym2(3,3),sym3(3,3)

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

!DEBUG
!write(std_out,*)'d3sym : enter'
!do i1dir = 1, 3
!do i2dir = 1, 3
!do i3dir = 1, 3
!write(std_out,*)i1dir,i2dir,i3dir,&
!&   blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
!end do
!end do
!end do
!stop
!ENDDEBUG

!First, take into account the permutations symmetry of
!(i1pert,i1dir) and (i3pert,i3dir)

 do i1pert = 1, mpert
   do i2pert = 1, mpert

     do i3pert = 1, mpert

       do i1dir = 1, 3
         do i2dir = 1, 3
           do i3dir = 1, 3

             if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. &
&             (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=1)) then

               d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = &
&              d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)

               blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = 1

             end if

           end do
         end do
       end do

     end do
   end do
 end do

!Big Big Loop : symmetrize three times, because
!of some cases in which one element is not yet available
!at the first pass, and even at the second one !

 do ithree=1,3

!  Loop over perturbations
   do i1pert = 1, mpert
     do i2pert = 1, mpert
       do i3pert = 1, mpert

         do i1dir = 1, 3
           do i2dir = 1, 3
             do i3dir = 1, 3

!              Will get element (idir1,ipert1,idir2,ipert2)
!              so this element should not yet be present ...
               if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then

                 d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp

                 do isym = 1, nsym

                   found = 1

                   if (i1pert <= natom) then
                     ipesy1 = indsym(4,isym,i1pert)
                     sym1(:,:) = symrec(:,:,isym)
                   else if (i1pert == natom + 2) then
                     ipesy1 = i1pert
                     sym1(:,:) = symrel(:,:,isym)
                   else
                     found = 0
                   end if

                   if (i2pert <= natom) then
                     ipesy2 = indsym(4,isym,i2pert)
                     sym2(:,:) = symrec(:,:,isym)
                   else if (i2pert == natom + 2) then
                     ipesy2 = i2pert
                     sym2(:,:) = symrel(:,:,isym)
                   else
                     found = 0
                   end if

                   if (i3pert <= natom) then
                     ipesy3 = indsym(4,isym,i3pert)
                     sym3(:,:) = symrec(:,:,isym)
                   else if (i3pert == natom + 2) then
                     ipesy3 = i3pert
                     sym3(:,:) = symrel(:,:,isym)
                   else
                     found = 0
                   end if

                   sumr = 0_dp ; sumi = 0_dp;
                   do idisy1 = 1, 3
                     do idisy2 = 1, 3
                       do idisy3 = 1, 3

                         if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) &
&                         .and.(sym3(i3dir,idisy3) /=0)) then

                           if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then

                             sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
&                             sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
                             sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
&                             sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)

                           else

                             found = 0

                           end if

                         end if

                       end do
                     end do
                   end do

                   if (found == 1) then
                     d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
                     d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
                     blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
                   end if

                 end do  ! isym

               end if  ! blkflg

!              Close loop over perturbations
             end do
           end do
         end do
       end do
     end do
   end do

 end do  ! close loop over ithree

end subroutine d3sym
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/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.
!!
!! 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
!! [qmissing]=Allocatable array with the indices of the q-points in the BZ that could not be obtained
!!    by symmetry. If qmissing is present, the routine does not stop if the full BZ cannot be reconstructed.
!!    The caller is responsible for filling the missing entries.
!!
!! TODO
!!   * A full description of the inputs should be included
!!   
!! NOTES
!!   Time-reversal symmetry is always assumed
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

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


!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
!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)
 integer,allocatable,optional,intent(out) :: qmissing(:)
 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,nqmiss
 real(dp),parameter :: tol=2.d-8
 real(dp) :: arg1,arg2,im,re,sumi,sumr
 logical :: allow_qmiss
 character(len=500) :: message
!arrays
 integer,allocatable :: qtest(:,:)
 integer :: qmiss_(nqpt)
 real(dp) :: qq(3),qsym(6),ss(3,3)
 real(dp),allocatable :: ddd(:,:,:,:,:)

! *********************************************************************
 ! Initialize output (some q-points might not be reconstructed if qmissing is present)
 dynmat = huge(one)
 allow_qmiss = (present(qmissing))

 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
! qtest(iqpt,1)=iblok
! qtest(iqpt,2)=isym
! qtest(iqpt,3)=time_reversal
 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 do ! End of the loop on the q points of the sampling
     end do ! End of the loop on the symmetries

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

! Check if all the informations relatives to the q points sampling are found in the DDB; 
! if not => stop message
 nqmiss = 0
 do iqpt=1,nqpt
   if (qtest(iqpt,1)==0) then
     nqmiss = nqmiss + 1
     qmiss_(nqmiss) = iqpt
     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,i6,a,a,a,3es16.6,a,a,a,a)' )&
&     '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.'
     if (.not.allow_qmiss) then
       MSG_ERROR(message)
     else
       !continue
       MSG_COMMENT(message)
     end if
   end if
 end do

 ! Will return a list with the index of the q-points that could not be symmetrized.
 if (allow_qmiss) then
   ABI_ALLOCATE(qmissing, (nqmiss))
   if (nqmiss > 0) qmissing = qmiss_(1:nqmiss) 
 end if

!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)
   ! Skip this q-point if don't have enough info and allow_qmiss
   if (allow_qmiss .and. q1==0) cycle

!  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
!    write(std_out,*)"ss",ss
   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,i6,a,a,a,4i4,a,a,a,a)' )&
&             'Informations are missing in the DDB.',ch10,&
&             'In block',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.'
             MSG_ERROR(message)
           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, '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 do ! 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

 ABI_DEALLOCATE(ddd)
 ABI_DEALLOCATE(qtest)

end subroutine symdm9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/axial9
!!
!! NAME
!! axial9
!!
!! FUNCTION
!! Generates the local coordinates system from the
!! knowledge of the first vector (longitudinal) and
!! the ifc matrix in cartesian coordinates
!!
!! INPUTS
!! ifccar(3,3)= matrix of interatomic force constants in cartesian coordinates
!! vect1(3)= cartesian coordinates of the first local vector
!!
!! OUTPUT
!! vect2(3)= cartesian coordinates of the second local vector
!! vect3(3)= cartesian coordinates of the third local vector
!!
!! PARENTS
!!      m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine axial9(ifccar,vect1,vect2,vect3)


!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 'axial9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!arrays
 real(dp),intent(in) :: ifccar(3,3),vect1(3)
 real(dp),intent(out) :: vect2(3),vect3(3)

!Local variables -------------------------
!scalars
 integer :: flag,ii,itrial,jj
 real(dp) :: innorm,scprod
!arrays
 real(dp) :: work(3)

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

 do jj=1,3
   work(jj)=0.0_dp
   do ii=1,3
     work(jj)=work(jj)+ifccar(jj,ii)*vect1(ii)
   end do
 end do

 flag=0
 do itrial=1,4
   scprod=0.0_dp
   do ii=1,3
     scprod=scprod+work(ii)*vect1(ii)
   end do

   do ii=1,3
     work(ii)=work(ii)-vect1(ii)*scprod
   end do

   scprod=0.0_dp
   do ii=1,3
     scprod=scprod+work(ii)**2
   end do

   if(scprod<1.0d-10)then
     work(1:3)=0.0_dp
     if(itrial>1)work(itrial-1)=1.0_dp
   else
     flag=1
   end if

   if(flag==1)exit
 end do

 innorm=scprod**(-0.5_dp)
 do ii=1,3
   vect2(ii)=work(ii)*innorm
 end do

 vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2)
 vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3)
 vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1)

end subroutine axial9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/dymfz9
!!
!! NAME
!! dymfz9
!!
!! FUNCTION
!! As the subroutine canatm has transformed the coordinates of the
!! atoms in normalized canonical coordinates, the corresponding
!! dynamical matrix should be multiplied by a phase shift corresponding
!! to the translation between New and Old coordinates of its two
!! corresponding atoms.
!!
!! INPUTS
!! dynmat = non-phase shifted dynamical matrices
!! natom = number of atoms
!! nqpt = number of qpoints
!! gprim = reciprocal lattice vectors (cartesian but dimensionless)
!! option=1 : the matrices are transformed from the old (tn)
!!  coordinate system to the new (normalized canonical)
!!        2 : the matrices are restored from the normalized
!!  canonical coordinate system to the usual (tn) one...
!! rcan = canonical coordinates of atoms
!! spqpt = qpoint coordinates (reduced reciprocal)
!! trans = Atomic translations : xred = rcan + trans
!!
!! OUTPUT
!! dynmat = phase shifted dynamical matrices
!!
!! PARENTS
!!      m_dynmat,m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine dymfz9(dynmat,natom,nqpt,gprim,option,spqpt,trans)


!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 'dymfz9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: natom,nqpt,option
!arrays
 real(dp),intent(in) :: gprim(3,3),spqpt(3,nqpt),trans(3,natom)
 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,iqpt,mu,nu
 real(dp) :: im,ktrans,re
!arrays
 real(dp) :: kk(3)

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

 DBG_ENTER("COLL")

 do iqpt=1,nqpt
   !  Definition of q in normalized reciprocal space
   kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
   kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
   kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)

   if(option==1)then
     kk(1)=-kk(1)
     kk(2)=-kk(2)
     kk(3)=-kk(3)
   end if

   do ia=1,natom
     do ib=1,natom
!      Product of q with the differences between the two atomic translations
       ktrans=kk(1)*(trans(1,ia)-trans(1,ib))+kk(2)*(trans(2,ia)-&
&       trans(2,ib))+kk(3)*(trans(3,ia)-trans(3,ib))

       do mu=1,3
         do nu=1,3
           re=dynmat(1,mu,ia,nu,ib,iqpt)
           im=dynmat(2,mu,ia,nu,ib,iqpt)
!          Transformation of the Old dynamical matrices by New ones by multiplication by a phase shift
           dynmat(1,mu,ia,nu,ib,iqpt)=re*cos(two_pi*ktrans)-im*sin(two_pi*ktrans)
           dynmat(2,mu,ia,nu,ib,iqpt)=re*sin(two_pi*ktrans)+im*cos(two_pi*ktrans)
         end do
       end do
     end do
   end do
 end do

 DBG_EXIT("COLL")

end subroutine dymfz9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/nanal9
!!
!! NAME
!! nanal9
!!
!! FUNCTION
!! If plus=0 then substracts the non-analytical part from one dynamical
!!           matrices, with number iqpt.
!! If plus=1 then adds the non-analytical part to the dynamical
!!           matrices, with number iqpt.
!!
!! INPUTS
!! dyew(2,3,natom,3,natom)= Non-analytical part
!! natom= Number of atoms in the unit cell
!! iqpt= Referenced q point for the dynamical matrix
!! nqpt= Number of q points
!! plus= (see above)
!!
!! OUTPUT
!! dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
!!
!! PARENTS
!!      m_dynmat,m_ifc
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine nanal9(dyew,dynmat,iqpt,natom,nqpt,plus)


!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 'nanal9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: iqpt,natom,nqpt,plus
!arrays
 real(dp),intent(in) :: dyew(2,3,natom,3,natom)
 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)

!Local variables -------------------------
!scalars
 integer :: ia,ib,mu,nu
 character(len=500) :: message

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

 if (plus==0) then

   do ia=1,natom
     do ib=1,natom
       do mu=1,3
         do nu=1,3
!          The following four lines are OK
           dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)&
&           -dyew(1,mu,ia,nu,ib)
           dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)&
&           -dyew(2,mu,ia,nu,ib)
         end do
       end do
     end do
   end do

 else if (plus==1) then
!  write(std_out,*)' nanal9 : now, only the analytic part '
!  write(std_out,*)' nanal9 : now, only the ewald part '
   do ia=1,natom
     do ib=1,natom
       do mu=1,3
         do nu=1,3
!          The following four lines arethe good ones
           dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)&
&           +dyew(1,mu,ia,nu,ib)
           dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)&
&           +dyew(2,mu,ia,nu,ib)
!          DEBUG
!          dynmat(1,mu,ia,nu,ib,iqpt)=dyew(1,mu,ia,nu,ib)
!          dynmat(2,mu,ia,nu,ib,iqpt)=dyew(2,mu,ia,nu,ib)
!          dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)
!          dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)
!          ENDDEBUG
         end do
       end do
     end do
   end do

 else
   write(message,'(a,a,a,i0,a)' )&
&   'The argument "plus" must be equal to 0 or 1.',ch10,&
&   'The value ',plus,' is not available.'
   MSG_BUG(message)
 end if

end subroutine nanal9
!!***

!----------------------------------------------------------------------

!!****f* m_dynmat/gtdyn9
!!
!! NAME
!! gtdyn9
!!
!! FUNCTION
!! Generates a dynamical matrix from interatomic force
!! constants and long-range electrostatic interactions.
!!
!! INPUTS
!! acell(3)=length scales by which rprim is to be multiplied
!! atmfrc(2,3,natom,3,natom,nrpt) = Interatomic Forces in real space
!!  (imaginary part only for debugging)
!! dielt(3,3) = dielectric tensor
!! dipdip= if 0, no dipole-dipole interaction was subtracted in atmfrc
!!  if 1, atmfrc has been build without dipole-dipole part
!! dyewq0(3,3,natom)= Ewald part of the dynamical matrix, at q=0.
!! gmet(3,3)= metric tensor in reciprocal space.
!! gprim(3,3)= Normalized coordinates in reciprocal space
!! mpert =maximum number of ipert
!! natom= Number of atoms in the unit cell
!! nrpt= Number of R points in the Big Box
!! qphnrm= Normalisation coefficient for qpt
!! qpt(3)= Reduced coordinates of the q vectors in reciprocal space
!! rmet(3,3)= Metric tensor in real space.
!! rprim(3,3)= dimensionless primitive translations in real space
!! rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
!!  These coordinates are normalized (=> * acell(3)!!)
!! trans(3,natom)= Atomic translations : xred = rcan + trans
!! ucvol= unit cell volume
!! wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
!! xred(3,natom)= relative coords of atoms in unit cell (dimensionless)
!! zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
!!
!! OUTPUT
!! d2cart(2,3,mpert,3,mpert)=dynamical matrix obtained for the wavevector qpt (normalized using qphnrm)
!!
!! PARENTS
!!      anaddb,inpphon,interpolate_phfrq,m_ifc,m_phdos,mkphbs,thm9
!!
!! CHILDREN
!!      dymfz9,ewald9,ftifc_r2q,nanal9,q0dy3_apply
!!
!! SOURCE

subroutine gtdyn9(acell,atmfrc,dielt,dipdip,&
& dyewq0,d2cart,gmet,gprim,mpert,natom,&
& nrpt,qphnrm,qpt,rmet,rprim,rpt,&
& trans,ucvol,wghatm,xred,zeff)


!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 'gtdyn9'
!End of the abilint section

 implicit none

!Arguments -------------------------------
!scalars
 integer,intent(in) :: dipdip,mpert,natom,nrpt
 real(dp),intent(in) :: qphnrm,ucvol
!arrays
 real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qpt(3)
 real(dp),intent(in) :: rmet(3,3),rprim(3,3),rpt(3,nrpt)
 real(dp),intent(in) :: trans(3,natom),wghatm(natom,natom,nrpt),xred(3,natom)
 real(dp),intent(in) :: zeff(3,3,natom)
 real(dp),intent(in) :: atmfrc(2,3,natom,3,natom,nrpt)
 real(dp),intent(in) :: dyewq0(3,3,natom)
 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)

!Local variables -------------------------
!scalars
 integer :: i1,i2,ib,iqpt,nqpt,nsize,option,plus
 integer :: sumg0
 !character(len=500) :: message
!arrays
 real(dp) :: qphon(3)
 real(dp),allocatable :: dq(:,:,:,:,:),dyew(:,:,:,:,:)

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

 !write(message, '(a,3es15.5)' )'  gtdyn9 : enter ftiaf9 with q =',qphon(1:3)
 !call wrtout(std_out,message,'COLL')

 ABI_ALLOCATE(dq,(2,3,natom,3,natom))

!Get the normalized wavevector
!write(std_out,*)' gtdyn9 : wavevector entirely blocked on zero'
 if(abs(qphnrm)<1.0d-7)then
   qphon(1:3)=zero
 else
   qphon(1:3)=qpt(1:3)/qphnrm
 end if

!Generate the analytical part from the interatomic forces
 nqpt=1
 call ftifc_r2q (atmfrc,dq,gprim,natom,nqpt,nrpt,rpt,qphon,wghatm)

!The analytical dynamical matrix dq has been generated
!in the normalized canonical coordinate system. Now, the
!phase is modified, in order to recover the usual (xred)
!coordinate of atoms.

 option=2; nqpt=1
 call dymfz9(dq,natom,nqpt,gprim,option,qphon,trans)

 if (dipdip==1) then
!  Add the non-analytical part
!  Compute dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
!  second energy derivative wrt xred(3,natom) in Hartrees
!  (Denoted A-bar in the notes)
   ABI_ALLOCATE(dyew,(2,3,natom,3,natom))

   sumg0=0
   call ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff)
   call q0dy3_apply(natom,dyewq0,dyew)

   plus=1; iqpt=1
   call nanal9(dyew,dq,iqpt,natom,nqpt,plus)
   ABI_DEALLOCATE(dyew)
 end if

!Copy the dynamical matrix in the proper location

!First zero all the elements
 nsize=2*(3*mpert)**2
 d2cart(:,:,:,:,:)=zero

!Copy the elements from dq to d2cart
 d2cart(:,:,1:natom,:,1:natom)=dq(:,:,1:natom,:,1:natom)

!In case we have the gamma point,
 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)then
!  Copy the effective charge and dielectric constant in the final array
   do i1=1,3
     do i2=1,3
       d2cart(1,i1,natom+2,i2,natom+2)=dielt(i1,i2)
       do ib=1,natom
         d2cart(1,i1,natom+2,i2,ib)=zeff(i1,i2,ib)
         d2cart(1,i2,ib,i1,natom+2)=zeff(i1,i2,ib)
       end do
     end do
   end do
 end if

 ABI_DEALLOCATE(dq)

end subroutine gtdyn9
!!***

!----------------------------------------------------------------------

end module m_dynmat
