TABLE OF CONTENTS


ABINIT/m_dynmat [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dynmat

FUNCTION

  This module provides low-level tools to operate on the dynamical matrix

COPYRIGHT

  Copyright (C) 2014-2018 ABINIT group (XG, JCC, MJV, NH, RC, MVeithen, MM, MG, MT, DCA)
  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

 24 #if defined HAVE_CONFIG_H
 25 #include "config.h"
 26 #endif
 27 
 28 #include "abi_common.h"
 29 
 30 module m_dynmat
 31 
 32  use defs_basis
 33  use m_abicore
 34  use m_errors
 35  use m_linalg_interfaces
 36 
 37  use m_fstrings,        only : itoa, sjoin
 38  use m_numeric_tools,   only : wrap2_pmhalf, mkherm
 39  use m_symtk,           only : mati3inv, matr3inv, littlegroup_q
 40  use m_cgtools,         only : fxphas_seq
 41  use m_ewald,           only : ewald9
 42 
 43  implicit none
 44 
 45  private
 46 
 47  public :: asria_calc           ! Calculate the correction for the Acoustic sum rule on
 48                                 !   the InterAtomic Forces or on the dynamical matrix directly
 49  public :: asria_corr           ! Imposition of the Acoustic sum rule on the InterAtomic Forces
 50                                 !   or on the dynamical matrix directly from the previously calculated d2asr
 51  public :: asrprs               ! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry
 52  public :: cart29               ! Transform a second-derivative matrix from reduced  coordinates to cartesian coordinates, and also
 53                                 !   1) add the ionic part of the effective charges,
 54                                 !   2) normalize the electronic dielectric tensor, and add the vacuum polarisation
 55  public :: cart39               ! Transform a vector from reduced coordinates to cartesian coordinates,
 56                                 !   taking into account the perturbation from which it was derived,
 57                                 !   and also check the existence of the new values.
 58  public :: d2cart_to_red        ! Transform a second-derivative matrix
 59                                 ! from cartesian to reduced coordinate.
 60  public :: chkph3               ! Check the completeness of the dynamical matrix
 61  public :: chneu9               ! Imposition of the Acoustic sum rule on the Effective charges
 62  public :: d2sym3               ! Build (nearly) all the other matrix elements that can be build using symmetries.
 63  public :: q0dy3_apply          ! Takes care of the inclusion of the ewald q=0 term in the dynamical matrix
 64  public :: q0dy3_calc           ! Calculate the q=0 correction term to the dynamical matrix
 65  ! TODO: 3 routines to symmetrize. Clarify different use cases
 66  public :: symdyma              ! Symmetrize the dynamical matrices
 67  public :: dfpt_sygra           ! Symmetrize derivatives of energy with respect to coordinates,
 68  public :: dfpt_sydy            ! Symmetrize dynamical matrix (eventually diagonal wrt to the atoms)
 69  public :: wings3               ! Suppress the wings of the cartesian 2DTE for which the diagonal element is not known
 70  public :: asrif9               ! Imposes the Acoustic Sum Rule to Interatomic Forces
 71  public :: make_bigbox          ! Generates a Big Box of R points for the Fourier Transforms the dynamical matrix
 72  public :: bigbx9               ! Helper functions that faciliates the generation  of a Big Box containing
 73  public :: canat9               ! From reduced to canonical coordinates
 74  public :: canct9               ! Convert from canonical coordinates to cartesian coordinates
 75  public :: chkrp9               ! Check if the rprim used for the definition of the unit cell (in the
 76                                 ! inputs) are consistent with the rprim used in the routine generating  the Big Box
 77  public :: dist9                ! Compute the distance between atoms in the big box
 78  public :: ftifc_q2r            ! Fourier transform of the dynamical matrices to obtain interatomic forces (real space).
 79  public :: ftifc_r2q            ! Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space).
 80  public :: dynmat_dq            ! Compute the derivative D(q)/dq via Fourier transform of the interatomic forces
 81  public :: ifclo9               ! Convert from cartesian coordinates to local coordinates
 82  public :: wght9                ! Generates a weight to each R points of the Big Box and for each pair of atoms
 83  public :: d3sym                ! Given a set of calculated elements of the 3DTE matrix,
 84                                 ! build (nearly) all the other matrix elements that can be build using symmetries.
 85  public :: sytens               ! Determines the set of irreductible elements of the non-linear optical susceptibility
 86                                 ! and Raman tensors
 87  public :: symdm9               ! Generate the dynamical matrix in the full BZ from the irreducible q-points.
 88  public :: axial9               ! Generates the local coordinates system from the  knowledge of the first vector (longitudinal) and
 89                                 !   the ifc matrix in cartesian coordinates
 90  public :: dymfz9               ! Multiply the dynamical matrix by a phase shift to account for normalized canonical coordinates.
 91  public :: nanal9               ! Subtract/Add the non-analytical part from one dynamical matrix with number iqpt.
 92  public :: gtdyn9               ! Generates a dynamical matrix from interatomic force constants and
 93                                 ! long-range electrostatic interactions.
 94  public :: dfpt_phfrq           ! Diagonalize IFC(q), return phonon frequencies and eigenvectors.
 95                                 ! If q is Gamma, the non-analytical behaviour can be included.
 96  public :: dfpt_prtph           ! Print phonon frequencies
 97  public :: massmult_and_breaksym  ! Multiply IFC(q) by atomic masses.
 98 
 99  ! TODO: Change name,
100  public :: ftgam
101  public :: ftgam_init
102 
103 
104 ! *************************************************************************
105 
106 contains

m_dynmat/asria_calc [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      dfpt_gatherdy,m_ddb,m_effective_potential_file

CHILDREN

SOURCE

135 subroutine asria_calc(asr,d2asr,d2cart,mpert,natom)
136 
137 
138 !This section has been created automatically by the script Abilint (TD).
139 !Do not modify the following lines by hand.
140 #undef ABI_FUNC
141 #define ABI_FUNC 'asria_calc'
142 !End of the abilint section
143 
144  implicit none
145 
146 !Arguments -------------------------------
147 !scalars
148  integer,intent(in) :: asr,mpert,natom
149 !arrays
150  real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
151  real(dp),intent(out) :: d2asr(2,3,natom,3,natom)
152 
153 !Local variables-------------------------------
154 !scalars
155  integer :: idir1,idir2,ii,ipert1,ipert2
156  integer :: constrank, imatelem, iconst, nconst, nd2_packed, info
157  !character(len=500) :: message
158 !arrays
159  integer, allocatable :: packingindex(:,:,:,:)
160  real(dp), allocatable :: constraints(:,:,:)
161  real(dp), allocatable :: d2cart_packed(:,:)
162  real(dp), allocatable :: singvals(:)
163  real(dp), allocatable :: constr_rhs(:,:)
164  real(dp), allocatable :: work(:,:),rwork(:)
165 
166 ! *********************************************************************
167 
168  d2asr = zero
169 
170  if (asr==0) return
171 
172  !call wrtout(std_out,' asria_calc: calculation of the correction to the ASR for the interatomic forces.')
173  do ipert1=1,natom
174    do idir1=1,3
175      do idir2=1,3
176 
177 !      Compute d2asr
178        do ipert2=1,natom
179          d2asr(:,idir1,ipert1,idir2,ipert1)=&
180 &         d2asr(:,idir1,ipert1,idir2,ipert1)+&
181 &         d2cart(:,idir1,ipert1,idir2,ipert2)
182        end do
183      end do
184    end do
185  end do
186 
187 !holistic method: overwrite d2asr with hermitian solution
188  if (asr == 5) then
189    nconst = 9*natom
190    nd2_packed = 3*natom*(3*natom+1)/2
191    ABI_ALLOCATE(constraints,(2,nconst, nd2_packed))
192    ABI_ALLOCATE(d2cart_packed,(2,nd2_packed))
193    ABI_ALLOCATE(constr_rhs,(2,nd2_packed))
194    ABI_ALLOCATE(singvals,(nconst))
195    ABI_ALLOCATE(work,(2,3*nd2_packed))
196    ABI_ALLOCATE(rwork,(5*nd2_packed))
197    ABI_ALLOCATE(packingindex,(3,natom,3,natom))
198    ii=1
199    packingindex=-1
200    do ipert2=1,natom
201      do idir2=1,3
202        do ipert1=1,ipert2-1
203          do idir1=1,3
204            packingindex(idir1,ipert1,idir2,ipert2) = ii
205            ii = ii+1
206          end do
207        end do
208        do idir1=1,idir2
209          packingindex(idir1,ipert2,idir2,ipert2) = ii
210          ii = ii+1
211        end do
212      end do
213    end do
214 !  setup constraint matrix
215    constraints = zero
216    do ipert1=1,natom
217      do idir1=1,3
218        do idir2=1,3
219          iconst = idir2+3*(idir1-1 + 3*(ipert1-1))
220 !        set all atom forces, this component
221          do ipert2=1,natom
222            imatelem = packingindex(idir1,ipert1,idir2,ipert2)
223            if (imatelem == -1) then
224              imatelem = packingindex(idir2,ipert2,idir1,ipert1)
225            end if
226            constraints(1,iconst,imatelem) = one
227          end do
228        end do
229      end do
230    end do
231 
232    d2cart_packed = -999.0d0
233    do ipert2=1,natom
234      do idir2=1,3
235        do ipert1=1,natom
236          do idir1=1,3
237            imatelem = packingindex(idir1,ipert1,idir2,ipert2)
238            if (imatelem == -1) cycle
239            d2cart_packed(:,imatelem) = d2cart(:,idir1,ipert1,idir2,ipert2)
240          end do
241        end do
242      end do
243    end do
244    constr_rhs = zero
245    constr_rhs(1,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(1,:))
246    constr_rhs(2,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(2,:))
247 
248 !  lwork = 3*nd2_packed
249    call zgelss (nconst,nd2_packed,1,constraints,nconst,constr_rhs,nd2_packed,&
250 &   singvals,-one,constrank,work,3*nd2_packed,rwork,info)
251    ABI_CHECK(info == 0, sjoin('zgelss returned:', itoa(info)))
252 
253 !  unpack
254    do ipert2=1,natom
255      do idir2=1,3
256        do ipert1=1,natom
257          do idir1=1,3
258            imatelem = packingindex(idir1,ipert1,idir2,ipert2)
259            if (imatelem == -1) then
260              imatelem = packingindex(idir2,ipert2,idir1,ipert1)
261 !            NOTE: should complex conjugate the correction below.
262            end if
263            d2asr(:,idir1,ipert1,idir2,ipert2) = constr_rhs(:,imatelem)
264          end do
265        end do
266      end do
267    end do
268 
269    ABI_DEALLOCATE(constraints)
270    ABI_DEALLOCATE(d2cart_packed)
271    ABI_DEALLOCATE(singvals)
272    ABI_DEALLOCATE(constr_rhs)
273    ABI_DEALLOCATE(work)
274    ABI_DEALLOCATE(rwork)
275    ABI_DEALLOCATE(packingindex)
276  end if
277 
278 end subroutine asria_calc

m_dynmat/asria_corr [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

PARENTS

      ddb_elast,ddb_internalstr,dfpt_gatherdy,m_ddb,thmeig

CHILDREN

SOURCE

310 subroutine asria_corr(asr,d2asr,d2cart,mpert,natom)
311 
312 
313 !This section has been created automatically by the script Abilint (TD).
314 !Do not modify the following lines by hand.
315 #undef ABI_FUNC
316 #define ABI_FUNC 'asria_corr'
317 !End of the abilint section
318 
319  implicit none
320 
321 !Arguments -------------------------------
322 !scalars
323  integer,intent(in) :: asr,mpert,natom
324 !arrays
325  real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
326  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
327 
328 !Local variables-------------------------------
329 !scalars
330  integer :: idir1,idir2,ipert1,ipert2
331 
332 ! *********************************************************************
333 
334  if (asr==0) return
335  !call wrtout(std_out,' asria_corr: imposition of the ASR for the interatomic forces.','COLL')
336 
337 ! Remove d2asr
338  do ipert2=1,natom
339    do idir2=1,3
340      do ipert1=1,natom
341        do idir1=1,3
342          d2cart(:,idir1,ipert1,idir2,ipert2)= d2cart(:,idir1,ipert1,idir2,ipert2) - d2asr(:,idir1,ipert1,idir2,ipert2)
343        end do
344      end do
345    end do
346  end do
347 
348 end subroutine asria_corr

m_dynmat/asrif9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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(3,natom,3,natom,nrpt)= Interatomic Forces

OUTPUT

 atmfrc(3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces

TODO

 List of ouput should be included.

PARENTS

      ddb_hybrid,m_ifc,m_tdep_abitypes

CHILDREN

SOURCE

2794 subroutine asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm)
2795 
2796 
2797 !This section has been created automatically by the script Abilint (TD).
2798 !Do not modify the following lines by hand.
2799 #undef ABI_FUNC
2800 #define ABI_FUNC 'asrif9'
2801 !End of the abilint section
2802 
2803  implicit none
2804 
2805 !Arguments -------------------------------
2806 !scalars
2807  integer,intent(in) :: asr,natom,nrpt
2808 !arrays
2809  real(dp),intent(in) :: rpt(3,nrpt),wghatm(natom,natom,nrpt)
2810  real(dp),intent(inout) :: atmfrc(3,natom,3,natom,nrpt)
2811 
2812 !Local variables -------------------------
2813 !scalars
2814  integer :: found,ia,ib,irpt,izero,mu,nu
2815  real(dp) :: sumifc
2816 
2817 ! *********************************************************************
2818 
2819  DBG_ENTER("COLL")
2820 
2821  if(asr==1.or.asr==2)then
2822 
2823    found=0
2824 !  Search for the R vector which is equal to ( 0 , 0 , 0 )
2825 !  This vector leaves the atom a on itself !
2826    do irpt=1,nrpt
2827      if (all(abs(rpt(:,irpt))<=1.0d-10)) then
2828        found=1
2829        izero=irpt
2830      end if
2831      if (found==1) exit
2832    end do
2833 
2834    if(found==0)then
2835      MSG_BUG('Not able to find the vector R=(0,0,0).')
2836    end if
2837 
2838    do mu=1,3
2839      do nu=1,3
2840        do ia=1,natom
2841          sumifc=zero
2842          do ib=1,natom
2843 
2844 !          Get the sumifc of interatomic forces acting on the atom ia,
2845 !          either in a symmetrical manner, or an unsymmetrical one.
2846            if(asr==1)then
2847              do irpt=1,nrpt
2848                sumifc=sumifc+wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)
2849              end do
2850            else if(asr==2)then
2851              do irpt=1,nrpt
2852                sumifc=sumifc+&
2853 &               (wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)+&
2854 &               wghatm(ia,ib,irpt)*atmfrc(nu,ia,mu,ib,irpt))/2
2855              end do
2856            end if
2857          end do
2858 
2859 !        Correct the self-interaction in order to fulfill the ASR
2860          atmfrc(mu,ia,nu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)-sumifc
2861          if(asr==2)then
2862            atmfrc(nu,ia,mu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)
2863          end if
2864 
2865        end do
2866      end do
2867    end do
2868  end if
2869 
2870  DBG_EXIT("COLL")
2871 
2872 end subroutine asrif9

m_dynmat/asrprs [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      m_ddb

CHILDREN

SOURCE

390 subroutine asrprs(asr,asrflag,rotinv,uinvers,vtinvers,singular,d2cart,mpert,natom,xcart)
391 
392 
393 !This section has been created automatically by the script Abilint (TD).
394 !Do not modify the following lines by hand.
395 #undef ABI_FUNC
396 #define ABI_FUNC 'asrprs'
397 !End of the abilint section
398 
399  implicit none
400 
401 !Arguments ------------------------------------
402 !scalars
403  integer,intent(in) :: asr,asrflag,mpert,natom,rotinv
404 !arrays
405  real(dp),intent(in) :: xcart(3,natom)
406  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
407  real(dp),intent(inout) :: singular(1:3*natom*(3*natom-1)/2)
408  real(dp),intent(inout) :: uinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)
409  real(dp),intent(inout) :: vtinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)
410 
411 !Local variables-------------------------------
412 !scalars
413  integer :: column,idir1,idir2,ii,info,ipert1,ipert2,jj,n3,row,superdim
414  real(dp) :: rcond,test
415 ! real(dp) :: tau ! tau is present but commented out in this routine
416  character(len=500) :: message
417 !arrays
418  integer :: check(3,natom,3)
419  real(dp) :: tmp(natom,3,3),weightf(1:natom,1:natom)
420  real(dp),allocatable :: d2cartold(:,:,:,:,:),d2vecc(:),d2veccnew(:),d2vecr(:)
421  real(dp),allocatable :: d2vecrnew(:),superm(:,:),umatrix(:,:),vtmatrix(:)
422  real(dp),allocatable :: work(:)
423 
424 ! *********************************************************************
425 
426  if(asr/=3 .and. asr/=4)then
427    write(message,'(3a,i0)')&
428 &   'The argument asr should be 3 or 4,',ch10,&
429 &   'however, asr = ',asr
430    MSG_BUG(message)
431  end if
432 
433  if (asr==3.or.asr==4)then
434    write(message, '(a,a)' ) ch10, &
435 &   'asrprs: imposition of the ASR for the interatomic forces and rotational invariance'
436    call wrtout(std_out,message,'COLL')
437  end if
438 
439  write(message,'(a,i6)')' asrflag is', asrflag
440  call wrtout(std_out,message,'COLL')
441  call wrtout(ab_out,message,'COLL')
442 
443 !variables for the dimensions of the matrices
444 
445 !n1=3*natom*(3*natom-1)/2
446 !n2=9*natom
447  n3=3*natom
448 
449  superdim=9*natom*(natom-1)/2+n3
450 
451  ABI_ALLOCATE(d2vecr,(1:superdim))
452  ABI_ALLOCATE(d2vecc,(1:superdim))
453  d2vecr=0d0
454  d2vecc=0d0
455 
456 !should be changed set to delta function for debugging
457  weightf=1d0
458 !tau=1d-10
459  do ii=1, natom
460 !  do jj=1, ii-1
461 !  weightf(ii,jj)= &
462 !  &     ((xcart(1,ii)-xcart(1,jj))**2+(xcart(2,ii)-xcart(2,jj))**2+(xcart(3,ii)-xcart(3,jj))**2)**tau
463 !  enddo
464    weightf(ii,ii)=0d0
465  end do
466 
467  ABI_ALLOCATE(d2cartold,(2,3,mpert,3,mpert))
468 
469  d2cartold=d2cart
470 
471 !setup vector with uncorrected derivatives
472 
473  do ipert1=1, natom
474    do ipert2=1, ipert1-1
475      do idir1=1,3
476        do idir2=1,3
477          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
478          if(abs(d2cart(1,idir1,ipert1,idir2,ipert2))<1d-6)then
479            d2cart(1,idir1,ipert1,idir2,ipert2)=0d0
480          else
481            d2vecr(row)=4*weightf(ipert1,ipert2)*d2cart(1,idir1,ipert1,idir2,ipert2)
482          end if
483          if(abs(d2cart(2,idir1,ipert1,idir2,ipert2))<1d-6) then
484            d2cart(2,idir1,ipert1,idir2,ipert2)=0d0
485          else
486            d2vecc(row)=4*weightf(ipert1,ipert2)*d2cart(2,idir1,ipert1,idir2,ipert2)
487          end if
488        end do
489      end do
490    end do
491  end do
492 
493  if(asrflag==1) then !calculate the pseudo-inverse of the supermatrix
494    ABI_ALLOCATE(superm,(1:superdim,1:superdim))
495 
496    superm=0d0
497 
498 !  Setting up the supermatrix containing G, A, D
499 
500    do ipert1=1, natom
501      do idir1=1, 3
502 !      Setting up G
503        idir2=mod(idir1,3)+1
504        row=3*(ipert1-1)+idir1
505        do ipert2=1, ipert1-1
506          column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
507          superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
508          column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir2
509          superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
510        end do
511        do ipert2=ipert1+1, natom
512          column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir1-1)+rotinv
513          superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
514          column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir2-1)+rotinv
515          superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
516        end do
517      end do
518      do idir1=1, 3
519 !      Setting up D
520        idir2=mod(idir1,3)+1
521        ii=mod(idir1+1,3)+1
522        do ipert2=1, ipert1-1
523          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
524          column=9*natom*(natom-1)/2+3*(ipert1-1)+idir1
525          superm(column,row)=superm(column,row)+xcart(idir2,ipert2)-xcart(idir2,ipert1)
526          column=9*natom*(natom-1)/2+3*(ipert1-1)+ii
527          superm(column,row)=superm(column,row)+xcart(ii,ipert1)-xcart(ii,ipert2)
528          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+rotinv
529          column=9*natom*(natom-1)/2+3*(ipert2-1)+idir1
530          superm(column,row)=superm(column,row)+xcart(idir2,ipert1)-xcart(idir2,ipert2)
531          column=9*natom*(natom-1)/2+3*(ipert2-1)+ii
532          superm(column,row)=superm(column,row)+xcart(ii,ipert2)-xcart(ii,ipert1)
533        end do
534 !      Setting up A
535        do idir2=1, 3
536          do ipert2=1, ipert1-1
537            column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
538            row=n3+column
539            superm(column,row)=4*weightf(ipert1,ipert2)
540          end do
541        end do
542      end do
543    end do
544 
545 !  calculate the pseudo-inverse of the supermatrix
546 
547    ABI_ALLOCATE(work,(1:6*superdim))
548    ABI_ALLOCATE(vtmatrix,(1:superdim))
549    ABI_ALLOCATE(umatrix,(1:superdim,1:superdim))
550 
551 !  singular value decomposition of superm
552 
553    call dgesvd('A','O',superdim,superdim,superm,superdim,singular,umatrix,superdim, &
554 &   vtmatrix, 1, work,6*superdim,info)
555    ABI_CHECK(info == 0, sjoin('dgesvd returned:', itoa(info)))
556 
557    ABI_DEALLOCATE(vtmatrix)
558    ABI_DEALLOCATE(work)
559 
560    write(message, '(a,es16.8,es16.8)' )&
561 &   ' Largest and smallest values from svd', singular(1), singular(superdim)
562    call wrtout(std_out,message,'COLL')
563    call wrtout(ab_out,message,'COLL')
564 
565 !  Invert U and V**T, orthogonal matrices
566 
567    do ii=1, superdim
568      do jj=1, superdim
569        uinvers(ii,jj)=umatrix(jj,ii)
570        vtinvers(ii,jj)=superm(jj,ii)
571      end do
572    end do
573 
574    ABI_DEALLOCATE(umatrix)
575    ABI_DEALLOCATE(superm)
576 
577    write(message,'(a,a)')' asrprs: done with asrflag 1', ch10
578    call wrtout(std_out,message,'COLL')
579 
580  end if !asrflag=1
581 
582  if(asrflag==2) then
583 
584    ABI_ALLOCATE(d2vecrnew,(1:superdim))
585    ABI_ALLOCATE(d2veccnew,(1:superdim))
586 
587 !  Calculate V**T**-1 Sigma**-1 U**-1 *rhs
588 
589    do ii=1, superdim
590      d2vecrnew(ii)=0d0
591      d2veccnew(ii)=0d0
592      do jj=1, superdim
593        d2vecrnew(ii)=d2vecrnew(ii)+uinvers(ii,jj)*d2vecr(jj)
594        d2veccnew(ii)=d2veccnew(ii)+uinvers(ii,jj)*d2vecc(jj)
595      end do
596    end do
597 
598    rcond=1d-10*singular(1)
599    do ii=1, superdim
600      if(singular(ii)>rcond) then
601        d2vecrnew(ii)=d2vecrnew(ii)/singular(ii)
602        d2veccnew(ii)=d2veccnew(ii)/singular(ii)
603      else
604        d2vecrnew(ii)=0d0
605        d2veccnew(ii)=0d0
606      end if
607    end do
608 
609    do ii=1, superdim
610      d2vecr(ii)=0d0
611      d2vecc(ii)=0d0
612      do jj=1, superdim
613        d2vecr(ii)=d2vecr(ii)+vtinvers(ii,jj)*d2vecrnew(jj)
614        d2vecc(ii)=d2vecc(ii)+vtinvers(ii,jj)*d2veccnew(jj)
615      end do
616    end do
617 
618 !  Store vector back into the matrix of 2nd order derivates
619 
620    do ipert1=1, natom
621      do ipert2=1, ipert1-1
622        do idir1=1,3
623          do idir2=1,3
624            row=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
625            d2cart(1,idir1,ipert1,idir2,ipert2)=d2vecr(row)
626            d2cart(2,idir1,ipert1,idir2,ipert2)=d2vecc(row)
627            d2cart(1,idir2,ipert2,idir1,ipert1)=d2vecr(row)
628            d2cart(2,idir2,ipert2,idir1,ipert1)=d2vecc(row)
629          end do
630        end do
631      end do
632    end do
633 
634    ABI_DEALLOCATE(d2vecrnew)
635    ABI_DEALLOCATE(d2veccnew)
636 
637    check=0
638 
639    do ipert1=1, natom
640      do idir1=1, 3
641        do idir2=1, 3
642          d2cart(1,idir1,ipert1,idir2,ipert1)=0d0
643          d2cart(2,idir1,ipert1,idir2,ipert1)=0d0
644          tmp(ipert1,idir1,idir2)=0d0
645          do ipert2=1, natom
646            if(ipert2/=ipert1) then
647              tmp(ipert1,idir1,idir2)=tmp(ipert1,idir1,idir2) &
648 &             -d2cart(1,idir1,ipert1,idir2,ipert2) &
649 &             -d2cart(1,idir2,ipert2,idir1,ipert1)
650            end if
651          end do
652        end do
653      end do
654    end do
655 
656    do ipert1=1, natom
657      do idir1=1, 3
658        do idir2=1, 3
659          d2cart(1,idir1,ipert1,idir2,ipert1)=tmp(ipert1,idir1,idir2)/2
660          d2cart(1,idir2,ipert1,idir1,ipert1)=d2cart(1,idir1,ipert1,idir2,ipert1)
661        end do
662      end do
663    end do
664 
665    write(std_out,*) 'this should all be zero'
666 
667    do ipert1=1, natom
668      do idir1=1, 3
669        do idir2=1, 3
670          test=0d0
671          do ipert2=1, natom
672            test=test+d2cart(1,idir1,ipert1,idir2,ipert2)+d2cart(1,idir2,ipert2,idir1,ipert1)
673          end do
674          write(std_out,'(i3,i3,i3,es11.3)') idir1,ipert1,idir2,test
675 
676          write(message, '(i3,i3,i3,es11.3)' ) idir1,ipert1,idir2,test
677          call wrtout(ab_out,message,'COLL')
678 
679        end do
680      end do
681    end do
682 
683    write(std_out,*) 'these as well'
684    do ipert2=1, natom
685      do idir1=1, 3
686        do idir2=1, 3
687          test=0d0
688          do ipert1=1, natom
689            test=test+d2cart(1,idir1,ipert1,idir2,ipert2)
690          end do
691          write(std_out,'(i3,i3,i3,i3,es11.3)') idir1,ipert1,idir2,ipert2,test
692        end do
693      end do
694    end do
695 
696    write(message,'(a,a)')' asrprs: done with asrflag 2', ch10
697    call wrtout(std_out,message,'COLL')
698 
699  end if !ends asrflag=2
700 
701  ABI_DEALLOCATE(d2vecr)
702  ABI_DEALLOCATE(d2vecc)
703  ABI_DEALLOCATE(d2cartold)
704 
705 end subroutine asrprs

m_dynmat/axial9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

5367 subroutine axial9(ifccar,vect1,vect2,vect3)
5368 
5369 
5370 !This section has been created automatically by the script Abilint (TD).
5371 !Do not modify the following lines by hand.
5372 #undef ABI_FUNC
5373 #define ABI_FUNC 'axial9'
5374 !End of the abilint section
5375 
5376  implicit none
5377 
5378 !Arguments -------------------------------
5379 !arrays
5380  real(dp),intent(in) :: ifccar(3,3),vect1(3)
5381  real(dp),intent(out) :: vect2(3),vect3(3)
5382 
5383 !Local variables -------------------------
5384 !scalars
5385  integer :: flag,ii,itrial,jj
5386  real(dp) :: innorm,scprod
5387 !arrays
5388  real(dp) :: work(3)
5389 
5390 ! *********************************************************************
5391 
5392  do jj=1,3
5393    work(jj)=zero
5394    do ii=1,3
5395      work(jj)=work(jj)+ifccar(jj,ii)*vect1(ii)
5396    end do
5397  end do
5398 
5399  flag=0
5400  do itrial=1,4
5401    scprod=zero
5402    do ii=1,3
5403      scprod=scprod+work(ii)*vect1(ii)
5404    end do
5405 
5406    do ii=1,3
5407      work(ii)=work(ii)-vect1(ii)*scprod
5408    end do
5409 
5410    scprod=zero
5411    do ii=1,3
5412      scprod=scprod+work(ii)**2
5413    end do
5414 
5415    if(scprod<1.0d-10)then
5416      work(1:3)=zero
5417      if(itrial>1)work(itrial-1)=1.0_dp
5418    else
5419      flag=1
5420    end if
5421 
5422    if(flag==1)exit
5423  end do
5424 
5425  innorm=scprod**(-0.5_dp)
5426  do ii=1,3
5427    vect2(ii)=work(ii)*innorm
5428  end do
5429 
5430  vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2)
5431  vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3)
5432  vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1)
5433 
5434 end subroutine axial9

m_dynmat/bigbx9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 or -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

 cell= (nrpt,3) Give the index of the the cell and irpt
 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

SOURCE

2993 subroutine bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)
2994 
2995 
2996 !This section has been created automatically by the script Abilint (TD).
2997 !Do not modify the following lines by hand.
2998 #undef ABI_FUNC
2999 #define ABI_FUNC 'bigbx9'
3000 !End of the abilint section
3001 
3002  implicit none
3003 
3004 !Arguments -------------------------------
3005 !scalars
3006  integer,intent(in) :: brav,choice,mrpt,nqshft
3007  integer,intent(out) :: nrpt
3008 !arrays
3009  integer,intent(in) :: ngqpt(3)
3010  real(dp),intent(in) :: rprim(3,3)
3011  real(dp),intent(out) :: rpt(3,mrpt)
3012  integer,intent(out) :: cell(3,mrpt)
3013 
3014 !Local variables -------------------------
3015 !In some cases, the atoms coordinates are not packed in the
3016 ! [0,1]^3 cube. Then, the parameter "buffer" might be increased,
3017 !to search relevant pairs of atoms in bigger boxes than usual.
3018 !scalars
3019  integer,parameter :: buffer=1
3020  integer :: irpt,lim1,lim2,lim3,lqshft,r1,r2,r3
3021  character(len=500) :: msg
3022 
3023 ! *********************************************************************
3024 
3025  lqshft=1
3026  if(nqshft/=1)lqshft=2
3027 
3028 
3029 !Simple Cubic Lattice
3030  if (abs(brav)==1) then
3031    lim1=((ngqpt(1))+1)*lqshft+buffer
3032    lim2=((ngqpt(2))+1)*lqshft+buffer
3033    lim3=((ngqpt(3))+1)*lqshft+buffer
3034    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
3035    if(choice/=0)then
3036      if (nrpt/=mrpt) then
3037        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
3038        MSG_BUG(msg)
3039      end if
3040      irpt=0
3041      do r1=-lim1,lim1
3042        do r2=-lim2,lim2
3043          do r3=-lim3,lim3
3044            irpt=irpt+1
3045            rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
3046            rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
3047            rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
3048            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
3049          end do
3050        end do
3051      end do
3052    end if
3053 
3054 !  Face Centered Cubic Lattice
3055  else if (brav==2) then
3056    lim1=((ngqpt(1)+3)/4)*lqshft+buffer
3057    lim2=((ngqpt(2)+3)/4)*lqshft+buffer
3058    lim3=((ngqpt(3)+3)/4)*lqshft+buffer
3059    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*4
3060    if(choice/=0)then
3061      if (nrpt/=mrpt) then
3062        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
3063        MSG_BUG(msg)
3064      end if
3065      irpt=0
3066      do r1=-lim1,lim1
3067        do r2=-lim2,lim2
3068          do r3=-lim3,lim3
3069            irpt=irpt+4
3070            rpt(1,irpt-3)=r1
3071            rpt(2,irpt-3)=r2
3072            rpt(3,irpt-3)=r3
3073            rpt(1,irpt-2)=r1
3074            rpt(2,irpt-2)=r2+0.5
3075            rpt(3,irpt-2)=r3+0.5
3076            rpt(1,irpt-1)=r1+0.5
3077            rpt(2,irpt-1)=r2
3078            rpt(3,irpt-1)=r3+0.5
3079            rpt(1,irpt)=r1+0.5
3080            rpt(2,irpt)=r2+0.5
3081            rpt(3,irpt)=r3
3082 !TEST_AM
3083 !           cell(irpt-3,1)=r1;cell(irpt-3,2)=r2;cell(irpt-3,3)=r3
3084            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
3085          end do
3086        end do
3087      end do
3088    end if
3089 
3090 !  Body Centered Cubic Lattice
3091  else if (brav==3) then
3092    lim1=((ngqpt(1)+3)/4)*lqshft+buffer
3093    lim2=((ngqpt(2)+3)/4)*lqshft+buffer
3094    lim3=((ngqpt(3)+3)/4)*lqshft+buffer
3095    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*2
3096    if(choice/=0)then
3097      if(nrpt/=mrpt) then
3098        write(msg,'(2(a,i0))')' nrpt= ',nrpt,' is not equal to mrpt= ',mrpt
3099        MSG_BUG(msg)
3100      end if
3101      irpt=0
3102      do r1=-lim1,lim1
3103        do r2=-lim2,lim2
3104          do r3=-lim3,lim3
3105            irpt=irpt+2
3106            rpt(1,irpt-1)=r1
3107            rpt(2,irpt-1)=r2
3108            rpt(3,irpt-1)=r3
3109            rpt(1,irpt)=r1+0.5
3110            rpt(2,irpt)=r2+0.5
3111            rpt(3,irpt)=r3+0.5
3112 !TEST_AM
3113 !           cell(irpt-1,1)=r1;cell(irpt-1,2)=r2;cell(irpt-1,3)=r3
3114            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
3115          end do
3116        end do
3117      end do
3118    end if
3119 
3120 !  Hexagonal Lattice
3121  else if (brav==4) then
3122    lim1=(ngqpt(1)+1)*lqshft+buffer
3123    lim2=(ngqpt(2)+1)*lqshft+buffer
3124    lim3=((ngqpt(3)/2)+1)*lqshft+buffer
3125    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
3126    if(choice/=0)then
3127      if(nrpt/=mrpt)then
3128        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt=',mrpt
3129        MSG_BUG(msg)
3130      end if
3131      irpt=0
3132      do r1=-lim1,lim1
3133        do r2=-lim2,lim2
3134          do r3=-lim3,lim3
3135            irpt=irpt+1
3136            rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
3137            rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
3138            rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
3139            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
3140          end do
3141        end do
3142      end do
3143    end if
3144 
3145  else
3146    write(msg,'(a,i0,a)')' The value of brav= ',brav,' is not allowed (should be -1, 1, 2 or 4).'
3147    MSG_BUG(msg)
3148  end if
3149 
3150 end subroutine bigbx9

m_dynmat/canat9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 or -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

SOURCE

3181 subroutine canat9(brav,natom,rcan,rprim,trans,xred)
3182 
3183 
3184 !This section has been created automatically by the script Abilint (TD).
3185 !Do not modify the following lines by hand.
3186 #undef ABI_FUNC
3187 #define ABI_FUNC 'canat9'
3188 !End of the abilint section
3189 
3190  implicit none
3191 
3192 !Arguments -------------------------------
3193 !scalars
3194  integer,intent(in) :: brav,natom
3195 !arrays
3196  real(dp),intent(in) :: rprim(3,3),xred(3,natom)
3197  real(dp),intent(out) :: rcan(3,natom),trans(3,natom)
3198 
3199 !Local variables -------------------------
3200 !scalars
3201  integer :: found,iatom,ii
3202  character(len=500) :: message
3203 !arrays
3204  real(dp) :: dontno(3,4),rec(3),rok(3),shift(3),tt(3)
3205 
3206 ! *********************************************************************
3207 
3208  DBG_ENTER("COLL")
3209 
3210 !Normalization of the cartesian atomic coordinates
3211 !If not normalized : rcan(i) <- rcan(i) * acell(i)
3212  do iatom=1,natom
3213    rcan(:,iatom)=xred(1,iatom)*rprim(:,1)+&
3214 &   xred(2,iatom)*rprim(:,2)+&
3215 &   xred(3,iatom)*rprim(:,3)
3216  end do
3217 
3218 !Study of the different cases for the Bravais lattice:
3219 
3220 !Simple Cubic Lattice
3221  if (abs(brav)==1) then
3222 
3223    do iatom=1,natom
3224 
3225 !    Canon will produces these coordinate transformations
3226 !    (Note : here we still use reduced coordinates )
3227      call wrap2_pmhalf(xred(1,iatom),rok(1),shift(1))
3228      call wrap2_pmhalf(xred(2,iatom),rok(2),shift(2))
3229      call wrap2_pmhalf(xred(3,iatom),rok(3),shift(3))
3230 
3231 !    New coordinates : rcan
3232      rcan(:,iatom)=rok(1)*rprim(:,1)+&
3233 &     rok(2)*rprim(:,2)+&
3234 &     rok(3)*rprim(:,3)
3235 !    Translations between New and Old coordinates
3236      tt(:)=xred(1,iatom)*rprim(:,1)+&
3237 &     xred(2,iatom)*rprim(:,2)+&
3238 &     xred(3,iatom)*rprim(:,3)
3239      trans(:,iatom)=tt(:)-rcan(:,iatom)
3240    end do
3241 
3242  else if (brav==2) then
3243 !  Face Centered Lattice
3244 !  Special possible translations in the F.C.C. case
3245    dontno(:,:)=zero
3246    dontno(2,2)=0.5_dp
3247    dontno(3,2)=0.5_dp
3248    dontno(1,3)=0.5_dp
3249    dontno(3,3)=0.5_dp
3250    dontno(1,4)=0.5_dp
3251    dontno(2,4)=0.5_dp
3252    do iatom=1,natom
3253      found=0
3254      do ii=1,4
3255        if (found==1) exit
3256 !      Canon will produces these coordinate transformations
3257        call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
3258        call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
3259        call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
3260 !      In the F.C.C., ABS[ Ri ] + ABS[ Rj ] < or = 1/2
3261 !      The equal signs has been treated using a tolerance parameter
3262 !      not to have twice the same point in the unit cell !
3263        rok(1)=rok(1)-1.0d-10
3264        rok(2)=rok(2)-2.0d-10
3265        rok(3)=rok(3)-5.0d-10
3266        if (abs(rok(1))+abs(rok(2))<=0.5_dp) then
3267          if (abs(rok(1))+abs(rok(3))<=0.5_dp) then
3268            if (abs(rok(2))+abs(rok(3))<=0.5_dp) then
3269              tt(:)=rcan(:,iatom)
3270 !            New coordinates : rcan
3271              rcan(1,iatom)=rok(1)+1.0d-10
3272              rcan(2,iatom)=rok(2)+2.0d-10
3273              rcan(3,iatom)=rok(3)+5.0d-10
3274 !            Translations between New and Old coordinates
3275              trans(:,iatom)=tt(:)-rcan(:,iatom)
3276              found=1
3277            end if
3278          end if
3279        end if
3280      end do
3281    end do
3282 
3283  else if (brav==3) then
3284 !  Body Centered Cubic Lattice
3285 !  Special possible translations in the B.C.C. case
3286 
3287    dontno(:,1)=zero
3288    dontno(:,2)=0.5_dp
3289    do iatom=1,natom
3290      found=0
3291      do ii=1,2
3292        if (found==1) exit
3293 !      Canon will produce these coordinate transformations
3294        call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
3295        call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
3296        call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
3297 !      In the F.C.C., ABS[ Ri ] < or = 1/2
3298 !      and    ABS[ R1 ] + ABS[ R2 ] + ABS[ R3 ] < or = 3/4
3299 !      The equal signs has been treated using a tolerance parameter
3300 !      not to have twice the same point in the unit cell !
3301        rok(1)=rok(1)-1.0d-10
3302        rok(2)=rok(2)-2.0d-10
3303        rok(3)=rok(3)-5.0d-10
3304        if(abs(rok(1))+abs(rok(2))+abs(rok(3))<=0.75_dp)then
3305          if ( abs(rok(1))<=0.5_dp .and.&
3306 &         abs(rok(2))<=0.5_dp .and.&
3307 &         abs(rok(3))<=0.5_dp       ) then
3308            tt(:)=rcan(:,iatom)
3309 !          New coordinates : rcan
3310            rcan(1,iatom)=rok(1)+1.0d-10
3311            rcan(2,iatom)=rok(2)+2.0d-10
3312            rcan(3,iatom)=rok(3)+5.0d-10
3313 !          Translations between New and Old coordinates
3314            trans(:,iatom)=tt(:)-rcan(:,iatom)
3315            found=1
3316          end if
3317        end if
3318      end do
3319    end do
3320 
3321  else if (brav==4) then
3322 !  Hexagonal Lattice
3323 !  In this case, it is easier first to work in reduced coordinates space !
3324    do iatom=1,natom
3325 !    Passage from the reduced space to the "lozenge" cell
3326      rec(1)=xred(1,iatom)-0.5_dp
3327      rec(2)=xred(2,iatom)-0.5_dp
3328      rec(3)=xred(3,iatom)
3329 !    Canon will produces these coordinate transformations
3330      call wrap2_pmhalf(rec(1),rok(1),shift(1))
3331      call wrap2_pmhalf(rec(2),rok(2),shift(2))
3332      call wrap2_pmhalf(rec(3),rok(3),shift(3))
3333      rec(1)=rok(1)+0.5_dp
3334      rec(2)=rok(2)+0.5_dp
3335      rec(3)=rok(3)
3336 !    Passage in Cartesian Normalized Coordinates
3337      rcan(:,iatom)=rec(1)*rprim(:,1)+&
3338 &     rec(2)*rprim(:,2)+&
3339 &     rec(3)*rprim(:,3)
3340 !    Use of a tolerance parameter not to have twice the same point
3341 !    in the unit cell !
3342      rcan(1,iatom)=rcan(1,iatom)-1.0d-10
3343      rcan(2,iatom)=rcan(2,iatom)-2.0d-10
3344 !    Passage to the honey-com hexagonal unit cell !
3345      if (rcan(1,iatom)>0.5_dp) then
3346        rcan(1,iatom)=rcan(1,iatom)-1.0_dp
3347      end if
3348      if (rcan(1,iatom)>zero.and.rcan(1,iatom)+sqrt(3.0_dp)*rcan&
3349 &     (2,iatom)>1.0_dp) then
3350        rcan(1,iatom)=rcan(1,iatom)-0.5_dp
3351        rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
3352      end if
3353      if (rcan(1,iatom)<=zero.and.sqrt(3.0_dp)*rcan(2,iatom)-rcan&
3354 &     (1,iatom)>1.0_dp) then
3355        rcan(1,iatom)=rcan(1,iatom)+0.5_dp
3356        rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
3357      end if
3358 !    Translations between New and Old coordinates
3359      tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)&
3360 &     +xred(3,iatom)*rprim(:,3)
3361      trans(:,iatom)=tt(:)-rcan(:,iatom)
3362    end do
3363 
3364 !  End of the possible cases for brav : -1, 1, 2, 4.
3365  else
3366 
3367    write(message, '(a,i0,a,a,a)' )&
3368 &   'The required value of brav=',brav,' is not available.',ch10,&
3369 &   'It should be -1, 1,2 or 4 .'
3370    MSG_BUG(message)
3371  end if
3372 
3373  call wrtout(std_out,' Canonical Atomic Coordinates ',"COLL")
3374 
3375  do iatom=1,natom
3376    write(message, '(a,i5,3es18.8)' )' atom',iatom,rcan(1,iatom),rcan(2,iatom),rcan(3,iatom)
3377    call wrtout(std_out,message,'COLL')
3378  end do
3379 
3380  DBG_EXIT("COLL")
3381 
3382 end subroutine canat9

m_dynmat/canct9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      ddb_hybrid,m_ifc

CHILDREN

SOURCE

3418 subroutine canct9(acell,gprim,ib,index,irpt,natom,nrpt,rcan,rcart,rprim,rpt)
3419 
3420 
3421 !This section has been created automatically by the script Abilint (TD).
3422 !Do not modify the following lines by hand.
3423 #undef ABI_FUNC
3424 #define ABI_FUNC 'canct9'
3425 !End of the abilint section
3426 
3427  implicit none
3428 
3429 !Arguments -------------------------------
3430 !scalars
3431  integer,intent(in) :: index,natom,nrpt
3432  integer,intent(out) :: ib,irpt
3433 !arrays
3434  real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3)
3435  real(dp),intent(in) :: rpt(3,nrpt)
3436  real(dp),intent(out) :: rcart(3)
3437 
3438 !Local variables -------------------------
3439 !scalars
3440  integer :: jj
3441 !arrays
3442  real(dp) :: xred(3)
3443 
3444 ! *********************************************************************
3445 
3446  irpt=(index-1)/natom+1
3447  ib=index-natom*(irpt-1)
3448 
3449 !Transform the canonical coordinates to reduced coord.
3450  do jj=1,3
3451    xred(jj)=gprim(1,jj)*(rpt(1,irpt)+rcan(1,ib))&
3452 &   +gprim(2,jj)*(rpt(2,irpt)+rcan(2,ib))&
3453 &   +gprim(3,jj)*(rpt(3,irpt)+rcan(3,ib))
3454  end do
3455 
3456 !Then to cartesian coordinates (here the position of the atom b)
3457  do jj=1,3
3458    rcart(jj)=xred(1)*acell(1)*rprim(jj,1)+&
3459 &   xred(2)*acell(2)*rprim(jj,2)+&
3460 &   xred(3)*acell(3)*rprim(jj,3)
3461  end do
3462 
3463 end subroutine canct9

m_dynmat/cart29 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      dfpt_gatherdy,m_ddb

CHILDREN

SOURCE

752 subroutine cart29(blkflg,blkval,carflg,d2cart,&
753 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
754 
755 
756 !This section has been created automatically by the script Abilint (TD).
757 !Do not modify the following lines by hand.
758 #undef ABI_FUNC
759 #define ABI_FUNC 'cart29'
760 !End of the abilint section
761 
762  implicit none
763 
764 !Arguments -------------------------------
765 !scalars
766  integer,intent(in) :: iblok,mpert,natom,nblok,ntypat
767  real(dp),intent(in) :: ucvol
768 !arrays
769  integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),typat(natom)
770  integer,intent(out) :: carflg(3,mpert,3,mpert)
771  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3)
772  real(dp),intent(in) :: zion(ntypat)
773  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
774 
775 !Local variables -------------------------
776 !scalars
777  integer :: idir1,idir2,ii,ipert1,ipert2
778 !arrays
779  integer :: flg1(3),flg2(3)
780  real(dp) :: vec1(3),vec2(3)
781 
782 ! *********************************************************************
783 
784 !First, copy the data blok in place.
785  d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok)
786 
787 !Cartesian coordinates transformation (in two steps)
788 !First step
789  do ipert1=1,mpert
790    do ipert2=1,mpert
791      do ii=1,2
792        do idir1=1,3
793          do idir2=1,3
794            vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2)
795 !          Note here blkflg
796            flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok)
797          end do
798          call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2)
799          do idir2=1,3
800            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2)
801 !          And here carflg
802            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2)
803          end do
804        end do
805      end do
806    end do
807  end do
808 
809 !Second step
810  do ipert1=1,mpert
811    do ipert2=1,mpert
812      do ii=1,2
813        do idir2=1,3
814          do idir1=1,3
815            vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2)
816 !          Note here carflg
817            flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2)
818          end do
819          call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2)
820          do idir1=1,3
821            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1)
822 !          And here carflg again
823            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1)
824          end do
825        end do
826      end do
827    end do
828  end do
829 
830 !For the dielectric tensor, takes into account the volume
831 !of the unit cell, and add the unit matrix (polarization of the vacuum)
832  do idir1=1,3
833    do idir2=1,3
834      do ii=1,2
835        d2cart(ii,idir1,natom+2,idir2,natom+2)=&
836 &       -four_pi/ucvol*d2cart(ii,idir1,natom+2,idir2,natom+2)
837      end do
838    end do
839  end do
840 
841  do idir1=1,3
842    d2cart(1,idir1,natom+2,idir1,natom+2)=&
843 &   1.0_dp+d2cart(1,idir1,natom+2,idir1,natom+2)
844  end do
845 
846 !Add the ionic charges to delta z to get the effective charges
847  do ipert1=1,natom
848    do idir1=1,3
849      d2cart(1,idir1,ipert1,idir1,natom+2)=&
850 &     zion(typat(ipert1))+d2cart(1,idir1,ipert1,idir1,natom+2)
851    end do
852  end do
853  do ipert2=1,natom
854    do idir2=1,3
855      d2cart(1,idir2,natom+2,idir2,ipert2)=&
856 &     zion(typat(ipert2))+d2cart(1,idir2,natom+2,idir2,ipert2)
857    end do
858  end do
859 
860 !For the piezoelectric tensor, takes into account the volume of the unit cell
861  do ipert2=natom+3,natom+4
862    do idir1=1,3
863      do idir2=1,3
864        do ii=1,2
865          d2cart(ii,idir1,natom+2,idir2,ipert2)=&
866 &         (1.0_dp/ucvol)*d2cart(ii,idir1,natom+2,idir2,ipert2)
867          d2cart(ii,idir2,ipert2,idir1,natom+2)=&
868 &         (1.0_dp/ucvol)*d2cart(ii,idir2,ipert2,idir1,natom+2)
869        end do
870      end do
871    end do
872  end do
873 
874 end subroutine cart29

m_dynmat/cart39 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      dfpt_gatherdy,m_ddb,m_dynmat

CHILDREN

SOURCE

908 subroutine cart39(flg1,flg2,gprimd,ipert,natom,rprimd,vec1,vec2)
909 
910 
911 !This section has been created automatically by the script Abilint (TD).
912 !Do not modify the following lines by hand.
913 #undef ABI_FUNC
914 #define ABI_FUNC 'cart39'
915 !End of the abilint section
916 
917  implicit none
918 
919 !Arguments -------------------------------
920 !scalars
921  integer,intent(in) :: ipert,natom
922 !arrays
923  integer,intent(in) :: flg1(3)
924  integer,intent(out) :: flg2(3)
925  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3),vec1(3)
926  real(dp),intent(out) :: vec2(3)
927 
928 !Local variables -------------------------
929 !scalars
930  integer :: idir,ii
931 
932 ! *********************************************************************
933 
934 !Treat phonon-type perturbation
935  if(ipert>=1.and.ipert<=natom)then
936 
937    do idir=1,3
938      vec2(idir)=zero
939      flg2(idir)=1
940      do ii=1,3
941        if(abs(gprimd(idir,ii))>1.0d-10)then
942          if(flg1(ii)==1)then
943            vec2(idir)=vec2(idir)+gprimd(idir,ii)*vec1(ii)
944          else
945            flg2(idir)=0
946          end if
947        end if
948      end do
949      if(flg2(idir)==0)vec2(idir)=zero
950    end do
951 
952 !  Treat electric field perturbation
953  else if(ipert==natom+2) then
954 !  OCL SCALAR
955    do idir=1,3
956      vec2(idir)=zero
957      flg2(idir)=1
958 !    OCL SCALAR
959      do ii=1,3
960        if(abs(rprimd(idir,ii))>1.0d-10)then
961          if(flg1(ii)==1)then
962            vec2(idir)=vec2(idir)+rprimd(idir,ii)*vec1(ii)/two_pi
963          else
964            flg2(idir)=0
965          end if
966        end if
967      end do
968      if(flg2(idir)==0)vec2(idir)=zero
969    end do
970 
971 !  Treat other perturbations
972  else
973    do idir=1,3
974      vec2(idir)=vec1(idir)
975      flg2(idir)=flg1(idir)
976    end do
977  end if
978 
979 end subroutine cart39

m_dynmat/chkph3 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

1172 subroutine chkph3(carflg,idir,mpert,natom)
1173 
1174 
1175 !This section has been created automatically by the script Abilint (TD).
1176 !Do not modify the following lines by hand.
1177 #undef ABI_FUNC
1178 #define ABI_FUNC 'chkph3'
1179 !End of the abilint section
1180 
1181  implicit none
1182 
1183 !Arguments -------------------------------
1184 !scalars
1185  integer,intent(in) :: idir,mpert,natom
1186 !arrays
1187  integer,intent(in) :: carflg(3,mpert,3,mpert)
1188 
1189 !Local variables -------------------------
1190 !scalars
1191  integer :: idir1,idir2,ipert1,ipert2,send
1192  character(len=500) :: message
1193 
1194 ! *********************************************************************
1195 
1196  send=0
1197 
1198 !Check the elements of the analytical part of the dynamical matrix
1199  do ipert2=1,natom
1200    do idir2=1,3
1201      do ipert1=1,natom
1202        do idir1=1,3
1203          if(carflg(idir1,ipert1,idir2,ipert2)==0)then
1204            send=1
1205          end if
1206        end do
1207      end do
1208    end do
1209  end do
1210 
1211 !If some electric field is present
1212  if(idir/=0)then
1213 
1214 !  Check the dielectric constant
1215    if(carflg(idir,natom+2,idir,natom+2)==0)then
1216      send=1
1217    end if
1218 
1219 !  Check the effective charges
1220    do ipert1=1,natom
1221      do idir1=1,3
1222        if(carflg(idir1,ipert1,idir,natom+2)==0)then
1223          send=1
1224        end if
1225      end do
1226    end do
1227 
1228  end if
1229 
1230 !If needed, send the message
1231  if(send==1)then
1232    write(message, '(a,a,a,a)' )&
1233 &   ' chkph3 : WARNING -',ch10,&
1234 &   '  The dynamical matrix was incomplete :',&
1235 &   ' phonon frequencies may be wrong ...'
1236    call wrtout(std_out,message,'COLL')
1237    call wrtout(ab_out,message,'COLL')
1238  end if
1239 
1240 end subroutine chkph3

m_dynmat/chkrp9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 or -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

SOURCE

3492 subroutine chkrp9(brav,rprim)
3493 
3494 
3495 !This section has been created automatically by the script Abilint (TD).
3496 !Do not modify the following lines by hand.
3497 #undef ABI_FUNC
3498 #define ABI_FUNC 'chkrp9'
3499 !End of the abilint section
3500 
3501  implicit none
3502 
3503 !Arguments -------------------------------
3504 !scalars
3505  integer,intent(in) :: brav
3506 !arrays
3507  real(dp),intent(in) :: rprim(3,3)
3508 
3509 !Local variables -------------------------
3510 !scalars
3511  integer :: ii,jj
3512  character(len=500) :: message
3513 
3514 ! *********************************************************************
3515 
3516  if (abs(brav)==1) then
3517 !  Simple Cubic Lattice No condition in this case !
3518    continue
3519 
3520  else if (brav==2) then
3521 !  Face Centered Lattice
3522    do ii=1,3
3523      do jj=1,3
3524        if (  ( ii==jj .and. abs(rprim(ii,jj))>tol10) .or.&
3525 &       ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
3526          write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3527 &         'The input variable rprim does not correspond to the',ch10,&
3528 &         'fixed rprim to be used with brav=2 and ifcflag=1 :',ch10,&
3529 &         '   0  1/2  1/2',ch10,&
3530 &         '  1/2  0   1/2',ch10,&
3531 &         '  1/2 1/2   0 ',ch10,&
3532 &         'Action: rebuild your DDB by using the latter rprim.'
3533          MSG_ERROR(message)
3534        end if
3535      end do
3536    end do
3537 
3538  else if (brav==3) then
3539 !  Body Centered Cubic Lattice
3540    do ii=1,3
3541      do jj=1,3
3542        if (  ( ii==jj .and. abs(rprim(ii,jj)+.5_dp)>tol10) .or.&
3543 &       ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
3544          write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3545 &         'The input variable rprim does not correspond to the',ch10,&
3546 &         'fixed rprim to be used with brav=3 and ifcflag=1 :',ch10,&
3547 &         '  -1/2  1/2  1/2',ch10,&
3548 &         '   1/2 -1/2  1/2',ch10,&
3549 &         '   1/2  1/2 -1/2',ch10,&
3550 &         'Action: rebuild your DDB by using the latter rprim.'
3551          MSG_ERROR(message)
3552        end if
3553      end do
3554    end do
3555 
3556  else if (brav==4) then
3557 !  Hexagonal Lattice
3558    if (abs(rprim(1,1)-1.0_dp)>tol10 .or. &
3559 &   abs(rprim(3,3)-1.0_dp)>tol10 .or. &
3560 &   abs(rprim(2,1)      )>tol10 .or. &
3561 &   abs(rprim(3,1)      )>tol10 .or. &
3562 &   abs(rprim(1,3)      )>tol10 .or. &
3563 &   abs(rprim(2,3)      )>tol10 .or. &
3564 &   abs(rprim(3,2)      )>tol10 .or. &
3565 &   abs(rprim(1,2)+0.5_dp)>tol10 .or. &
3566 &   abs(rprim(2,2)-0.5_dp*sqrt(3.0_dp))>tol10 ) then
3567      write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3568 &     'The input variable rprim does not correspond to the',ch10,&
3569 &     'fixed rprim to be used with brav=4 and ifcflag=1 :',ch10,&
3570 &     '   1      0      0',ch10,&
3571 &     '  -1/2 sqrt[3]/2 0',ch10,&
3572 &     '   0      0      1',ch10,&
3573 &     'Action: rebuild your DDB by using the latter rprim.'
3574      MSG_ERROR(message)
3575    end if
3576 
3577  else
3578 
3579    write(message, '(a,i4,a,a,a,a,a)' )&
3580 &   'The value of brav=',brav,' is not allowed.',ch10,&
3581 &   'Only  -1, 1,2,3 or 4 are allowed.',ch10,&
3582 &   'Action: change the value of brav in your input file.'
3583    MSG_ERROR(message)
3584  end if
3585 
3586 end subroutine chkrp9

m_dynmat/chneu9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

      dfpt_gatherdy,m_ddb

CHILDREN

SOURCE

1275 subroutine chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion)
1276 
1277 
1278 !This section has been created automatically by the script Abilint (TD).
1279 !Do not modify the following lines by hand.
1280 #undef ABI_FUNC
1281 #define ABI_FUNC 'chneu9'
1282 !End of the abilint section
1283 
1284  implicit none
1285 
1286 !Arguments -------------------------------
1287 !scalars
1288  integer,intent(in) :: chneut,mpert,natom,ntypat,selectz
1289 !arrays
1290  integer,intent(in) :: typat(natom)
1291  real(dp),intent(in) :: zion(ntypat)
1292  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
1293 
1294 !Local variables -------------------------
1295 !scalars
1296  integer :: idir1,idir2,ii,ipert1,ipert2
1297  character(len=500) :: message
1298 !arrays
1299  real(dp) :: sumwght(2)
1300  real(dp),allocatable :: wghtat(:)
1301 
1302 ! *********************************************************************
1303 
1304  ABI_ALLOCATE(wghtat,(natom))
1305 
1306 !In case of acoustic sum rule imposition, compute the weights on each atom.
1307  if (chneut==1)then
1308 
1309 !  The weight is the same for all atom
1310    do ipert1=1,natom
1311      wghtat(ipert1)=1./natom
1312    end do
1313 
1314  else if (chneut==2) then
1315 
1316 !  The weight is proportional to the diagonal electronic screening charge of the atom
1317    sumwght(1)=zero
1318    do ipert1=1,natom
1319      wghtat(ipert1)=zero
1320      do idir1=1,3
1321        wghtat(ipert1)=wghtat(ipert1)+&
1322 &       d2cart(1,idir1,ipert1,idir1,natom+2)+&
1323 &       d2cart(1,idir1,natom+2,idir1,ipert1)-2*zion(typat(ipert1))
1324      end do
1325      sumwght(1)=sumwght(1)+wghtat(ipert1)
1326    end do
1327 
1328 !  Normalize the weights to unity
1329    do ipert1=1,natom
1330      wghtat(ipert1)=wghtat(ipert1)/sumwght(1)
1331    end do
1332  end if
1333 
1334 !Calculation of the violation of the charge neutrality
1335 !and imposition of the charge neutrality condition
1336  if (chneut/=0)then
1337    write(message, '(a,a,a,a,a,a,a)' )&
1338 &   ' The violation of the charge neutrality conditions',ch10,&
1339 &   ' by the effective charges is as follows :',ch10,&
1340 &   '    atom        electric field',ch10,&
1341 &   ' displacement     direction   '
1342    call wrtout(ab_out,message,'COLL')
1343    do idir1=1,3
1344      do idir2=1,3
1345        do ii=1,2
1346          sumwght(ii)=zero
1347          do ipert1=1,natom
1348            sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir2,natom+2)
1349          end do
1350          do ipert1=1,natom
1351            d2cart(ii,idir1,ipert1,idir2,natom+2)=&
1352 &           d2cart(ii,idir1,ipert1,idir2,natom+2)-sumwght(ii)*wghtat(ipert1)
1353          end do
1354        end do
1355        write(message, '(i8,i16,2f16.6)' ) idir1,idir2,sumwght(1),sumwght(2)
1356        call wrtout(ab_out,message,'COLL')
1357      end do
1358    end do
1359    write(message, '(a)' )' '
1360    call wrtout(ab_out,message,'COLL')
1361 
1362 !  The same for the symmetrical part
1363    do idir1=1,3
1364      do idir2=1,3
1365        do ii=1,2
1366          sumwght(ii)=zero
1367          do ipert2=1,natom
1368            sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir2,ipert2)
1369          end do
1370          do ipert2=1,natom
1371            d2cart(ii,idir1,natom+2,idir2,ipert2)=&
1372 &           d2cart(ii,idir1,natom+2,idir2,ipert2)-sumwght(ii)*wghtat(ipert2)
1373          end do
1374        end do
1375      end do
1376    end do
1377  end if
1378 
1379 !Selection of the trace of the effective charge tensor attached to each atom
1380  if(selectz==1)then
1381    do ipert1=1,natom
1382      do ii=1,2
1383        sumwght(ii)=zero
1384        do idir1=1,3
1385          sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir1,natom+2)
1386        end do
1387        do idir1=1,3
1388          do idir2=1,3
1389            d2cart(ii,idir1,ipert1,idir2,natom+2)=zero
1390          end do
1391        end do
1392        do idir1=1,3
1393          d2cart(ii,idir1,ipert1,idir1,natom+2)=sumwght(ii)/3.0_dp
1394        end do
1395      end do
1396    end do
1397 !  Do the same for the symmetrical part of d2cart
1398    do ipert2=1,natom
1399      do ii=1,2
1400        sumwght(ii)=zero
1401        do idir1=1,3
1402          sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir1,ipert2)
1403        end do
1404        do idir1=1,3
1405          do idir2=1,3
1406            d2cart(ii,idir1,natom+2,idir2,ipert2)=zero
1407          end do
1408        end do
1409        do idir1=1,3
1410          d2cart(ii,idir1,natom+2,idir1,ipert2)=sumwght(ii)/3.0_dp
1411        end do
1412      end do
1413    end do
1414  end if
1415 
1416 !Selection of the symmetric part of the effective charge tensor attached to each atom
1417  if(selectz==2)then
1418    do ipert1=1,natom
1419      do ii=1,2
1420        do idir1=1,3
1421          do idir2=1,3
1422            sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
1423 &           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
1424            d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii)
1425            d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii)
1426          end do
1427        end do
1428      end do
1429    end do
1430 !  Do the same for the symmetrical part of d2cart
1431    do ipert1=1,natom
1432      do ii=1,2
1433        do idir1=1,3
1434          do idir2=1,3
1435            sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
1436 &           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
1437            d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii)
1438            d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii)
1439          end do
1440        end do
1441      end do
1442    end do
1443  end if
1444 
1445 !Write the effective charge tensor
1446  write(message, '(a,a,a,a,a,a,a)' )&
1447 & ' Effective charge tensors after ',ch10,&
1448 & ' imposition of the charge neutrality,',ch10,&
1449 & ' and eventual restriction to some part :',ch10,&
1450 & '   atom    displacement  '
1451  call wrtout(ab_out,message,'COLL')
1452 
1453  do ipert1=1,natom
1454    do idir1=1,3
1455      write(message, '(2i10,3es16.6)' )ipert1,idir1,&
1456 &     (d2cart(1,idir1,ipert1,idir2,natom+2),idir2=1,3)
1457      call wrtout(ab_out,message,'COLL')
1458    end do
1459  end do
1460 
1461 !Zero the imaginary part of the dynamical matrix
1462  write(message, '(a)' )' Now, the imaginary part of the dynamical matrix is zeroed '
1463  call wrtout(ab_out,message,'COLL')
1464  call wrtout(std_out,message,'COLL')
1465 
1466  do ipert1=1,natom
1467    do ipert2=1,natom
1468      do idir1=1,3
1469        do idir2=1,3
1470          d2cart(2,idir1,ipert1,idir2,ipert2)=zero
1471        end do
1472      end do
1473    end do
1474  end do
1475 
1476  ABI_DEALLOCATE(wghtat)
1477 
1478 end subroutine chneu9

m_dynmat/d2cart_to_red [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 d2cart_to_red

FUNCTION

 Transform a second-derivative matrix from cartesian
 coordinates to reduced coordinate. Also,
 1) remove the ionic part of the effective charges,
 2) remove the vacuum polarisation from the dielectric tensor
    and scale it with the unit cell volume
 In short, does the inverse operation of cart29.

INPUTS

  d2cart(2,3,mpert,3,mpert)=
    second-derivative matrix in cartesian coordinates
  gprimd(3,3)=basis vector in the reciprocal space
  rprimd(3,3)=basis vector in the real space
  mpert =maximum number of ipert
  natom=number of atom

OUTPUT

  d2red(2,3,mpert,3,mpert)=
    second-derivative matrix in reduced coordinates

PARENTS

      ddb_interpolate

CHILDREN

SOURCE

1016 subroutine d2cart_to_red(d2cart, d2red, gprimd, rprimd, mpert, natom, &
1017 &                        ntypat,typat,ucvol,zion)
1018 
1019 
1020 !This section has been created automatically by the script Abilint (TD).
1021 !Do not modify the following lines by hand.
1022 #undef ABI_FUNC
1023 #define ABI_FUNC 'd2cart_to_red'
1024 !End of the abilint section
1025 
1026  implicit none
1027 
1028 !Arguments -------------------------------
1029 !scalars
1030  integer,intent(in) :: mpert,natom,ntypat
1031  real(dp),intent(in) :: ucvol
1032 !arrays
1033  integer,intent(in) :: typat(natom)
1034  real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
1035  real(dp),intent(out) :: d2red(2,3,mpert,3,mpert)
1036  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
1037  real(dp),intent(in) :: zion(ntypat)
1038 
1039 !Local variables -------------------------
1040 !scalars
1041  integer :: idir1,idir2,ii,ipert1,ipert2
1042  real(dp) :: fac
1043 !arrays
1044  integer :: flg1(3),flg2(3)
1045  real(dp) :: vec1(3),vec2(3)
1046 
1047 ! *********************************************************************
1048 
1049  flg1 = one
1050  flg2 = one
1051 
1052  d2red = d2cart
1053 
1054 !Remove the ionic charges to z to get the change in effective charges
1055  do ipert1=1,natom
1056    do idir1=1,3
1057      d2red(1,idir1,ipert1,idir1,natom+2)=&
1058 &     d2red(1,idir1,ipert1,idir1,natom+2) - zion(typat(ipert1))
1059    end do
1060  end do
1061  do ipert2=1,natom
1062    do idir2=1,3
1063      d2red(1,idir2,natom+2,idir2,ipert2)=&
1064 &     d2red(1,idir2,natom+2,idir2,ipert2) - zion(typat(ipert2))
1065    end do
1066  end do
1067 
1068  ! Remove the vacuum polarizability from the dielectric tensor
1069  do idir1=1,3
1070    d2red(1,idir1,natom+2,idir1,natom+2)=&
1071 &   d2red(1,idir1,natom+2,idir1,natom+2) - 1.0_dp
1072  end do
1073 
1074 ! Scale the dielectric tensor with the volue of the unit cell
1075  do idir1=1,3
1076    do idir2=1,3
1077      do ii=1,2
1078        d2red(ii,idir1,natom+2,idir2,natom+2)=&
1079 &       - (ucvol / four_pi) * d2red(ii,idir1,natom+2,idir2,natom+2)
1080      end do
1081    end do
1082  end do
1083 
1084 !For the piezoelectric tensor, takes into account the volume of the unit cell
1085  do ipert2=natom+3,natom+4
1086    do idir1=1,3
1087      do idir2=1,3
1088        do ii=1,2
1089          d2red(ii,idir1,natom+2,idir2,ipert2)=&
1090 &         (ucvol)*d2red(ii,idir1,natom+2,idir2,ipert2)
1091          d2red(ii,idir2,ipert2,idir1,natom+2)=&
1092 &         (ucvol)*d2red(ii,idir2,ipert2,idir1,natom+2)
1093        end do
1094      end do
1095    end do
1096  end do
1097 
1098 ! Reduced coordinates transformation (in two steps)
1099 ! Note that rprimd and gprimd are swapped, compared to what cart39 expects
1100 ! A factor of (2pi) ** 2 is added to transform the electric field perturbations
1101 
1102 !First step
1103  do ipert1=1,mpert
1104    fac = one; if (ipert1==natom+2) fac = two_pi ** 2
1105 
1106    do ipert2=1,mpert
1107      do ii=1,2
1108        do idir1=1,3
1109          do idir2=1,3
1110            vec1(idir2)=d2red(ii,idir1,ipert1,idir2,ipert2)
1111          end do
1112          ! Transform vector from cartesian to reduced coordinates
1113          call cart39(flg1,flg2,transpose(rprimd),ipert1,natom,transpose(gprimd),vec1,vec2)
1114          do idir2=1,3
1115            d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) * fac
1116          end do
1117        end do
1118      end do
1119    end do
1120  end do
1121 
1122 !Second step
1123  do ipert1=1,mpert
1124    do ipert2=1,mpert
1125      fac = one; if (ipert2==natom+2) fac = two_pi ** 2
1126 
1127      do ii=1,2
1128        do idir2=1,3
1129          do idir1=1,3
1130            vec1(idir1)=d2red(ii,idir1,ipert1,idir2,ipert2)
1131          end do
1132          ! Transform vector from cartesian to reduced coordinates
1133          call cart39(flg1,flg2,transpose(rprimd),ipert2,natom,transpose(gprimd),vec1,vec2)
1134          do idir1=1,3
1135            d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) * fac
1136          end do
1137        end do
1138      end do
1139    end do
1140  end do
1141 
1142 
1143 end subroutine d2cart_to_red

m_dynmat/d2sym3 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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
  zero_by_symm= if 1, set blkflg to 1 for the elements that must be zero by symmetry, and zero them.
    This has the indirect effect of being able to resymmetrize the whole matrix, thus
    enforcing better the symmetry for the 2DTE.

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,respfn

CHILDREN

SOURCE

1546 subroutine d2sym3(blkflg,d2,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev,zero_by_symm)
1547 
1548 
1549 !This section has been created automatically by the script Abilint (TD).
1550 !Do not modify the following lines by hand.
1551 #undef ABI_FUNC
1552 #define ABI_FUNC 'd2sym3'
1553 !End of the abilint section
1554 
1555  implicit none
1556 
1557 !Arguments -------------------------------
1558 !scalars
1559  integer,intent(in) :: mpert,natom,nsym,timrev,zero_by_symm
1560 !arrays
1561  integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym)
1562  integer,intent(in),target :: symrec(3,3,nsym),symrel(3,3,nsym)
1563  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
1564  real(dp),intent(in) :: qpt(3)
1565  real(dp),intent(inout) :: d2(2,3,mpert,3,mpert)
1566 
1567 !Local variables -------------------------
1568 !scalars
1569  logical, parameter :: do_final_sym=.true.
1570  logical :: qzero
1571  integer :: exch12,found,idir1,idir2,idisy1,idisy2,ipert1,ipert2
1572  integer :: ipesy1,ipesy2,isgn,isym,ithree,itirev,nblkflg_is_one,noccur,nsym_used,quit,quit1
1573  real(dp) :: arg1,arg2,im,norm,re,sumi,sumr,xi,xr
1574 !arrays
1575  integer,pointer :: sym1_(:,:,:),sym2_(:,:,:)
1576  real(dp),allocatable :: d2tmp1(:,:,:),d2tmp2(:,:,:),d2work(:,:,:,:,:)
1577 
1578 ! *********************************************************************
1579 
1580  qzero=(qpt(1)**2+qpt(2)**2+qpt(3)**2<tol16)
1581 
1582 !Here look after exchange of 1 and 2 axis,
1583 !for electric field in diamond symmetry
1584  exch12=0
1585  if (qzero) then
1586    do isym=1,nsym
1587      exch12=1
1588      if(symrel(1,1,isym)/=0)exch12=0
1589      if(symrel(1,2,isym)/=1)exch12=0
1590      if(symrel(1,3,isym)/=0)exch12=0
1591      if(symrel(2,1,isym)/=1)exch12=0
1592      if(symrel(2,2,isym)/=0)exch12=0
1593      if(symrel(2,3,isym)/=0)exch12=0
1594      if(symrel(3,1,isym)/=0)exch12=0
1595      if(symrel(3,2,isym)/=0)exch12=0
1596      if(symrel(3,3,isym)/=1)exch12=0
1597 !    if(exch12==1) write(std_out,*)' d2sym3 : found exchange 1 2 =',isym
1598      if(exch12==1)exit
1599    end do
1600  end if
1601 
1602 !Exchange of perturbations
1603 
1604 !Consider two cases : either time-reversal symmetry
1605 !conserves the wavevector, or not
1606  if(timrev==0)then
1607 
1608 !  do ipert1=1,mpert  See notes
1609    do ipert1=1,min(natom+2,mpert)
1610      do idir1=1,3
1611 
1612 !      Since the matrix is hermitian, the diagonal elements are real
1613        d2(2,idir1,ipert1,idir1,ipert1)=zero
1614 
1615 !      do ipert2=1,mpert See notes
1616        do ipert2=1,min(natom+2,mpert)
1617          do idir2=1,3
1618 
1619 !          If an element exists
1620            if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
1621 
1622 !            Either complete the symmetric missing element
1623              if(blkflg(idir2,ipert2,idir1,ipert1)==0)then
1624 
1625                d2(1,idir2,ipert2,idir1,ipert1)= d2(1,idir1,ipert1,idir2,ipert2)
1626                d2(2,idir2,ipert2,idir1,ipert1)=-d2(2,idir1,ipert1,idir2,ipert2)
1627 
1628                blkflg(idir2,ipert2,idir1,ipert1)=1
1629 
1630 !              Or symmetrize (the matrix is hermitian) in case both exists
1631 !              (Note : this opportunity has been disabled for more
1632 !              obvious search for bugs in the code )
1633 !              else
1634 !              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
1635 !              sumi=d2(1,idir2,ipert2,idir1,ipert1)-d2(1,idir1,ipert1,idir2,ipert2)
1636 !              d2(1,idir2,ipert2,idir1,ipert1)=half*sumr
1637 !              d2(1,idir1,ipert1,idir2,ipert2)=half*sumr
1638 !              d2(2,idir2,ipert2,idir1,ipert1)=half*sumi
1639 !              d2(2,idir1,ipert1,idir2,ipert2)=-half*sumi
1640              end if
1641            end if
1642 
1643          end do
1644        end do
1645 
1646      end do
1647    end do
1648 
1649 !  Here, case with time-reversal symmetry
1650  else
1651 
1652 !  do ipert1=1,mpert See notes
1653    do ipert1=1,min(natom+2,mpert)
1654      do idir1=1,3
1655 !      do ipert2=1,mpert See notes
1656        do ipert2=1,min(natom+2,mpert)
1657          do idir2=1,3
1658            d2(2,idir1,ipert1,idir2,ipert2)=zero
1659 
1660 !          If an element exists
1661            if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
1662 
1663 !            Either complete the symmetric missing element
1664              if(blkflg(idir2,ipert2,idir1,ipert1)==0)then
1665 
1666                d2(1,idir2,ipert2,idir1,ipert1)=d2(1,idir1,ipert1,idir2,ipert2)
1667                blkflg(idir2,ipert2,idir1,ipert1)=1
1668 
1669 !              Or symmetrize (the matrix is hermitian) in case both exists
1670 !              (Note : this opportunity has been disabled for more
1671 !              obvious search for bugs in the code )
1672 !              else
1673 !              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
1674 !              d2(1,idir2,ipert2,idir1,ipert1)=half*sumr
1675 !              d2(1,idir1,ipert1,idir2,ipert2)=half*sumr
1676              end if
1677 
1678            end if
1679          end do
1680        end do
1681      end do
1682    end do
1683  end if
1684 
1685 !Big Big Loop : symmetrize three times, because
1686 !of some cases in which one element is not yet available
1687 !at the first pass, and even at the second one !
1688  do ithree=1,3
1689 
1690 !  Big loop on all elements
1691 !  do ipert1=1,mpert See notes
1692    do ipert1=1,min(natom+2,mpert)
1693 
1694 !    Select the symmetries according to pertubation 1
1695      if (ipert1<=natom)then
1696        sym1_ => symrec
1697      else
1698        sym1_ => symrel
1699      end if
1700 
1701      do idir1=1,3
1702 !      do ipert2=1,mpert See notes
1703        do ipert2=1,min(natom+2,mpert)
1704 
1705     !    Select the symmetries according to pertubation 2
1706          if (ipert2<=natom)then
1707            sym2_ => symrec
1708          else
1709            sym2_ => symrel
1710          end if
1711 
1712          do idir2=1,3
1713 
1714 !          Will get element (idir1,ipert1,idir2,ipert2)
1715 !          so this element should not yet be present ...
1716            if(blkflg(idir1,ipert1,idir2,ipert2)/=1)then
1717 
1718              d2(1,idir1,ipert1,idir2,ipert2)=zero
1719              d2(2,idir1,ipert1,idir2,ipert2)=zero
1720 
1721 !            Loop on all symmetries, including time-reversal
1722              quit1=0
1723              do isym=1,nsym
1724                do itirev=1,2
1725                  isgn=3-2*itirev
1726 
1727                  if(symq(4,itirev,isym)/=0)then
1728                    found=1
1729 
1730 !                  Here select the symmetric of ipert1
1731                    if(ipert1<=natom)then
1732                      ipesy1=indsym(4,isym,ipert1)
1733                    else if(ipert1==(natom+2).and.qzero)then
1734                      ipesy1=ipert1
1735                    else
1736                      found=0
1737                    end if
1738 
1739 !                  Here select the symmetric of ipert2
1740                    if(ipert2<=natom)then
1741                      ipesy2=indsym(4,isym,ipert2)
1742                    else if(ipert2==(natom+2).and.qzero)then
1743                      ipesy2=ipert2
1744                    else
1745                      found=0
1746                    end if
1747 
1748 !                  Now that a symmetric perturbation has been obtained,
1749 !                  including the expression of the symmetry matrix, see
1750 !                  if the symmetric values are available
1751                    if( found==1 ) then
1752 
1753                      sumr=zero
1754                      sumi=zero
1755                      noccur=0
1756                      nblkflg_is_one=0
1757                      quit=0
1758                      do idisy1=1,3
1759                        do idisy2=1,3
1760                          if(sym1_(idir1,idisy1,isym)/=0 .and. sym2_(idir2,idisy2,isym)/=0 )then
1761                            if(blkflg(idisy1,ipesy1,idisy2,ipesy2)==1)then
1762                              nblkflg_is_one=nblkflg_is_one+1
1763                              sumr=sumr+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*&
1764 &                             d2(1,idisy1,ipesy1,idisy2,ipesy2)
1765                              sumi=sumi+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*&
1766 &                             d2(2,idisy1,ipesy1,idisy2,ipesy2)
1767 
1768 !                            Here, in case the symmetric of the element
1769 !                            is the element, or the symmetric with
1770 !                            respect to permutation of perturbations
1771 !                            (some more conditions on the time-reversal
1772 !                            symmetry must be fulfilled although)
1773                            else if(  idisy1==idir1 .and. ipesy1==ipert1&
1774 &                             .and. idisy2==idir2 .and. ipesy2==ipert2&
1775 &                             .and.(isgn==1 .or. timrev==1 &
1776 &                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
1777 &                             then
1778                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1779                            else if(  idisy1==idir2 .and. ipesy1==ipert2&
1780 &                             .and. idisy2==idir1 .and. ipesy2==ipert1&
1781 &                             .and.(isgn==-1 .or. timrev==1&
1782 &                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
1783 &                             then
1784                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1785 
1786 !                            Here, electric field case
1787                            else if( exch12==1 .and. &
1788 &                             ipert1==natom+2 .and. ipert2==natom+2&
1789 &                             .and.(( idisy1+idir1 ==3 &
1790 &                             .and. idisy2==3 .and. idir2==3)&
1791 &                             .or. ( idisy1+idir2 ==3&
1792 &                             .and. idisy2==3 .and. idir1==3)&
1793 &                             .or. ( idisy2+idir2 ==3&
1794 &                             .and. idisy1==3 .and. idir1==3)&
1795 &                             .or. ( idisy2+idir1 ==3&
1796 &                             .and. idisy1==3 .and. idir2==3)))&
1797 &                             then
1798                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1799 
1800                            else
1801 !                            Not found
1802                              found=0
1803                              quit=1
1804                              exit
1805                            end if
1806 
1807                          end if
1808                        end do
1809                        if(quit==1)exit
1810                      end do
1811                    end if
1812 
1813 !                  In case zero_by_symm==0, the computed materix element must be associated to at least one really computed matrix element
1814                    if(zero_by_symm==0 .and. nblkflg_is_one==0)then
1815                      found=0
1816                    endif
1817 
1818 !                  Now, if still found and associated to at least one really computed matrix element, put the correct value into array d2
1819                    if(found==1)then
1820 
1821 !                    In case of phonons, need to take into account the
1822 !                    time-reversal symmetry, and the shift back to the unit cell
1823 !
1824 !                    XG990712 : I am not sure this must be kept for electric field ...
1825 !                    1) Consider time-reversal symmetry
1826                      sumi=isgn*sumi
1827 
1828                      if(ipert1<=natom .and. ipert2<=natom)then
1829 !                      2) Shift the atoms back to the unit cell.
1830                        arg1=two_pi*( qpt(1)*indsym(1,isym,ipert1)&
1831 &                       +qpt(2)*indsym(2,isym,ipert1)&
1832 &                       +qpt(3)*indsym(3,isym,ipert1) )
1833                        arg2=two_pi*( qpt(1)*indsym(1,isym,ipert2)&
1834 &                       +qpt(2)*indsym(2,isym,ipert2)&
1835 &                       +qpt(3)*indsym(3,isym,ipert2) )
1836                        re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
1837 !                      XG010117 Must use isgn
1838                        im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))
1839                      else
1840                        re=one
1841                        im=zero
1842                      end if
1843 
1844 !                    Final check, could still fail if the
1845 !                    element was its own symmetric
1846                      if (abs(1.0_dp-re*noccur)< 1.0d-6.and.abs(im*noccur) <1.0d-6) then
1847                        found=0
1848                      end if
1849 
1850                    end if
1851 
1852                    if(found==1)then
1853 
1854                      if(noccur==0)then
1855                        d2(1,idir1,ipert1,idir2,ipert2)=re*sumr-im*sumi
1856                        d2(2,idir1,ipert1,idir2,ipert2)=re*sumi+im*sumr
1857                      else
1858 !                      See page July 2, 1994 in computer codes notebook
1859                        xr=re*sumr-im*sumi
1860                        xi=re*sumi+im*sumr
1861                        norm=one+noccur**2-two*re*noccur
1862                        xr=xr/norm
1863                        xi=xi/norm
1864                        d2(1,idir1,ipert1,idir2,ipert2)=&
1865 &                       (one-re*noccur)*xr-im*noccur*xi
1866                        d2(2,idir1,ipert1,idir2,ipert2)=&
1867 &                       (one-re*noccur)*xi+im*noccur*xr
1868                      end if
1869 
1870 !                    The element has been constructed !
1871                      blkflg(idir1,ipert1,idir2,ipert2)=1
1872 
1873                      quit1=1
1874                      exit ! Exit loop on symmetry operations
1875                    end if
1876 
1877 !                  End loop on all symmetries + time-reversal
1878                  end if
1879                end do
1880                if(quit1==1)exit
1881              end do
1882 
1883            end if
1884          end do ! End big loop on all elements
1885        end do
1886      end do
1887    end do
1888 
1889  end do !  End Big Big Loop
1890 
1891 !MT oct. 20, 2014:
1892 !Once the matrix has been built, it does not necessarily fulfill the correct symmetries.
1893 !It has just been filled up from rows or columns that only fulfill symmetries preserving
1894 !one particular perturbation.
1895 !An additional symmetrization might solve this (do not consider TR-symmetry)
1896  if (do_final_sym) then
1897    ABI_ALLOCATE(d2tmp1,(2,3,3))
1898    ABI_ALLOCATE(d2tmp2,(2,3,3))
1899    ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert))
1900    d2work(:,:,:,:,:)=d2(:,:,:,:,:)
1901    do ipert1=1,min(natom+2,mpert)
1902      if ((ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).or.(ipert1==natom+2.and.(.not.qzero))) cycle
1903      if (ipert1<=natom)then
1904        sym1_ => symrec
1905      else
1906        sym1_ => symrel
1907      end if
1908      do ipert2=1,min(natom+2,mpert)
1909 !      if (any(blkflg(:,ipert1,:,ipert2)==0)) cycle
1910        if ((ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11).or.(ipert2==natom+2.and.(.not.qzero))) cycle
1911        if (ipert2<=natom)then
1912          sym2_ => symrec
1913        else
1914          sym2_ => symrel
1915        end if
1916        nsym_used=0
1917        d2tmp2(:,:,:)=zero
1918        do isym=1,nsym
1919          if (symq(4,1,isym)==1) then
1920            ipesy1=ipert1;if (ipert1<=natom) ipesy1=indsym(4,isym,ipert1)
1921            ipesy2=ipert2;if (ipert2<=natom) ipesy2=indsym(4,isym,ipert2)
1922 !          The condition on next line is too severe, since some elements of sym1_ or sym2_ might be zero,
1923 !          which means not all blkflg(:,ipesy1,:,ipesy2) would need to be 1 to symmetrize the matrix.
1924 !          However, coding something more refined is really more difficult.
1925 !          This condition then has the side effect that more symmetries can be applied when zero_by_symm==1,
1926 !          since blkflg can be set to 1 when the symmetries guarantee the matrix element to be zero.
1927            if (all(blkflg(:,ipesy1,:,ipesy2)==1)) then
1928              nsym_used=nsym_used+1
1929              re=one;im=zero
1930              if (ipert1<=natom.and.ipert2<=natom.and.(.not.qzero)) then
1931                arg1=two_pi*(qpt(1)*(indsym(1,isym,ipert1)-indsym(1,isym,ipert2)) &
1932 &                          +qpt(2)*(indsym(2,isym,ipert1)-indsym(2,isym,ipert2)) &
1933 &                          +qpt(3)*(indsym(3,isym,ipert1)-indsym(3,isym,ipert2)))
1934                re=cos(arg1);im=sin(arg1)
1935              end if
1936              d2tmp1(:,:,:)=zero
1937              do idir2=1,3 !kappa
1938                do idir1=1,3 !mu
1939                  do idisy1=1,3 !nu
1940                    d2tmp1(:,idir1,idir2)=d2tmp1(:,idir1,idir2) &
1941 &                     +sym1_(idir1,idisy1,isym)*d2(:,idisy1,ipesy1,idir2,ipesy2)
1942                  end do
1943                end do
1944              end do
1945              do idir2=1,3 !mu
1946                do idir1=1,3 !kappa
1947                  do idisy2=1,3 !nu
1948                    d2tmp2(1,idir1,idir2)=d2tmp2(1,idir1,idir2) &
1949 &                  +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*re-d2tmp1(2,idir1,idisy2)*im)
1950                    d2tmp2(2,idir1,idir2)=d2tmp2(2,idir1,idir2) &
1951 &                  +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*im+d2tmp1(2,idir1,idisy2)*re)
1952                  end do
1953                end do
1954              end do
1955            end if
1956          end if
1957        end do ! isym
1958        if (nsym_used>0) d2work(:,1:3,ipert1,1:3,ipert2)=d2tmp2(:,1:3,1:3)/dble(nsym_used)
1959      end do !ipert2
1960    end do !ipert1
1961    if (mpert>=natom)   d2(:,1:3,1:natom,1:3,1:natom)=d2work(:,1:3,1:natom,1:3,1:natom)
1962    if (mpert>=natom+2) then
1963      d2(:,1:3,natom+2,1:3,1:natom)=d2work(:,1:3,natom+2,1:3,1:natom)
1964      d2(:,1:3,1:natom,1:3,natom+2)=d2work(:,1:3,1:natom,1:3,natom+2)
1965      d2(:,1:3,natom+2,1:3,natom+2)=d2work(:,1:3,natom+2,1:3,natom+2)
1966    end if
1967    ABI_DEALLOCATE(d2tmp1)
1968    ABI_DEALLOCATE(d2tmp2)
1969    ABI_DEALLOCATE(d2work)
1970  end if
1971 
1972 end subroutine d2sym3

m_dynmat/d3sym [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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,nonlinear

CHILDREN

SOURCE

4568 subroutine d3sym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel)
4569 
4570 
4571 !This section has been created automatically by the script Abilint (TD).
4572 !Do not modify the following lines by hand.
4573 #undef ABI_FUNC
4574 #define ABI_FUNC 'd3sym'
4575 !End of the abilint section
4576 
4577  implicit none
4578 
4579 !Arguments -------------------------------
4580 !scalars
4581  integer,intent(in) :: mpert,natom,nsym
4582 !arrays
4583  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4584  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
4585  real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert)
4586 
4587 !Local variables -------------------------
4588 !scalars
4589  integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3
4590  integer :: ipesy1,ipesy2,ipesy3,isym,ithree
4591  real(dp) :: sumi,sumr
4592 !arrays
4593  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
4594 
4595 ! *********************************************************************
4596 
4597 !DEBUG
4598 !write(std_out,*)'d3sym : enter'
4599 !do i1dir = 1, 3
4600 !do i2dir = 1, 3
4601 !do i3dir = 1, 3
4602 !write(std_out,*)i1dir,i2dir,i3dir,blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
4603 !end do
4604 !end do
4605 !end do
4606 !stop
4607 !ENDDEBUG
4608 
4609 !First, take into account the permutations symmetry of
4610 !(i1pert,i1dir) and (i3pert,i3dir)
4611  do i1pert = 1, mpert
4612    do i2pert = 1, mpert
4613      do i3pert = 1, mpert
4614 
4615        do i1dir = 1, 3
4616          do i2dir = 1, 3
4617            do i3dir = 1, 3
4618 
4619              if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. &
4620 &             (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=1)) then
4621 
4622                d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = &
4623 &              d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
4624 
4625                blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = 1
4626 
4627              end if
4628 
4629            end do
4630          end do
4631        end do
4632 
4633      end do
4634    end do
4635  end do
4636 
4637 !Big Big Loop : symmetrize three times, because
4638 !of some cases in which one element is not yet available
4639 !at the first pass, and even at the second one !
4640 
4641  do ithree=1,3
4642 
4643 !  Loop over perturbations
4644    do i1pert = 1, mpert
4645      do i2pert = 1, mpert
4646        do i3pert = 1, mpert
4647 
4648          do i1dir = 1, 3
4649            do i2dir = 1, 3
4650              do i3dir = 1, 3
4651 
4652 !              Will get element (idir1,ipert1,idir2,ipert2)
4653 !              so this element should not yet be present ...
4654                if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then
4655 
4656                  d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp
4657 
4658                  do isym = 1, nsym
4659 
4660                    found = 1
4661 
4662                    if (i1pert <= natom) then
4663                      ipesy1 = indsym(4,isym,i1pert)
4664                      sym1(:,:) = symrec(:,:,isym)
4665                    else if (i1pert == natom + 2) then
4666                      ipesy1 = i1pert
4667                      sym1(:,:) = symrel(:,:,isym)
4668                    else
4669                      found = 0
4670                    end if
4671 
4672                    if (i2pert <= natom) then
4673                      ipesy2 = indsym(4,isym,i2pert)
4674                      sym2(:,:) = symrec(:,:,isym)
4675                    else if (i2pert == natom + 2) then
4676                      ipesy2 = i2pert
4677                      sym2(:,:) = symrel(:,:,isym)
4678                    else
4679                      found = 0
4680                    end if
4681 
4682                    if (i3pert <= natom) then
4683                      ipesy3 = indsym(4,isym,i3pert)
4684                      sym3(:,:) = symrec(:,:,isym)
4685                    else if (i3pert == natom + 2) then
4686                      ipesy3 = i3pert
4687                      sym3(:,:) = symrel(:,:,isym)
4688                    else
4689                      found = 0
4690                    end if
4691 
4692                    sumr = 0_dp ; sumi = 0_dp;
4693                    do idisy1 = 1, 3
4694                      do idisy2 = 1, 3
4695                        do idisy3 = 1, 3
4696 
4697                          if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) &
4698 &                         .and.(sym3(i3dir,idisy3) /=0)) then
4699 
4700                            if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then
4701 
4702                              sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4703 &                             sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4704                              sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4705 &                             sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4706 
4707                            else
4708 
4709                              found = 0
4710 
4711                            end if
4712 
4713                          end if
4714 
4715                        end do
4716                      end do
4717                    end do
4718 
4719                    if (found == 1) then
4720                      d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
4721                      d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
4722                      blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4723                    end if
4724 
4725                  end do  ! isym
4726 
4727                end if  ! blkflg
4728 
4729 !              Close loop over perturbations
4730              end do
4731            end do
4732          end do
4733        end do
4734      end do
4735    end do
4736 
4737  end do  ! close loop over ithree
4738 
4739 end subroutine d3sym

m_dynmat/dfpt_phfrq [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 dfpt_phfrq

FUNCTION

 Get the phonon frequencies and eigenvectors (as well as the corresponding displacements)
 If q is at Gamma, the non-analytical behaviour can be included.
 Then, the effective dielectric tensor, the effective charges
 and oscillator strengths for the limiting direction are also returned

INPUTS

  amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms)
  d2cart(2,3,mpert,3,mpert)=dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates
  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.
  mpert =maximum number of ipert
  msym=maximum number of symmetries
  natom=number of atoms in unit cell
  nsym=number of space group symmetries
  ntypat=number of atom types
  qphnrm=(described above)
  qphon(3)= to be divided by qphnrm, give the phonon wavevector;
     if qphnrm==0.0_dp, then the wavevector is zero (Gamma point)
     and qphon gives the direction of the induced electric field in **CARTESIAN** coordinates.
     in the latter case, if qphon is zero, no non-analytical contribution is included.
  rprimd(3,3)=dimensional primitive translations (bohr)
  symdynmat=if 1, (re)symmetrize the dynamical matrix, except if Gamma wavevector with electric field added.
  symrel(3,3,nsym)=matrices of the group symmetries (real space)
  typat(natom)=integer label of each type of atom (1,2,...)
  ucvol=unit cell volume

OUTPUT

  displ(2*3*natom*3*natom)= at the end, contains the displacements of atoms in cartesian coordinates.
    The first index means either the real or the imaginary part,
    The second index runs on the direction and the atoms displaced
    The third index runs on the modes.
  eigval(3*natom)=contains the eigenvalues of the dynamical matrix
  eigvec(2*3*natom*3*natom)= at the end, contains the eigenvectors of the dynamical matrix in cartesian coordinates.
  phfrq(3*natom)=phonon frequencies (square root of the dynamical matrix eigenvalues,
    except if these are negative, and in this case, give minus the square root of the absolute value
    of the matrix eigenvalues). Hartree units.

NOTES

   1) One makes the dynamical matrix hermitian...
   2) In case of q=Gamma, only the real part is used.

PARENTS

      anaddb,m_ddb,m_effective_potential_file,m_ifc,m_phonons,respfn,thmeig

CHILDREN

SOURCE

5831 subroutine dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,&
5832 & mpert,msym,natom,nsym,ntypat,phfrq,qphnrm,qphon,rprimd,&
5833 & symdynmat,symrel,symafm,typat,ucvol)
5834 
5835 
5836 !This section has been created automatically by the script Abilint (TD).
5837 !Do not modify the following lines by hand.
5838 #undef ABI_FUNC
5839 #define ABI_FUNC 'dfpt_phfrq'
5840 !End of the abilint section
5841 
5842  implicit none
5843 
5844 !Arguments -------------------------------
5845 !scalars
5846  integer,intent(in) :: mpert,msym,natom,nsym,ntypat,symdynmat
5847  real(dp),intent(in) :: qphnrm,ucvol
5848 !arrays
5849  integer,intent(in) :: indsym(4,msym*natom),symrel(3,3,nsym),typat(natom)
5850  integer,intent(in) :: symafm(nsym)
5851  real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),rprimd(3,3)
5852  real(dp),intent(inout) :: qphon(3)
5853  real(dp),intent(out) :: displ(2*3*natom*3*natom),eigval(3*natom)
5854  real(dp),intent(out) :: eigvec(2*3*natom*3*natom),phfrq(3*natom)
5855 
5856 !Local variables -------------------------
5857 !scalars
5858  integer :: analyt,i1,i2,idir1,idir2,ier,ii,imode,ipert1,ipert2
5859  integer :: jmode,indexi,indexj,index
5860  real(dp) :: epsq,norm,qphon2
5861  logical,parameter :: debug = .False.
5862  real(dp) :: sc_prod
5863 !arrays
5864  real(dp) :: qptn(3),dum(2,0)
5865  real(dp),allocatable :: matrx(:,:),zeff(:,:),zhpev1(:,:),zhpev2(:)
5866 
5867 ! *********************************************************************
5868 
5869 !Prepare the diagonalisation: analytical part.
5870 !Note: displ is used as work space here
5871  i1=0
5872  do ipert1=1,natom
5873    do idir1=1,3
5874      i1=i1+1
5875      i2=0
5876      do ipert2=1,natom
5877        do idir2=1,3
5878          i2=i2+1
5879          index=i1+3*natom*(i2-1)
5880          displ(2*index-1)=d2cart(1,idir1,ipert1,idir2,ipert2)
5881          displ(2*index  )=d2cart(2,idir1,ipert1,idir2,ipert2)
5882        end do
5883      end do
5884    end do
5885  end do
5886 
5887 !Determine the analyticity of the matrix.
5888  analyt=1; if(abs(qphnrm)<tol8) analyt=0
5889  if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=2
5890 
5891 !In case of q=Gamma, only the real part is used
5892  if(analyt==0 .or. analyt==2)then
5893    do i1=1,3*natom
5894      do i2=1,3*natom
5895        index=i1+3*natom*(i2-1)
5896        displ(2*index)=zero
5897      end do
5898    end do
5899  end if
5900 
5901 !In the case the non-analyticity is required :
5902 ! MG: the tensor is in cartesian coordinates and this means that qphon must be in
5903 !     given in Cartesian coordinates.
5904  if(analyt==0)then
5905 
5906 !  Normalize the limiting direction
5907    qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2
5908    qphon(:)=qphon(:)/sqrt(qphon2)
5909 
5910 !  Get the dielectric constant for the limiting direction
5911    epsq=zero
5912    do idir1=1,3
5913      do idir2=1,3
5914        epsq= epsq + qphon(idir1)*qphon(idir2) * d2cart(1,idir1,natom+2,idir2,natom+2)
5915      end do
5916    end do
5917 
5918    ABI_ALLOCATE(zeff,(3,natom))
5919 
5920 !  Get the effective charges for the limiting direction
5921    do idir1=1,3
5922      do ipert1=1,natom
5923        zeff(idir1,ipert1)=zero
5924        do idir2=1,3
5925          zeff(idir1,ipert1) = zeff(idir1,ipert1) + qphon(idir2)* d2cart(1,idir1,ipert1,idir2,natom+2)
5926        end do
5927      end do
5928    end do
5929 
5930 !  Get the non-analytical part of the dynamical matrix, and suppress its imaginary part.
5931    i1=0
5932    do ipert1=1,natom
5933      do idir1=1,3
5934        i1=i1+1
5935        i2=0
5936        do ipert2=1,natom
5937          do idir2=1,3
5938            i2=i2+1
5939            index=i1+3*natom*(i2-1)
5940            displ(2*index-1)=displ(2*index-1)+four_pi/ucvol*zeff(idir1,ipert1)*zeff(idir2,ipert2)/epsq
5941            displ(2*index  )=zero
5942          end do
5943        end do
5944      end do
5945    end do
5946 
5947    ABI_DEALLOCATE(zeff)
5948  end if !  End of the non-analyticity treatment :
5949 
5950  ! Multiply IFC(q) by masses
5951  call massmult_and_breaksym(natom, ntypat, typat, amu, displ)
5952 
5953 !***********************************************************************
5954 !Diagonalize the dynamical matrix
5955 
5956 !Symmetrize the dynamical matrix
5957 !FIXME: swap the next 2 lines and update test files to include symmetrization for Gamma point too (except in non-analytic case)
5958 !if (symdynmat==1 .and. analyt > 0) then
5959  if (symdynmat==1 .and. analyt == 1) then
5960    qptn(:)=qphon(:)
5961    if (analyt==1) qptn(:)=qphon(:)/qphnrm
5962    call symdyma(displ,indsym,natom,nsym,qptn,rprimd,symrel,symafm)
5963  end if
5964 
5965  ier=0; ii=1
5966  ABI_ALLOCATE(matrx,(2,(3*natom*(3*natom+1))/2))
5967  do i2=1,3*natom
5968    do i1=1,i2
5969      matrx(1,ii)=displ(1+2*(i1-1)+2*(i2-1)*3*natom)
5970      matrx(2,ii)=displ(2+2*(i1-1)+2*(i2-1)*3*natom)
5971      ii=ii+1
5972    end do
5973  end do
5974 
5975  ABI_ALLOCATE(zhpev1,(2,2*3*natom-1))
5976  ABI_ALLOCATE(zhpev2,(3*3*natom-2))
5977 
5978  call ZHPEV ('V','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,zhpev2,ier)
5979  ABI_CHECK(ier == 0, sjoin('zhpev returned:', itoa(ier)))
5980 
5981  ABI_DEALLOCATE(matrx)
5982  ABI_DEALLOCATE(zhpev1)
5983  ABI_DEALLOCATE(zhpev2)
5984 
5985  if (debug) then
5986    ! Check the orthonormality of the eigenvectors
5987    do imode=1,3*natom
5988      do jmode=imode,3*natom
5989        indexi=2*3*natom*(imode-1)
5990        indexj=2*3*natom*(jmode-1)
5991        sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom))
5992        write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod
5993      end do
5994    end do
5995  end if
5996 
5997 !***********************************************************************
5998 
5999 !Get the phonon frequencies (negative by convention, if
6000 !the eigenvalue of the dynamical matrix is negative)
6001  do imode=1,3*natom
6002    if(eigval(imode)>=1.0d-16)then
6003      phfrq(imode)=sqrt(eigval(imode))
6004    else if(eigval(imode)>=-1.0d-16)then
6005      phfrq(imode)=zero
6006    else
6007      phfrq(imode)=-sqrt(-eigval(imode))
6008    end if
6009  end do
6010 
6011 !Fix the phase of the eigenvectors
6012  call fxphas_seq(eigvec,dum,0,0,1,3*natom*3*natom,0,3*natom,3*natom,0)
6013 
6014 !Normalise the eigenvectors
6015  do imode=1,3*natom
6016    norm=zero
6017    do idir1=1,3
6018      do ipert1=1,natom
6019        i1=idir1+(ipert1-1)*3
6020        index=i1+3*natom*(imode-1)
6021        norm=norm+eigvec(2*index-1)**2+eigvec(2*index)**2
6022      end do
6023    end do
6024    norm=sqrt(norm)
6025    do idir1=1,3
6026      do ipert1=1,natom
6027        i1=idir1+(ipert1-1)*3
6028        index=i1+3*natom*(imode-1)
6029        eigvec(2*index-1)=eigvec(2*index-1)/norm
6030        eigvec(2*index)=eigvec(2*index)/norm
6031      end do
6032    end do
6033  end do
6034 
6035 !Get the phonon displacements
6036  do imode=1,3*natom
6037    do idir1=1,3
6038      do ipert1=1,natom
6039        i1=idir1+(ipert1-1)*3
6040        index=i1+3*natom*(imode-1)
6041        displ(2*index-1)=eigvec(2*index-1) / sqrt(amu(typat(ipert1))*amu_emass)
6042        displ(2*index  )=eigvec(2*index  ) / sqrt(amu(typat(ipert1))*amu_emass)
6043      end do
6044    end do
6045  end do
6046 
6047  if (debug) then
6048    write(std_out,'(a)')' Phonon eigenvectors and displacements '
6049    do imode=1,3*natom
6050      indexi=2*3*natom*(imode-1)
6051      write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' eigvec(1:6*natom)=',eigvec(indexi+1:indexi+6*natom)
6052      write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' displ(1:6*natom)=',displ(indexi+1:indexi+6*natom)
6053    end do
6054 
6055    ! Check the orthonormality of the eigenvectors
6056    do imode=1,3*natom
6057      do jmode=imode,3*natom
6058        indexi=2*3*natom*(imode-1)
6059        indexj=2*3*natom*(jmode-1)
6060        sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom))
6061        write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod
6062      end do
6063    end do
6064  end if
6065 
6066 end subroutine dfpt_phfrq

m_dynmat/dfpt_prtph [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 dfpt_prtph

FUNCTION

 Print the phonon frequencies, on unit 6 as well as the printing
 unit (except if the associated number -iout- is negative),
 and for the latter, in Hartree, meV, Thz, Kelvin or cm-1.
 If eivec==1,2, also print the eigenmodes : displacements in cartesian coordinates.
 If eivec==4, generate output files for band2eps (drawing tool for the phonon band structure

INPUTS

  displ(2,3*natom,3*natom)= contains the displacements of atoms in cartesian coordinates.
  The first index means either the real or the imaginary part,
  The second index runs on the direction and the atoms displaced
  The third index runs on the modes.
  eivec=(if eivec==0, the eigendisplacements are not printed,
    if eivec==1,2, the eigendisplacements are printed,
    if eivec==4, files for band2eps
  enunit=units for output of the phonon frequencies :
    0=> Hartree and cm-1, 1=> eV and Thz, other=> Ha,Thz,eV,cm-1 and K
  iout= unit for long print (if negative, the routine only print on unit 6, and in Hartree only).
  natom= number of atom
  phfreq(3*natom)= phonon frequencies in Hartree
  qphnrm=phonon wavevector normalisation factor
  qphon(3)=phonon wavevector

OUTPUT

  Only printing

NOTES

 called by one processor only

PARENTS

      anaddb,m_effective_potential_file,m_ifc,m_phonons,respfn

CHILDREN

      wrtout

SOURCE

6110 subroutine dfpt_prtph(displ,eivec,enunit,iout,natom,phfrq,qphnrm,qphon)
6111 
6112 
6113 !This section has been created automatically by the script Abilint (TD).
6114 !Do not modify the following lines by hand.
6115 #undef ABI_FUNC
6116 #define ABI_FUNC 'dfpt_prtph'
6117 !End of the abilint section
6118 
6119  implicit none
6120 
6121 !Arguments -------------------------------
6122 !scalars
6123  integer,intent(in) :: eivec,enunit,iout,natom
6124  real(dp),intent(in) :: qphnrm
6125 !arrays
6126  real(dp),intent(in) :: displ(2,3*natom,3*natom),phfrq(3*natom),qphon(3)
6127 
6128 !Local variables -------------------------
6129 !scalars
6130  integer :: i,idir,ii,imode,jj
6131  real(dp) :: tolerance
6132  logical :: t_degenerate
6133  character(len=500) :: message
6134 !arrays
6135  real(dp) :: vecti(3),vectr(3)
6136  character(len=1) :: metacharacter(3*natom)
6137 
6138 ! *********************************************************************
6139 
6140 !Check the value of eivec
6141  if (all(eivec /= [0,1,2,4])) then
6142    write(message, '(a,i0,a,a)' )&
6143 &   'In the calling subroutine, eivec is',eivec,ch10,&
6144 &   'but allowed values are between 0 and 4.'
6145    MSG_BUG(message)
6146  end if
6147 
6148 !write the phonon frequencies on unit std_out
6149  write(message,'(4a)' )' ',ch10,&
6150 & ' phonon wavelength (reduced coordinates) , ','norm, and energies in hartree'
6151  call wrtout(std_out,message,'COLL')
6152 
6153 !The next format should be rewritten
6154  write(message,'(a,4f5.2)' )' ',(qphon(i),i=1,3),qphnrm
6155  call wrtout(std_out,message,'COLL')
6156  do jj=1,3*natom,5
6157    if (3*natom-jj<5) then
6158      write(message,'(5es17.9)') (phfrq(ii),ii=jj,3*natom)
6159    else
6160      write(message,'(5es17.9)') (phfrq(ii),ii=jj,jj+4)
6161    end if
6162    call wrtout(std_out,message,'COLL')
6163  end do
6164  write(message,'(a,es17.9)')' Zero Point Motion energy (sum of freqs/2)=',sum(phfrq(1:3*natom))/2
6165  call wrtout(std_out,message,'COLL')
6166 
6167 !Put the wavevector in nice format
6168  if(iout>=0)then
6169    call wrtout(iout,' ','COLL')
6170    if(qphnrm/=0.0_dp)then
6171      write(message, '(a,3f9.5)' )&
6172 &     '  Phonon wavevector (reduced coordinates) :',(qphon(i)/qphnrm+tol10,i=1,3)
6173    else
6174      write(message, '(3a,3f9.5)' )&
6175 &     '  Phonon at Gamma, with non-analyticity in the',ch10,&
6176 &     '  direction (cartesian coordinates)',qphon(1:3)+tol10
6177    end if
6178    call wrtout(iout,message,'COLL')
6179 
6180 !  Write it, in different units.
6181    if(enunit/=1)then
6182      write(iout, '(a)' )' Phonon energies in Hartree :'
6183      do jj=1,3*natom,5
6184        if (3*natom-jj<5) then
6185          write(message, '(1x,5es14.6)') (phfrq(ii),ii=jj,3*natom)
6186        else
6187          write(message, '(1x,5es14.6)') (phfrq(ii),ii=jj,jj+4)
6188        end if
6189        call wrtout(iout,message,'COLL')
6190      end do
6191    end if
6192    if(enunit/=0)then
6193      write(iout, '(a)' )' Phonon energies in meV     :'
6194      do jj=1,3*natom,5
6195        if (3*natom-jj<5) then
6196          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,3*natom)
6197        else
6198          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,jj+4)
6199        end if
6200        call wrtout(iout,message,'COLL')
6201      end do
6202    end if
6203    if(enunit/=1)then
6204      write(iout, '(a)' )' Phonon frequencies in cm-1    :'
6205      do jj=1,3*natom,5
6206        if (3*natom-jj<5) then
6207          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,3*natom)
6208        else
6209          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,jj+4)
6210        end if
6211        call wrtout(iout,message,'COLL')
6212      end do
6213    end if
6214    if(enunit/=0)then
6215      write(iout, '(a)' )' Phonon frequencies in Thz     :'
6216      do jj=1,3*natom,5
6217        if (3*natom-jj<5) then
6218          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,3*natom)
6219        else
6220          write(message, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,jj+4)
6221        end if
6222        call wrtout(iout,message,'COLL')
6223      end do
6224    end if
6225    if(enunit/=0.and.enunit/=1)then
6226      write(iout, '(a)' )' Phonon energies in Kelvin  :'
6227      do jj=1,3*natom,5
6228        if (3*natom-jj<5) then
6229          write(message, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,3*natom)
6230        else
6231          write(message, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,jj+4)
6232        end if
6233        call wrtout(iout,message,'COLL')
6234      end do
6235    end if
6236  end if
6237 
6238 !Take care of the eigendisplacements
6239  if(eivec==1 .or. eivec==2)then
6240    write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
6241 &   ' Eigendisplacements ',ch10,&
6242 &   ' (will be given, for each mode : in cartesian coordinates',ch10,&
6243 &   '   for each atom the real part of the displacement vector,',ch10,&
6244 &   '   then the imaginary part of the displacement vector - absolute values smaller than 1.0d-7 are set to zero)'
6245    call wrtout(std_out,message,'COLL')
6246    if(iout>=0) then
6247      call wrtout(iout,message,'COLL')
6248    end if
6249 
6250 !  Examine the degeneracy of each mode. The portability of the echo of the eigendisplacements
6251 !  is very hard to obtain, and has not been attempted.
6252    do imode=1,3*natom
6253 !    The degenerate modes are not portable
6254      t_degenerate=.false.
6255      if(imode>1)then
6256        if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
6257      end if
6258      if(imode<3*natom)then
6259        if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
6260      end if
6261      metacharacter(imode)=';'; if(t_degenerate)metacharacter(imode)='-'
6262    end do
6263 
6264    do imode=1,3*natom
6265      write(message,'(a,i4,a,es16.6)' )'  Mode number ',imode,'   Energy',phfrq(imode)
6266      call wrtout(std_out,message,'COLL')
6267      if(iout>=0)then
6268        write(message, '(a,i4,a,es16.6)' )'  Mode number ',imode,'   Energy',phfrq(imode)
6269        call wrtout(iout,message,'COLL')
6270      end if
6271      tolerance=1.0d-7
6272      if(abs(phfrq(imode))<1.0d-5)tolerance=2.0d-7
6273      if(phfrq(imode)<1.0d-5)then
6274        write(message,'(3a)' )' Attention : low frequency mode.',ch10,&
6275 &       '   (Could be unstable or acoustic mode)'
6276        call wrtout(std_out,message,'COLL')
6277        if(iout>=0)then
6278          write(iout, '(3a)' )' Attention : low frequency mode.',ch10,&
6279 &         '   (Could be unstable or acoustic mode)'
6280        end if
6281      end if
6282      do ii=1,natom
6283        do idir=1,3
6284          vectr(idir)=displ(1,idir+(ii-1)*3,imode)
6285          if(abs(vectr(idir))<tolerance)vectr(idir)=0.0_dp
6286          vecti(idir)=displ(2,idir+(ii-1)*3,imode)
6287          if(abs(vecti(idir))<tolerance)vecti(idir)=0.0_dp
6288        end do
6289        write(message,'(i4,3es16.8,a,4x,3es16.8)' ) ii,vectr(:),ch10,vecti(:)
6290        call wrtout(std_out,message,'COLL')
6291        if(iout>=0)then
6292          write(message,'(a,i3,3es16.8,2a,3x,3es16.8)') metacharacter(imode),ii,vectr(:),ch10,&
6293 &         metacharacter(imode),   vecti(:)
6294          call wrtout(iout,message,'COLL')
6295        end if
6296      end do
6297    end do
6298  end if
6299 
6300 end subroutine dfpt_prtph

m_dynmat/dfpt_sydy [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 dfpt_sydy

FUNCTION

 Symmetrize dynamical matrix (eventually diagonal wrt to the atoms)
 Unsymmetrized dynamical matrix is input as dyfrow;
 symmetrized dynamical matrix is then placed in sdyfro.
 If nsym=1 simply copy dyfrow into sdyfro.

INPUTS

  cplex=1 if dynamical matrix is real, 2 if it is complex
  dyfrow(3,3,natom,1+(natom-1)*nondiag)=unsymmetrized dynamical matrix
  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 cell.
  nondiag=0 if dynamical matrix is     diagonal with respect to atoms
  nsym=number of symmetry operators in group.
  qphon(3)=wavevector of the phonon
  symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q
  symrec(3,3,nsym)=symmetries of group in terms of operations on real
    space primitive translations (integers).

OUTPUT

  sdyfro(3,3,natom,1+(natom-1)*nondiag)=symmetrized dynamical matrix

NOTES

 Symmetrization of gradients with respect to reduced
 coordinates tn is conducted according to the expression
 $[d(e)/d(t(n,a))]_{symmetrized} = (1/Nsym)*Sum(S)*symrec(n,m,S)*
              [d(e)/d(t(m,b))]_{unsymmetrized}$
 where $t(m,b)= (symrel^{-1})(m,n)*(t(n,a)-tnons(n))$ and tnons
 is a possible nonsymmorphic translation.  The label "b" here
 refers to the atom which gets rotated into "a" under symmetry "S".
 symrel is the symmetry matrix in real space, which is the inverse
 transpose of symrec.  symrec is the symmetry matrix in reciprocal
 space.  $sym_{cartesian} = R * symrel * R^{-1} = G * symrec * G^{-1}$
 where the columns of R and G are the dimensional primitive translations
 in real and reciprocal space respectively.
 Note the use of "symrec" in the symmetrization expression above.

PARENTS

      dfpt_dyfro

CHILDREN

SOURCE

2554 subroutine dfpt_sydy(cplex,dyfrow,indsym,natom,nondiag,nsym,qphon,sdyfro,symq,symrec)
2555 
2556 
2557 !This section has been created automatically by the script Abilint (TD).
2558 !Do not modify the following lines by hand.
2559 #undef ABI_FUNC
2560 #define ABI_FUNC 'dfpt_sydy'
2561 !End of the abilint section
2562 
2563  implicit none
2564 
2565 !Arguments -------------------------------
2566 !scalars
2567  integer,intent(in) :: cplex,natom,nondiag,nsym
2568 !arrays
2569  integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym),symrec(3,3,nsym)
2570  real(dp),intent(in) :: dyfrow(cplex,3,3,natom,1+(natom-1)*nondiag),qphon(3)
2571  real(dp),intent(out) :: sdyfro(cplex,3,3,natom,1+(natom-1)*nondiag)
2572 
2573 !Local variables -------------------------
2574 !scalars
2575  integer :: ia,indi,indj,isym,ja,kappa,mu,natom_nondiag,nsym_used,nu
2576  logical :: qeq0
2577  real(dp) :: arg,div,phasei,phaser
2578 !arrays
2579  real(dp) :: work(cplex,3,3)
2580 
2581 ! *********************************************************************
2582 
2583  if (nsym==1) then
2584 
2585 !  Only symmetry is identity so simply copy
2586    sdyfro(:,:,:,:,:)=dyfrow(:,:,:,:,:)
2587 
2588  else
2589 
2590 !  Actually carry out symmetrization
2591    sdyfro(:,:,:,:,:)=zero
2592    qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)
2593 !  === Diagonal dyn. matrix OR q=0
2594    if (nondiag==0.or.qeq0) then
2595      natom_nondiag=1;if (nondiag==1) natom_nondiag=natom
2596      do ja=1,natom_nondiag
2597        do ia=1,natom
2598          do isym=1,nsym
2599            indi=indsym(4,isym,ia)
2600            indj=1;if (nondiag==1) indj=indsym(4,isym,ja)
2601            work(:,:,:)=zero
2602            do mu=1,3
2603              do nu=1,3
2604                do kappa=1,3
2605                  work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj)
2606                end do
2607              end do
2608            end do
2609            do mu=1,3
2610              do nu=1,3
2611                do kappa=1,3
2612                  sdyfro(:,kappa,mu,ia,ja)=sdyfro(:,kappa,mu,ia,ja)+symrec(mu,nu,isym)*work(:,kappa,nu)
2613                end do
2614              end do
2615            end do
2616          end do
2617        end do
2618      end do
2619      div=one/dble(nsym)
2620      sdyfro(:,:,:,:,:)=div*sdyfro(:,:,:,:,:)
2621 !    === Non diagonal dyn. matrix AND q<>0
2622    else
2623      do ja=1,natom
2624        do ia=1,natom
2625          nsym_used=0
2626          do isym=1,nsym
2627            if (symq(4,1,isym)==1) then
2628              arg=two_pi*(qphon(1)*(indsym(1,isym,ia)-indsym(1,isym,ja)) &
2629 &             +qphon(2)*(indsym(2,isym,ia)-indsym(2,isym,ja)) &
2630 &             +qphon(3)*(indsym(3,isym,ia)-indsym(3,isym,ja)))
2631              phaser=cos(arg);phasei=sin(arg)
2632              nsym_used=nsym_used+1
2633              indi=indsym(4,isym,ia)
2634              indj=indsym(4,isym,ja)
2635              work(:,:,:)=zero
2636              do mu=1,3
2637                do nu=1,3
2638                  do kappa=1,3
2639                    work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj)
2640                  end do
2641                end do
2642              end do
2643              do mu=1,3
2644                do nu=1,3
2645                  do kappa=1,3
2646                    sdyfro(1,kappa,mu,ia,ja)=sdyfro(1,kappa,mu,ia,ja) &
2647 &                   +symrec(mu,nu,isym)*(work(1,kappa,nu)*phaser-work(2,kappa,nu)*phasei)
2648                  end do
2649                end do
2650              end do
2651              if (cplex==2) then
2652                do mu=1,3
2653                  do nu=1,3
2654                    do kappa=1,3
2655                      sdyfro(2,kappa,mu,ia,ja)=sdyfro(2,kappa,mu,ia,ja) &
2656 &                     +symrec(mu,nu,isym)*(work(1,kappa,nu)*phasei+work(2,kappa,nu)*phaser)
2657                    end do
2658                  end do
2659                end do
2660              end if
2661            end if
2662          end do
2663          div=one/dble(nsym_used)
2664          sdyfro(:,:,:,ia,ja)=div*sdyfro(:,:,:,ia,ja)
2665        end do
2666      end do
2667    end if
2668 
2669  end if
2670 
2671 end subroutine dfpt_sydy

m_dynmat/dfpt_sygra [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 dfpt_sygra

FUNCTION

 Symmetrize derivatives of energy with respect to coordinates,
 as appearing in phonon calculations.
 Unsymmetrized gradients are input as deunsy; symmetrized grads are then placed in desym.
 If nsym=1 simply copy deunsy into desym (only symmetry is identity).
 The index of the initial perturbation is needed, in case there is a change
 of atom position (moved in another cell) due to the symmetry operation.

INPUTS

  natom=number of atoms in cell
  deunsy(2,3,natom)=unsymmetrized gradients wrt dimensionless tn (hartree)
  note: there is a real and a imaginary part ...
  indsym(4,nsym,natom)=label given by subroutine symatm, indicating atom
   label which gets rotated into given atom by given symmetry
   (first three elements are related primitive translation--
   see symatm where this is computed)
  nsym=number of symmetry operators in group
  ipert=index of the initial perturbation
  qpt(3)= wavevector of the phonon, in reduced coordinates
  symrec(3,3,nsym)=symmetries of group in terms of operations on
    reciprocal space primitive translations--see comments below

OUTPUT

 desym(2,3,natom)=symmetrized gradients wrt dimensionless tn (hartree)

NOTES

 Written by X. Gonze starting from sygrad, written by D.C. Allan:
    introduction of the q vector for phonon symmetrization
 This routine should once be merged with sygrad...

PARENTS

      dfpt_nstdy,dfpt_nstpaw

CHILDREN

SOURCE

2413 subroutine dfpt_sygra(natom,desym,deunsy,indsym,ipert,nsym,qpt,symrec)
2414 
2415 
2416 !This section has been created automatically by the script Abilint (TD).
2417 !Do not modify the following lines by hand.
2418 #undef ABI_FUNC
2419 #define ABI_FUNC 'dfpt_sygra'
2420 !End of the abilint section
2421 
2422  implicit none
2423 
2424 !Arguments -------------------------------
2425 !scalars
2426  integer,intent(in) :: ipert,natom,nsym
2427 !arrays
2428  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym)
2429  real(dp),intent(in) :: deunsy(2,3,natom),qpt(3)
2430  real(dp),intent(out) :: desym(2,3,natom)
2431 
2432 !Local variables -------------------------
2433 !scalars
2434  integer :: ia,ind,isym,mu
2435  real(dp) :: arg,im,re,sumi,sumr
2436 
2437 ! *********************************************************************
2438 
2439 !DEBUG
2440 !write(std_out,*)' dfpt_sygra : enter '
2441 !write(std_out,*)' dfpt_sygra : qpt(:)',qpt(:)
2442 !do ia=1,natom
2443 !do mu=1,3
2444 !write(std_out,*)' dfpt_sygra : deunsy(:2,mu,ia)',deunsy(:2,mu,ia)
2445 !enddo
2446 !enddo
2447 !ENDDEBUG
2448 
2449  if (nsym==1) then
2450 
2451 !  Only symmetry is identity so simply copy
2452    desym(:,:,:)=deunsy(:,:,:)
2453 
2454  else
2455 
2456 !  Actually conduct symmetrization
2457 !  write(std_out,*)' dfpt_sygra : desym(:2,:3,:natom),qpt(:)',desym(:2,:3,:natom),qpt(:)
2458    do ia=1,natom
2459 !    write(std_out,*)' dfpt_sygra : ia=',ia
2460      do mu=1,3
2461        sumr=zero
2462        sumi=zero
2463 !      write(std_out,*)' dfpt_sygra : mu=',mu
2464        do isym=1,nsym
2465          ind=indsym(4,isym,ia)
2466 !        Must shift the atoms back to the unit cell.
2467 !        arg=two_pi*( qpt(1)*indsym(1,isym,ia)&
2468 !        &         +qpt(2)*indsym(2,isym,ia)&
2469 !        &         +qpt(3)*indsym(3,isym,ia) )
2470 !        Selection of non-zero q point, to avoid ipert being outside the 1 ... natom range
2471          if(qpt(1)**2+qpt(2)**2+qpt(3)**2 > tol16)then
2472            arg=two_pi*( qpt(1)*(indsym(1,isym,ia)-indsym(1,isym,ipert))&
2473 &           +qpt(2)* (indsym(2,isym,ia)-indsym(2,isym,ipert))&
2474 &           +qpt(3)* (indsym(3,isym,ia)-indsym(3,isym,ipert)))
2475          else
2476            arg=zero
2477          end if
2478 
2479          re=dble(symrec(mu,1,isym))*deunsy(1,1,ind)+&
2480 &         dble(symrec(mu,2,isym))*deunsy(1,2,ind)+&
2481 &         dble(symrec(mu,3,isym))*deunsy(1,3,ind)
2482          im=dble(symrec(mu,1,isym))*deunsy(2,1,ind)+&
2483 &         dble(symrec(mu,2,isym))*deunsy(2,2,ind)+&
2484 &         dble(symrec(mu,3,isym))*deunsy(2,3,ind)
2485          sumr=sumr+re*cos(arg)-im*sin(arg)
2486          sumi=sumi+re*sin(arg)+im*cos(arg)
2487 !        sumr=sumr+re
2488 !        sumi=sumi+im
2489 !        write(std_out,*)' dfpt_sygra : isym,indsym(4,isym,ia),arg,re,im,sumr,sumi',&
2490 !        &      isym,indsym(4,isym,ia),arg,re,im,sumr,sumi
2491        end do
2492        desym(1,mu,ia)=sumr/dble(nsym)
2493        desym(2,mu,ia)=sumi/dble(nsym)
2494 !      write(std_out,*)' dfpt_sygra : desym(:,mu,ia)',desym(:,mu,ia)
2495      end do
2496    end do
2497  end if
2498 
2499 end subroutine dfpt_sygra

m_dynmat/dist9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

3618 subroutine dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt)
3619 
3620 
3621 !This section has been created automatically by the script Abilint (TD).
3622 !Do not modify the following lines by hand.
3623 #undef ABI_FUNC
3624 #define ABI_FUNC 'dist9'
3625 !End of the abilint section
3626 
3627  implicit none
3628 
3629 !Arguments -------------------------------
3630 !scalars
3631  integer,intent(in) :: natom,nrpt
3632 !arrays
3633  real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3)
3634  real(dp),intent(in) :: rpt(3,nrpt)
3635  real(dp),intent(out) :: dist(natom,natom,nrpt)
3636 
3637 !Local variables -------------------------
3638 !scalars
3639  integer :: ia,ib,ii,irpt
3640 !arrays
3641  real(dp) :: ra(3),rb(3),rdiff(3),red(3),rptcar(3),xred(3)
3642 
3643 ! *********************************************************************
3644 
3645 !BIG loop on all generic atoms
3646  do ia=1,natom
3647 
3648 !  First transform canonical coordinates to reduced coordinates
3649    do ii=1,3
3650      xred(ii)=gprim(1,ii)*rcan(1,ia)+gprim(2,ii)*rcan(2,ia)&
3651 &     +gprim(3,ii)*rcan(3,ia)
3652    end do
3653 !  Then to cartesian coordinates
3654    ra(:)=xred(1)*acell(1)*rprim(:,1)+xred(2)*acell(2)*rprim(:,2)+&
3655 &   xred(3)*acell(3)*rprim(:,3)
3656    do ib=1,natom
3657      do ii=1,3
3658        xred(ii)=gprim(1,ii)*rcan(1,ib)+gprim(2,ii)*rcan(2,ib)&
3659 &       +gprim(3,ii)*rcan(3,ib)
3660      end do
3661      do ii=1,3
3662        rb(ii)=xred(1)*acell(1)*rprim(ii,1)+&
3663 &       xred(2)*acell(2)*rprim(ii,2)+&
3664 &       xred(3)*acell(3)*rprim(ii,3)
3665      end do
3666      do irpt=1,nrpt
3667 !      First transform it to reduced coordinates
3668        do ii=1,3
3669          red(ii)=gprim(1,ii)*rpt(1,irpt)+gprim(2,ii)*rpt(2,irpt)&
3670 &         +gprim(3,ii)*rpt(3,irpt)
3671        end do
3672 !      Then to cartesian coordinates
3673        do ii=1,3
3674          rptcar(ii)=red(1)*acell(1)*rprim(ii,1)+&
3675 &         red(2)*acell(2)*rprim(ii,2)+&
3676 &         red(3)*acell(3)*rprim(ii,3)
3677        end do
3678        do ii=1,3
3679          rdiff(ii)=-rptcar(ii)+ra(ii)-rb(ii)
3680        end do
3681        dist(ia,ib,irpt)=(rdiff(1)**2+rdiff(2)**2+rdiff(3)**2)**0.5
3682 
3683      end do
3684 
3685    end do
3686  end do
3687 
3688 end subroutine dist9

m_dynmat/dymfz9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

5474 subroutine dymfz9(dynmat,natom,nqpt,gprim,option,spqpt,trans)
5475 
5476 
5477 !This section has been created automatically by the script Abilint (TD).
5478 !Do not modify the following lines by hand.
5479 #undef ABI_FUNC
5480 #define ABI_FUNC 'dymfz9'
5481 !End of the abilint section
5482 
5483  implicit none
5484 
5485 !Arguments -------------------------------
5486 !scalars
5487  integer,intent(in) :: natom,nqpt,option
5488 !arrays
5489  real(dp),intent(in) :: gprim(3,3),spqpt(3,nqpt),trans(3,natom)
5490  real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)
5491 
5492 !Local variables -------------------------
5493 !scalars
5494  integer :: ia,ib,iqpt,mu,nu
5495  real(dp) :: im,ktrans,re
5496 !arrays
5497  real(dp) :: kk(3)
5498 
5499 ! *********************************************************************
5500 
5501  DBG_ENTER("COLL")
5502 
5503  do iqpt=1,nqpt
5504    !  Definition of q in normalized reciprocal space
5505    kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
5506    kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
5507    kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
5508 
5509    if(option==1)then
5510      kk(1)=-kk(1)
5511      kk(2)=-kk(2)
5512      kk(3)=-kk(3)
5513    end if
5514 
5515    do ia=1,natom
5516      do ib=1,natom
5517 !      Product of q with the differences between the two atomic translations
5518        ktrans=kk(1)*(trans(1,ia)-trans(1,ib))+kk(2)*(trans(2,ia)-&
5519 &       trans(2,ib))+kk(3)*(trans(3,ia)-trans(3,ib))
5520        do mu=1,3
5521          do nu=1,3
5522            re=dynmat(1,mu,ia,nu,ib,iqpt)
5523            im=dynmat(2,mu,ia,nu,ib,iqpt)
5524 !          Transformation of the Old dynamical matrices by New ones by multiplication by a phase shift
5525            dynmat(1,mu,ia,nu,ib,iqpt)=re*cos(two_pi*ktrans)-im*sin(two_pi*ktrans)
5526            dynmat(2,mu,ia,nu,ib,iqpt)=re*sin(two_pi*ktrans)+im*cos(two_pi*ktrans)
5527          end do
5528        end do
5529      end do
5530    end do
5531  end do
5532 
5533  DBG_EXIT("COLL")
5534 
5535 end subroutine dymfz9

m_dynmat/dynmat_dq [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 dynmat_dq

FUNCTION

   Compute the derivative D(q)/dq of the dynamical matrix via Fourier transform
   of the interatomic forces

INPUTS

 qpt(3)= Reduced coordinates of the q vector in reciprocal space
 natom= Number of atoms in the unit cell
 gprim(3,3)= Normalized coordinates in reciprocal space
 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)!!)
 atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space
 wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector

OUTPUT

 dddq(2,3,natom,3,natom,3)= Derivate of the dynamical matrix in cartesian coordinates.
  The tree directions are stored in the last dimension.
  These coordinates are normalized (=> * acell(3)!!)

PARENTS

      m_ifc

CHILDREN

SOURCE

3941 subroutine dynmat_dq(qpt,natom,gprim,nrpt,rpt,atmfrc,wghatm,dddq)
3942 
3943 
3944 !This section has been created automatically by the script Abilint (TD).
3945 !Do not modify the following lines by hand.
3946 #undef ABI_FUNC
3947 #define ABI_FUNC 'dynmat_dq'
3948 !End of the abilint section
3949 
3950  implicit none
3951 
3952 !Arguments -------------------------------
3953 !scalars
3954  integer,intent(in) :: natom,nrpt
3955 !arrays
3956  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt(3)
3957  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
3958  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
3959  real(dp),intent(out) :: dddq(2,3,natom,3,natom,3)
3960 
3961 !Local variables -------------------------
3962 !scalars
3963  integer :: ia,ib,irpt,mu,nu,ii
3964  real(dp) :: im,kr,re
3965 !arrays
3966  real(dp) :: kk(3),fact(2,3)
3967 
3968 ! *********************************************************************
3969 
3970  dddq = zero
3971 
3972  do irpt=1,nrpt
3973    ! Calculation of the k coordinates in Normalized Reciprocal coordinates
3974    kk(1) = qpt(1)*gprim(1,1)+ qpt(2)*gprim(1,2) + qpt(3)*gprim(1,3)
3975    kk(2) = qpt(1)*gprim(2,1)+ qpt(2)*gprim(2,2) + qpt(3)*gprim(2,3)
3976    kk(3) = qpt(1)*gprim(3,1)+ qpt(2)*gprim(3,2) + qpt(3)*gprim(3,3)
3977 
3978    ! Product of k and r
3979    kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3980 
3981    ! Get phase factor
3982    re=cos(two_pi*kr); im=sin(two_pi*kr)
3983 
3984    ! Inner loop on atoms and directions
3985    do ib=1,natom
3986      do ia=1,natom
3987        if (abs(wghatm(ia,ib,irpt))>1.0d-10) then
3988          ! take into account rotation due to i.
3989          fact(1,:) = -im * wghatm(ia,ib,irpt) * rpt(:,irpt)
3990          fact(2,:) =  re * wghatm(ia,ib,irpt) * rpt(:,irpt)
3991          do nu=1,3
3992            do mu=1,3
3993              ! Real and imaginary part of the dynamical matrices
3994              ! Atmfrc should be real
3995              do ii=1,3
3996                dddq(1,mu,ia,nu,ib,ii) = dddq(1,mu,ia,nu,ib,ii) + fact(1,ii) * atmfrc(mu,ia,nu,ib,irpt)
3997                dddq(2,mu,ia,nu,ib,ii) = dddq(2,mu,ia,nu,ib,ii) + fact(2,ii) * atmfrc(mu,ia,nu,ib,irpt)
3998              end do
3999            end do
4000          end do
4001        end if
4002      end do
4003    end do
4004  end do
4005 
4006 end subroutine dynmat_dq

m_dynmat/ftgam [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 ftgam

FUNCTION

 If qtor=1 (q->r):
  Generates the Fourier transform of the recip space gkk matrices
  to obtain the real space ones.
 If qtor=0 (r->q):
  Generates the Fourier transform of the real space gkk matrices
  to obtain the reciprocal space ones.

INPUTS

 natom= Number of atoms in the unit cell
 nqpt= Number of q points in the Brillouin zone
           if qtor=0 this number is read in the input file
 nrpt= Number of R points in the Big Box
 qtor= ( q to r : see above )
 rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
           These coordinates are normalized (=> * acell(3)!!)
 qpt_full(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
           if qtor=0 these vectors are read in the input file
 wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/output
 gam_qpt(2,3*natom*3*natom,nqpt)
  = gamma matrices in recip space coming from the Derivative Data Base
 gam_rpt(2,3*natom*3*natom,nrpt)
  = gamma matrices in real space stored in file unit_gkk_rpt

PARENTS

      elphon,get_tau_k,integrate_gamma_alt,m_phgamma,mka2f,mka2f_tr
      mka2f_tr_lova,mkph_linwid

CHILDREN

NOTES

   copied from ftiaf9.f
   recip to real space: real space is forced to disk file unit_gkk_rpt
                        recip space depends on gkqwrite and unitgkq3
   real to recip space: real space is forced to disk file unit_gkk_rpt
                        recip space is necessarily in memory in gkk_qpt

    real space elements are complex, but could be reduced, as (-r) = (+r)*

SOURCE

6442 subroutine ftgam (wghatm,gam_qpt,gam_rpt,natom,nqpt,nrpt,qtor,coskr, sinkr)
6443 
6444 
6445 !This section has been created automatically by the script Abilint (TD).
6446 !Do not modify the following lines by hand.
6447 #undef ABI_FUNC
6448 #define ABI_FUNC 'ftgam'
6449 !End of the abilint section
6450 
6451  implicit none
6452 
6453 !Arguments -------------------------------
6454 !scalars
6455  integer,intent(in) :: natom,nqpt,nrpt,qtor
6456 !arrays
6457  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
6458  real(dp),intent(inout) :: gam_qpt(2,3*natom*3*natom,nqpt)
6459  real(dp),intent(inout) :: gam_rpt(2,3*natom*3*natom,nrpt)
6460  real(dp),intent(in) :: coskr(nqpt,nrpt)
6461  real(dp),intent(in) :: sinkr(nqpt,nrpt)
6462 
6463 !Local variables -------------------------
6464 !scalars
6465  integer :: iatom,idir,ip,iqpt,irpt,jatom,jdir
6466  real(dp) :: im,re
6467  character(len=500) :: message
6468 
6469 ! *********************************************************************
6470 
6471  select case (qtor)
6472 !
6473    case (1)  !Recip to real space
6474      gam_rpt(:,:,:) = zero
6475      do irpt=1,nrpt
6476        do iqpt=1,nqpt
6477 !        Get the phase factor with normalization!
6478          re=coskr(iqpt,irpt)
6479          im=sinkr(iqpt,irpt)
6480          do ip=1,3*natom*3*natom
6481 !          Real and imaginary part of the real-space gam matrices
6482            gam_rpt(1,ip,irpt) = gam_rpt(1,ip,irpt)&
6483 &           +re*gam_qpt(1,ip,iqpt) &
6484 &           +im*gam_qpt(2,ip,iqpt)
6485            gam_rpt(2,ip,irpt) = gam_rpt(2,ip,irpt)&
6486 &           +re*gam_qpt(2,ip,iqpt) &
6487 &           -im*gam_qpt(1,ip,iqpt)
6488          end do
6489        end do
6490      end do
6491      gam_rpt = gam_rpt/nqpt
6492 !
6493    case (0) ! Recip space from real space
6494 
6495      gam_qpt(:,:,:)=zero
6496 
6497      do irpt=1,nrpt
6498        do iqpt=1,nqpt
6499 
6500          do iatom=1,natom
6501            do jatom=1,natom
6502              re = coskr(iqpt,irpt)*wghatm(iatom,jatom,irpt)
6503              im = sinkr(iqpt,irpt)*wghatm(iatom,jatom,irpt)
6504 
6505              do idir=1,3
6506                do jdir=1,3
6507 !                Get phase factor
6508 
6509                  ip= jdir + (jatom-1)*3 + (idir-1)*3*natom + (iatom-1)*9*natom
6510 !                Real and imaginary part of the interatomic forces
6511                  gam_qpt(1,ip,iqpt)=&
6512 &                 gam_qpt(1,ip,iqpt)&
6513 &                 +re*gam_rpt(1,ip,irpt)&
6514 &                 -im*gam_rpt(2,ip,irpt)
6515 !                !DEBUG
6516                  gam_qpt(2,ip,iqpt)=&
6517 &                 gam_qpt(2,ip,iqpt)&
6518 &                 +im*gam_rpt(1,ip,irpt)&
6519 &                 +re*gam_rpt(2,ip,irpt)
6520 !                !ENDDEBUG
6521 
6522                end do ! end jdir
6523              end do ! end idir
6524            end do
6525          end do ! end iatom
6526 
6527        end do ! end iqpt
6528      end do ! end irpt
6529 
6530    case default ! There is no other space to Fourier transform from
6531      write(message,'(a,i0,a)' )'  The only allowed values for qtor are 0 or 1, while  qtor=',qtor,' has been required.'
6532      MSG_BUG(message)
6533  end select
6534 
6535 end subroutine ftgam

m_dynmat/ftgam_init [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 ftgam_init

FUNCTION

  Generates the sin and cos phases for the Fourier transform of the gkk matrices

INPUTS

 gprim = reciprocal space vectors to get cartesian coord for qpt
 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)!!)
 qpt_full(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
           if qtor=0 these vectors are read in the input file

OUTPUT

 coskr, sinkr = cosine and sine of phase factors for given r and q points

PARENTS

      elphon,get_tau_k,integrate_gamma_alt,m_phgamma,mka2f,mka2f_tr
      mka2f_tr_lova,mkph_linwid

CHILDREN

SOURCE

6566 subroutine ftgam_init (gprim,nqpt,nrpt,qpt_full,rpt,coskr, sinkr)
6567 
6568 
6569 !This section has been created automatically by the script Abilint (TD).
6570 !Do not modify the following lines by hand.
6571 #undef ABI_FUNC
6572 #define ABI_FUNC 'ftgam_init'
6573 !End of the abilint section
6574 
6575  implicit none
6576 
6577 !Arguments -------------------------------
6578 !scalars
6579  integer,intent(in) :: nqpt,nrpt
6580 !arrays
6581  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt_full(3,nqpt)
6582  real(dp),intent(out) :: coskr(nqpt,nrpt)
6583  real(dp),intent(out) :: sinkr(nqpt,nrpt)
6584 
6585 !Local variables -------------------------
6586 !scalars
6587  integer :: iqpt,irpt
6588  real(dp) :: kr
6589 !arrays
6590  real(dp) :: kk(3)
6591 
6592 ! *********************************************************************
6593 
6594 ! Prepare the phase factors
6595  do iqpt=1,nqpt
6596    ! Calculation of the k coordinates in Normalized Reciprocal coordinates
6597    kk(1)=   qpt_full(1,iqpt)*gprim(1,1)+&
6598 &   qpt_full(2,iqpt)*gprim(1,2)+&
6599 &   qpt_full(3,iqpt)*gprim(1,3)
6600 
6601    kk(2)=   qpt_full(1,iqpt)*gprim(2,1)+&
6602 &   qpt_full(2,iqpt)*gprim(2,2)+&
6603 &   qpt_full(3,iqpt)*gprim(2,3)
6604 
6605    kk(3)=   qpt_full(1,iqpt)*gprim(3,1)+&
6606 &   qpt_full(2,iqpt)*gprim(3,2)+&
6607 &   qpt_full(3,iqpt)*gprim(3,3)
6608    do irpt=1,nrpt
6609      ! Product of k and r
6610      kr = kk(1)*rpt(1,irpt)+ kk(2)*rpt(2,irpt)+ kk(3)*rpt(3,irpt)
6611      coskr(iqpt,irpt)=cos(two_pi*kr)
6612      sinkr(iqpt,irpt)=sin(two_pi*kr)
6613    end do
6614  end do
6615 
6616 end subroutine ftgam_init

m_dynmat/ftifc_q2r [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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(3,natom,3,natom,nrpt)= Interatomic Forces in real space !!
  We used the imaginary part just for debugging !

PARENTS

      m_ifc

CHILDREN

SOURCE

3723 subroutine ftifc_q2r(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt)
3724 
3725 
3726 !This section has been created automatically by the script Abilint (TD).
3727 !Do not modify the following lines by hand.
3728 #undef ABI_FUNC
3729 #define ABI_FUNC 'ftifc_q2r'
3730 !End of the abilint section
3731 
3732  implicit none
3733 
3734 !Arguments -------------------------------
3735 !scalars
3736  integer,intent(in) :: natom,nqpt,nrpt
3737 !arrays
3738  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
3739  real(dp),intent(out) :: atmfrc(3,natom,3,natom,nrpt)
3740  real(dp),intent(in) :: dynmat(2,3,natom,3,natom,nqpt)
3741 
3742 !Local variables -------------------------
3743 !scalars
3744  integer :: ia,ib,iqpt,irpt,mu,nu
3745  real(dp) :: im,kr,re
3746 !arrays
3747  real(dp) :: kk(3)
3748 
3749 ! *********************************************************************
3750 
3751  DBG_ENTER("COLL")
3752 
3753 !Interatomic Forces from Dynamical Matrices
3754  atmfrc = zero
3755  do irpt=1,nrpt
3756    do iqpt=1,nqpt
3757 
3758 !    Calculation of the k coordinates in Normalized Reciprocal coordinates
3759      kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
3760      kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
3761      kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
3762 
3763 !    Product of k and r
3764      kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3765 
3766 !    Get the phase factor
3767      re=cos(two_pi*kr)
3768      im=sin(two_pi*kr)
3769 
3770 !    Now, big inner loops on atoms and directions
3771 
3772 !    The indices are ordered to give better speed
3773      do ib=1,natom
3774        do nu=1,3
3775          do ia=1,natom
3776            do mu=1,3
3777 !            Real and imaginary part of the interatomic forces
3778              atmfrc(mu,ia,nu,ib,irpt)=atmfrc(mu,ia,nu,ib,irpt)&
3779 &             +re*dynmat(1,mu,ia,nu,ib,iqpt)&
3780 &             +im*dynmat(2,mu,ia,nu,ib,iqpt)
3781 !            The imaginary part should be equal to zero !!!!!!
3782 !            atmfrc(2,mu,ia,nu,ib,irpt)=atmfrc(2,mu,ia,nu,ib,irpt)
3783 !            &          +re*dynmat(2,mu,ia,nu,ib,iqpt)
3784 !            &          -im*dynmat(1,mu,ia,nu,ib,iqpt)
3785            end do
3786          end do
3787        end do
3788      end do
3789 
3790    end do
3791  end do
3792 
3793 !The sumifc has to be weighted by a normalization factor of 1/nqpt
3794  atmfrc = atmfrc/nqpt
3795 
3796  DBG_EXIT("COLL")
3797 
3798 end subroutine ftifc_q2r

m_dynmat/ftifc_r2q [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 ftifc_r2q

FUNCTION

 (r->q):
   Generates the Fourier transform of the interatomic forces
   to obtain dynamical matrices (reciprocal space).

INPUTS

 atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space
 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

SOURCE

3834 subroutine ftifc_r2q(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt,wghatm)
3835 
3836 
3837 !This section has been created automatically by the script Abilint (TD).
3838 !Do not modify the following lines by hand.
3839 #undef ABI_FUNC
3840 #define ABI_FUNC 'ftifc_r2q'
3841 !End of the abilint section
3842 
3843  implicit none
3844 
3845 !Arguments -------------------------------
3846 !scalars
3847  integer,intent(in) :: natom,nqpt,nrpt
3848 !arrays
3849  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
3850  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
3851  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
3852  real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt)
3853 
3854 !Local variables -------------------------
3855 !scalars
3856  integer :: ia,ib,iqpt,irpt,mu,nu
3857  real(dp) :: facti,factr,im,kr,re
3858 !arrays
3859  real(dp) :: kk(3)
3860 
3861 ! *********************************************************************
3862 
3863  dynmat = zero
3864 
3865  do iqpt=1,nqpt
3866    do irpt=1,nrpt
3867 
3868 !    Calculation of the k coordinates in Normalized Reciprocal coordinates
3869      kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)* gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
3870      kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)* gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
3871      kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)* gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
3872 
3873 !    Product of k and r
3874      kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3875 
3876 !    Get phase factor
3877      re=cos(two_pi*kr)
3878      im=sin(two_pi*kr)
3879 
3880 !    Inner loop on atoms and directions
3881      do ib=1,natom
3882        do ia=1,natom
3883          if(abs(wghatm(ia,ib,irpt))>1.0d-10)then
3884            factr=re*wghatm(ia,ib,irpt)
3885            facti=im*wghatm(ia,ib,irpt)
3886            do nu=1,3
3887              do mu=1,3
3888 !              Real and imaginary part of the dynamical matrices
3889                dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)&
3890 &               +factr*atmfrc(mu,ia,nu,ib,irpt)
3891 !              Atmfrc should be real
3892 !              &       -im*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt)
3893                dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)&
3894 &               +facti*atmfrc(mu,ia,nu,ib,irpt)
3895 !              Atmfrc should be real
3896 !              &        +re*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt)
3897              end do
3898            end do
3899          end if
3900        end do
3901      end do
3902    end do
3903  end do
3904 
3905 end subroutine ftifc_r2q

m_dynmat/gtdyn9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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(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,ddb_interpolate,m_effective_potential_file,m_ifc,m_phonons

CHILDREN

SOURCE

5677 subroutine gtdyn9(acell,atmfrc,dielt,dipdip,&
5678 & dyewq0,d2cart,gmet,gprim,mpert,natom,&
5679 & nrpt,qphnrm,qpt,rmet,rprim,rpt,&
5680 & trans,ucvol,wghatm,xred,zeff)
5681 
5682 
5683 !This section has been created automatically by the script Abilint (TD).
5684 !Do not modify the following lines by hand.
5685 #undef ABI_FUNC
5686 #define ABI_FUNC 'gtdyn9'
5687 !End of the abilint section
5688 
5689  implicit none
5690 
5691 !Arguments -------------------------------
5692 !scalars
5693  integer,intent(in) :: dipdip,mpert,natom,nrpt
5694  real(dp),intent(in) :: qphnrm,ucvol
5695 !arrays
5696  real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qpt(3)
5697  real(dp),intent(in) :: rmet(3,3),rprim(3,3),rpt(3,nrpt)
5698  real(dp),intent(in) :: trans(3,natom),wghatm(natom,natom,nrpt),xred(3,natom)
5699  real(dp),intent(in) :: zeff(3,3,natom)
5700  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
5701  real(dp),intent(in) :: dyewq0(3,3,natom)
5702  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
5703 
5704 !Local variables -------------------------
5705 !scalars
5706  integer,parameter :: nqpt1=1,option2=2,sumg0=0,plus1=1,iqpt1=1
5707  integer :: i1,i2,ib,nsize
5708 !arrays
5709  real(dp) :: qphon(3)
5710  real(dp),allocatable :: dq(:,:,:,:,:),dyew(:,:,:,:,:)
5711 
5712 ! *********************************************************************
5713 
5714  ABI_ALLOCATE(dq,(2,3,natom,3,natom))
5715 
5716 !Get the normalized wavevector
5717  if(abs(qphnrm)<1.0d-7)then
5718    qphon(1:3)=zero
5719  else
5720    qphon(1:3)=qpt(1:3)/qphnrm
5721  end if
5722 
5723 !Generate the analytical part from the interatomic forces
5724  call ftifc_r2q (atmfrc,dq,gprim,natom,nqpt1,nrpt,rpt,qphon,wghatm)
5725 
5726 !The analytical dynamical matrix dq has been generated
5727 !in the normalized canonical coordinate system. Now, the
5728 !phase is modified, in order to recover the usual (xred) coordinate of atoms.
5729  call dymfz9(dq,natom,nqpt1,gprim,option2,qphon,trans)
5730 
5731  if (dipdip==1) then
5732 !  Add the non-analytical part
5733 !  Compute dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
5734 !  second energy derivative wrt xred(3,natom) in Hartrees
5735 !  (Denoted A-bar in the notes)
5736    ABI_ALLOCATE(dyew,(2,3,natom,3,natom))
5737 
5738    call ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff)
5739    call q0dy3_apply(natom,dyewq0,dyew)
5740    call nanal9(dyew,dq,iqpt1,natom,nqpt1,plus1)
5741 
5742    ABI_DEALLOCATE(dyew)
5743  end if
5744 
5745 !Copy the dynamical matrix in the proper location
5746 
5747 !First zero all the elements
5748  nsize=2*(3*mpert)**2
5749  d2cart = zero
5750 
5751 !Copy the elements from dq to d2cart
5752  d2cart(:,:,1:natom,:,1:natom)=dq(:,:,1:natom,:,1:natom)
5753 
5754 !In case we have the gamma point,
5755  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)then
5756 !  Copy the effective charge and dielectric constant in the final array
5757    do i1=1,3
5758      do i2=1,3
5759        d2cart(1,i1,natom+2,i2,natom+2)=dielt(i1,i2)
5760        do ib=1,natom
5761          d2cart(1,i1,natom+2,i2,ib)=zeff(i1,i2,ib)
5762          d2cart(1,i2,ib,i1,natom+2)=zeff(i1,i2,ib)
5763        end do
5764      end do
5765    end do
5766  end if
5767 
5768  ABI_DEALLOCATE(dq)
5769 
5770 end subroutine gtdyn9

m_dynmat/ifclo9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

4036 subroutine ifclo9(ifccar,ifcloc,vect1,vect2,vect3)
4037 
4038 
4039 !This section has been created automatically by the script Abilint (TD).
4040 !Do not modify the following lines by hand.
4041 #undef ABI_FUNC
4042 #define ABI_FUNC 'ifclo9'
4043 !End of the abilint section
4044 
4045  implicit none
4046 
4047 !Arguments -------------------------------
4048 !arrays
4049  real(dp),intent(in) :: ifccar(3,3),vect1(3),vect2(3),vect3(3)
4050  real(dp),intent(out) :: ifcloc(3,3)
4051 
4052 !Local variables -------------------------
4053 !scalars
4054  integer :: ii,jj
4055 !arrays
4056  real(dp) :: work(3,3)
4057 
4058 ! *********************************************************************
4059 
4060  do jj=1,3
4061    do ii=1,3
4062      work(jj,ii)=zero
4063    end do
4064    do ii=1,3
4065      work(jj,1)=work(jj,1)+ifccar(jj,ii)*vect1(ii)
4066      work(jj,2)=work(jj,2)+ifccar(jj,ii)*vect2(ii)
4067      work(jj,3)=work(jj,3)+ifccar(jj,ii)*vect3(ii)
4068    end do
4069  end do
4070 
4071  do jj=1,3
4072    do ii=1,3
4073      ifcloc(ii,jj)=zero
4074    end do
4075    do ii=1,3
4076      ifcloc(1,jj)=ifcloc(1,jj)+vect1(ii)*work(ii,jj)
4077      ifcloc(2,jj)=ifcloc(2,jj)+vect2(ii)*work(ii,jj)
4078      ifcloc(3,jj)=ifcloc(3,jj)+vect3(ii)*work(ii,jj)
4079    end do
4080  end do
4081 
4082 end subroutine ifclo9

m_dynmat/make_bigbox [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 or -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

 cell= (nrpt,3) Give the index of the the cell and irpt
 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 with the proper dimension. Client code is responsible
  for the deallocation.

PARENTS

      m_ifc

CHILDREN

SOURCE

2911 subroutine make_bigbox(brav,cell,ngqpt,nqshft,rprim,nrpt,rpt)
2912 
2913 
2914 !This section has been created automatically by the script Abilint (TD).
2915 !Do not modify the following lines by hand.
2916 #undef ABI_FUNC
2917 #define ABI_FUNC 'make_bigbox'
2918 !End of the abilint section
2919 
2920  implicit none
2921 
2922 !Arguments -------------------------------
2923 !scalars
2924  integer,intent(in) :: brav,nqshft
2925  integer,intent(out) :: nrpt
2926 !arrays
2927  integer,intent(in) :: ngqpt(3)
2928  real(dp),intent(in) :: rprim(3,3)
2929  real(dp),allocatable,intent(out) :: rpt(:,:)
2930  integer,allocatable,intent(out) :: cell(:,:)
2931 
2932 !Local variables -------------------------
2933 !scalars
2934  integer :: choice,mrpt
2935 !arrays
2936  real(dp) :: dummy_rpt(3,1)
2937  integer:: dummy_cell(1,3)
2938 
2939 ! *********************************************************************
2940 
2941  ! Compute the number of points (cells) in real space
2942  choice=0
2943  call bigbx9(brav,dummy_cell,choice,1,ngqpt,nqshft,mrpt,rprim,dummy_rpt)
2944 
2945  ! Now we can allocate and calculate the points and the weights.
2946  nrpt = mrpt
2947  ABI_ALLOCATE(rpt,(3,nrpt))
2948  ABI_ALLOCATE(cell,(3,nrpt))
2949 
2950  choice=1
2951  call bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)
2952 
2953 end subroutine make_bigbox

m_dynmat/massmult_and_breaksym [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

  mult_masses_and_break_symms

FUNCTION

  Multiply the IFC(q) by the atomic masses, slightly break symmetry to make tests more
  portable and make the matrix hermitian before returning.

INPUTS

  amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms)
  natom=number of atoms in unit cell
  ntypat=number of atom types
  typat(natom)=integer label of each type of atom (1,2,...)

SIDE EFFECTS

  mat(2*3*natom*3*natom)=Multiplies by atomic masses in output.

PARENTS

CHILDREN

SOURCE

6327 subroutine massmult_and_breaksym(natom, ntypat, typat, amu, mat)
6328 
6329 
6330 !This section has been created automatically by the script Abilint (TD).
6331 !Do not modify the following lines by hand.
6332 #undef ABI_FUNC
6333 #define ABI_FUNC 'massmult_and_breaksym'
6334 !End of the abilint section
6335 
6336  implicit none
6337 
6338 !Arguments -------------------------------
6339 !scalars
6340  integer,intent(in) :: natom,ntypat
6341 !arrays
6342  integer,intent(in) :: typat(natom)
6343  real(dp),intent(in) :: amu(ntypat)
6344  real(dp),intent(inout) :: mat(2*3*natom*3*natom)
6345 
6346 !Local variables -------------------------
6347 !scalars
6348  integer :: i1,i2,idir1,idir2,index,ipert1,ipert2
6349  real(dp),parameter :: break_symm=1.0d-12
6350  !real(dp),parameter :: break_symm=zero
6351  real(dp) :: fac
6352 !arrays
6353  real(dp) :: nearidentity(3,3)
6354 
6355 ! *********************************************************************
6356 
6357 !This slight breaking of the symmetry allows the
6358 !results to be more portable between machines
6359  nearidentity(:,:)=one
6360  nearidentity(1,1)=one+break_symm
6361  nearidentity(3,3)=one-break_symm
6362 
6363 !Include the masses in the dynamical matrix
6364  do ipert1=1,natom
6365    do ipert2=1,natom
6366      fac=1.0_dp/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
6367      do idir1=1,3
6368        do idir2=1,3
6369          i1=idir1+(ipert1-1)*3
6370          i2=idir2+(ipert2-1)*3
6371          index=i1+3*natom*(i2-1)
6372          mat(2*index-1)=mat(2*index-1)*fac*nearidentity(idir1,idir2)
6373          mat(2*index  )=mat(2*index  )*fac*nearidentity(idir1,idir2)
6374 !        This is to break slightly the translation invariance, and make
6375 !        the automatic tests more portable
6376          if(ipert1==ipert2 .and. idir1==idir2)then
6377            mat(2*index-1)=mat(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp
6378          end if
6379        end do
6380      end do
6381    end do
6382  end do
6383 
6384  ! Make the dynamical matrix hermitian
6385  call mkherm(mat,3*natom)
6386 
6387 end subroutine massmult_and_breaksym

m_dynmat/nanal9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

5568 subroutine nanal9(dyew,dynmat,iqpt,natom,nqpt,plus)
5569 
5570 
5571 !This section has been created automatically by the script Abilint (TD).
5572 !Do not modify the following lines by hand.
5573 #undef ABI_FUNC
5574 #define ABI_FUNC 'nanal9'
5575 !End of the abilint section
5576 
5577  implicit none
5578 
5579 !Arguments -------------------------------
5580 !scalars
5581  integer,intent(in) :: iqpt,natom,nqpt,plus
5582 !arrays
5583  real(dp),intent(in) :: dyew(2,3,natom,3,natom)
5584  real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)
5585 
5586 !Local variables -------------------------
5587 !scalars
5588  integer :: ia,ib,mu,nu
5589  character(len=500) :: message
5590 
5591 ! *********************************************************************
5592 
5593  if (plus==0) then
5594 
5595    do ia=1,natom
5596      do ib=1,natom
5597        do mu=1,3
5598          do nu=1,3
5599 !          The following four lines are OK
5600            dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) - dyew(1,mu,ia,nu,ib)
5601            dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) - dyew(2,mu,ia,nu,ib)
5602          end do
5603        end do
5604      end do
5605    end do
5606 
5607  else if (plus==1) then
5608 !  write(std_out,*)' nanal9 : now, only the analytic part '
5609 !  write(std_out,*)' nanal9 : now, only the ewald part '
5610    do ia=1,natom
5611      do ib=1,natom
5612        do mu=1,3
5613          do nu=1,3
5614            dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) + dyew(1,mu,ia,nu,ib)
5615            dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) + dyew(2,mu,ia,nu,ib)
5616          end do
5617        end do
5618      end do
5619    end do
5620 
5621  else
5622    write(message,'(a,a,a,i0,a)' )&
5623 &   'The argument "plus" must be equal to 0 or 1.',ch10,&
5624 &   'The value ',plus,' is not available.'
5625    MSG_BUG(message)
5626  end if
5627 
5628 end subroutine nanal9

m_dynmat/q0dy3_apply [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 dfpt_ewald, for both
 q==0 and the real wavelength.

 The q0dy3_apply should be used in conjunction with the subroutine
 dfpt_ewald (or ewald9):
 First, the call of dfpt_ewald 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 dfpt_ewald 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

SOURCE

2017 subroutine q0dy3_apply(natom,dyewq0,dyew)
2018 
2019 
2020 !This section has been created automatically by the script Abilint (TD).
2021 !Do not modify the following lines by hand.
2022 #undef ABI_FUNC
2023 #define ABI_FUNC 'q0dy3_apply'
2024 !End of the abilint section
2025 
2026  implicit none
2027 
2028 !Arguments -------------------------------
2029 !scalars
2030  integer,intent(in) :: natom
2031 !arrays
2032  real(dp),intent(in) :: dyewq0(3,3,natom)
2033  real(dp),intent(inout) :: dyew(2,3,natom,3,natom)
2034 
2035 !Local variables -------------------------
2036 !scalars
2037  integer :: ia,mu,nu
2038 
2039 ! *********************************************************************
2040 
2041  do mu=1,3
2042    do nu=1,3
2043      do ia=1,natom
2044        dyew(1,mu,ia,nu,ia)=dyew(1,mu,ia,nu,ia)-dyewq0(mu,nu,ia)
2045      end do
2046    end do
2047  end do
2048 
2049 end subroutine q0dy3_apply

m_dynmat/q0dy3_calc [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 dfpt_ewald, 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
 dfpt_ewald (or ewald9).
 First, the call of dfpt_ewald 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 dfpt_ewald 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

      ddb_hybrid,m_ifc,respfn

CHILDREN

SOURCE

2101 subroutine q0dy3_calc(natom,dyewq0,dyew,option)
2102 
2103 
2104 !This section has been created automatically by the script Abilint (TD).
2105 !Do not modify the following lines by hand.
2106 #undef ABI_FUNC
2107 #define ABI_FUNC 'q0dy3_calc'
2108 !End of the abilint section
2109 
2110  implicit none
2111 
2112 !Arguments -------------------------------
2113 !scalars
2114  integer,intent(in) :: natom,option
2115 !arrays
2116  real(dp),intent(in) :: dyew(2,3,natom,3,natom)
2117  real(dp),intent(out) :: dyewq0(3,3,natom)
2118 
2119 !Local variables -------------------------
2120 !scalars
2121  integer :: ia,ib,mu,nu
2122  character(len=500) :: message
2123 
2124 ! *********************************************************************
2125 
2126  if(option==1.or.option==2)then
2127    do mu=1,3
2128      do nu=1,3
2129        do ia=1,natom
2130          dyewq0(mu,nu,ia)=zero
2131          do ib=1,natom
2132            dyewq0(mu,nu,ia)=dyewq0(mu,nu,ia)+dyew(1,mu,ia,nu,ib)
2133          end do
2134        end do
2135      end do
2136    end do
2137  else
2138    write (message, '(3a)')&
2139 &   'option should be 1 or 2.',ch10,&
2140 &   'action: correct calling routine'
2141    MSG_BUG(message)
2142  end if
2143 
2144  if(option==2)then
2145    do ia=1,natom
2146      do mu=1,3
2147        do nu=mu,3
2148          dyewq0(mu,nu,ia)=(dyewq0(mu,nu,ia)+dyewq0(nu,mu,ia))/2
2149          dyewq0(nu,mu,ia)=dyewq0(mu,nu,ia)
2150        end do
2151      end do
2152    end do
2153  end if
2154 
2155 end subroutine q0dy3_calc

m_dynmat/symdm9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 symdm9

FUNCTION

 Use the set of special k points calculated by the Monkhorst & Pack Technique.
 Check if all the information 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 or -1 if non-stationary block
   2 or -2 if stationary block
   3 or -3 if third order derivatives
   positive if symmetries are used to set elements to zero whenever possible, negative to prevent this to happen.
 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

SOURCE

5058 subroutine symdm9(blkflg,blknrm,blkqpt,blktyp,blkval,&
5059 & dynmat,gprim,indsym,mpert,natom,nblok,nqpt,nsym,rfmeth,&
5060 & rprim,spqpt,symrec,symrel,qmissing)
5061 
5062 
5063 !This section has been created automatically by the script Abilint (TD).
5064 !Do not modify the following lines by hand.
5065 #undef ABI_FUNC
5066 #define ABI_FUNC 'symdm9'
5067 !End of the abilint section
5068 
5069  implicit none
5070 
5071 !Arguments -------------------------------
5072 !scalars
5073  integer,intent(in) :: mpert,natom,nblok,nqpt,nsym,rfmeth
5074 !arrays
5075  integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),blktyp(nblok)
5076  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
5077  integer,allocatable,optional,intent(out) :: qmissing(:)
5078  real(dp),intent(in) :: blknrm(3,nblok),blkqpt(9,nblok)
5079  real(dp),intent(in) :: blkval(2,3*mpert*3*mpert,nblok),gprim(3,3),rprim(3,3)
5080  real(dp),intent(in) :: spqpt(3,nqpt)
5081  real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt)
5082 
5083 !Local variables -------------------------
5084 !tol sets tolerance for equality of q points between those of
5085 !the DDB and those of the sampling grid
5086 !scalars
5087  integer :: ia,ib,iblok,idir1,idir2,ii,ipert1,ipert2,iqpt,isym,jj,kk,ll
5088  integer :: mu,nu,q1,q2,nqmiss
5089  real(dp),parameter :: tol=2.d-8
5090  real(dp) :: arg1,arg2,im,re,sumi,sumr
5091  logical :: allow_qmiss
5092  character(len=500) :: message
5093 !arrays
5094  integer,allocatable :: qtest(:,:)
5095  integer :: qmiss_(nqpt)
5096  real(dp) :: qq(3),qsym(6),ss(3,3)
5097  real(dp),allocatable :: ddd(:,:,:,:,:)
5098 
5099 ! *********************************************************************
5100  ! Initialize output (some q-points might not be reconstructed if qmissing is present)
5101  dynmat = huge(one)
5102  allow_qmiss = (present(qmissing))
5103 
5104  ABI_ALLOCATE(ddd,(2,3,natom,3,natom))
5105 !Check if the blkqpt points and their symmetrics are sufficient
5106 !in the DDB to retrieve all the q points of the B.Z. sampling
5107 
5108 !Initialization of a test variable
5109 ! qtest(iqpt,1)=iblok
5110 ! qtest(iqpt,2)=isym
5111 ! qtest(iqpt,3)=time_reversal
5112  ABI_ALLOCATE(qtest,(nqpt,3))
5113  do iqpt=1,nqpt
5114    qtest(iqpt,1)=0
5115  end do
5116 
5117 !Q points coming from the DDB
5118 !write(std_out,*)' Nbr. of Blocks -> ',nblok
5119 
5120  do iblok=1,nblok
5121 
5122    if (abs(blktyp(iblok))==abs(rfmeth)) then
5123      qq(1)=blkqpt(1,iblok)/blknrm(1,iblok)
5124      qq(2)=blkqpt(2,iblok)/blknrm(1,iblok)
5125      qq(3)=blkqpt(3,iblok)/blknrm(1,iblok)
5126 
5127 !    Calculation of the symmetric points (including Time Reversal)
5128      do isym=1,nsym
5129        qsym(1)=qq(1)*symrec(1,1,isym)+qq(2)*symrec(1,2,isym)+qq(3)*symrec(1,3,isym)
5130        qsym(2)=qq(1)*symrec(2,1,isym)+qq(2)*symrec(2,2,isym)+qq(3)*symrec(2,3,isym)
5131        qsym(3)=qq(1)*symrec(3,1,isym)+qq(2)*symrec(3,2,isym)+qq(3)*symrec(3,3,isym)
5132 
5133 !      Dont forget the Time Reversal symmetry
5134        qsym(4)=-qq(1)*symrec(1,1,isym)-qq(2)*symrec(1,2,isym)-qq(3)*symrec(1,3,isym)
5135        qsym(5)=-qq(1)*symrec(2,1,isym)-qq(2)*symrec(2,2,isym)-qq(3)*symrec(2,3,isym)
5136        qsym(6)=-qq(1)*symrec(3,1,isym)-qq(2)*symrec(3,2,isym)-qq(3)*symrec(3,3,isym)
5137 
5138 !      Comparison between the q points and their symmetric points
5139 !      and the set of q points which samples the entire Brillouin zone
5140        do iqpt=1,nqpt
5141 
5142          if (mod(abs(spqpt(1,iqpt)-qsym(1))+tol,1._dp)<2*tol)then
5143            if (mod(abs(spqpt(2,iqpt)-qsym(2))+tol,1._dp)<2*tol)then
5144              if (mod(abs(spqpt(3,iqpt)-qsym(3))+tol,1._dp)<2*tol)then
5145 
5146 !              write(std_out,*)' q point from the DDB ! '
5147 !              write(std_out,*)' block -> ',iblok
5148 !              write(std_out,*)' sym.  -> ',isym
5149 !              write(std_out,*)' No Time Reversal '
5150 !              write(std_out,*)'(',qsym(1),',',qsym(2),',',qsym(3),')'
5151 !              write(std_out,*)' '
5152 
5153                qtest(iqpt,1)=iblok
5154                qtest(iqpt,2)=isym
5155                qtest(iqpt,3)=0
5156              end if
5157            end if
5158          end if
5159 
5160          if (mod(abs(spqpt(1,iqpt)-qsym(4))+tol,1._dp)<2*tol)then
5161            if (mod(abs(spqpt(2,iqpt)-qsym(5))+tol,1._dp)<2*tol)then
5162              if (mod(abs(spqpt(3,iqpt)-qsym(6))+tol,1._dp)<2*tol)then
5163 
5164 !              write(std_out,*)' q point from the DDB ! '
5165 !              write(std_out,*)' block -> ',iblok
5166 !              write(std_out,*)' sym.  -> ',isym
5167 !              write(std_out,*)' Time Reversal '
5168 !              write(std_out,*)'(',qsym(4),',',qsym(5),',',qsym(6),')'
5169 !              write(std_out,*)' '
5170 
5171                qtest(iqpt,1)=iblok
5172                qtest(iqpt,2)=isym
5173                qtest(iqpt,3)=1
5174              end if
5175            end if
5176          end if
5177 
5178        end do ! End of the loop on the q points of the sampling
5179      end do ! End of the loop on the symmetries
5180 
5181    end if
5182  end do !  End of the loop on the q points of the DDB
5183 
5184 ! Check if all the information relatives to the q points sampling are found in the DDB;
5185 ! if not => stop message
5186  nqmiss = 0
5187  do iqpt=1,nqpt
5188    if (qtest(iqpt,1)==0) then
5189      nqmiss = nqmiss + 1
5190      qmiss_(nqmiss) = iqpt
5191      write(message, '(a,a,a)' )&
5192 &     ' symdm9 : the bloks found in the DDB are characterized',ch10,&
5193 &     '  by the following wavevectors :'
5194      call wrtout(std_out,message,'COLL')
5195      do iblok=1,nblok
5196        write(message, '(a,4d20.12)')&
5197 &       ' ',blkqpt(1,iblok),blkqpt(2,iblok),blkqpt(3,iblok),blknrm(1,iblok)
5198        call wrtout(std_out,message,'COLL')
5199      end do
5200      write(message, '(a,a,a,i0,a,a,a,3es16.6,a,a,a,a)' )&
5201 &     'Information is missing in the DDB.',ch10,&
5202 &     'The dynamical matrix number ',iqpt,' cannot be built,',ch10,&
5203 &     'since no block with wavevector',spqpt(1:3,iqpt),ch10,&
5204 &     'has been found.',ch10,&
5205 &     'Action: add the required blok in the DDB, or modify your input file.'
5206      if (.not.allow_qmiss) then
5207        MSG_ERROR(message)
5208      else
5209        !continue
5210        MSG_COMMENT(message)
5211      end if
5212    end if
5213  end do
5214 
5215  ! Will return a list with the index of the q-points that could not be symmetrized.
5216  if (allow_qmiss) then
5217    ABI_ALLOCATE(qmissing, (nqmiss))
5218    if (nqmiss > 0) qmissing = qmiss_(1:nqmiss)
5219  end if
5220 
5221 !Generation of the dynamical matrices relative to the q points
5222 !of the set which samples the entire Brillouin zone
5223  do iqpt=1,nqpt
5224    q1=qtest(iqpt,1)
5225    q2=qtest(iqpt,2)
5226    ! Skip this q-point if don't have enough info and allow_qmiss
5227    if (allow_qmiss .and. q1==0) cycle
5228 
5229 !  Check if the symmetry accompagnied with time reversal : q <- -q
5230    if (qtest(iqpt,3)==0) then
5231      do ii=1,3
5232        qq(ii)=blkqpt(ii,q1)/blknrm(1,q1)
5233        do jj=1,3
5234          ss(ii,jj)=zero
5235          do kk=1,3
5236            do ll=1,3
5237              ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2)
5238            end do
5239          end do
5240        end do
5241      end do
5242 !    write(std_out,*)"ss",ss
5243    else
5244      do ii=1,3
5245        qq(ii)=-blkqpt(ii,q1)/blknrm(1,q1)
5246        do jj=1,3
5247          ss(ii,jj)=zero
5248          do kk=1,3
5249            do ll=1,3
5250              ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2)
5251            end do
5252          end do
5253        end do
5254      end do
5255    end if
5256 !
5257 !  Check whether all the information is contained in the DDB
5258    do ipert2=1,natom
5259      do idir2=1,3
5260        do ipert1=1,natom
5261          do idir1=1,3
5262            if(blkflg(idir1,ipert1,idir2,ipert2,q1)/=1)then
5263              write(message, '(a,a,a,i6,a,a,a,4i4,a,a,a,a)' )&
5264 &             'Information are missing in the DDB.',ch10,&
5265 &             'In block',q1,' the following element is missing :',ch10,&
5266 &             'idir1,ipert1,idir2,ipert2=',idir1,ipert1,idir2,ipert2,ch10,&
5267 &             'Action: add the required information in the DDB,',ch10,&
5268 &             'or modify your input file.'
5269              MSG_ERROR(message)
5270            end if
5271          end do
5272        end do
5273      end do
5274    end do
5275 
5276 !  Read the dynamical matrices in the DDB
5277    do ipert2=1,natom
5278      do idir2=1,3
5279        do ipert1=1,natom
5280          do idir1=1,3
5281            ddd(:,idir1,ipert1,idir2,ipert2)=blkval(:,idir1+3*(ipert1-1+mpert*(idir2-1+3*(ipert2-1))),q1)
5282          end do
5283        end do
5284      end do
5285    end do
5286 
5287 !  Calculation of the dynamical matrix of a symmetrical q point
5288    do ia=1,natom
5289      do ib=1,natom
5290 !      write(std_out,*)'atom-> ',ia,indsym(4,q2,ia)
5291 !      write(std_out,*)'atom-> ',ib,indsym(4,q2,ib)
5292        arg1=two_pi*(qq(1)*indsym(1,q2,ia)+qq(2)*indsym(2,q2,ia)+qq(3)*indsym(3,q2,ia))
5293        arg2=two_pi*(qq(1)*indsym(1,q2,ib)+qq(2)*indsym(2,q2,ib)+qq(3)*indsym(3,q2,ib))
5294        re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
5295        im=cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)
5296 !      write(std_out,*)'re : ',re, 'im : ',im
5297        do mu=1,3
5298          do nu=1,3
5299 !          write(std_out,*)' '
5300            sumr=zero
5301            sumi=zero
5302            do ii=1,3
5303              do jj=1,3
5304 !              If there is Time Reversal : D.M. <- Complex Conjugate D.M.
5305                if (qtest(iqpt,3)==0) then
5306                  sumr=sumr+ss(mu,ii)*ss(nu,jj)*ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
5307                  sumi=sumi+ss(mu,ii)*ss(nu,jj)*ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
5308                else
5309                  sumr=sumr+ss(mu,ii)*ss(nu,jj)*ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
5310                  sumi=sumi-ss(mu,ii)*ss(nu,jj)*ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib))
5311                end if
5312              end do
5313            end do
5314 
5315 !          Dynmat -> Dynamical Matrix for the q point of the sampling
5316 !          write(std_out,*)' Sumr -> ',mu,nu,sumr
5317 !          write(std_out,*)' Sumi -> ',mu,nu,sumi
5318            dynmat(1,mu,ia,nu,ib,iqpt)=re*sumr-im*sumi
5319            dynmat(2,mu,ia,nu,ib,iqpt)=re*sumi+im*sumr
5320 
5321 !          DEBUG
5322 !          if((ia==2 .or. ia==3) .and. ib==1)then
5323 !          write(std_out,'(5i3,2es16.8)' )mu,ia,nu,ib,iqpt,dynmat(1:2,mu,ia,nu,ib,iqpt)
5324 !          end if
5325 !          ENDDEBUG
5326 
5327          end do ! End loop on the coordinates
5328        end do
5329 
5330      end do ! End loop on the ia atoms
5331    end do ! End loop on the ib atoms
5332  end do ! End loop on the q points of the sampling
5333 
5334  ABI_DEALLOCATE(ddd)
5335  ABI_DEALLOCATE(qtest)
5336 
5337 end subroutine symdm9

m_dynmat/symdyma [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 symdyma

FUNCTION

 Symmetrize the dynamical matrices

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.
 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 in cartesian coordinates 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) [[cite:Hendrikse1995]]

TODO

 A full description of the equations should be included

PARENTS

      m_dynmat,m_phgamma,relaxpol

CHILDREN

SOURCE

2199 subroutine symdyma(dmati,indsym,natom,nsym,qptn,rprimd,symrel,symafm)
2200 
2201 
2202 !This section has been created automatically by the script Abilint (TD).
2203 !Do not modify the following lines by hand.
2204 #undef ABI_FUNC
2205 #define ABI_FUNC 'symdyma'
2206 !End of the abilint section
2207 
2208  implicit none
2209 
2210 !Arguments -------------------------------
2211 !scalars
2212  integer,intent(in) :: natom,nsym
2213 !arrays
2214  integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
2215  integer,intent(in) :: symafm(nsym)
2216  real(dp),intent(in) :: qptn(3),rprimd(3,3)
2217  real(dp),intent(inout) :: dmati(2*3*natom*3*natom)
2218 
2219 !Local variables -------------------------
2220 !scalars
2221  integer :: i1,i2,iat,idir,ii,index,isgn,isym,itirev,jat,jdir,jj,kk,ll
2222  integer :: niat,njat,timrev
2223  real(dp) :: arg1,arg2,dmint,im,re,sumi,sumr
2224 !arrays
2225  integer :: indij(natom,natom),symq(4,2,nsym),symrec(3,3,nsym)
2226  real(dp) :: TqR(3,3),TqS_(3,3),dynmat(2,3,natom,3,natom)
2227  real(dp) :: dynmatint(2*nsym,2,3,natom,3,natom),gprimd(3,3)
2228  real(dp) :: symcart(3,3,nsym)
2229 
2230 ! *********************************************************************
2231 
2232 !0) initializations
2233  call matr3inv(rprimd,gprimd)
2234  do isym=1,nsym
2235    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2236  end do
2237 
2238  TqR=zero
2239  TqS_=zero
2240  dynmat=zero
2241 
2242 !Note: dynmat is used as work space here
2243  i1=0
2244  do iat=1,natom
2245    do idir=1,3
2246      i1=i1+1
2247      i2=0
2248      do jat=1,natom
2249        do jdir=1,3
2250          i2=i2+1
2251          index=i1+3*natom*(i2-1)
2252          dynmat(1,idir,iat,jdir,jat)=dmati(2*index-1)
2253          dynmat(2,idir,iat,jdir,jat)=dmati(2*index  )
2254        end do
2255      end do
2256    end do
2257  end do
2258 
2259 !Transform symrel to cartesian coordinates (RC coding)
2260 !do isym=1,nsym
2261 !symcart(:,:,isym)=matmul(rprimd,matmul(dble(symrel(:,:,isym)),gprimd))
2262 !end do
2263 
2264 !Coding from symdm9
2265  do isym=1,nsym
2266    do jj=1,3
2267      symcart(:,jj,isym)=zero
2268      do kk=1,3
2269        do ll=1,3
2270          symcart(:,jj,isym)=symcart(:,jj,isym)+rprimd(:,kk)*gprimd(jj,ll)*symrel(kk,ll,isym)
2271        end do
2272      end do
2273    end do
2274  end do
2275 
2276 
2277 !Get the symq of the CURRENT Q POINT
2278 !mjv: set prtvol=0 for production runs.
2279  call littlegroup_q(nsym,qptn,symq,symrec,symafm,timrev,prtvol=0)
2280 
2281  indij(:,:)=0
2282  dynmatint=zero
2283 
2284  do isym=1,nsym                                 ! loop over all the symmetries
2285 !  write(std_out,*) 'current symmetry',isym
2286    do itirev=1,2                                 ! loop over the time-reversal symmetry
2287      isgn=3-2*itirev                             ! to take into accont the time-reversal
2288 !    write(std_out,*) 'timereversal',isgn
2289      if (symq(4,itirev,isym)==1) then             ! isym belongs to the wave vector point group
2290 !      write(std_out,*) 'isym belongs to the wave vector point group'
2291        do iat=1,natom                              ! loop over the atoms
2292          do jat=1,natom                            ! loop over the atoms
2293            niat=indsym(4,isym,iat)                   ! niat={R|t}iat
2294            njat=indsym(4,isym,jat)                   ! njat={R|t}jat
2295            indij(niat,njat)=indij(niat,njat)+1
2296 !          write(std_out,'(a,5i5)') 'current status:',iat,jat,niat,njat,indij(niat,njat)
2297 !          phase calculation, arg1 and arg2 because of two-atom derivative
2298            arg1=two_pi*( qptn(1)*indsym(1,isym,iat)+&
2299 &           qptn(2)*indsym(2,isym,iat)+&
2300 &           qptn(3)*indsym(3,isym,iat) )
2301            arg2=two_pi*( qptn(1)*indsym(1,isym,jat)+&
2302 &           qptn(2)*indsym(2,isym,jat)+&
2303 &           qptn(3)*indsym(3,isym,jat) )
2304            re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
2305            im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))
2306 
2307            do idir=1,3                               ! loop over displacements
2308              do jdir=1,3                              ! loop over displacements
2309 !              we pick the (iat,jat) (3x3) block of the dyn.mat.
2310                sumr=zero
2311                sumi=zero
2312                do ii=1,3
2313                  do jj=1,3
2314                    sumr=sumr+symcart(idir,ii,isym)*dynmat(1,ii,niat,jj,njat)*symcart(jdir,jj,isym)
2315                    sumi=sumi+symcart(idir,ii,isym)*dynmat(2,ii,niat,jj,njat)*symcart(jdir,jj,isym)
2316                  end do
2317                end do
2318                sumi=isgn*sumi
2319 
2320                dynmatint(nsym*(itirev-1)+isym,1,idir,iat,jdir,jat)=re*sumr-im*sumi
2321                dynmatint(nsym*(itirev-1)+isym,2,idir,iat,jdir,jat)=re*sumi+im*sumr
2322 
2323              end do
2324            end do
2325          end do
2326        end do ! end treatment of the (iat,jat) (3x3) block of dynmat
2327      end if ! symmetry check
2328    end do ! time-reversal
2329  end do ! symmetries
2330 
2331 
2332 !4) make the average, get the final symmetric dynamical matrix
2333  do iat=1,natom
2334    do jat=1,natom
2335      do idir=1,3
2336        do jdir=1,3
2337          dmint=zero
2338          do isym=1,2*nsym
2339            dmint=dmint+dynmatint(isym,1,idir,iat,jdir,jat)
2340          end do
2341          dynmat(1,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
2342          dmint=zero
2343          do isym=1,2*nsym
2344            dmint=dmint+dynmatint(isym,2,idir,iat,jdir,jat)
2345          end do
2346          dynmat(2,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
2347        end do
2348      end do
2349    end do
2350  end do
2351 
2352  i1=0
2353  do iat=1,natom
2354    do idir=1,3
2355      i1=i1+1
2356      i2=0
2357      do jat=1,natom
2358        do jdir=1,3
2359          i2=i2+1
2360          index=i1+3*natom*(i2-1)
2361          dmati(2*index-1)=dynmat(1,idir,iat,jdir,jat)
2362          dmati(2*index  )=dynmat(2,idir,iat,jdir,jat)
2363        end do
2364      end do
2365    end do
2366  end do
2367 
2368 end subroutine symdyma

m_dynmat/sytens [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 sytens

FUNCTION

 Determines the set of irreductible elements of the non-linear
 optical susceptibility and Raman tensors

INPUTS

  indsym(4,nsym,natom)=indirect indexing array described above: 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)

OUTPUT

  (see side effects)

SIDE EFFECTS

  rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
       that have to be computed
    At the input :
       1   ->   element has to be computed explicitely
    At the output :
       1   ->   element has to be computed explicitely
      -1   ->   use symmetry operations to obtain the corresponding element
      -2   ->   element is zero by symmetry

PARENTS

      m_ddb,nonlinear,respfn

CHILDREN

SOURCE

4785 subroutine sytens(indsym,mpert,natom,nsym,rfpert,symrec,symrel)
4786 
4787  use defs_basis
4788  use m_abicore
4789 
4790 !This section has been created automatically by the script Abilint (TD).
4791 !Do not modify the following lines by hand.
4792 #undef ABI_FUNC
4793 #define ABI_FUNC 'sytens'
4794 !End of the abilint section
4795 
4796  implicit none
4797 
4798 !Arguments -------------------------------
4799 !scalars
4800  integer,intent(in) :: mpert,natom,nsym
4801 !arrays
4802  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4803  integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert)
4804 
4805 !Local variables -------------------------
4806 !scalars
4807  integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_
4808  integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2
4809  integer :: ipesy3,isym
4810 !arrays
4811  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
4812  integer,allocatable :: pertsy(:,:,:,:,:,:)
4813 
4814 !***********************************************************************
4815 
4816 !DEBUG
4817 !write(std_out,*)'sytens : enter'
4818 !write(std_out,*)'indsym = '
4819 !write(std_out,*)indsym
4820 !stop
4821 !ENDDEBUG
4822 
4823  ABI_ALLOCATE(pertsy,(3,mpert,3,mpert,3,mpert))
4824  pertsy(:,:,:,:,:,:) = 0
4825 
4826 !Loop over perturbations
4827 
4828  do i1pert_ = 1, mpert
4829    do i2pert_ = 1, mpert
4830      do i3pert_ = 1, mpert
4831 
4832        do i1dir_ = 1, 3
4833          do i2dir_ = 1, 3
4834            do i3dir_ = 1, 3
4835 
4836              i1pert = (mpert - i1pert_ + 1)
4837              if (i1pert <= natom) i1pert = natom + 1 - i1pert
4838              i2pert = (mpert - i2pert_ + 1)
4839              if (i2pert <= natom) i2pert = natom + 1 - i2pert
4840              i3pert = (mpert - i3pert_ + 1)
4841              if (i3pert <= natom) i3pert = natom + 1 - i3pert
4842 
4843              if (i1pert <= natom) then
4844                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
4845              else if (i2pert <= natom) then
4846                i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_
4847              else if (i3pert <= natom) then
4848                i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_
4849              else
4850                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
4851              end if
4852 
4853              if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then
4854 
4855 !              Loop over all symmetries
4856 
4857                flag = 0
4858                do isym = 1, nsym
4859 
4860                  found = 1
4861 
4862 !                Select the symmetric element of i1pert,i2pert,i3pert
4863 
4864                  if (i1pert <= natom) then
4865                    ipesy1 = indsym(4,isym,i1pert)
4866                    sym1(:,:) = symrec(:,:,isym)
4867                  else if (i1pert == natom + 2) then
4868                    ipesy1 = i1pert
4869                    sym1(:,:) = symrel(:,:,isym)
4870                  else
4871                    found = 0
4872                  end if
4873 
4874                  if (i2pert <= natom) then
4875                    ipesy2 = indsym(4,isym,i2pert)
4876                    sym2(:,:) = symrec(:,:,isym)
4877                  else if (i2pert == natom + 2) then
4878                    ipesy2 = i2pert
4879                    sym2(:,:) = symrel(:,:,isym)
4880                  else
4881                    found = 0
4882                  end if
4883 
4884                  if (i3pert <= natom) then
4885                    ipesy3 = indsym(4,isym,i3pert)
4886                    sym3(:,:) = symrec(:,:,isym)
4887                  else if (i3pert == natom + 2) then
4888                    ipesy3 = i3pert
4889                    sym3(:,:) = symrel(:,:,isym)
4890                  else
4891                    found = 0
4892                  end if
4893 
4894 !                See if the symmetric element is available and check if some
4895 !                of the elements may be zeor. In the latter case, they do not need
4896 !                to be computed.
4897 
4898 
4899                  if ((flag /= -1).and.&
4900 &                 (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
4901                    flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir)
4902                  end if
4903 
4904 
4905                  do idisy1 = 1, 3
4906                    do idisy2 = 1, 3
4907                      do idisy3 = 1, 3
4908 
4909                        if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.&
4910 &                       (sym3(i3dir,idisy3) /= 0)) then
4911                          if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
4912                            found = 0
4913 !                          exit      ! exit loop over symmetries
4914                          end if
4915                        end if
4916 
4917 
4918                        if ((flag == -1).and.&
4919 &                       ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then
4920                          if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.&
4921 &                         (sym3(i3dir,idisy3)/=0)) then
4922                            flag = 0
4923                          end if
4924                        end if
4925 
4926 
4927 
4928                      end do
4929                    end do
4930                  end do
4931 
4932 
4933                  if (found == 1) then
4934                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1
4935                  end if
4936 
4937 
4938 !                In case a symmetry operation only changes the sign of an
4939 !                element, this element has to be equal to zero
4940 
4941                  if (flag == -1) then
4942                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2
4943                    exit
4944                  end if
4945 
4946                end do    ! close loop on symmetries
4947 
4948 
4949 
4950 !              If the elemetn i1pert,i2pert,i3pert is not symmetric
4951 !              to a basis element, it is a basis element
4952 
4953                if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then
4954                  pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4955                end if
4956 
4957              end if ! rfpert /= 0
4958 
4959            end do        ! close loop over perturbations
4960          end do
4961        end do
4962      end do
4963    end do
4964  end do
4965 
4966 !Now, take into account the permutation of (i1pert,i1dir)
4967 !and (i3pert,i3dir)
4968 
4969  do i1pert = 1, mpert
4970    do i2pert = 1, mpert
4971      do i3pert = 1, mpert
4972 
4973        do i1dir = 1, 3
4974          do i2dir = 1, 3
4975            do i3dir = 1, 3
4976 
4977              if ((i1pert /= i3pert).or.(i1dir /= i3dir)) then
4978 
4979                if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.&
4980 &               (pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) == 1)) then
4981                  pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = -1
4982                end if
4983 
4984              end if
4985 
4986            end do
4987          end do
4988        end do
4989 
4990      end do
4991    end do
4992  end do
4993 
4994  rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:)
4995 
4996  ABI_DEALLOCATE(pertsy)
4997 
4998 
4999 end subroutine sytens

m_dynmat/wght9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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 or -1=S.C.;2=F.C.C.;4=Hex. -1 is for old algo to find weights, =1 is for Wigner-Seitz algo)
 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))
 rprimd(3,3)=dimensional primitive translations for real space (bohr)

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

SOURCE

4126 subroutine wght9(brav,gprim,natom,ngqpt,nqpt,nqshft,nrpt,qshft,rcan,rpt,rprimd,r_inscribed_sphere,wghatm)
4127 
4128 
4129 !This section has been created automatically by the script Abilint (TD).
4130 !Do not modify the following lines by hand.
4131 #undef ABI_FUNC
4132 #define ABI_FUNC 'wght9'
4133 !End of the abilint section
4134 
4135  implicit none
4136 
4137 !Arguments -------------------------------
4138 !scalars
4139  integer,intent(in) :: brav,natom,nqpt,nqshft,nrpt
4140  real(dp),intent(out) :: r_inscribed_sphere
4141 !arrays
4142  integer,intent(inout) :: ngqpt(9)
4143  real(dp),intent(in) :: gprim(3,3),qshft(3,4),rcan(3,natom),rpt(3,nrpt),rprimd(3,3)
4144  real(dp),intent(out) :: wghatm(natom,natom,nrpt)
4145 
4146 !Local variables -------------------------
4147 !scalars
4148  integer :: ia,ib,ii,jj,kk,iqshft,irpt,jqshft,nbordh,tok,new_wght,nptws,nreq
4149  integer :: idir
4150  real(dp), parameter :: tolsym=tol8
4151  real(dp) :: factor,sumwght,normsq,proj
4152  character(len=500) :: message
4153 !arrays
4154  integer :: nbord(9)
4155  real(dp) :: rdiff(9),red(3,3),ptws(4, 729),pp(3),rdiff_tmp(3)
4156 
4157 ! *********************************************************************
4158 
4159  DBG_ENTER("COLL")
4160 
4161 !First analyze the vectors qshft
4162  if(nqshft/=1)then
4163 
4164    if(brav==4)then
4165      write(message,'(a,a,a,i5,a,a,a)' )&
4166 &     'For the time being, only nqshft=1',ch10,&
4167 &     'is allowed with brav=4, while it is nqshft=',nqshft,'.',ch10,&
4168 &     'Action: in the input file, correct either brav or nqshft.'
4169      MSG_ERROR(message)
4170    end if
4171 
4172    if(nqshft==2)then
4173 !    Make sure that the q vectors form a BCC lattice
4174      do ii=1,3
4175        if(abs(abs(qshft(ii,1)-qshft(ii,2))-.5_dp)>1.d-10)then
4176          write(message, '(a,a,a,a,a,a,a)' )&
4177 &         'The test of the q1shft vectors shows that they',ch10,&
4178 &         'do not generate a body-centered lattice, which',ch10,&
4179 &         'is mandatory for nqshft=2.',ch10,&
4180 &         'Action: change the q1shft vectors in your input file.'
4181          MSG_ERROR(message)
4182        end if
4183      end do
4184    else if(nqshft==4)then
4185 !    Make sure that the q vectors form a FCC lattice
4186      do iqshft=1,3
4187        do jqshft=iqshft+1,4
4188          tok=0
4189          do ii=1,3
4190 !          Test on the presence of a +-0.5 difference
4191            if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-.5_dp) <1.d-10) tok=tok+1
4192 !          Test on the presence of a 0 or +-1.0 difference
4193            if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-1._dp) <1.d-10  .or.&
4194 &             abs(qshft(ii,iqshft)-qshft(ii,jqshft)) < 1.d-10) tok=tok+4
4195          end do
4196 !        Test 1 should be satisfied twice, and test 2 once
4197          if(tok/=6)then
4198            write(message, '(7a)' )&
4199 &           'The test of the q1shft vectors shows that they',ch10,&
4200 &           'do not generate a face-centered lattice, which',ch10,&
4201 &           'is mandatory for nqshft=4.',ch10,&
4202 &           'Action: change the q1shft vectors in your input file.'
4203            MSG_ERROR(message)
4204          end if
4205        end do
4206      end do
4207    else
4208      write(message, '(a,i4,3a)' )&
4209 &     'nqshft must be 1, 2 or 4. It is nqshft=',nqshft,'.',ch10,&
4210 &     'Action: change nqshft in your input file.'
4211      MSG_ERROR(message)
4212    end if
4213 
4214  end if
4215 
4216  factor=0.5_dp
4217  if(brav==2 .or. brav==3)then
4218    factor=0.25_dp
4219  end if
4220  if(nqshft/=1)factor=factor*2
4221 
4222  if (brav==1) then
4223 
4224    ! Does not support multiple shifts
4225    if (nqshft/=1) then
4226      write(message, '(a)' ) 'This version of the weights does not support nqshft/=1.'
4227      MSG_ERROR(message)
4228    end if
4229 
4230    ! Find the points of the lattice given by ngqpt*acell. These are used to define
4231    ! a Wigner-Seitz cell around the origin. The origin is excluded from the list.
4232    ! TODO : in principle this should be only -1 to +1 for ii jj kk!
4233    nptws=0
4234    do ii=-2,2
4235      do jj=-2,2
4236        do kk=-2,2
4237          do idir=1,3
4238            pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3)
4239          end do
4240          normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3)
4241          if (normsq > tol6) then
4242            nptws = nptws + 1
4243            ptws(:3,nptws) = pp(:)
4244            ptws(4,nptws) = half*normsq
4245          end if
4246        end do
4247      end do
4248    end do
4249  end if ! end new_wght
4250 
4251 !DEBUG
4252 !write(std_out,*)'factor,ngqpt',factor,ngqpt(1:3)
4253 !ENDDEBUG
4254 
4255  r_inscribed_sphere = sum((matmul(rprimd(:,:),ngqpt(1:3)))**2)
4256  do ii=-1,1
4257    do jj=-1,1
4258      do kk=-1,1
4259        if (ii==0 .and. jj==0 .and. kk==0) cycle
4260        do idir=1,3
4261          pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3)
4262        end do
4263        normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3)
4264        r_inscribed_sphere = min(r_inscribed_sphere, normsq)
4265      end do
4266    end do
4267  end do
4268  r_inscribed_sphere = sqrt(r_inscribed_sphere)
4269 
4270 
4271 !Begin the big loop on ia and ib
4272  do ia=1,natom
4273    do ib=1,natom
4274 
4275 !    Simple Lattice
4276      if (abs(brav)==1) then
4277 !      In this case, it is better to work in reduced coordinates
4278 !      As rcan is in canonical coordinates, => multiplication by gprim
4279        do ii=1,3
4280          red(1,ii)=  rcan(1,ia)*gprim(1,ii) +rcan(2,ia)*gprim(2,ii) +rcan(3,ia)*gprim(3,ii)
4281          red(2,ii)=  rcan(1,ib)*gprim(1,ii) +rcan(2,ib)*gprim(2,ii) +rcan(3,ib)*gprim(3,ii)
4282        end do
4283      end if
4284 
4285      do irpt=1,nrpt
4286 
4287 !      Initialization of the weights to 1.0
4288        wghatm(ia,ib,irpt)=1.0_dp
4289 
4290 !      Compute the difference vector
4291 
4292 !      Simple Cubic Lattice
4293        if (abs(brav)==1) then
4294 !        Change of rpt to reduced coordinates
4295          do ii=1,3
4296            red(3,ii)=  rpt(1,irpt)*gprim(1,ii) +rpt(2,irpt)*gprim(2,ii) +rpt(3,irpt)*gprim(3,ii)
4297            rdiff(ii)=red(2,ii)-red(1,ii)+red(3,ii)
4298          end do
4299          if (brav==1) then
4300            ! rdiff in cartesian coordinates
4301            do ii=1,3
4302              rdiff_tmp(ii)=rdiff(1)*rprimd(ii,1)+rdiff(2)*rprimd(ii,2)+rdiff(3)*rprimd(ii,3)
4303            end do
4304            rdiff(1:3)=rdiff_tmp(1:3)
4305          end if
4306 !        Other lattices
4307        else
4308          do ii=1,3
4309            rdiff(ii)=rcan(ii,ib)-rcan(ii,ia)+rpt(ii,irpt)
4310          end do
4311        end if
4312 
4313 !      ***************************************************************
4314 
4315 !      Assignement of weights
4316 
4317        if(nqshft==1 .and. brav/=4)then
4318 
4319          if (brav/=1) then
4320            do ii=1,3
4321 !            If the rpt vector is greater than
4322 !            the allowed space => weight = 0.0
4323              if (abs(rdiff(ii))-tol10>factor*ngqpt(ii)) then
4324                wghatm(ia,ib,irpt)=zero
4325              else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4326 !              If the point is in a boundary position => weight/2
4327                wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4328              end if
4329            end do
4330          else
4331            ! new weights
4332            wghatm(ia,ib,irpt)=zero
4333            nreq = 1
4334            do ii=1,nptws
4335              proj = rdiff(1)*ptws(1,ii)+rdiff(2)*ptws(2,ii)+rdiff(3)*ptws(3,ii)
4336              ! if rdiff closer to ptws than the origin the weight is zero
4337              ! if rdiff close to the origin with respect to all the other ptws the weight is 1
4338              ! if rdiff is equidistant from the origin and N other ptws the weight is 1/(N+1)
4339              if (proj-ptws(4,ii)>tolsym) then
4340                nreq = 0
4341                EXIT
4342              else if(abs(proj-ptws(4,ii))<=tolsym) then
4343                nreq=nreq+1
4344              end if
4345            end do
4346            if (nreq>0) then
4347              wghatm(ia,ib,irpt)=one/DBLE(nreq)
4348            end if
4349          end if
4350 
4351        else if(brav==4)then
4352 !        Hexagonal
4353 
4354 !        Examination of the X and Y boundaries in order to form an hexagon
4355 !        First generate the relevant boundaries
4356          rdiff(4)=0.5_dp*( rdiff(1)+sqrt(3.0_dp)*rdiff(2) )
4357          ngqpt(4)=ngqpt(1)
4358          rdiff(5)=0.5_dp*( rdiff(1)-sqrt(3.0_dp)*rdiff(2) )
4359          ngqpt(5)=ngqpt(1)
4360 
4361 !        Test the four inequalities
4362          do ii=1,5
4363            if(ii/=2)then
4364 
4365              nbord(ii)=0
4366 !            If the rpt vector is greater than
4367 !            the allowed space => weight = 0.0
4368              if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4369                wghatm(ia,ib,irpt)=zero
4370              else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4371 !              If the point is in a boundary position increment nbord(ii)
4372                nbord(ii)=1
4373              end if
4374 
4375            end if
4376          end do
4377 
4378 !        Computation of weights
4379          nbordh=nbord(1)+nbord(4)+nbord(5)
4380          if (nbordh==1) then
4381            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4382          else if (nbordh==2) then
4383            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4384          else if (nbordh/=0) then
4385            message = 'There is a problem of borders and weights (hex).'
4386            MSG_BUG(message)
4387          end if
4388          if (nbord(3)==1)then
4389            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4390          end if
4391 
4392 !        BCC packing of k-points
4393        else if(nqshft==2 .and. brav/=4)then
4394 
4395 !        First, generate the relevant boundaries
4396          rdiff(4)= rdiff(1)+rdiff(2)
4397          rdiff(5)= rdiff(1)-rdiff(2)
4398          rdiff(6)= rdiff(1)+rdiff(3)
4399          rdiff(7)= rdiff(1)-rdiff(3)
4400          rdiff(8)= rdiff(3)+rdiff(2)
4401          rdiff(9)= rdiff(3)-rdiff(2)
4402          if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
4403            write(message, '(a,a,a,3i6,a,a,a,a)' )&
4404 &           'In the BCC case, the three ngqpt numbers ',ch10,&
4405 &           '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
4406 &           'should be equal.',ch10,&
4407 &           'Action: use identical ngqpt(1:3) in your input file.'
4408            MSG_ERROR(message)
4409          end if
4410          do ii=4,9
4411            ngqpt(ii)=ngqpt(1)
4412          end do
4413 
4414 !        Test the relevant inequalities
4415          nbord(1)=0
4416          do ii=4,9
4417 !          If the rpt vector is greater than the allowed space => weight = 0.0
4418            if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4419              wghatm(ia,ib,irpt)=zero
4420            else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4421 !            If the point is in a boundary position increment nbord(1)
4422              nbord(1)=nbord(1)+1
4423            end if
4424          end do
4425 
4426 !        Computation of weights
4427          if (nbord(1)==1) then
4428            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4429          else if (nbord(1)==2) then
4430            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4431          else if (nbord(1)==3) then
4432            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
4433          else if (nbord(1)==4) then
4434            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/6
4435          else if (nbord(1)/=0) then
4436            message = ' There is a problem of borders and weights (BCC).'
4437            MSG_ERROR(message)
4438          end if
4439 
4440 !        FCC packing of k-points
4441        else if(nqshft==4 .and. brav/=4)then
4442 
4443 !        First, generate the relevant boundaries
4444          rdiff(4)= (rdiff(1)+rdiff(2)+rdiff(3))*2._dp/3._dp
4445          rdiff(5)= (rdiff(1)-rdiff(2)+rdiff(3))*2._dp/3._dp
4446          rdiff(6)= (rdiff(1)+rdiff(2)-rdiff(3))*2._dp/3._dp
4447          rdiff(7)= (rdiff(1)-rdiff(2)-rdiff(3))*2._dp/3._dp
4448          if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
4449            write(message, '(a,a,a,3i6,a,a,a,a)' )&
4450 &           'In the FCC case, the three ngqpt numbers ',ch10,&
4451 &           '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
4452 &           'should be equal.',ch10,&
4453 &           'Action: use identical ngqpt(1:3) in your input file.'
4454            MSG_ERROR(message)
4455          end if
4456          do ii=4,7
4457            ngqpt(ii)=ngqpt(1)
4458          end do
4459 
4460 !        Test the relevant inequalities
4461          nbord(1)=0
4462          do ii=1,7
4463 
4464 !          If the rpt vector is greater than the allowed space => weight = 0.0
4465            if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4466              wghatm(ia,ib,irpt)=zero
4467 !            If the point is in a boundary position increment nbord(1)
4468            else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4469              nbord(1)=nbord(1)+1
4470            end if
4471          end do
4472 
4473 !        Computation of weights
4474          if (nbord(1)==1) then
4475            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4476          else if (nbord(1)==2) then
4477            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4478          else if (nbord(1)==3) then
4479            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
4480          else if (nbord(1)/=0 .and. wghatm(ia,ib,irpt)>1.d-10) then
4481 !          Interestingly nbord(1)==4 happens for some points outside of the volume
4482            message = ' There is a problem of borders and weights (FCC).'
4483            MSG_BUG(message)
4484          end if
4485 
4486        else
4487          write(message, '(3a,i0,a)' )&
4488 &         'One should not arrive here ... ',ch10,&
4489 &         'The value nqshft ',nqshft,' is not available'
4490          MSG_BUG(message)
4491        end if
4492      end do ! Assignement of weights is done
4493 
4494    end do ! End of the double loop on ia and ib
4495  end do
4496 
4497 ! Check the results
4498  do ia=1,natom
4499    do ib=1,natom
4500      sumwght=zero
4501      do irpt=1,nrpt
4502 !      Check if the sum of the weights is equal to the number of q points
4503        sumwght=sumwght+wghatm(ia,ib,irpt)
4504 !      write(std_out,'(a,3i5)' )' atom1, atom2, irpt ; rpt ; wghatm ',ia,ib,irpt
4505 !      write(std_out,'(3es16.6,es18.6)' )&
4506 !      &    rpt(1,irpt),rpt(2,irpt),rpt(3,irpt),wghatm(ia,ib,irpt)
4507      end do
4508      if (abs(sumwght-nqpt)>tol10) then
4509        write(message, '(a,a,a,2i4,a,a,es14.4,a,a,i4)' )&
4510 &       'The sum of the weight is not equal to nqpt.',ch10,&
4511 &       'atoms :',ia,ib,ch10,&
4512 &       'The sum of the weights is : ',sumwght,ch10,&
4513 &       'The number of q points is : ',nqpt
4514        call wrtout(std_out,message,'COLL')
4515        write(message, '(13a)')&
4516 &       'This might have several sources.',ch10,&
4517 &       'If tolsym is larger than 1.0e-8, the atom positions might be loose',ch10,&
4518 &       'and the q point weights not computed properly.',ch10,&
4519 &       'Action: make input atomic positions more symmetric.',ch10,&
4520 &       'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine, and recompile.',ch10,&
4521 &       'Actually, this can also happen when ngqpt is 0 0 0,',ch10,&
4522 &       'if abs(brav)/=1, in which case you should change brav to 1.'
4523        MSG_BUG(message)
4524      end if
4525    end do
4526  end do
4527 
4528  DBG_EXIT("COLL")
4529 
4530 end subroutine wght9

m_dynmat/wings3 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

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

SOURCE

2715 subroutine wings3(carflg,d2cart,mpert)
2716 
2717 
2718 !This section has been created automatically by the script Abilint (TD).
2719 !Do not modify the following lines by hand.
2720 #undef ABI_FUNC
2721 #define ABI_FUNC 'wings3'
2722 !End of the abilint section
2723 
2724  implicit none
2725 
2726 !Arguments -------------------------------
2727 !scalars
2728  integer,intent(in) :: mpert
2729 !arrays
2730  integer,intent(inout) :: carflg(3,mpert,3,mpert)
2731  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
2732 
2733 !Local variables -------------------------
2734 !scalars
2735  integer :: idir,idir1,ipert,ipert1
2736 
2737 ! *********************************************************************
2738 
2739  do ipert=1,mpert
2740    do idir=1,3
2741      if(carflg(idir,ipert,idir,ipert)==0)then
2742        do ipert1=1,mpert
2743          do idir1=1,3
2744            carflg(idir,ipert,idir1,ipert1)=0
2745            carflg(idir1,ipert1,idir,ipert)=0
2746            d2cart(1,idir,ipert,idir1,ipert1)=zero
2747            d2cart(2,idir,ipert,idir1,ipert1)=zero
2748            d2cart(1,idir1,ipert1,idir,ipert)=zero
2749            d2cart(2,idir1,ipert1,idir,ipert)=zero
2750          end do
2751        end do
2752      end if
2753    end do
2754  end do
2755 
2756 end subroutine wings3