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-2022 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!

SOURCE

 19 #if defined HAVE_CONFIG_H
 20 #include "config.h"
 21 #endif
 22 
 23 #include "abi_common.h"
 24 
 25 module m_dynmat
 26 
 27  use, intrinsic :: iso_c_binding
 28  use defs_basis
 29  use m_abicore
 30  use m_errors
 31  use m_linalg_interfaces
 32  use m_xmpi
 33 
 34  use m_fstrings,        only : itoa, sjoin
 35  use m_numeric_tools,   only : wrap2_pmhalf, mkherm
 36  use m_symtk,           only : mati3inv, matr3inv, littlegroup_q
 37  use m_cgtools,         only : fxphas_seq
 38  use m_ewald,           only : ewald9
 39  use m_time,            only : timab
 40 
 41  implicit none
 42 
 43  private
 44 
 45  public :: asria_calc           ! Calculate the correction for the Acoustic sum rule on
 46                                 !   the InterAtomic Forces or on the dynamical matrix directly
 47  public :: asria_corr           ! Imposition of the Acoustic sum rule on the InterAtomic Forces
 48                                 !   or on the dynamical matrix directly from the previously calculated d2asr
 49  public :: asrprs               ! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry
 50  public :: cart29               ! Transform a second-derivative matrix from reduced  coordinates to cartesian coordinates, and also
 51                                 !   1) add the ionic part of the effective charges,
 52                                 !   2) normalize the electronic dielectric tensor, and add the vacuum polarisation
 53  public :: cart39               ! Transform a vector from reduced coordinates to cartesian coordinates,
 54                                 !   taking into account the perturbation from which it was derived,
 55                                 !   and also check the existence of the new values.
 56  public :: d2cart_to_red        ! Transform a second-derivative matrix
 57                                 ! from cartesian to reduced coordinate.
 58  public :: chkph3               ! Check the completeness of the dynamical matrix
 59  public :: chneu9               ! Imposition of the charge neutrality sum rule on the Effective charges
 60  public :: d2sym3               ! Build (nearly) all the other matrix elements that can be build using symmetries.
 61  public :: q0dy3_apply          ! Takes care of the inclusion of the ewald q=0 term in the dynamical matrix
 62  public :: q0dy3_calc           ! Calculate the q=0 correction term to the dynamical matrix
 63  ! TODO: 3 routines to symmetrize. Clarify different use cases
 64  public :: symdyma              ! Symmetrize the dynamical matrices
 65  public :: dfpt_sygra           ! Symmetrize derivatives of energy with respect to coordinates,
 66  public :: dfpt_sydy            ! Symmetrize dynamical matrix (eventually diagonal wrt to the atoms)
 67  public :: wings3               ! Suppress the wings of the cartesian 2DTE for which the diagonal element is not known
 68  public :: asrif9               ! Imposes the Acoustic Sum Rule to Interatomic Forces
 69  public :: get_bigbox_and_weights ! Compute
 70  public :: bigbx9               ! Generates a Big Box of R points for the Fourier Transforms the dynamical matrix
 71  public :: make_bigbox          ! Helper functions that faciliates the generation  of a Big Box containing
 72  public :: canat9               ! From reduced to canonical coordinates
 73  public :: canct9               ! Convert from canonical coordinates to cartesian coordinates
 74  public :: chkrp9               ! Check if the rprim used for the definition of the unit cell (in the
 75                                 ! inputs) are consistent with the rprim used in the routine generating  the Big Box
 76  public :: dist9                ! Compute the distance between atoms in the big box
 77  public :: ftifc_q2r            ! Fourier transform of the dynamical matrices to obtain interatomic forces (real space).
 78  private :: ftifc_r2q            ! Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space).
 79  public :: dynmat_dq            ! Compute the derivative D(q)/dq via Fourier transform of the interatomic forces
 80  public :: ifclo9               ! Convert from cartesian coordinates to local coordinates
 81  public :: wght9                ! Generates a weight to each R points of the Big Box and for each pair of atoms
 82  public :: d3sym                ! Given a set of calculated elements of the 3DTE matrix,
 83                                 ! build (nearly) all the other matrix elements that can be build using symmetries.
 84  public :: d3lwsym                                
 85  public :: sylwtens             ! Determines the set of irreductible elements of the spatial-dispersion tensors
 86  public :: sytens               ! Determines the set of irreductible elements of the nonlinear optical susceptibility
 87                                 ! and Raman tensors
 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 :: pheigvec_normalize   ! Normalize input eigenvectors in cartesian coordinates.
 97  public :: phdispl_from_eigvec  ! Phonon displacements from eigenvectors
 98  public :: dfpt_prtph           ! Print phonon frequencies
 99  public :: massmult_and_breaksym  ! Multiply IFC(q) by atomic masses.
100  public :: massmult_and_breaksym_cplx  ! Version for complex array
101 
102  ! TODO: Change name,
103  public :: ftgam
104  public :: ftgam_init
105 
106 
107 ! *************************************************************************
108 
109 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.

SOURCE

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

SOURCE

294 subroutine asria_corr(asr,d2asr,d2cart,mpert,natom)
295 
296 !Arguments -------------------------------
297 !scalars
298  integer,intent(in) :: asr,mpert,natom
299 !arrays
300  real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
301  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
302 
303 !Local variables-------------------------------
304 !scalars
305  integer :: idir1,idir2,ipert1,ipert2
306 
307 ! *********************************************************************
308 
309  if (asr==0) return
310  !call wrtout(std_out,' asria_corr: imposition of the ASR for the interatomic forces.')
311 
312  ! Remove d2asr
313  do ipert2=1,natom
314    do idir2=1,3
315      do ipert1=1,natom
316        do idir1=1,3
317          d2cart(:,idir1,ipert1,idir2,ipert2)= d2cart(:,idir1,ipert1,idir2,ipert2) - d2asr(:,idir1,ipert1,idir2,ipert2)
318        end do
319      end do
320    end do
321  end do
322 
323 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.

SOURCE

2571 subroutine asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm)
2572 
2573 !Arguments -------------------------------
2574 !scalars
2575  integer,intent(in) :: asr,natom,nrpt
2576 !arrays
2577  real(dp),intent(in) :: rpt(3,nrpt),wghatm(natom,natom,nrpt)
2578  real(dp),intent(inout) :: atmfrc(3,natom,3,natom,nrpt)
2579 
2580 !Local variables -------------------------
2581 !scalars
2582  integer :: found,ia,ib,irpt,izero,mu,nu
2583  real(dp) :: sumifc
2584 
2585 ! *********************************************************************
2586 
2587  if(asr==1.or.asr==2)then
2588    found=0
2589    ! Search for the R vector which is equal to ( 0 , 0 , 0 )
2590    ! This vector leaves the atom a on itself !
2591    do irpt=1,nrpt
2592      if (all(abs(rpt(:,irpt))<=1.0d-10)) then
2593        found=1
2594        izero=irpt
2595      end if
2596      if (found==1) exit
2597    end do
2598 
2599    if(found==0)then
2600      ABI_BUG('Not able to find the vector R=(0,0,0).')
2601    end if
2602 
2603    do mu=1,3
2604      do nu=1,3
2605        do ia=1,natom
2606          sumifc=zero
2607          do ib=1,natom
2608 
2609            ! Get the sumifc of interatomic forces acting on the atom ia,
2610            ! either in a symmetrical manner, or an unsymmetrical one.
2611            if(asr==1)then
2612              do irpt=1,nrpt
2613                sumifc=sumifc+wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)
2614              end do
2615            else if(asr==2)then
2616              do irpt=1,nrpt
2617                sumifc=sumifc+&
2618                 (wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)+&
2619                  wghatm(ia,ib,irpt)*atmfrc(nu,ia,mu,ib,irpt))/2
2620              end do
2621            end if
2622          end do
2623 
2624          ! Correct the self-interaction in order to fulfill the ASR
2625          atmfrc(mu,ia,nu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)-sumifc
2626          if (asr==2) atmfrc(nu,ia,mu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)
2627        end do
2628      end do
2629    end do
2630  end if
2631 
2632 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

SOURCE

360 subroutine asrprs(asr,asrflag,rotinv,uinvers,vtinvers,singular,d2cart,mpert,natom,xcart)
361 
362 !Arguments ------------------------------------
363 !scalars
364  integer,intent(in) :: asr,asrflag,mpert,natom,rotinv
365 !arrays
366  real(dp),intent(in) :: xcart(3,natom)
367  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
368  real(dp),intent(inout) :: singular(1:3*natom*(3*natom-1)/2)
369  real(dp),intent(inout) :: uinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)
370  real(dp),intent(inout) :: vtinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2)
371 
372 !Local variables-------------------------------
373 !scalars
374  integer :: column,idir1,idir2,ii,info,ipert1,ipert2,jj,n3,row,superdim
375  real(dp) :: rcond,test
376 ! real(dp) :: tau ! tau is present but commented out in this routine
377  character(len=500) :: msg
378 !arrays
379  integer :: check(3,natom,3)
380  real(dp) :: tmp(natom,3,3),weightf(1:natom,1:natom)
381  real(dp),allocatable :: d2cartold(:,:,:,:,:),d2vecc(:),d2veccnew(:),d2vecr(:)
382  real(dp),allocatable :: d2vecrnew(:),superm(:,:),umatrix(:,:),vtmatrix(:)
383  real(dp),allocatable :: work(:)
384 
385 ! *********************************************************************
386 
387  if(asr/=3 .and. asr/=4)then
388    write(msg,'(3a,i0)')&
389    'The argument asr should be 3 or 4,',ch10, 'however, asr = ',asr
390    ABI_BUG(msg)
391  end if
392 
393  if (asr==3.or.asr==4)then
394    write(msg, '(a,a)' ) ch10, &
395    'asrprs: imposition of the ASR for the interatomic forces and rotational invariance'
396    call wrtout(std_out,msg)
397  end if
398 
399  write(msg,'(a,i0)')' asrflag is ', asrflag
400  call wrtout([std_out, ab_out], msg)
401 
402 !variables for the dimensions of the matrices
403 
404 !n1=3*natom*(3*natom-1)/2
405 !n2=9*natom
406  n3=3*natom
407 
408  superdim=9*natom*(natom-1)/2+n3
409 
410  ABI_MALLOC(d2vecr,(1:superdim))
411  ABI_MALLOC(d2vecc,(1:superdim))
412  d2vecr=0d0
413  d2vecc=0d0
414 
415 !should be changed set to delta function for debugging
416  weightf=1d0
417 !tau=1d-10
418  do ii=1, natom
419 !  do jj=1, ii-1
420 !  weightf(ii,jj)= &
421 !  &     ((xcart(1,ii)-xcart(1,jj))**2+(xcart(2,ii)-xcart(2,jj))**2+(xcart(3,ii)-xcart(3,jj))**2)**tau
422 !  enddo
423    weightf(ii,ii)=0d0
424  end do
425 
426  ABI_MALLOC(d2cartold,(2,3,mpert,3,mpert))
427 
428  d2cartold=d2cart
429 
430 !setup vector with uncorrected derivatives
431 
432  do ipert1=1, natom
433    do ipert2=1, ipert1-1
434      do idir1=1,3
435        do idir2=1,3
436          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
437          if(abs(d2cart(1,idir1,ipert1,idir2,ipert2))<1d-6)then
438            d2cart(1,idir1,ipert1,idir2,ipert2)=0d0
439          else
440            d2vecr(row)=4*weightf(ipert1,ipert2)*d2cart(1,idir1,ipert1,idir2,ipert2)
441          end if
442          if(abs(d2cart(2,idir1,ipert1,idir2,ipert2))<1d-6) then
443            d2cart(2,idir1,ipert1,idir2,ipert2)=0d0
444          else
445            d2vecc(row)=4*weightf(ipert1,ipert2)*d2cart(2,idir1,ipert1,idir2,ipert2)
446          end if
447        end do
448      end do
449    end do
450  end do
451 
452  if(asrflag==1) then !calculate the pseudo-inverse of the supermatrix
453    ABI_MALLOC(superm,(1:superdim,1:superdim))
454 
455    superm=0d0
456 
457 !  Setting up the supermatrix containing G, A, D
458 
459    do ipert1=1, natom
460      do idir1=1, 3
461 !      Setting up G
462        idir2=mod(idir1,3)+1
463        row=3*(ipert1-1)+idir1
464        do ipert2=1, ipert1-1
465          column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
466          superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
467          column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir2
468          superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
469        end do
470        do ipert2=ipert1+1, natom
471          column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir1-1)+rotinv
472          superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1)
473          column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir2-1)+rotinv
474          superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2)
475        end do
476      end do
477      do idir1=1, 3
478 !      Setting up D
479        idir2=mod(idir1,3)+1
480        ii=mod(idir1+1,3)+1
481        do ipert2=1, ipert1-1
482          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1
483          column=9*natom*(natom-1)/2+3*(ipert1-1)+idir1
484          superm(column,row)=superm(column,row)+xcart(idir2,ipert2)-xcart(idir2,ipert1)
485          column=9*natom*(natom-1)/2+3*(ipert1-1)+ii
486          superm(column,row)=superm(column,row)+xcart(ii,ipert1)-xcart(ii,ipert2)
487          row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+rotinv
488          column=9*natom*(natom-1)/2+3*(ipert2-1)+idir1
489          superm(column,row)=superm(column,row)+xcart(idir2,ipert1)-xcart(idir2,ipert2)
490          column=9*natom*(natom-1)/2+3*(ipert2-1)+ii
491          superm(column,row)=superm(column,row)+xcart(ii,ipert2)-xcart(ii,ipert1)
492        end do
493 !      Setting up A
494        do idir2=1, 3
495          do ipert2=1, ipert1-1
496            column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
497            row=n3+column
498            superm(column,row)=4*weightf(ipert1,ipert2)
499          end do
500        end do
501      end do
502    end do
503 
504 !  calculate the pseudo-inverse of the supermatrix
505 
506    ABI_MALLOC(work,(1:6*superdim))
507    ABI_MALLOC(vtmatrix,(1:superdim))
508    ABI_MALLOC(umatrix,(1:superdim,1:superdim))
509 
510 !  singular value decomposition of superm
511 
512    call dgesvd('A','O',superdim,superdim,superm,superdim,singular,umatrix,superdim, &
513                vtmatrix, 1, work,6*superdim,info)
514    ABI_CHECK(info == 0, sjoin('dgesvd returned:', itoa(info)))
515 
516    ABI_FREE(vtmatrix)
517    ABI_FREE(work)
518 
519    write(msg, '(a,es16.8,es16.8)' )' Largest and smallest values from svd', singular(1), singular(superdim)
520    call wrtout([std_out, ab_out], msg)
521 
522 !  Invert U and V**T, orthogonal matrices
523 
524    do ii=1, superdim
525      do jj=1, superdim
526        uinvers(ii,jj)=umatrix(jj,ii)
527        vtinvers(ii,jj)=superm(jj,ii)
528      end do
529    end do
530 
531    ABI_FREE(umatrix)
532    ABI_FREE(superm)
533 
534    write(msg,'(a,a)')' asrprs: done with asrflag 1', ch10
535    call wrtout(std_out,msg)
536 
537  end if !asrflag=1
538 
539  if(asrflag==2) then
540 
541    ABI_MALLOC(d2vecrnew,(1:superdim))
542    ABI_MALLOC(d2veccnew,(1:superdim))
543 
544 !  Calculate V**T**-1 Sigma**-1 U**-1 *rhs
545 
546    do ii=1, superdim
547      d2vecrnew(ii)=0d0
548      d2veccnew(ii)=0d0
549      do jj=1, superdim
550        d2vecrnew(ii)=d2vecrnew(ii)+uinvers(ii,jj)*d2vecr(jj)
551        d2veccnew(ii)=d2veccnew(ii)+uinvers(ii,jj)*d2vecc(jj)
552      end do
553    end do
554 
555    rcond=1d-10*singular(1)
556    do ii=1, superdim
557      if(singular(ii)>rcond) then
558        d2vecrnew(ii)=d2vecrnew(ii)/singular(ii)
559        d2veccnew(ii)=d2veccnew(ii)/singular(ii)
560      else
561        d2vecrnew(ii)=0d0
562        d2veccnew(ii)=0d0
563      end if
564    end do
565 
566    do ii=1, superdim
567      d2vecr(ii)=0d0
568      d2vecc(ii)=0d0
569      do jj=1, superdim
570        d2vecr(ii)=d2vecr(ii)+vtinvers(ii,jj)*d2vecrnew(jj)
571        d2vecc(ii)=d2vecc(ii)+vtinvers(ii,jj)*d2veccnew(jj)
572      end do
573    end do
574 
575 !  Store vector back into the matrix of 2nd order derivates
576 
577    do ipert1=1, natom
578      do ipert2=1, ipert1-1
579        do idir1=1,3
580          do idir2=1,3
581            row=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2
582            d2cart(1,idir1,ipert1,idir2,ipert2)=d2vecr(row)
583            d2cart(2,idir1,ipert1,idir2,ipert2)=d2vecc(row)
584            d2cart(1,idir2,ipert2,idir1,ipert1)=d2vecr(row)
585            d2cart(2,idir2,ipert2,idir1,ipert1)=d2vecc(row)
586          end do
587        end do
588      end do
589    end do
590 
591    ABI_FREE(d2vecrnew)
592    ABI_FREE(d2veccnew)
593 
594    check=0
595 
596    do ipert1=1, natom
597      do idir1=1, 3
598        do idir2=1, 3
599          d2cart(1,idir1,ipert1,idir2,ipert1)=0d0
600          d2cart(2,idir1,ipert1,idir2,ipert1)=0d0
601          tmp(ipert1,idir1,idir2)=0d0
602          do ipert2=1, natom
603            if(ipert2/=ipert1) then
604              tmp(ipert1,idir1,idir2)=tmp(ipert1,idir1,idir2) &
605 &             -d2cart(1,idir1,ipert1,idir2,ipert2) &
606 &             -d2cart(1,idir2,ipert2,idir1,ipert1)
607            end if
608          end do
609        end do
610      end do
611    end do
612 
613    do ipert1=1, natom
614      do idir1=1, 3
615        do idir2=1, 3
616          d2cart(1,idir1,ipert1,idir2,ipert1)=tmp(ipert1,idir1,idir2)/2
617          d2cart(1,idir2,ipert1,idir1,ipert1)=d2cart(1,idir1,ipert1,idir2,ipert1)
618        end do
619      end do
620    end do
621 
622    write(std_out,*) 'this should all be zero'
623 
624    do ipert1=1, natom
625      do idir1=1, 3
626        do idir2=1, 3
627          test=0d0
628          do ipert2=1, natom
629            test=test+d2cart(1,idir1,ipert1,idir2,ipert2)+d2cart(1,idir2,ipert2,idir1,ipert1)
630          end do
631          write(std_out,'(i3,i3,i3,es11.3)') idir1,ipert1,idir2,test
632 
633          write(msg, '(i3,i3,i3,es11.3)' ) idir1,ipert1,idir2,test
634          call wrtout(ab_out,msg)
635        end do
636      end do
637    end do
638 
639    write(std_out,*) 'these as well'
640    do ipert2=1, natom
641      do idir1=1, 3
642        do idir2=1, 3
643          test=0d0
644          do ipert1=1, natom
645            test=test+d2cart(1,idir1,ipert1,idir2,ipert2)
646          end do
647          write(std_out,'(i3,i3,i3,i3,es11.3)') idir1,ipert1,idir2,ipert2,test
648        end do
649      end do
650    end do
651 
652    write(msg,'(a,a)')' asrprs: done with asrflag 2', ch10
653    call wrtout(std_out,msg)
654 
655  end if !ends asrflag=2
656 
657  ABI_FREE(d2vecr)
658  ABI_FREE(d2vecc)
659  ABI_FREE(d2cartold)
660 
661 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

SOURCE

5244 subroutine axial9(ifccar,vect1,vect2,vect3)
5245 
5246 !Arguments -------------------------------
5247 !arrays
5248  real(dp),intent(in) :: ifccar(3,3),vect1(3)
5249  real(dp),intent(out) :: vect2(3),vect3(3)
5250 
5251 !Local variables -------------------------
5252 !scalars
5253  integer :: flag,ii,itrial,jj
5254  real(dp) :: innorm,scprod
5255 !arrays
5256  real(dp) :: work(3)
5257 
5258 ! *********************************************************************
5259 
5260  do jj=1,3
5261    work(jj)=zero
5262    do ii=1,3
5263      work(jj)=work(jj)+ifccar(jj,ii)*vect1(ii)
5264    end do
5265  end do
5266 
5267  flag=0
5268  do itrial=1,4
5269    scprod=zero
5270    do ii=1,3
5271      scprod=scprod+work(ii)*vect1(ii)
5272    end do
5273 
5274    do ii=1,3
5275      work(ii)=work(ii)-vect1(ii)*scprod
5276    end do
5277 
5278    scprod=zero
5279    do ii=1,3
5280      scprod=scprod+work(ii)**2
5281    end do
5282 
5283    if(scprod<1.0d-10)then
5284      work(1:3)=zero
5285      if(itrial>1)work(itrial-1)=1.0_dp
5286    else
5287      flag=1
5288    end if
5289 
5290    if(flag==1)exit
5291  end do
5292 
5293  innorm=scprod**(-0.5_dp)
5294  do ii=1,3
5295    vect2(ii)=work(ii)*innorm
5296  end do
5297 
5298  vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2)
5299  vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3)
5300  vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1)
5301 
5302 end subroutine axial9

m_dynmat/bigbx9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 bigbx9

FUNCTION

 Generation of a Big Box containing all the R points (cells) 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(3,nrpt)= integer coordinates of the cells (R points) in the rprim basis
 nprt= Number of cells (R points) in the Big Box
 rpt(3,mrpt)= canonical coordinates of the cells (R points)
  These coordinates are normalized (=> * acell(3)!!)
  (output only if choice=1)

SOURCE

2885 subroutine bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)
2886 
2887 !Arguments -------------------------------
2888 !scalars
2889  integer,intent(in) :: brav,choice,mrpt,nqshft
2890  integer,intent(out) :: nrpt
2891 !arrays
2892  integer,intent(in) :: ngqpt(3)
2893  real(dp),intent(in) :: rprim(3,3)
2894  real(dp),intent(out) :: rpt(3,mrpt)
2895  integer,intent(out) :: cell(3,mrpt)
2896 
2897 !Local variables -------------------------
2898 !In some cases, the atoms coordinates are not packed in the
2899 ! [0,1]^3 cube. Then, the parameter "buffer" might be increased,
2900 !to search relevant pairs of atoms in bigger boxes than usual.
2901 !scalars
2902  integer,parameter :: buffer=1
2903  integer :: irpt,lim1,lim2,lim3,lqshft,r1,r2,r3
2904  character(len=500) :: msg
2905 
2906 ! *********************************************************************
2907 
2908  lqshft=1
2909  if(nqshft/=1)lqshft=2
2910 
2911 
2912 !Simple Cubic Lattice
2913  if (abs(brav)==1) then
2914    lim1=((ngqpt(1))+1)*lqshft+buffer
2915    lim2=((ngqpt(2))+1)*lqshft+buffer
2916    lim3=((ngqpt(3))+1)*lqshft+buffer
2917    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
2918    if(choice/=0)then
2919      if (nrpt/=mrpt) then
2920        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
2921        ABI_BUG(msg)
2922      end if
2923      irpt=0
2924      do r1=-lim1,lim1
2925        do r2=-lim2,lim2
2926          do r3=-lim3,lim3
2927            irpt=irpt+1
2928            rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
2929            rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
2930            rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
2931            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
2932          end do
2933        end do
2934      end do
2935    end if
2936 
2937 !  Face Centered Cubic Lattice
2938  else if (brav==2) then
2939    lim1=((ngqpt(1)+3)/4)*lqshft+buffer
2940    lim2=((ngqpt(2)+3)/4)*lqshft+buffer
2941    lim3=((ngqpt(3)+3)/4)*lqshft+buffer
2942    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*4
2943    if(choice/=0)then
2944      if (nrpt/=mrpt) then
2945        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt
2946        ABI_BUG(msg)
2947      end if
2948      irpt=0
2949      do r1=-lim1,lim1
2950        do r2=-lim2,lim2
2951          do r3=-lim3,lim3
2952            irpt=irpt+4
2953            rpt(1,irpt-3)=r1
2954            rpt(2,irpt-3)=r2
2955            rpt(3,irpt-3)=r3
2956            rpt(1,irpt-2)=r1
2957            rpt(2,irpt-2)=r2+0.5
2958            rpt(3,irpt-2)=r3+0.5
2959            rpt(1,irpt-1)=r1+0.5
2960            rpt(2,irpt-1)=r2
2961            rpt(3,irpt-1)=r3+0.5
2962            rpt(1,irpt)=r1+0.5
2963            rpt(2,irpt)=r2+0.5
2964            rpt(3,irpt)=r3
2965 !TEST_AM
2966 !           cell(irpt-3,1)=r1;cell(irpt-3,2)=r2;cell(irpt-3,3)=r3
2967            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
2968          end do
2969        end do
2970      end do
2971    end if
2972 
2973 !  Body Centered Cubic Lattice
2974  else if (brav==3) then
2975    lim1=((ngqpt(1)+3)/4)*lqshft+buffer
2976    lim2=((ngqpt(2)+3)/4)*lqshft+buffer
2977    lim3=((ngqpt(3)+3)/4)*lqshft+buffer
2978    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*2
2979    if(choice/=0)then
2980      if(nrpt/=mrpt) then
2981        write(msg,'(2(a,i0))')' nrpt= ',nrpt,' is not equal to mrpt= ',mrpt
2982        ABI_BUG(msg)
2983      end if
2984      irpt=0
2985      do r1=-lim1,lim1
2986        do r2=-lim2,lim2
2987          do r3=-lim3,lim3
2988            irpt=irpt+2
2989            rpt(1,irpt-1)=r1
2990            rpt(2,irpt-1)=r2
2991            rpt(3,irpt-1)=r3
2992            rpt(1,irpt)=r1+0.5
2993            rpt(2,irpt)=r2+0.5
2994            rpt(3,irpt)=r3+0.5
2995 !TEST_AM
2996 !           cell(irpt-1,1)=r1;cell(irpt-1,2)=r2;cell(irpt-1,3)=r3
2997            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
2998          end do
2999        end do
3000      end do
3001    end if
3002 
3003 !  Hexagonal Lattice
3004  else if (brav==4) then
3005    lim1=(ngqpt(1)+1)*lqshft+buffer
3006    lim2=(ngqpt(2)+1)*lqshft+buffer
3007    lim3=((ngqpt(3)/2)+1)*lqshft+buffer
3008    nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
3009    if(choice/=0)then
3010      if(nrpt/=mrpt)then
3011        write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt=',mrpt
3012        ABI_BUG(msg)
3013      end if
3014      irpt=0
3015      do r1=-lim1,lim1
3016        do r2=-lim2,lim2
3017          do r3=-lim3,lim3
3018            irpt=irpt+1
3019            rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3)
3020            rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3)
3021            rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3)
3022            cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
3023          end do
3024        end do
3025      end do
3026    end if
3027 
3028  else
3029    write(msg,'(a,i0,a)')' The value of brav= ',brav,' is not allowed (should be -1, 1, 2 or 4).'
3030    ABI_BUG(msg)
3031  end if
3032 
3033 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

SOURCE

3059 subroutine canat9(brav,natom,rcan,rprim,trans,xred)
3060 
3061 !Arguments -------------------------------
3062 !scalars
3063  integer,intent(in) :: brav,natom
3064 !arrays
3065  real(dp),intent(in) :: rprim(3,3),xred(3,natom)
3066  real(dp),intent(out) :: rcan(3,natom),trans(3,natom)
3067 
3068 !Local variables -------------------------
3069 !scalars
3070  integer :: found,iatom,ii
3071  character(len=500) :: msg
3072 !arrays
3073  real(dp) :: dontno(3,4),rec(3),rok(3),shift(3),tt(3)
3074 
3075 ! *********************************************************************
3076 
3077 !Normalization of the cartesian atomic coordinates
3078 !If not normalized : rcan(i) <- rcan(i) * acell(i)
3079  do iatom=1,natom
3080    rcan(:,iatom)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3)
3081  end do
3082 
3083  !Study of the different cases for the Bravais lattice:
3084  if (abs(brav)==1) then
3085    !Simple Cubic Lattice
3086 
3087    do iatom=1,natom
3088      ! Canon will produces these coordinate transformations
3089      ! (Note: here we still use reduced coordinates )
3090      call wrap2_pmhalf(xred(1,iatom),rok(1),shift(1))
3091      call wrap2_pmhalf(xred(2,iatom),rok(2),shift(2))
3092      call wrap2_pmhalf(xred(3,iatom),rok(3),shift(3))
3093 
3094 !    New coordinates : rcan
3095      rcan(:,iatom)=rok(1)*rprim(:,1)+rok(2)*rprim(:,2)+rok(3)*rprim(:,3)
3096 !    Translations between New and Old coordinates
3097      tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3)
3098      trans(:,iatom)=tt(:)-rcan(:,iatom)
3099    end do
3100 
3101  else if (brav==2) then
3102    ! Face Centered Lattice
3103    ! Special possible translations in the F.C.C. case
3104    dontno(:,:)=zero
3105    dontno(2,2)=0.5_dp
3106    dontno(3,2)=0.5_dp
3107    dontno(1,3)=0.5_dp
3108    dontno(3,3)=0.5_dp
3109    dontno(1,4)=0.5_dp
3110    dontno(2,4)=0.5_dp
3111    do iatom=1,natom
3112      found=0
3113      do ii=1,4
3114        if (found==1) exit
3115        ! Canon will produce these coordinate transformations
3116        call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
3117        call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
3118        call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
3119        ! In the F.C.C., ABS[ Ri ] + ABS[ Rj ] < or = 1/2
3120        ! The equal sign hase been treated using a tolerance parameter
3121        ! not to have twice the same point in the unit cell !
3122        rok(1)=rok(1)-1.0d-10
3123        rok(2)=rok(2)-2.0d-10
3124        rok(3)=rok(3)-5.0d-10
3125        if (abs(rok(1))+abs(rok(2))<=0.5_dp) then
3126          if (abs(rok(1))+abs(rok(3))<=0.5_dp) then
3127            if (abs(rok(2))+abs(rok(3))<=0.5_dp) then
3128              tt(:)=rcan(:,iatom)
3129              ! New coordinates : rcan
3130              rcan(1,iatom)=rok(1)+1.0d-10
3131              rcan(2,iatom)=rok(2)+2.0d-10
3132              rcan(3,iatom)=rok(3)+5.0d-10
3133              ! Translations between New and Old coordinates
3134              trans(:,iatom)=tt(:)-rcan(:,iatom)
3135              found=1
3136            end if
3137          end if
3138        end if
3139      end do
3140    end do
3141 
3142  else if (brav==3) then
3143    ! Body Centered Cubic Lattice
3144    ! Special possible translations in the B.C.C. case
3145    dontno(:,1)=zero
3146    dontno(:,2)=0.5_dp
3147    do iatom=1,natom
3148      found=0
3149      do ii=1,2
3150        if (found==1) exit
3151        ! Canon will produce these coordinate transformations
3152        call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1))
3153        call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2))
3154        call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3))
3155        ! In the F.C.C., ABS[ Ri ] < or = 1/2
3156        ! and    ABS[ R1 ] + ABS[ R2 ] + ABS[ R3 ] < or = 3/4
3157        ! The equal signs have been treated using a tolerance parameter
3158        ! not to have twice the same point in the unit cell !
3159        rok(1)=rok(1)-1.0d-10
3160        rok(2)=rok(2)-2.0d-10
3161        rok(3)=rok(3)-5.0d-10
3162        if(abs(rok(1))+abs(rok(2))+abs(rok(3))<=0.75_dp)then
3163          if ( abs(rok(1))<=0.5_dp .and. abs(rok(2))<=0.5_dp .and. abs(rok(3))<=0.5_dp) then
3164            tt(:)=rcan(:,iatom)
3165            ! New coordinates : rcan
3166            rcan(1,iatom)=rok(1)+1.0d-10
3167            rcan(2,iatom)=rok(2)+2.0d-10
3168            rcan(3,iatom)=rok(3)+5.0d-10
3169            ! Translations between New and Old coordinates
3170            trans(:,iatom)=tt(:)-rcan(:,iatom)
3171            found=1
3172          end if
3173        end if
3174      end do
3175    end do
3176 
3177  else if (brav==4) then
3178    ! Hexagonal Lattice
3179    ! In this case, it is easier first to work in reduced coordinates space !
3180    do iatom=1,natom
3181      ! Passage from the reduced space to the "lozenge" cell
3182      rec(1)=xred(1,iatom)-0.5_dp
3183      rec(2)=xred(2,iatom)-0.5_dp
3184      rec(3)=xred(3,iatom)
3185      ! Canon will produces these coordinate transformations
3186      call wrap2_pmhalf(rec(1),rok(1),shift(1))
3187      call wrap2_pmhalf(rec(2),rok(2),shift(2))
3188      call wrap2_pmhalf(rec(3),rok(3),shift(3))
3189      rec(1)=rok(1)+0.5_dp
3190      rec(2)=rok(2)+0.5_dp
3191      rec(3)=rok(3)
3192      ! Passage in Cartesian Normalized Coordinates
3193      rcan(:,iatom)=rec(1)*rprim(:,1)+rec(2)*rprim(:,2)+rec(3)*rprim(:,3)
3194      ! Use of a tolerance parameter not to have twice the same point in the unit cell !
3195      rcan(1,iatom)=rcan(1,iatom)-1.0d-10
3196      rcan(2,iatom)=rcan(2,iatom)-2.0d-10
3197      ! Passage to the honeycomb hexagonal unit cell !
3198      if (rcan(1,iatom)>0.5_dp) then
3199        rcan(1,iatom)=rcan(1,iatom)-1.0_dp
3200      end if
3201      if (rcan(1,iatom)>zero.and.rcan(1,iatom)+sqrt(3.0_dp)*rcan(2,iatom)>1.0_dp) then
3202        rcan(1,iatom)=rcan(1,iatom)-0.5_dp
3203        rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
3204      end if
3205      if (rcan(1,iatom)<=zero.and.sqrt(3.0_dp)*rcan(2,iatom)-rcan(1,iatom)>1.0_dp) then
3206        rcan(1,iatom)=rcan(1,iatom)+0.5_dp
3207        rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp
3208      end if
3209      ! Translations between New and Old coordinates
3210      tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3)
3211      trans(:,iatom)=tt(:)-rcan(:,iatom)
3212    end do
3213 
3214    ! End of the possible cases for brav : -1, 1, 2, 4.
3215  else
3216    write(msg, '(a,i0,a,a,a)' )&
3217    'The required value of brav=',brav,' is not available.',ch10,&
3218    'It should be -1, 1,2 or 4 .'
3219    ABI_BUG(msg)
3220  end if
3221 
3222  call wrtout(std_out,' Canonical Atomic Coordinates ')
3223  do iatom=1,natom
3224    write(msg, '(a,i5,3es18.8)' )' atom',iatom,rcan(1,iatom),rcan(2,iatom),rcan(3,iatom)
3225    call wrtout(std_out,msg)
3226  end do
3227 
3228 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.

SOURCE

3259 subroutine canct9(acell,gprim,ib,index,irpt,natom,nrpt,rcan,rcart,rprim,rpt)
3260 
3261 !Arguments -------------------------------
3262 !scalars
3263  integer,intent(in) :: index,natom,nrpt
3264  integer,intent(out) :: ib,irpt
3265 !arrays
3266  real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3)
3267  real(dp),intent(in) :: rpt(3,nrpt)
3268  real(dp),intent(out) :: rcart(3)
3269 
3270 !Local variables -------------------------
3271 !scalars
3272  integer :: jj
3273 !arrays
3274  real(dp) :: xred(3)
3275 
3276 ! *********************************************************************
3277 
3278  irpt=(index-1)/natom+1
3279  ib=index-natom*(irpt-1)
3280 
3281 !Transform the canonical coordinates to reduced coord.
3282  do jj=1,3
3283    xred(jj)=gprim(1,jj)*(rpt(1,irpt)+rcan(1,ib))&
3284 &   +gprim(2,jj)*(rpt(2,irpt)+rcan(2,ib))&
3285 &   +gprim(3,jj)*(rpt(3,irpt)+rcan(3,ib))
3286  end do
3287 
3288 !Then to cartesian coordinates (here the position of the atom b)
3289  do jj=1,3
3290    rcart(jj)=xred(1)*acell(1)*rprim(jj,1)+&
3291 &   xred(2)*acell(2)*rprim(jj,2)+&
3292 &   xred(3)*acell(3)*rprim(jj,3)
3293  end do
3294 
3295 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

SOURCE

703 subroutine cart29(blkflg,blkval,carflg,d2cart,&
704 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
705 
706 !Arguments -------------------------------
707 !scalars
708  integer,intent(in) :: iblok,mpert,natom,nblok,ntypat
709  real(dp),intent(in) :: ucvol
710 !arrays
711  integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),typat(natom)
712  integer,intent(out) :: carflg(3,mpert,3,mpert)
713  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3)
714  real(dp),intent(in) :: zion(ntypat)
715  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
716 
717 !Local variables -------------------------
718 !scalars
719  integer :: idir1,idir2,ii,ipert1,ipert2
720 !arrays
721  integer :: flg1(3),flg2(3)
722  real(dp) :: vec1(3),vec2(3)
723 
724 ! *********************************************************************
725 
726 !First, copy the data blok in place.
727  d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok)
728 
729 !Cartesian coordinates transformation (in two steps)
730 !First step
731  do ipert1=1,mpert
732    do ipert2=1,mpert
733      do ii=1,2
734        do idir1=1,3
735          do idir2=1,3
736            vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2)
737 !          Note here blkflg
738            flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok)
739          end do
740          call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2)
741          do idir2=1,3
742            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2)
743 !          And here carflg
744            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2)
745          end do
746        end do
747      end do
748    end do
749  end do
750 
751 !Second step
752  do ipert1=1,mpert
753    do ipert2=1,mpert
754      do ii=1,2
755        do idir2=1,3
756          do idir1=1,3
757            vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2)
758 !          Note here carflg
759            flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2)
760          end do
761          call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2)
762          do idir1=1,3
763            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1)
764 !          And here carflg again
765            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1)
766          end do
767        end do
768      end do
769    end do
770  end do
771 
772 !For the dielectric tensor, takes into account the volume
773 !of the unit cell, and add the unit matrix (polarization of the vacuum)
774  do idir1=1,3
775    do idir2=1,3
776      do ii=1,2
777        d2cart(ii,idir1,natom+2,idir2,natom+2)=&
778 &       -four_pi/ucvol*d2cart(ii,idir1,natom+2,idir2,natom+2)
779      end do
780    end do
781  end do
782 
783  do idir1=1,3
784    d2cart(1,idir1,natom+2,idir1,natom+2)=&
785 &   1.0_dp+d2cart(1,idir1,natom+2,idir1,natom+2)
786  end do
787 
788 !Add the ionic charges to delta z to get the effective charges
789  do ipert1=1,natom
790    do idir1=1,3
791      d2cart(1,idir1,ipert1,idir1,natom+2)=&
792 &     zion(typat(ipert1))+d2cart(1,idir1,ipert1,idir1,natom+2)
793    end do
794  end do
795  do ipert2=1,natom
796    do idir2=1,3
797      d2cart(1,idir2,natom+2,idir2,ipert2)=&
798 &     zion(typat(ipert2))+d2cart(1,idir2,natom+2,idir2,ipert2)
799    end do
800  end do
801 
802 !For the piezoelectric tensor, takes into account the volume of the unit cell
803  do ipert2=natom+3,natom+4
804    do idir1=1,3
805      do idir2=1,3
806        do ii=1,2
807          d2cart(ii,idir1,natom+2,idir2,ipert2)=&
808 &         (1.0_dp/ucvol)*d2cart(ii,idir1,natom+2,idir2,ipert2)
809          d2cart(ii,idir2,ipert2,idir1,natom+2)=&
810 &         (1.0_dp/ucvol)*d2cart(ii,idir2,ipert2,idir1,natom+2)
811        end do
812      end do
813    end do
814  end do
815 
816 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

SOURCE

845 subroutine cart39(flg1,flg2,gprimd,ipert,natom,rprimd,vec1,vec2)
846 
847 !Arguments -------------------------------
848 !scalars
849  integer,intent(in) :: ipert,natom
850 !arrays
851  integer,intent(in) :: flg1(3)
852  integer,intent(out) :: flg2(3)
853  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3),vec1(3)
854  real(dp),intent(out) :: vec2(3)
855 
856 !Local variables -------------------------
857 !scalars
858  integer :: idir,ii
859 
860 ! *********************************************************************
861 
862 !Treat phonon-type perturbation
863  if(ipert>=1.and.ipert<=natom)then
864 
865    do idir=1,3
866      vec2(idir)=zero
867      flg2(idir)=1
868      do ii=1,3
869        if(abs(gprimd(idir,ii))>1.0d-10)then
870          if(flg1(ii)==1)then
871            vec2(idir)=vec2(idir)+gprimd(idir,ii)*vec1(ii)
872          else
873            flg2(idir)=0
874          end if
875        end if
876      end do
877      if(flg2(idir)==0)vec2(idir)=zero
878    end do
879 
880 !  Treat electric field and qvec perturbations
881  else if(ipert==natom+2.or.ipert==natom+8) then
882 !  OCL SCALAR
883    do idir=1,3
884      vec2(idir)=zero
885      flg2(idir)=1
886 !    OCL SCALAR
887      do ii=1,3
888        if(abs(rprimd(idir,ii))>1.0d-10)then
889          if(flg1(ii)==1)then
890            vec2(idir)=vec2(idir)+rprimd(idir,ii)*vec1(ii)/two_pi
891          else
892            flg2(idir)=0
893          end if
894        end if
895      end do
896      if(flg2(idir)==0)vec2(idir)=zero
897    end do
898 
899 !  Treat other perturbations
900  else
901    do idir=1,3
902      vec2(idir)=vec1(idir)
903      flg2(idir)=flg1(idir)
904    end do
905  end if
906 
907 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

SOURCE

1081 subroutine chkph3(carflg,idir,mpert,natom)
1082 
1083 !Arguments -------------------------------
1084 !scalars
1085  integer,intent(in) :: idir,mpert,natom
1086 !arrays
1087  integer,intent(in) :: carflg(3,mpert,3,mpert)
1088 
1089 !Local variables -------------------------
1090 !scalars
1091  integer :: idir1,idir2,ipert1,ipert2,send
1092  character(len=500) :: msg
1093 
1094 ! *********************************************************************
1095 
1096  send=0
1097 
1098 !Check the elements of the analytical part of the dynamical matrix
1099  do ipert2=1,natom
1100    do idir2=1,3
1101      do ipert1=1,natom
1102        do idir1=1,3
1103          if(carflg(idir1,ipert1,idir2,ipert2)==0)then
1104            send=1
1105          end if
1106        end do
1107      end do
1108    end do
1109  end do
1110 
1111 !If some electric field is present
1112  if(idir/=0)then
1113 
1114 !  Check the dielectric constant
1115    if(carflg(idir,natom+2,idir,natom+2)==0)then
1116      send=1
1117    end if
1118 
1119 !  Check the effective charges
1120    do ipert1=1,natom
1121      do idir1=1,3
1122        if(carflg(idir1,ipert1,idir,natom+2)==0)then
1123          send=1
1124        end if
1125      end do
1126    end do
1127 
1128  end if
1129 
1130  ! If needed, send the message
1131  if(send==1)then
1132    write(msg, '(a,a,a,a)' )&
1133    ' chkph3 : WARNING -',ch10,&
1134    '  Dynamical matrix incomplete, phonon frequencies may be wrong, check input variables rfatpol and rfdir.'
1135    call wrtout([std_out, ab_out],msg)
1136  end if
1137 
1138 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)

SOURCE

3319 subroutine chkrp9(brav,rprim)
3320 
3321 !Arguments -------------------------------
3322 !scalars
3323  integer,intent(in) :: brav
3324 !arrays
3325  real(dp),intent(in) :: rprim(3,3)
3326 
3327 !Local variables -------------------------
3328 !scalars
3329  integer :: ii,jj
3330  character(len=500) :: msg
3331 
3332 ! *********************************************************************
3333 
3334  if (abs(brav)==1) then
3335 !  Simple Cubic Lattice No condition in this case !
3336    continue
3337 
3338  else if (brav==2) then
3339 !  Face Centered Lattice
3340    do ii=1,3
3341      do jj=1,3
3342        if (  ( ii==jj .and. abs(rprim(ii,jj))>tol10) .or. (ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
3343          write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3344          'The input variable rprim does not correspond to the',ch10,&
3345          'fixed rprim to be used with brav=2 and ifcflag=1 :',ch10,&
3346          '   0  1/2  1/2',ch10,&
3347          '  1/2  0   1/2',ch10,&
3348          '  1/2 1/2   0 ',ch10,&
3349          'Action: rebuild your DDB by using the latter rprim.'
3350          ABI_ERROR(msg)
3351        end if
3352      end do
3353    end do
3354 
3355  else if (brav==3) then
3356 !  Body Centered Cubic Lattice
3357    do ii=1,3
3358      do jj=1,3
3359        if (  ( ii==jj .and. abs(rprim(ii,jj)+.5_dp)>tol10) .or. (ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then
3360          write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3361          'The input variable rprim does not correspond to the',ch10,&
3362          'fixed rprim to be used with brav=3 and ifcflag=1 :',ch10,&
3363          '  -1/2  1/2  1/2',ch10,&
3364          '   1/2 -1/2  1/2',ch10,&
3365          '   1/2  1/2 -1/2',ch10,&
3366          'Action: rebuild your DDB by using the latter rprim.'
3367          ABI_ERROR(msg)
3368        end if
3369      end do
3370    end do
3371 
3372  else if (brav==4) then
3373 !  Hexagonal Lattice
3374    if (abs(rprim(1,1)-1.0_dp)>tol10 .or. &
3375        abs(rprim(3,3)-1.0_dp)>tol10 .or. &
3376        abs(rprim(2,1)      )>tol10 .or. &
3377        abs(rprim(3,1)      )>tol10 .or. &
3378        abs(rprim(1,3)      )>tol10 .or. &
3379        abs(rprim(2,3)      )>tol10 .or. &
3380        abs(rprim(3,2)      )>tol10 .or. &
3381        abs(rprim(1,2)+0.5_dp)>tol10 .or. &
3382        abs(rprim(2,2)-0.5_dp*sqrt(3.0_dp))>tol10 ) then
3383      write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )&
3384       'The input variable rprim does not correspond to the',ch10,&
3385       'fixed rprim to be used with brav=4 and ifcflag=1 :',ch10,&
3386       '   1      0      0',ch10,&
3387       '  -1/2 sqrt[3]/2 0',ch10,&
3388       '   0      0      1',ch10,&
3389       'Action: rebuild your DDB by using the latter rprim.'
3390      ABI_ERROR(msg)
3391    end if
3392 
3393  else
3394    write(msg, '(a,i4,a,a,a,a,a)' )&
3395    'The value of brav=',brav,' is not allowed.',ch10,&
3396    'Only  -1, 1,2,3 or 4 are allowed.',ch10,&
3397    'Action: change the value of brav in your input file.'
3398    ABI_ERROR(msg)
3399  end if
3400 
3401 end subroutine chkrp9

m_dynmat/chneu9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 chneu9

FUNCTION

 Imposition of the charge neutrality 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

SOURCE

1168 subroutine chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion)
1169 
1170 !Arguments -------------------------------
1171 !scalars
1172  integer,intent(in) :: chneut,mpert,natom,ntypat,selectz
1173 !arrays
1174  integer,intent(in) :: typat(natom)
1175  real(dp),intent(in) :: zion(ntypat)
1176  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
1177 
1178 !Local variables -------------------------
1179 !scalars
1180  integer :: idir1,idir2,ii,ipert1,ipert2
1181  character(len=500) :: msg
1182 !arrays
1183  real(dp) :: sumwght(2)
1184  real(dp),allocatable :: wghtat(:)
1185 
1186 ! *********************************************************************
1187 
1188  ABI_MALLOC(wghtat,(natom))
1189 
1190 !In case of acoustic sum rule imposition, compute the weights on each atom.
1191  if (chneut==1)then
1192 
1193 !  The weight is the same for all atom
1194    do ipert1=1,natom
1195      wghtat(ipert1)=1./natom
1196    end do
1197 
1198  else if (chneut==2) then
1199 
1200 !  The weight is proportional to the diagonal electronic screening charge of the atom
1201    sumwght(1)=zero
1202    do ipert1=1,natom
1203      wghtat(ipert1)=zero
1204      do idir1=1,3
1205        wghtat(ipert1)=wghtat(ipert1)+&
1206 &       d2cart(1,idir1,ipert1,idir1,natom+2)+&
1207 &       d2cart(1,idir1,natom+2,idir1,ipert1)-2*zion(typat(ipert1))
1208      end do
1209      sumwght(1)=sumwght(1)+wghtat(ipert1)
1210    end do
1211 
1212 !  Normalize the weights to unity
1213    do ipert1=1,natom
1214      wghtat(ipert1)=wghtat(ipert1)/sumwght(1)
1215    end do
1216  end if
1217 
1218 !Calculation of the violation of the charge neutrality
1219 !and imposition of the charge neutrality condition
1220  if (chneut/=0)then
1221    write(msg, '(a,a,a,a,a,a,a)' )&
1222     ' The violation of the charge neutrality conditions',ch10,&
1223     ' by the effective charges is as follows :',ch10,&
1224     '    atom        electric field',ch10,&
1225     ' displacement     direction   '
1226    call wrtout(ab_out,msg)
1227    do idir1=1,3
1228      do idir2=1,3
1229        do ii=1,2
1230          sumwght(ii)=zero
1231          do ipert1=1,natom
1232            sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir2,natom+2)
1233          end do
1234          do ipert1=1,natom
1235            d2cart(ii,idir1,ipert1,idir2,natom+2)=&
1236            d2cart(ii,idir1,ipert1,idir2,natom+2)-sumwght(ii)*wghtat(ipert1)
1237          end do
1238        end do
1239        write(msg, '(i8,i16,2f16.6)' ) idir1,idir2,sumwght(1),sumwght(2)
1240        call wrtout(ab_out,msg)
1241      end do
1242    end do
1243    write(msg, '(a)' )' '
1244    call wrtout(ab_out,msg)
1245 
1246 !  The same for the symmetrical part
1247    do idir1=1,3
1248      do idir2=1,3
1249        do ii=1,2
1250          sumwght(ii)=zero
1251          do ipert2=1,natom
1252            sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir2,ipert2)
1253          end do
1254          do ipert2=1,natom
1255            d2cart(ii,idir1,natom+2,idir2,ipert2)=&
1256            d2cart(ii,idir1,natom+2,idir2,ipert2)-sumwght(ii)*wghtat(ipert2)
1257          end do
1258        end do
1259      end do
1260    end do
1261  end if
1262 
1263 !Selection of the trace of the effective charge tensor attached to each atom
1264  if(selectz==1)then
1265    do ipert1=1,natom
1266      do ii=1,2
1267        sumwght(ii)=zero
1268        do idir1=1,3
1269          sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir1,natom+2)
1270        end do
1271        do idir1=1,3
1272          do idir2=1,3
1273            d2cart(ii,idir1,ipert1,idir2,natom+2)=zero
1274          end do
1275        end do
1276        do idir1=1,3
1277          d2cart(ii,idir1,ipert1,idir1,natom+2)=sumwght(ii)/3.0_dp
1278        end do
1279      end do
1280    end do
1281 !  Do the same for the symmetrical part of d2cart
1282    do ipert2=1,natom
1283      do ii=1,2
1284        sumwght(ii)=zero
1285        do idir1=1,3
1286          sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir1,ipert2)
1287        end do
1288        do idir1=1,3
1289          do idir2=1,3
1290            d2cart(ii,idir1,natom+2,idir2,ipert2)=zero
1291          end do
1292        end do
1293        do idir1=1,3
1294          d2cart(ii,idir1,natom+2,idir1,ipert2)=sumwght(ii)/3.0_dp
1295        end do
1296      end do
1297    end do
1298  end if
1299 
1300 !Selection of the symmetric part of the effective charge tensor attached to each atom
1301  if(selectz==2)then
1302    do ipert1=1,natom
1303      do ii=1,2
1304        do idir1=1,3
1305          do idir2=1,3
1306            sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
1307 &           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
1308            d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii)
1309            d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii)
1310          end do
1311        end do
1312      end do
1313    end do
1314 !  Do the same for the symmetrical part of d2cart
1315    do ipert1=1,natom
1316      do ii=1,2
1317        do idir1=1,3
1318          do idir2=1,3
1319            sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)&
1320 &           +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp
1321            d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii)
1322            d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii)
1323          end do
1324        end do
1325      end do
1326    end do
1327  end if
1328 
1329 !Write the effective charge tensor
1330  write(msg, '(a,a,a,a,a,a,a)' )&
1331    ' Effective charge tensors after ',ch10,&
1332    ' imposition of the charge neutrality,',ch10,&
1333    ' and eventual restriction to some part :',ch10,&
1334   '   atom    displacement  '
1335  call wrtout(ab_out,msg)
1336 
1337  do ipert1=1,natom
1338    do idir1=1,3
1339      write(msg, '(2i10,3es16.6)' )ipert1,idir1,(d2cart(1,idir1,ipert1,idir2,natom+2),idir2=1,3)
1340      call wrtout(ab_out,msg)
1341    end do
1342  end do
1343 
1344 !Zero the imaginary part of the dynamical matrix
1345  write(msg, '(a)' )' Now, the imaginary part of the dynamical matrix is zeroed '
1346  call wrtout(ab_out,msg)
1347  call wrtout(std_out,msg)
1348 
1349  do ipert1=1,natom
1350    do ipert2=1,natom
1351      do idir1=1,3
1352        do idir2=1,3
1353          d2cart(2,idir1,ipert1,idir2,ipert2)=zero
1354        end do
1355      end do
1356    end do
1357  end do
1358 
1359  ABI_FREE(wghtat)
1360 
1361 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

SOURCE

 939 subroutine d2cart_to_red(d2cart, d2red, gprimd, rprimd, mpert, natom, &
 940 &                        ntypat,typat,ucvol,zion)
 941 
 942 !Arguments -------------------------------
 943 !scalars
 944  integer,intent(in) :: mpert,natom,ntypat
 945  real(dp),intent(in) :: ucvol
 946 !arrays
 947  integer,intent(in) :: typat(natom)
 948  real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
 949  real(dp),intent(out) :: d2red(2,3,mpert,3,mpert)
 950  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
 951  real(dp),intent(in) :: zion(ntypat)
 952 
 953 !Local variables -------------------------
 954 !scalars
 955  integer :: idir1,idir2,ii,ipert1,ipert2
 956  real(dp) :: fac
 957 !arrays
 958  integer :: flg1(3),flg2(3)
 959  real(dp) :: vec1(3),vec2(3)
 960 
 961 ! *********************************************************************
 962 
 963  flg1 = one
 964  flg2 = one
 965 
 966  d2red = d2cart
 967 
 968 !Remove the ionic charges to z to get the change in effective charges
 969  do ipert1=1,natom
 970    do idir1=1,3
 971      d2red(1,idir1,ipert1,idir1,natom+2)=&
 972 &     d2red(1,idir1,ipert1,idir1,natom+2) - zion(typat(ipert1))
 973    end do
 974  end do
 975  do ipert2=1,natom
 976    do idir2=1,3
 977      d2red(1,idir2,natom+2,idir2,ipert2)=&
 978 &     d2red(1,idir2,natom+2,idir2,ipert2) - zion(typat(ipert2))
 979    end do
 980  end do
 981 
 982  ! Remove the vacuum polarizability from the dielectric tensor
 983  do idir1=1,3
 984    d2red(1,idir1,natom+2,idir1,natom+2)=&
 985 &   d2red(1,idir1,natom+2,idir1,natom+2) - 1.0_dp
 986  end do
 987 
 988 ! Scale the dielectric tensor with the volue of the unit cell
 989  do idir1=1,3
 990    do idir2=1,3
 991      do ii=1,2
 992        d2red(ii,idir1,natom+2,idir2,natom+2)=&
 993 &       - (ucvol / four_pi) * d2red(ii,idir1,natom+2,idir2,natom+2)
 994      end do
 995    end do
 996  end do
 997 
 998 !For the piezoelectric tensor, takes into account the volume of the unit cell
 999  do ipert2=natom+3,natom+4
1000    do idir1=1,3
1001      do idir2=1,3
1002        do ii=1,2
1003          d2red(ii,idir1,natom+2,idir2,ipert2)=&
1004 &         (ucvol)*d2red(ii,idir1,natom+2,idir2,ipert2)
1005          d2red(ii,idir2,ipert2,idir1,natom+2)=&
1006 &         (ucvol)*d2red(ii,idir2,ipert2,idir1,natom+2)
1007        end do
1008      end do
1009    end do
1010  end do
1011 
1012 ! Reduced coordinates transformation (in two steps)
1013 ! Note that rprimd and gprimd are swapped, compared to what cart39 expects
1014 ! A factor of (2pi) ** 2 is added to transform the electric field perturbations
1015 
1016 !First step
1017  do ipert1=1,mpert
1018    fac = one; if (ipert1==natom+2) fac = two_pi ** 2
1019 
1020    do ipert2=1,mpert
1021      do ii=1,2
1022        do idir1=1,3
1023          do idir2=1,3
1024            vec1(idir2)=d2red(ii,idir1,ipert1,idir2,ipert2)
1025          end do
1026          ! Transform vector from cartesian to reduced coordinates
1027          call cart39(flg1,flg2,transpose(rprimd),ipert1,natom,transpose(gprimd),vec1,vec2)
1028          do idir2=1,3
1029            d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) * fac
1030          end do
1031        end do
1032      end do
1033    end do
1034  end do
1035 
1036 !Second step
1037  do ipert1=1,mpert
1038    do ipert2=1,mpert
1039      fac = one; if (ipert2==natom+2) fac = two_pi ** 2
1040 
1041      do ii=1,2
1042        do idir2=1,3
1043          do idir1=1,3
1044            vec1(idir1)=d2red(ii,idir1,ipert1,idir2,ipert2)
1045          end do
1046          ! Transform vector from cartesian to reduced coordinates
1047          call cart39(flg1,flg2,transpose(rprimd),ipert2,natom,transpose(gprimd),vec1,vec2)
1048          do idir1=1,3
1049            d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) * fac
1050          end do
1051        end do
1052      end do
1053    end do
1054  end do
1055 
1056 
1057 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.

SOURCE

1424 subroutine d2sym3(blkflg,d2,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev,zero_by_symm)
1425 
1426 !Arguments -------------------------------
1427 !scalars
1428  integer,intent(in) :: mpert,natom,nsym,timrev,zero_by_symm
1429 !arrays
1430  integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym)
1431  integer,intent(in),target :: symrec(3,3,nsym),symrel(3,3,nsym)
1432  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
1433  real(dp),intent(in) :: qpt(3)
1434  real(dp),intent(inout) :: d2(2,3,mpert,3,mpert)
1435 
1436 !Local variables -------------------------
1437 !scalars
1438  logical, parameter :: do_final_sym=.true.
1439  logical :: qzero
1440  integer :: exch12,found,idir1,idir2,idisy1,idisy2,ipert1,ipert2
1441  integer :: ipesy1,ipesy2,isgn,isym,ithree,itirev,nblkflg_is_one,noccur,nsym_used,quit,quit1
1442  real(dp) :: arg1,arg2,im,norm,re,sumi,sumr,xi,xr
1443 !arrays
1444  integer,pointer :: sym1_(:,:,:),sym2_(:,:,:)
1445  real(dp),allocatable :: d2tmp1(:,:,:),d2tmp2(:,:,:),d2work(:,:,:,:,:)
1446 
1447 ! *********************************************************************
1448 
1449  qzero=(qpt(1)**2+qpt(2)**2+qpt(3)**2<tol16)
1450 
1451 !Here look after exchange of 1 and 2 axis,
1452 !for electric field in diamond symmetry
1453  exch12=0
1454  if (qzero) then
1455    do isym=1,nsym
1456      exch12=1
1457      if(symrel(1,1,isym)/=0)exch12=0
1458      if(symrel(1,2,isym)/=1)exch12=0
1459      if(symrel(1,3,isym)/=0)exch12=0
1460      if(symrel(2,1,isym)/=1)exch12=0
1461      if(symrel(2,2,isym)/=0)exch12=0
1462      if(symrel(2,3,isym)/=0)exch12=0
1463      if(symrel(3,1,isym)/=0)exch12=0
1464      if(symrel(3,2,isym)/=0)exch12=0
1465      if(symrel(3,3,isym)/=1)exch12=0
1466 !    if(exch12==1) write(std_out,*)' d2sym3 : found exchange 1 2 =',isym
1467      if(exch12==1)exit
1468    end do
1469  end if
1470 
1471 !Exchange of perturbations
1472 
1473 !Consider two cases : either time-reversal symmetry
1474 !conserves the wavevector, or not
1475  if(timrev==0)then
1476 
1477 !  do ipert1=1,mpert  See notes
1478    do ipert1=1,min(natom+2,mpert)
1479      do idir1=1,3
1480 
1481 !      Since the matrix is hermitian, the diagonal elements are real
1482        d2(2,idir1,ipert1,idir1,ipert1)=zero
1483 
1484 !      do ipert2=1,mpert See notes
1485        do ipert2=1,min(natom+2,mpert)
1486          do idir2=1,3
1487 
1488 !          If an element exists
1489            if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
1490 
1491 !            Either complete the symmetric missing element
1492              if(blkflg(idir2,ipert2,idir1,ipert1)==0)then
1493 
1494                d2(1,idir2,ipert2,idir1,ipert1)= d2(1,idir1,ipert1,idir2,ipert2)
1495                d2(2,idir2,ipert2,idir1,ipert1)=-d2(2,idir1,ipert1,idir2,ipert2)
1496 
1497                blkflg(idir2,ipert2,idir1,ipert1)=1
1498 
1499 !              Or symmetrize (the matrix is hermitian) in case both exists
1500 !              (Note : this opportunity has been disabled for more
1501 !              obvious search for bugs in the code )
1502 !              else
1503 !              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
1504 !              sumi=d2(1,idir2,ipert2,idir1,ipert1)-d2(1,idir1,ipert1,idir2,ipert2)
1505 !              d2(1,idir2,ipert2,idir1,ipert1)=half*sumr
1506 !              d2(1,idir1,ipert1,idir2,ipert2)=half*sumr
1507 !              d2(2,idir2,ipert2,idir1,ipert1)=half*sumi
1508 !              d2(2,idir1,ipert1,idir2,ipert2)=-half*sumi
1509              end if
1510            end if
1511 
1512          end do
1513        end do
1514 
1515      end do
1516    end do
1517 
1518 !  Here, case with time-reversal symmetry
1519  else
1520 
1521 !  do ipert1=1,mpert See notes
1522    do ipert1=1,min(natom+2,mpert)
1523      do idir1=1,3
1524 !      do ipert2=1,mpert See notes
1525        do ipert2=1,min(natom+2,mpert)
1526          do idir2=1,3
1527            d2(2,idir1,ipert1,idir2,ipert2)=zero
1528 
1529 !          If an element exists
1530            if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
1531 
1532 !            Either complete the symmetric missing element
1533              if(blkflg(idir2,ipert2,idir1,ipert1)==0)then
1534 
1535                d2(1,idir2,ipert2,idir1,ipert1)=d2(1,idir1,ipert1,idir2,ipert2)
1536                blkflg(idir2,ipert2,idir1,ipert1)=1
1537 
1538 !              Or symmetrize (the matrix is hermitian) in case both exists
1539 !              (Note : this opportunity has been disabled for more
1540 !              obvious search for bugs in the code )
1541 !              else
1542 !              sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2)
1543 !              d2(1,idir2,ipert2,idir1,ipert1)=half*sumr
1544 !              d2(1,idir1,ipert1,idir2,ipert2)=half*sumr
1545              end if
1546 
1547            end if
1548          end do
1549        end do
1550      end do
1551    end do
1552  end if
1553 
1554 !Big Big Loop : symmetrize three times, because
1555 !of some cases in which one element is not yet available
1556 !at the first pass, and even at the second one !
1557  do ithree=1,3
1558 
1559 !  Big loop on all elements
1560 !  do ipert1=1,mpert See notes
1561    do ipert1=1,min(natom+2,mpert)
1562 
1563 !    Select the symmetries according to pertubation 1
1564      if (ipert1<=natom)then
1565        sym1_ => symrec
1566      else
1567        sym1_ => symrel
1568      end if
1569 
1570      do idir1=1,3
1571 !      do ipert2=1,mpert See notes
1572        do ipert2=1,min(natom+2,mpert)
1573 
1574     !    Select the symmetries according to pertubation 2
1575          if (ipert2<=natom)then
1576            sym2_ => symrec
1577          else
1578            sym2_ => symrel
1579          end if
1580 
1581          do idir2=1,3
1582 
1583 !          Will get element (idir1,ipert1,idir2,ipert2)
1584 !          so this element should not yet be present ...
1585            if(blkflg(idir1,ipert1,idir2,ipert2)/=1)then
1586 
1587              d2(1,idir1,ipert1,idir2,ipert2)=zero
1588              d2(2,idir1,ipert1,idir2,ipert2)=zero
1589 
1590 !            Loop on all symmetries, including time-reversal
1591              quit1=0
1592              do isym=1,nsym
1593                do itirev=1,2
1594                  isgn=3-2*itirev
1595 
1596                  if(symq(4,itirev,isym)/=0)then
1597                    found=1
1598 
1599 !                  Here select the symmetric of ipert1
1600                    if(ipert1<=natom)then
1601                      ipesy1=indsym(4,isym,ipert1)
1602                    else if(ipert1==(natom+2).and.qzero)then
1603                      ipesy1=ipert1
1604                    else
1605                      found=0
1606                    end if
1607 
1608 !                  Here select the symmetric of ipert2
1609                    if(ipert2<=natom)then
1610                      ipesy2=indsym(4,isym,ipert2)
1611                    else if(ipert2==(natom+2).and.qzero)then
1612                      ipesy2=ipert2
1613                    else
1614                      found=0
1615                    end if
1616 
1617 !                  Now that a symmetric perturbation has been obtained,
1618 !                  including the expression of the symmetry matrix, see
1619 !                  if the symmetric values are available
1620                    if( found==1 ) then
1621 
1622                      sumr=zero
1623                      sumi=zero
1624                      noccur=0
1625                      nblkflg_is_one=0
1626                      quit=0
1627                      do idisy1=1,3
1628                        do idisy2=1,3
1629                          if(sym1_(idir1,idisy1,isym)/=0 .and. sym2_(idir2,idisy2,isym)/=0 )then
1630                            if(blkflg(idisy1,ipesy1,idisy2,ipesy2)==1)then
1631                              nblkflg_is_one=nblkflg_is_one+1
1632                              sumr=sumr+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*&
1633 &                             d2(1,idisy1,ipesy1,idisy2,ipesy2)
1634                              sumi=sumi+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*&
1635 &                             d2(2,idisy1,ipesy1,idisy2,ipesy2)
1636 
1637 !                            Here, in case the symmetric of the element
1638 !                            is the element, or the symmetric with
1639 !                            respect to permutation of perturbations
1640 !                            (some more conditions on the time-reversal
1641 !                            symmetry must be fulfilled although)
1642                            else if(  idisy1==idir1 .and. ipesy1==ipert1&
1643 &                             .and. idisy2==idir2 .and. ipesy2==ipert2&
1644 &                             .and.(isgn==1 .or. timrev==1 &
1645 &                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
1646 &                             then
1647                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1648                            else if(  idisy1==idir2 .and. ipesy1==ipert2&
1649 &                             .and. idisy2==idir1 .and. ipesy2==ipert1&
1650 &                             .and.(isgn==-1 .or. timrev==1&
1651 &                             .or. (idir1==idir2 .and. ipert1==ipert2)))&
1652 &                             then
1653                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1654 
1655 !                            Here, electric field case
1656                            else if( exch12==1 .and. &
1657 &                             ipert1==natom+2 .and. ipert2==natom+2&
1658 &                             .and.(( idisy1+idir1 ==3 &
1659 &                             .and. idisy2==3 .and. idir2==3)&
1660 &                             .or. ( idisy1+idir2 ==3&
1661 &                             .and. idisy2==3 .and. idir1==3)&
1662 &                             .or. ( idisy2+idir2 ==3&
1663 &                             .and. idisy1==3 .and. idir1==3)&
1664 &                             .or. ( idisy2+idir1 ==3&
1665 &                             .and. idisy1==3 .and. idir2==3)))&
1666 &                             then
1667                              noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)
1668 
1669                            else
1670 !                            Not found
1671                              found=0
1672                              quit=1
1673                              exit
1674                            end if
1675 
1676                          end if
1677                        end do
1678                        if(quit==1)exit
1679                      end do
1680                    end if
1681 
1682 !                  In case zero_by_symm==0, the computed materix element must be associated to at least one really computed matrix element
1683                    if(zero_by_symm==0 .and. nblkflg_is_one==0)then
1684                      found=0
1685                    endif
1686 
1687 !                  Now, if still found and associated to at least one really computed matrix element, put the correct value into array d2
1688                    if(found==1)then
1689 
1690 !                    In case of phonons, need to take into account the
1691 !                    time-reversal symmetry, and the shift back to the unit cell
1692 !
1693 !                    XG990712 : I am not sure this must be kept for electric field ...
1694 !                    1) Consider time-reversal symmetry
1695                      sumi=isgn*sumi
1696 
1697                      if(ipert1<=natom .and. ipert2<=natom)then
1698 !                      2) Shift the atoms back to the unit cell.
1699                        arg1=two_pi*( qpt(1)*indsym(1,isym,ipert1)&
1700 &                       +qpt(2)*indsym(2,isym,ipert1)&
1701 &                       +qpt(3)*indsym(3,isym,ipert1) )
1702                        arg2=two_pi*( qpt(1)*indsym(1,isym,ipert2)&
1703 &                       +qpt(2)*indsym(2,isym,ipert2)&
1704 &                       +qpt(3)*indsym(3,isym,ipert2) )
1705                        re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
1706 !                      XG010117 Must use isgn
1707                        im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))
1708                      else
1709                        re=one
1710                        im=zero
1711                      end if
1712 
1713 !                    Final check, could still fail if the
1714 !                    element was its own symmetric
1715                      if (abs(1.0_dp-re*noccur)< 1.0d-6.and.abs(im*noccur) <1.0d-6) then
1716                        found=0
1717                      end if
1718 
1719                    end if
1720 
1721                    if(found==1)then
1722 
1723                      if(noccur==0)then
1724                        d2(1,idir1,ipert1,idir2,ipert2)=re*sumr-im*sumi
1725                        d2(2,idir1,ipert1,idir2,ipert2)=re*sumi+im*sumr
1726                      else
1727 !                      See page July 2, 1994 in computer codes notebook
1728                        xr=re*sumr-im*sumi
1729                        xi=re*sumi+im*sumr
1730                        norm=one+noccur**2-two*re*noccur
1731                        xr=xr/norm
1732                        xi=xi/norm
1733                        d2(1,idir1,ipert1,idir2,ipert2)=&
1734 &                       (one-re*noccur)*xr-im*noccur*xi
1735                        d2(2,idir1,ipert1,idir2,ipert2)=&
1736 &                       (one-re*noccur)*xi+im*noccur*xr
1737                      end if
1738 
1739 !                    The element has been constructed !
1740                      blkflg(idir1,ipert1,idir2,ipert2)=1
1741 
1742                      quit1=1
1743                      exit ! Exit loop on symmetry operations
1744                    end if
1745 
1746 !                  End loop on all symmetries + time-reversal
1747                  end if
1748                end do
1749                if(quit1==1)exit
1750              end do
1751 
1752            end if
1753          end do ! End big loop on all elements
1754        end do
1755      end do
1756    end do
1757 
1758  end do !  End Big Big Loop
1759 
1760 !MT oct. 20, 2014:
1761 !Once the matrix has been built, it does not necessarily fulfill the correct symmetries.
1762 !It has just been filled up from rows or columns that only fulfill symmetries preserving
1763 !one particular perturbation.
1764 !An additional symmetrization might solve this (do not consider TR-symmetry)
1765  if (do_final_sym) then
1766    ABI_MALLOC(d2tmp1,(2,3,3))
1767    ABI_MALLOC(d2tmp2,(2,3,3))
1768    ABI_MALLOC(d2work,(2,3,mpert,3,mpert))
1769    d2work(:,:,:,:,:)=d2(:,:,:,:,:)
1770    do ipert1=1,min(natom+2,mpert)
1771      if ((ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).or.(ipert1==natom+2.and.(.not.qzero))) cycle
1772      if (ipert1<=natom)then
1773        sym1_ => symrec
1774      else
1775        sym1_ => symrel
1776      end if
1777      do ipert2=1,min(natom+2,mpert)
1778 !      if (any(blkflg(:,ipert1,:,ipert2)==0)) cycle
1779        if ((ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11).or.(ipert2==natom+2.and.(.not.qzero))) cycle
1780        if (ipert2<=natom)then
1781          sym2_ => symrec
1782        else
1783          sym2_ => symrel
1784        end if
1785        nsym_used=0
1786        d2tmp2(:,:,:)=zero
1787        do isym=1,nsym
1788          if (symq(4,1,isym)==1) then
1789            ipesy1=ipert1;if (ipert1<=natom) ipesy1=indsym(4,isym,ipert1)
1790            ipesy2=ipert2;if (ipert2<=natom) ipesy2=indsym(4,isym,ipert2)
1791 !          The condition on next line is too severe, since some elements of sym1_ or sym2_ might be zero,
1792 !          which means not all blkflg(:,ipesy1,:,ipesy2) would need to be 1 to symmetrize the matrix.
1793 !          However, coding something more refined is really more difficult.
1794 !          This condition then has the side effect that more symmetries can be applied when zero_by_symm==1,
1795 !          since blkflg can be set to 1 when the symmetries guarantee the matrix element to be zero.
1796            if (all(blkflg(:,ipesy1,:,ipesy2)==1)) then
1797              nsym_used=nsym_used+1
1798              re=one;im=zero
1799              if (ipert1<=natom.and.ipert2<=natom.and.(.not.qzero)) then
1800                arg1=two_pi*(qpt(1)*(indsym(1,isym,ipert1)-indsym(1,isym,ipert2)) &
1801 &                          +qpt(2)*(indsym(2,isym,ipert1)-indsym(2,isym,ipert2)) &
1802 &                          +qpt(3)*(indsym(3,isym,ipert1)-indsym(3,isym,ipert2)))
1803                re=cos(arg1);im=sin(arg1)
1804              end if
1805              d2tmp1(:,:,:)=zero
1806              do idir2=1,3 !kappa
1807                do idir1=1,3 !mu
1808                  do idisy1=1,3 !nu
1809                    d2tmp1(:,idir1,idir2)=d2tmp1(:,idir1,idir2) &
1810 &                     +sym1_(idir1,idisy1,isym)*d2(:,idisy1,ipesy1,idir2,ipesy2)
1811                  end do
1812                end do
1813              end do
1814              do idir2=1,3 !mu
1815                do idir1=1,3 !kappa
1816                  do idisy2=1,3 !nu
1817                    d2tmp2(1,idir1,idir2)=d2tmp2(1,idir1,idir2) &
1818 &                  +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*re-d2tmp1(2,idir1,idisy2)*im)
1819                    d2tmp2(2,idir1,idir2)=d2tmp2(2,idir1,idir2) &
1820 &                  +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*im+d2tmp1(2,idir1,idisy2)*re)
1821                  end do
1822                end do
1823              end do
1824            end if
1825          end if
1826        end do ! isym
1827        if (nsym_used>0) d2work(:,1:3,ipert1,1:3,ipert2)=d2tmp2(:,1:3,1:3)/dble(nsym_used)
1828      end do !ipert2
1829    end do !ipert1
1830    if (mpert>=natom)   d2(:,1:3,1:natom,1:3,1:natom)=d2work(:,1:3,1:natom,1:3,1:natom)
1831    if (mpert>=natom+2) then
1832      d2(:,1:3,natom+2,1:3,1:natom)=d2work(:,1:3,natom+2,1:3,1:natom)
1833      d2(:,1:3,1:natom,1:3,natom+2)=d2work(:,1:3,1:natom,1:3,natom+2)
1834      d2(:,1:3,natom+2,1:3,natom+2)=d2work(:,1:3,natom+2,1:3,natom+2)
1835    end if
1836    ABI_FREE(d2tmp1)
1837    ABI_FREE(d2tmp2)
1838    ABI_FREE(d2work)
1839  end if
1840 
1841 end subroutine d2sym3

m_dynmat/d3lwsym [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 d3lwsym

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

  has_strain = if .true. i2pert includes strain perturbation
  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 reduced space)
  symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space)
  symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian 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,m_nonlinear

CHILDREN

SOURCE

4466 !subroutine d3lwsym(blkflg,d3,has_strain,indsym,mpert,natom,nsym,symrec,symrel,symrel_cart)
4467 subroutine d3lwsym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel)
4468 
4469 !Arguments -------------------------------
4470 !scalars
4471  integer,intent(in) :: mpert,natom,nsym
4472 ! logical,intent(in) :: has_strain
4473 !arrays
4474  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4475  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
4476  real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert)
4477 ! real(dp),intent(in) :: symrel_cart(3,3,nsym)
4478 
4479 !Local variables -------------------------
4480 !scalars
4481  integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3
4482  integer :: ipesy1,ipesy2,ipesy3,isym,ithree
4483 !integer :: istr,i2dir_a,i2dir_b,disy2_a,idisy2_b
4484  real(dp) :: sumi,sumr
4485  logical :: is_strain
4486 !arrays
4487 ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
4488  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
4489 ! integer :: strflg(3,mpert,3,3,3,mpert),strflg_car(3,mpert,3,3,3,mpert)
4490 ! real(dp) :: d3str(2,3,mpert,3,3,3,mpert)
4491 
4492 ! *********************************************************************
4493 
4494 !First, take into account the permutations symmetry of
4495 !(i1pert,i1dir) and (i2pert,i2dir)
4496  do i1pert = 1, mpert
4497    do i2pert = 1, mpert
4498      do i3pert = 1, mpert
4499 
4500        do i1dir = 1, 3
4501          do i2dir = 1, 3
4502            do i3dir = 1, 3
4503 
4504              if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. &
4505               (blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert)/=1)) then
4506 
4507                d3(1,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = &
4508                d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
4509                d3(2,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = &
4510               -d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
4511 
4512                blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = 1
4513 
4514              end if
4515 
4516            end do
4517          end do
4518        end do
4519 
4520      end do
4521    end do
4522  end do
4523 
4524 !For strain perturbation we need an array with the two strain indexes
4525 ! if (has_strain) then
4526 !   strflg=0
4527 !   d3str=zero
4528 !   do i3pert=1, mpert
4529 !     do i3dir=1,3
4530 !       do i2pert=natom+3,natom+4
4531 !         do i2dir=1,3
4532 !           if (i2pert==natom+3) istr=i2dir
4533 !           if (i2pert==natom+4) istr=3+i2dir
4534 !           i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
4535 !           do i1pert=1,mpert
4536 !             do i1dir=1,3
4537 !               if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
4538 !                 strflg(i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)=1
4539 !                 d3str(:,i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)= &
4540 !               & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)          
4541 !                 if (i2pert==natom+4) then
4542 !                   strflg(i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)=1
4543 !                   d3str(:,i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)= &
4544 !                 & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)          
4545 !                 end if
4546 !               end if
4547 !             end do
4548 !           end do
4549 !         end do
4550 !       end do
4551 !     end do
4552 !   end do
4553 ! end if
4554 
4555 !Big Big Loop : symmetrize three times, because
4556 !of some cases in which one element is not yet available
4557 !at the first pass, and even at the second one !
4558 
4559  do ithree=1,3
4560 
4561 !  Loop over perturbations
4562    do i1pert = 1, mpert
4563      do i2pert = 1, mpert
4564        is_strain=.false.    
4565        do i3pert = 1, mpert
4566 
4567          do i1dir = 1, 3
4568            do i2dir = 1, 3
4569              do i3dir = 1, 3
4570 
4571 !              Will get element (idir1,ipert1,idir2,ipert2)
4572 !              so this element should not yet be present ...
4573                if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then
4574 
4575                  d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp
4576 
4577                  do isym = 1, nsym
4578 
4579                    found = 1
4580 
4581                    if (i1pert <= natom) then
4582                      ipesy1 = indsym(4,isym,i1pert)
4583                      sym1(:,:) = symrec(:,:,isym)
4584                    else if (i1pert == natom + 2) then
4585                      ipesy1 = i1pert
4586                      sym1(:,:) = symrel(:,:,isym)
4587                    else
4588                      found = 0
4589                    end if
4590 
4591                    if (i2pert <= natom) then
4592                      ipesy2 = indsym(4,isym,i2pert)
4593                      sym2(:,:) = symrec(:,:,isym)
4594                    else if (i2pert == natom + 2) then
4595                      ipesy2 = i2pert
4596                      sym2(:,:) = symrel(:,:,isym)
4597                    else if (i2pert == natom + 3.or. i2pert == natom + 4) then
4598                      !TODO: Symmetries on strain perturbation do not work yet.
4599                      found = 0
4600                      is_strain=.true.
4601 
4602 !                     ipesy2 = i2pert
4603 !                     sym2(:,:) = NINT(symrel_cart(:,:,isym))
4604 !                     if (i2pert==natom+3) istr=i2dir
4605 !                     if (i2pert==natom+4) istr=3+i2dir
4606 !                     i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
4607                    else
4608                      found = 0
4609                    end if
4610 
4611                    if (i3pert <= natom) then
4612                      ipesy3 = indsym(4,isym,i3pert)
4613                      sym3(:,:) = symrec(:,:,isym)
4614                    else if (i3pert == natom + 2.or.i3pert == natom + 8) then
4615                      ipesy3 = i3pert
4616                      sym3(:,:) = symrel(:,:,isym)
4617                    else
4618                      found = 0
4619                    end if
4620 
4621                    sumr = 0_dp ; sumi = 0_dp;
4622                    if (.not.is_strain) then
4623                      do idisy1 = 1, 3
4624                        do idisy2 = 1, 3
4625                          do idisy3 = 1, 3
4626 
4627                            if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) &
4628 &                           .and.(sym3(i3dir,idisy3) /=0)) then
4629 
4630                              if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then
4631 
4632                                sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4633 &                               sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4634                                sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4635 &                               sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4636 
4637                              else
4638 
4639                                found = 0
4640 
4641                              end if
4642 
4643                            end if
4644 
4645                          end do
4646                        end do
4647                      end do
4648                    else 
4649 !                     do idisy1 = 1, 3
4650 !                       !do idisy2_a = 1, 3
4651 !                       !  do idisy2_b = 1, 3
4652 !                       do idisy2 = 1, 3
4653 !                         if (ipesy2==natom+3) istr=idisy2
4654 !                         if (ipesy2==natom+4) istr=3+idisy2
4655 !                         idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr)
4656 !                           do idisy3 = 1, 3
4657 !
4658 !                             if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir_a,idisy2_a) /=0) &
4659 !&                             .and.(sym2(i2dir_b,idisy2_b) /=0).and.(sym3(i3dir,idisy3) /=0)) then
4660 !
4661 !                               if (strflg(idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3) == 1) then
4662 !
4663 !                                 sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_a)* &
4664 !&                                sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*& 
4665 !&                                d3str(1,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3)
4666 !                                 sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_b)*&
4667 !&                                sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*& 
4668 !&                                d3str(2,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3)
4669 !
4670 !                               else
4671 !
4672 !                                 found = 0
4673 !
4674 !                               end if
4675 !
4676 !                             end if
4677 !
4678 !                           end do
4679 !                         !end do
4680 !                       end do
4681 !                     end do
4682                    end if        
4683 
4684                    if (found == 1) then
4685                      d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
4686                      d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
4687                      blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4688                    end if
4689 
4690                  end do  ! isym
4691 
4692                end if  ! blkflg
4693 
4694 !              Close loop over perturbations
4695              end do
4696            end do
4697          end do
4698        end do
4699      end do
4700    end do
4701 
4702  end do  ! close loop over ithree
4703 
4704 end subroutine d3lwsym

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

SOURCE

4266 subroutine d3sym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel)
4267 
4268 !Arguments -------------------------------
4269 !scalars
4270  integer,intent(in) :: mpert,natom,nsym
4271 !arrays
4272  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4273  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
4274  real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert)
4275 
4276 !Local variables -------------------------
4277 !scalars
4278  integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3
4279  integer :: ipesy1,ipesy2,ipesy3,isym,ithree
4280  real(dp) :: sumi,sumr
4281 !arrays
4282  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
4283 
4284 ! *********************************************************************
4285 
4286 !DEBUG
4287 !write(std_out,*)'d3sym : enter'
4288 !do i1dir = 1, 3
4289 !do i2dir = 1, 3
4290 !do i3dir = 1, 3
4291 !write(std_out,*)i1dir,i2dir,i3dir,blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
4292 !end do
4293 !end do
4294 !end do
4295 !stop
4296 !ENDDEBUG
4297 
4298 !First, take into account the permutations symmetry of
4299 !(i1pert,i1dir) and (i3pert,i3dir)
4300  do i1pert = 1, mpert
4301    do i2pert = 1, mpert
4302      do i3pert = 1, mpert
4303 
4304        do i1dir = 1, 3
4305          do i2dir = 1, 3
4306            do i3dir = 1, 3
4307 
4308              if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. &
4309 &             (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=1)) then
4310 
4311                d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = &
4312 &              d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
4313 
4314                blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = 1
4315 
4316              end if
4317 
4318            end do
4319          end do
4320        end do
4321 
4322      end do
4323    end do
4324  end do
4325 
4326 !Big Big Loop : symmetrize three times, because
4327 !of some cases in which one element is not yet available
4328 !at the first pass, and even at the second one !
4329 
4330  do ithree=1,3
4331 
4332 !  Loop over perturbations
4333    do i1pert = 1, mpert
4334      do i2pert = 1, mpert
4335        do i3pert = 1, mpert
4336 
4337          do i1dir = 1, 3
4338            do i2dir = 1, 3
4339              do i3dir = 1, 3
4340 
4341 !              Will get element (idir1,ipert1,idir2,ipert2)
4342 !              so this element should not yet be present ...
4343                if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then
4344 
4345                  d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp
4346 
4347                  do isym = 1, nsym
4348 
4349                    found = 1
4350 
4351                    if (i1pert <= natom) then
4352                      ipesy1 = indsym(4,isym,i1pert)
4353                      sym1(:,:) = symrec(:,:,isym)
4354                    else if (i1pert == natom + 2) then
4355                      ipesy1 = i1pert
4356                      sym1(:,:) = symrel(:,:,isym)
4357                    else
4358                      found = 0
4359                    end if
4360 
4361                    if (i2pert <= natom) then
4362                      ipesy2 = indsym(4,isym,i2pert)
4363                      sym2(:,:) = symrec(:,:,isym)
4364                    else if (i2pert == natom + 2) then
4365                      ipesy2 = i2pert
4366                      sym2(:,:) = symrel(:,:,isym)
4367                    else
4368                      found = 0
4369                    end if
4370 
4371                    if (i3pert <= natom) then
4372                      ipesy3 = indsym(4,isym,i3pert)
4373                      sym3(:,:) = symrec(:,:,isym)
4374                    else if (i3pert == natom + 2) then
4375                      ipesy3 = i3pert
4376                      sym3(:,:) = symrel(:,:,isym)
4377                    else
4378                      found = 0
4379                    end if
4380 
4381                    sumr = 0_dp ; sumi = 0_dp;
4382                    do idisy1 = 1, 3
4383                      do idisy2 = 1, 3
4384                        do idisy3 = 1, 3
4385 
4386                          if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) &
4387 &                         .and.(sym3(i3dir,idisy3) /=0)) then
4388 
4389                            if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then
4390 
4391                              sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4392 &                             sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4393                              sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
4394 &                             sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
4395 
4396                            else
4397 
4398                              found = 0
4399 
4400                            end if
4401 
4402                          end if
4403 
4404                        end do
4405                      end do
4406                    end do
4407 
4408                    if (found == 1) then
4409                      d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
4410                      d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
4411                      blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4412                    end if
4413 
4414                  end do  ! isym
4415 
4416                end if  ! blkflg
4417 
4418 !              Close loop over perturbations
4419              end do
4420            end do
4421          end do
4422        end do
4423      end do
4424    end do
4425 
4426  end do  ! close loop over ithree
4427 
4428 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.

SOURCE

5660 subroutine dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,&
5661 & mpert,msym,natom,nsym,ntypat,phfrq,qphnrm,qphon,rprimd,&
5662 & symdynmat,symrel,symafm,typat,ucvol)
5663 
5664 !Arguments -------------------------------
5665 !scalars
5666  integer,intent(in) :: mpert,msym,natom,nsym,ntypat,symdynmat
5667  real(dp),intent(in) :: qphnrm,ucvol
5668 !arrays
5669  integer,intent(in) :: indsym(4,msym*natom),symrel(3,3,nsym),typat(natom)
5670  integer,intent(in) :: symafm(nsym)
5671  real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),rprimd(3,3)
5672  real(dp),intent(inout) :: qphon(3)
5673  real(dp),intent(out) :: displ(2*3*natom*3*natom),eigval(3*natom)
5674  real(dp),intent(out) :: eigvec(2*3*natom*3*natom),phfrq(3*natom)
5675 
5676 !Local variables -------------------------
5677 !scalars
5678  integer :: analyt,i1,i2,idir1,idir2,ier,ii,imode,ipert1,ipert2
5679  integer :: jmode,indexi,indexj,index
5680  real(dp) :: epsq,qphon2
5681  logical,parameter :: debug = .False.
5682  real(dp) :: sc_prod
5683 !arrays
5684  real(dp) :: qptn(3),dum(2,0) !, tsec(2)
5685  real(dp),allocatable :: matrx(:,:),zeff(:,:),zhpev1(:,:),zhpev2(:)
5686 
5687 ! *********************************************************************
5688 
5689  ! Keep track of time spent in dfpt_phfrq
5690  !call timab(1751, 1, tsec)
5691 
5692  ! Prepare the diagonalisation: analytical part.
5693  ! Note: displ is used as work space here
5694  i1=0
5695  do ipert1=1,natom
5696    do idir1=1,3
5697      i1=i1+1
5698      i2=0
5699      do ipert2=1,natom
5700        do idir2=1,3
5701          i2=i2+1
5702          index=i1+3*natom*(i2-1)
5703          displ(2*index-1)=d2cart(1,idir1,ipert1,idir2,ipert2)
5704          displ(2*index  )=d2cart(2,idir1,ipert1,idir2,ipert2)
5705        end do
5706      end do
5707    end do
5708  end do
5709 
5710  ! Determine the analyticity of the matrix.
5711  analyt=1; if(abs(qphnrm)<tol8) analyt=0
5712  if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=2
5713 
5714  ! In case of q=Gamma, only the real part is used
5715  if(analyt==0 .or. analyt==2)then
5716    do i1=1,3*natom
5717      do i2=1,3*natom
5718        index=i1+3*natom*(i2-1)
5719        displ(2*index)=zero
5720      end do
5721    end do
5722  end if
5723 
5724  ! In the case the non-analyticity is required:
5725  ! the tensor is in cartesian coordinates and this means that qphon must be in given in Cartesian coordinates.
5726  if(analyt==0)then
5727 
5728    ! Normalize the limiting direction
5729    qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2
5730    qphon(:)=qphon(:)/sqrt(qphon2)
5731 
5732    ! Get the dielectric constant for the limiting direction
5733    epsq=zero
5734    do idir1=1,3
5735      do idir2=1,3
5736        epsq= epsq + qphon(idir1)*qphon(idir2) * d2cart(1,idir1,natom+2,idir2,natom+2)
5737      end do
5738    end do
5739 
5740    ABI_MALLOC(zeff,(3,natom))
5741 
5742    ! Get the effective charges for the limiting direction
5743    do idir1=1,3
5744      do ipert1=1,natom
5745        zeff(idir1,ipert1)=zero
5746        do idir2=1,3
5747          zeff(idir1,ipert1) = zeff(idir1,ipert1) + qphon(idir2)* d2cart(1,idir1,ipert1,idir2,natom+2)
5748        end do
5749      end do
5750    end do
5751 
5752    ! Get the non-analytical part of the dynamical matrix, and suppress its imaginary part.
5753    i1=0
5754    do ipert1=1,natom
5755      do idir1=1,3
5756        i1=i1+1
5757        i2=0
5758        do ipert2=1,natom
5759          do idir2=1,3
5760            i2=i2+1
5761            index=i1+3*natom*(i2-1)
5762            displ(2*index-1)=displ(2*index-1)+four_pi/ucvol*zeff(idir1,ipert1)*zeff(idir2,ipert2)/epsq
5763            displ(2*index  )=zero
5764          end do
5765        end do
5766      end do
5767    end do
5768 
5769    ABI_FREE(zeff)
5770  end if !  End of the non-analyticity treatment
5771 
5772  ! Multiply IFC(q) by masses
5773  call massmult_and_breaksym(natom, ntypat, typat, amu, displ)
5774 
5775  ! ***********************************************************************
5776  ! Diagonalize the dynamical matrix
5777 
5778  !Symmetrize the dynamical matrix
5779  !FIXME: swap the next 2 lines and update test files to include symmetrization
5780  !       for Gamma point too (except in non-analytic case)
5781  !if (symdynmat==1 .and. analyt > 0) then
5782  if (symdynmat==1 .and. analyt == 1) then
5783    qptn(:)=qphon(:)
5784    if (analyt==1) qptn(:)=qphon(:)/qphnrm
5785    call symdyma(displ,indsym,natom,nsym,qptn,rprimd,symrel,symafm)
5786  end if
5787 
5788  ii=1
5789  ABI_MALLOC(matrx,(2,(3*natom*(3*natom+1))/2))
5790  do i2=1,3*natom
5791    do i1=1,i2
5792      matrx(1,ii)=displ(1+2*(i1-1)+2*(i2-1)*3*natom)
5793      matrx(2,ii)=displ(2+2*(i1-1)+2*(i2-1)*3*natom)
5794      ii=ii+1
5795    end do
5796  end do
5797 
5798  ABI_MALLOC(zhpev1,(2,2*3*natom-1))
5799  ABI_MALLOC(zhpev2,(3*3*natom-2))
5800 
5801  call ZHPEV ('V','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,zhpev2,ier)
5802  ABI_CHECK(ier == 0, sjoin('zhpev returned:', itoa(ier)))
5803 
5804  ABI_FREE(matrx)
5805  ABI_FREE(zhpev1)
5806  ABI_FREE(zhpev2)
5807 
5808  if (debug) then
5809    ! Check the orthonormality of the eigenvectors
5810    do imode=1,3*natom
5811      do jmode=imode,3*natom
5812        indexi=2*3*natom*(imode-1)
5813        indexj=2*3*natom*(jmode-1)
5814        sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom))
5815        write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod
5816      end do
5817    end do
5818  end if
5819 
5820  !***********************************************************************
5821 
5822  ! Get the phonon frequencies (negative by convention, if the eigenvalue of the dynamical matrix is negative)
5823  do imode=1,3*natom
5824    if(eigval(imode)>=1.0d-16)then
5825      phfrq(imode)=sqrt(eigval(imode))
5826    else if(eigval(imode)>=-1.0d-16)then
5827      phfrq(imode)=zero
5828    else
5829      phfrq(imode)=-sqrt(-eigval(imode))
5830    end if
5831  end do
5832 
5833  ! Fix the phase of the eigenvectors
5834  call fxphas_seq(eigvec,dum, 0, 0, 1, 3*natom*3*natom, 0, 3*natom, 3*natom, 0)
5835 
5836  ! Normalise the eigenvectors
5837  call pheigvec_normalize(natom, eigvec)
5838 
5839  ! Get the phonon displacements
5840  call phdispl_from_eigvec(natom, ntypat, typat, amu, eigvec, displ)
5841 
5842  if (debug) then
5843    write(std_out,'(a)')' Phonon eigenvectors and displacements '
5844    do imode=1,3*natom
5845      indexi=2*3*natom*(imode-1)
5846      write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' eigvec(1:6*natom)=',eigvec(indexi+1:indexi+6*natom)
5847      write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' displ(1:6*natom)=',displ(indexi+1:indexi+6*natom)
5848    end do
5849 
5850    ! Check the orthonormality of the eigenvectors
5851    do imode=1,3*natom
5852      do jmode=imode,3*natom
5853        indexi=2*3*natom*(imode-1)
5854        indexj=2*3*natom*(jmode-1)
5855        sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom))
5856        write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod
5857      end do
5858    end do
5859  end if
5860 
5861  !call timab(1751, 2, tsec)
5862 
5863 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

SOURCE

6014 subroutine dfpt_prtph(displ,eivec,enunit,iout,natom,phfrq,qphnrm,qphon)
6015 
6016 !Arguments -------------------------------
6017 !scalars
6018  integer,intent(in) :: eivec,enunit,iout,natom
6019  real(dp),intent(in) :: qphnrm
6020 !arrays
6021  real(dp),intent(in) :: displ(2,3*natom,3*natom),phfrq(3*natom),qphon(3)
6022 
6023 !Local variables -------------------------
6024 !scalars
6025  integer :: i,idir,ii,imode,jj
6026  real(dp) :: tolerance
6027  logical :: t_degenerate
6028  character(len=500) :: msg
6029 !arrays
6030  real(dp) :: vecti(3),vectr(3)
6031  character(len=1) :: metacharacter(3*natom)
6032 
6033 ! *********************************************************************
6034 
6035 !Check the value of eivec
6036  if (all(eivec /= [0,1,2,4])) then
6037    write(msg, '(a,i0,a,a)' )&
6038    'In the calling subroutine, eivec is',eivec,ch10,&
6039    'but allowed values are between 0 and 4.'
6040    ABI_BUG(msg)
6041  end if
6042 
6043 !write the phonon frequencies on unit std_out
6044  write(msg,'(4a)' )' ',ch10,' phonon wavelength (reduced coordinates) , ','norm, and energies in hartree'
6045  call wrtout(std_out,msg)
6046 
6047 !The next format should be rewritten
6048  write(msg,'(a,4f5.2)' )' ',(qphon(i),i=1,3),qphnrm
6049  call wrtout(std_out,msg)
6050  do jj=1,3*natom,5
6051    if (3*natom-jj<5) then
6052      write(msg,'(5es17.9)') (phfrq(ii),ii=jj,3*natom)
6053    else
6054      write(msg,'(5es17.9)') (phfrq(ii),ii=jj,jj+4)
6055    end if
6056    call wrtout(std_out,msg)
6057  end do
6058  write(msg,'(a,a,es17.9)') ch10,' Zero Point Motion energy (sum of freqs/2)=',sum(phfrq(1:3*natom))/2
6059  call wrtout(std_out,msg)
6060 
6061 !Put the wavevector in nice format
6062  if(iout>=0)then
6063    call wrtout(iout,' ')
6064    if(qphnrm/=0.0_dp)then
6065      write(msg, '(a,3f9.5)' )&
6066      '  Phonon wavevector (reduced coordinates) :',(qphon(i)/qphnrm+tol10,i=1,3)
6067    else
6068      write(msg, '(3a,3f9.5)' )&
6069      '  Phonon at Gamma, with non-analyticity in the',ch10,&
6070      '  direction (cartesian coordinates)',qphon(1:3)+tol10
6071    end if
6072    call wrtout(iout,msg)
6073 
6074 !  Write it, in different units.
6075    if(enunit/=1)then
6076      write(iout, '(a)' )' Phonon energies in Hartree :'
6077      do jj=1,3*natom,5
6078        if (3*natom-jj<5) then
6079          write(msg, '(1x,5es14.6)') (phfrq(ii),ii=jj,3*natom)
6080        else
6081          write(msg, '(1x,5es14.6)') (phfrq(ii),ii=jj,jj+4)
6082        end if
6083        call wrtout(iout,msg)
6084      end do
6085    end if
6086    if(enunit/=0)then
6087      write(iout, '(a)' )' Phonon energies in meV     :'
6088      do jj=1,3*natom,5
6089        if (3*natom-jj<5) then
6090          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,3*natom)
6091        else
6092          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,jj+4)
6093        end if
6094        call wrtout(iout,msg)
6095      end do
6096    end if
6097    if(enunit/=1)then
6098      write(iout, '(a)' )' Phonon frequencies in cm-1    :'
6099      do jj=1,3*natom,5
6100        if (3*natom-jj<5) then
6101          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,3*natom)
6102        else
6103          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,jj+4)
6104        end if
6105        call wrtout(iout,msg)
6106      end do
6107    end if
6108    if(enunit/=0)then
6109      write(iout, '(a)' )' Phonon frequencies in Thz     :'
6110      do jj=1,3*natom,5
6111        if (3*natom-jj<5) then
6112          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,3*natom)
6113        else
6114          write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,jj+4)
6115        end if
6116        call wrtout(iout,msg)
6117      end do
6118    end if
6119    if(enunit/=0.and.enunit/=1)then
6120      write(iout, '(a)' )' Phonon energies in Kelvin  :'
6121      do jj=1,3*natom,5
6122        if (3*natom-jj<5) then
6123          write(msg, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,3*natom)
6124        else
6125          write(msg, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,jj+4)
6126        end if
6127        call wrtout(iout,msg)
6128      end do
6129    end if
6130  end if
6131 
6132 !Take care of the eigendisplacements
6133  if(eivec==1 .or. eivec==2)then
6134    write(msg, '(a,a,a,a,a,a,a,a)' ) ch10,&
6135    ' Eigendisplacements ',ch10,&
6136    ' (will be given, for each mode : in cartesian coordinates',ch10,&
6137    '   for each atom the real part of the displacement vector,',ch10,&
6138    '   then the imaginary part of the displacement vector - absolute values smaller than 1.0d-7 are set to zero)'
6139    call wrtout(std_out,msg)
6140    if(iout>=0) then
6141      call wrtout(iout,msg)
6142    end if
6143 
6144 !  Examine the degeneracy of each mode. The portability of the echo of the eigendisplacements
6145 !  is very hard to obtain, and has not been attempted.
6146    do imode=1,3*natom
6147 !    The degenerate modes are not portable
6148      t_degenerate=.false.
6149      if(imode>1)then
6150        if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
6151      end if
6152      if(imode<3*natom)then
6153        if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
6154      end if
6155      metacharacter(imode)=';'; if(t_degenerate)metacharacter(imode)='-'
6156    end do
6157 
6158    do imode=1,3*natom
6159      write(msg,'(a,i4,a,es16.6)' )'  Mode number ',imode,'   Energy',phfrq(imode)
6160      call wrtout(std_out,msg)
6161      if(iout>=0)then
6162        write(msg, '(a,i4,a,es16.6)' )'  Mode number ',imode,'   Energy',phfrq(imode)
6163        call wrtout(iout,msg)
6164      end if
6165      tolerance=1.0d-7
6166      if(abs(phfrq(imode))<1.0d-5)tolerance=2.0d-7
6167      if(phfrq(imode)<1.0d-5)then
6168        write(msg,'(3a)' )' Attention : low frequency mode.',ch10,&
6169        '   (Could be unstable or acoustic mode)'
6170        call wrtout(std_out,msg)
6171        if(iout>=0)then
6172          write(iout, '(3a)' )' Attention : low frequency mode.',ch10,&
6173          '   (Could be unstable or acoustic mode)'
6174        end if
6175      end if
6176      do ii=1,natom
6177        do idir=1,3
6178          vectr(idir)=displ(1,idir+(ii-1)*3,imode)
6179          if(abs(vectr(idir))<tolerance)vectr(idir)=0.0_dp
6180          vecti(idir)=displ(2,idir+(ii-1)*3,imode)
6181          if(abs(vecti(idir))<tolerance)vecti(idir)=0.0_dp
6182        end do
6183        write(msg,'(i4,3es16.8,a,4x,3es16.8)' ) ii,vectr(:),ch10,vecti(:)
6184        call wrtout(std_out,msg)
6185        if(iout>=0)then
6186          write(msg,'(a,i3,3es16.8,2a,3x,3es16.8)') metacharacter(imode),ii,vectr(:),ch10,&
6187            metacharacter(imode), vecti(:)
6188          call wrtout(iout,msg)
6189        end if
6190      end do
6191    end do
6192  end if
6193 
6194 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.

SOURCE

2359 subroutine dfpt_sydy(cplex,dyfrow,indsym,natom,nondiag,nsym,qphon,sdyfro,symq,symrec)
2360 
2361 !Arguments -------------------------------
2362 !scalars
2363  integer,intent(in) :: cplex,natom,nondiag,nsym
2364 !arrays
2365  integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym),symrec(3,3,nsym)
2366  real(dp),intent(in) :: dyfrow(cplex,3,3,natom,1+(natom-1)*nondiag),qphon(3)
2367  real(dp),intent(out) :: sdyfro(cplex,3,3,natom,1+(natom-1)*nondiag)
2368 
2369 !Local variables -------------------------
2370 !scalars
2371  integer :: ia,indi,indj,isym,ja,kappa,mu,natom_nondiag,nsym_used,nu
2372  logical :: qeq0
2373  real(dp) :: arg,div,phasei,phaser
2374 !arrays
2375  real(dp) :: work(cplex,3,3)
2376 
2377 ! *********************************************************************
2378 
2379  if (nsym==1) then
2380 
2381 !  Only symmetry is identity so simply copy
2382    sdyfro(:,:,:,:,:)=dyfrow(:,:,:,:,:)
2383 
2384  else
2385 
2386 !  Actually carry out symmetrization
2387    sdyfro(:,:,:,:,:)=zero
2388    qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)
2389 !  === Diagonal dyn. matrix OR q=0
2390    if (nondiag==0.or.qeq0) then
2391      natom_nondiag=1;if (nondiag==1) natom_nondiag=natom
2392      do ja=1,natom_nondiag
2393        do ia=1,natom
2394          do isym=1,nsym
2395            indi=indsym(4,isym,ia)
2396            indj=1;if (nondiag==1) indj=indsym(4,isym,ja)
2397            work(:,:,:)=zero
2398            do mu=1,3
2399              do nu=1,3
2400                do kappa=1,3
2401                  work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj)
2402                end do
2403              end do
2404            end do
2405            do mu=1,3
2406              do nu=1,3
2407                do kappa=1,3
2408                  sdyfro(:,kappa,mu,ia,ja)=sdyfro(:,kappa,mu,ia,ja)+symrec(mu,nu,isym)*work(:,kappa,nu)
2409                end do
2410              end do
2411            end do
2412          end do
2413        end do
2414      end do
2415      div=one/dble(nsym)
2416      sdyfro(:,:,:,:,:)=div*sdyfro(:,:,:,:,:)
2417 !    === Non diagonal dyn. matrix AND q<>0
2418    else
2419      do ja=1,natom
2420        do ia=1,natom
2421          nsym_used=0
2422          do isym=1,nsym
2423            if (symq(4,1,isym)==1) then
2424              arg=two_pi*(qphon(1)*(indsym(1,isym,ia)-indsym(1,isym,ja)) &
2425 &             +qphon(2)*(indsym(2,isym,ia)-indsym(2,isym,ja)) &
2426 &             +qphon(3)*(indsym(3,isym,ia)-indsym(3,isym,ja)))
2427              phaser=cos(arg);phasei=sin(arg)
2428              nsym_used=nsym_used+1
2429              indi=indsym(4,isym,ia)
2430              indj=indsym(4,isym,ja)
2431              work(:,:,:)=zero
2432              do mu=1,3
2433                do nu=1,3
2434                  do kappa=1,3
2435                    work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj)
2436                  end do
2437                end do
2438              end do
2439              do mu=1,3
2440                do nu=1,3
2441                  do kappa=1,3
2442                    sdyfro(1,kappa,mu,ia,ja)=sdyfro(1,kappa,mu,ia,ja) &
2443 &                   +symrec(mu,nu,isym)*(work(1,kappa,nu)*phaser-work(2,kappa,nu)*phasei)
2444                  end do
2445                end do
2446              end do
2447              if (cplex==2) then
2448                do mu=1,3
2449                  do nu=1,3
2450                    do kappa=1,3
2451                      sdyfro(2,kappa,mu,ia,ja)=sdyfro(2,kappa,mu,ia,ja) &
2452 &                     +symrec(mu,nu,isym)*(work(1,kappa,nu)*phasei+work(2,kappa,nu)*phaser)
2453                    end do
2454                  end do
2455                end do
2456              end if
2457            end if
2458          end do
2459          div=one/dble(nsym_used)
2460          sdyfro(:,:,:,ia,ja)=div*sdyfro(:,:,:,ia,ja)
2461        end do
2462      end do
2463    end if
2464 
2465  end if
2466 
2467 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...

SOURCE

2232 subroutine dfpt_sygra(natom,desym,deunsy,indsym,ipert,nsym,qpt,symrec)
2233 
2234 !Arguments -------------------------------
2235 !scalars
2236  integer,intent(in) :: ipert,natom,nsym
2237 !arrays
2238  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym)
2239  real(dp),intent(in) :: deunsy(2,3,natom),qpt(3)
2240  real(dp),intent(out) :: desym(2,3,natom)
2241 
2242 !Local variables -------------------------
2243 !scalars
2244  integer :: ia,ind,isym,mu
2245  real(dp) :: arg,im,re,sumi,sumr
2246 
2247 ! *********************************************************************
2248 
2249 !DEBUG
2250 !write(std_out,*)' dfpt_sygra : enter '
2251 !write(std_out,*)' dfpt_sygra : qpt(:)',qpt(:)
2252 !do ia=1,natom
2253 !do mu=1,3
2254 !write(std_out,*)' dfpt_sygra : deunsy(:2,mu,ia)',deunsy(:2,mu,ia)
2255 !enddo
2256 !enddo
2257 !ENDDEBUG
2258 
2259  if (nsym==1) then
2260 
2261 !  Only symmetry is identity so simply copy
2262    desym(:,:,:)=deunsy(:,:,:)
2263 
2264  else
2265 
2266 !  Actually conduct symmetrization
2267 !  write(std_out,*)' dfpt_sygra : desym(:2,:3,:natom),qpt(:)',desym(:2,:3,:natom),qpt(:)
2268    do ia=1,natom
2269 !    write(std_out,*)' dfpt_sygra : ia=',ia
2270      do mu=1,3
2271        sumr=zero
2272        sumi=zero
2273 !      write(std_out,*)' dfpt_sygra : mu=',mu
2274        do isym=1,nsym
2275          ind=indsym(4,isym,ia)
2276 !        Must shift the atoms back to the unit cell.
2277 !        arg=two_pi*( qpt(1)*indsym(1,isym,ia)&
2278 !        &         +qpt(2)*indsym(2,isym,ia)&
2279 !        &         +qpt(3)*indsym(3,isym,ia) )
2280 !        Selection of non-zero q point, to avoid ipert being outside the 1 ... natom range
2281          if(qpt(1)**2+qpt(2)**2+qpt(3)**2 > tol16)then
2282            arg=two_pi*( qpt(1)*(indsym(1,isym,ia)-indsym(1,isym,ipert))&
2283 &           +qpt(2)* (indsym(2,isym,ia)-indsym(2,isym,ipert))&
2284 &           +qpt(3)* (indsym(3,isym,ia)-indsym(3,isym,ipert)))
2285          else
2286            arg=zero
2287          end if
2288 
2289          re=dble(symrec(mu,1,isym))*deunsy(1,1,ind)+&
2290 &         dble(symrec(mu,2,isym))*deunsy(1,2,ind)+&
2291 &         dble(symrec(mu,3,isym))*deunsy(1,3,ind)
2292          im=dble(symrec(mu,1,isym))*deunsy(2,1,ind)+&
2293 &         dble(symrec(mu,2,isym))*deunsy(2,2,ind)+&
2294 &         dble(symrec(mu,3,isym))*deunsy(2,3,ind)
2295          sumr=sumr+re*cos(arg)-im*sin(arg)
2296          sumi=sumi+re*sin(arg)+im*cos(arg)
2297 !        sumr=sumr+re
2298 !        sumi=sumi+im
2299 !        write(std_out,*)' dfpt_sygra : isym,indsym(4,isym,ia),arg,re,im,sumr,sumi',&
2300 !        &      isym,indsym(4,isym,ia),arg,re,im,sumr,sumi
2301        end do
2302        desym(1,mu,ia)=sumr/dble(nsym)
2303        desym(2,mu,ia)=sumi/dble(nsym)
2304 !      write(std_out,*)' dfpt_sygra : desym(:,mu,ia)',desym(:,mu,ia)
2305      end do
2306    end do
2307  end if
2308 
2309 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

SOURCE

3428 subroutine dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt)
3429 
3430 !Arguments -------------------------------
3431 !scalars
3432  integer,intent(in) :: natom,nrpt
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) :: dist(natom,natom,nrpt)
3437 
3438 !Local variables -------------------------
3439 !scalars
3440  integer :: ia,ib,ii,irpt
3441 !arrays
3442  real(dp) :: ra(3),rb(3),rdiff(3),red(3),rptcar(3),xred(3)
3443 
3444 ! *********************************************************************
3445 
3446 !BIG loop on all generic atoms
3447  do ia=1,natom
3448    ! First transform canonical coordinates to reduced coordinates
3449    do ii=1,3
3450      xred(ii)=gprim(1,ii)*rcan(1,ia)+gprim(2,ii)*rcan(2,ia)+gprim(3,ii)*rcan(3,ia)
3451    end do
3452    ! Then to cartesian coordinates
3453    ra(:)=xred(1)*acell(1)*rprim(:,1)+xred(2)*acell(2)*rprim(:,2)+xred(3)*acell(3)*rprim(:,3)
3454    do ib=1,natom
3455      do ii=1,3
3456        xred(ii)=gprim(1,ii)*rcan(1,ib)+gprim(2,ii)*rcan(2,ib)+gprim(3,ii)*rcan(3,ib)
3457      end do
3458      do ii=1,3
3459        rb(ii)=xred(1)*acell(1)*rprim(ii,1)+xred(2)*acell(2)*rprim(ii,2)+xred(3)*acell(3)*rprim(ii,3)
3460      end do
3461      do irpt=1,nrpt
3462        ! First transform it to reduced coordinates
3463        do ii=1,3
3464          red(ii)=gprim(1,ii)*rpt(1,irpt)+gprim(2,ii)*rpt(2,irpt)+gprim(3,ii)*rpt(3,irpt)
3465        end do
3466        ! Then to cartesian coordinates
3467        do ii=1,3
3468          rptcar(ii)=red(1)*acell(1)*rprim(ii,1)+red(2)*acell(2)*rprim(ii,2)+red(3)*acell(3)*rprim(ii,3)
3469        end do
3470        do ii=1,3
3471          rdiff(ii)=-rptcar(ii)+ra(ii)-rb(ii)
3472        end do
3473        dist(ia,ib,irpt)=(rdiff(1)**2+rdiff(2)**2+rdiff(3)**2)**0.5
3474      end do
3475    end do
3476  end do
3477 
3478 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

SOURCE

5337 subroutine dymfz9(dynmat,natom,nqpt,gprim,option,spqpt,trans)
5338 
5339 !Arguments -------------------------------
5340 !scalars
5341  integer,intent(in) :: natom,nqpt,option
5342 !arrays
5343  real(dp),intent(in) :: gprim(3,3),spqpt(3,nqpt),trans(3,natom)
5344  real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)
5345 
5346 !Local variables -------------------------
5347 !scalars
5348  integer :: ia,ib,iqpt,mu,nu
5349  real(dp) :: im,ktrans,re
5350 !arrays
5351  real(dp) :: kk(3)
5352 
5353 ! *********************************************************************
5354 
5355  do iqpt=1,nqpt
5356    ! Definition of q in normalized reciprocal space
5357    kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
5358    kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
5359    kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
5360 
5361    if(option==1)then
5362      kk(1)=-kk(1)
5363      kk(2)=-kk(2)
5364      kk(3)=-kk(3)
5365    end if
5366 
5367    do ia=1,natom
5368      do ib=1,natom
5369        ! Product of q with the differences between the two atomic translations
5370        ktrans=kk(1)*(trans(1,ia)-trans(1,ib))+kk(2)*(trans(2,ia)-trans(2,ib))+kk(3)*(trans(3,ia)-trans(3,ib))
5371        do mu=1,3
5372          do nu=1,3
5373            re=dynmat(1,mu,ia,nu,ib,iqpt)
5374            im=dynmat(2,mu,ia,nu,ib,iqpt)
5375            ! Transformation of the Old dynamical matrices by New ones by multiplication by a phase shift
5376            dynmat(1,mu,ia,nu,ib,iqpt)=re*cos(two_pi*ktrans)-im*sin(two_pi*ktrans)
5377            dynmat(2,mu,ia,nu,ib,iqpt)=re*sin(two_pi*ktrans)+im*cos(two_pi*ktrans)
5378          end do
5379        end do
5380      end do
5381    end do
5382  end do
5383 
5384 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)!!)

SOURCE

3716 subroutine dynmat_dq(qpt,natom,gprim,nrpt,rpt,atmfrc,wghatm,dddq)
3717 
3718 !Arguments -------------------------------
3719 !scalars
3720  integer,intent(in) :: natom,nrpt
3721 !arrays
3722  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt(3)
3723  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
3724  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
3725  real(dp),intent(out) :: dddq(2,3,natom,3,natom,3)
3726 
3727 !Local variables -------------------------
3728 !scalars
3729  integer :: ia,ib,irpt,mu,nu,ii
3730  real(dp) :: im,kr,re
3731 !arrays
3732  real(dp) :: kk(3),fact(2,3)
3733 
3734 ! *********************************************************************
3735 
3736  dddq = zero
3737 
3738  do irpt=1,nrpt
3739    ! Calculation of the k coordinates in Normalized Reciprocal coordinates
3740    kk(1) = qpt(1)*gprim(1,1)+ qpt(2)*gprim(1,2) + qpt(3)*gprim(1,3)
3741    kk(2) = qpt(1)*gprim(2,1)+ qpt(2)*gprim(2,2) + qpt(3)*gprim(2,3)
3742    kk(3) = qpt(1)*gprim(3,1)+ qpt(2)*gprim(3,2) + qpt(3)*gprim(3,3)
3743 
3744    ! Product of k and r
3745    kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3746 
3747    ! Get phase factor
3748    re=cos(two_pi*kr); im=sin(two_pi*kr)
3749 
3750    ! Inner loop on atoms and directions
3751    do ib=1,natom
3752      do ia=1,natom
3753        if (abs(wghatm(ia,ib,irpt))>1.0d-10) then
3754          ! take into account rotation due to i.
3755          fact(1,:) = -im * wghatm(ia,ib,irpt) * rpt(:,irpt)
3756          fact(2,:) =  re * wghatm(ia,ib,irpt) * rpt(:,irpt)
3757          do nu=1,3
3758            do mu=1,3
3759              ! Real and imaginary part of the dynamical matrices
3760              ! Atmfrc should be real
3761              do ii=1,3
3762                dddq(1,mu,ia,nu,ib,ii) = dddq(1,mu,ia,nu,ib,ii) + fact(1,ii) * atmfrc(mu,ia,nu,ib,irpt)
3763                dddq(2,mu,ia,nu,ib,ii) = dddq(2,mu,ia,nu,ib,ii) + fact(2,ii) * atmfrc(mu,ia,nu,ib,irpt)
3764              end do
3765            end do
3766          end do
3767        end if
3768      end do
3769    end do
3770  end do
3771 
3772 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

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

6350 subroutine ftgam (wghatm,gam_qpt,gam_rpt,natom,nqpt,nrpt,qtor,coskr, sinkr)
6351 
6352 !Arguments -------------------------------
6353 !scalars
6354  integer,intent(in) :: natom,nqpt,nrpt,qtor
6355 !arrays
6356  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
6357  real(dp),intent(inout) :: gam_qpt(2,3*natom*3*natom,nqpt)
6358  real(dp),intent(inout) :: gam_rpt(2,3*natom*3*natom,nrpt)
6359  real(dp),intent(in) :: coskr(nqpt,nrpt)
6360  real(dp),intent(in) :: sinkr(nqpt,nrpt)
6361 
6362 !Local variables -------------------------
6363 !scalars
6364  integer :: iatom,idir,ip,iqpt,irpt,jatom,jdir
6365  real(dp) :: im,re
6366  character(len=500) :: msg
6367 
6368 ! *********************************************************************
6369 
6370  select case (qtor)
6371  case (1)
6372    ! Recip to real space
6373    gam_rpt(:,:,:) = zero
6374    do irpt=1,nrpt
6375      do iqpt=1,nqpt
6376        ! Get the phase factor with normalization!
6377        re=coskr(iqpt,irpt)
6378        im=sinkr(iqpt,irpt)
6379        do ip=1,3*natom*3*natom
6380          ! Real and imaginary part of the real-space gam matrices
6381          gam_rpt(1,ip,irpt) = gam_rpt(1,ip,irpt) + re*gam_qpt(1,ip,iqpt) + im*gam_qpt(2,ip,iqpt)
6382          gam_rpt(2,ip,irpt) = gam_rpt(2,ip,irpt) + re*gam_qpt(2,ip,iqpt) - im*gam_qpt(1,ip,iqpt)
6383        end do
6384      end do
6385    end do
6386    gam_rpt = gam_rpt/nqpt
6387 
6388  case (0)
6389    ! Recip space from real space
6390    gam_qpt(:,:,:)=zero
6391 
6392    do irpt=1,nrpt
6393      do iqpt=1,nqpt
6394 
6395        do iatom=1,natom
6396          do jatom=1,natom
6397            re = coskr(iqpt,irpt)*wghatm(iatom,jatom,irpt)
6398            im = sinkr(iqpt,irpt)*wghatm(iatom,jatom,irpt)
6399 
6400            do idir=1,3
6401              do jdir=1,3
6402                ! Get phase factor
6403 
6404                ip= jdir + (jatom-1)*3 + (idir-1)*3*natom + (iatom-1)*9*natom
6405                ! Real and imaginary part of the interatomic forces
6406                gam_qpt(1,ip,iqpt) = gam_qpt(1,ip,iqpt) + re*gam_rpt(1,ip,irpt) - im*gam_rpt(2,ip,irpt)
6407                !DEBUG
6408                gam_qpt(2,ip,iqpt) = gam_qpt(2,ip,iqpt) + im*gam_rpt(1,ip,irpt) + re*gam_rpt(2,ip,irpt)
6409                !ENDDEBUG
6410              end do ! end jdir
6411            end do ! end idir
6412          end do
6413        end do ! end iatom
6414 
6415      end do ! end iqpt
6416    end do ! end irpt
6417 
6418  case default
6419    write(msg,'(a,i0,a)' )'The only allowed values for qtor are 0 or 1, while qtor= ',qtor,' has been required.'
6420    ABI_BUG(msg)
6421  end select
6422 
6423 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

SOURCE

6448 subroutine ftgam_init (gprim,nqpt,nrpt,qpt_full,rpt,coskr, sinkr)
6449 
6450 !Arguments -------------------------------
6451 !scalars
6452  integer,intent(in) :: nqpt,nrpt
6453 !arrays
6454  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt_full(3,nqpt)
6455  real(dp),intent(out) :: coskr(nqpt,nrpt)
6456  real(dp),intent(out) :: sinkr(nqpt,nrpt)
6457 
6458 !Local variables -------------------------
6459 !scalars
6460  integer :: iqpt,irpt
6461  real(dp) :: kr
6462 !arrays
6463  real(dp) :: kk(3)
6464 
6465 ! *********************************************************************
6466 
6467 ! Prepare the phase factors
6468  do iqpt=1,nqpt
6469    ! Calculation of the k coordinates in Normalized Reciprocal coordinates
6470    kk(1) = qpt_full(1,iqpt)*gprim(1,1) + qpt_full(2,iqpt)*gprim(1,2) + qpt_full(3,iqpt)*gprim(1,3)
6471    kk(2) = qpt_full(1,iqpt)*gprim(2,1) + qpt_full(2,iqpt)*gprim(2,2) + qpt_full(3,iqpt)*gprim(2,3)
6472    kk(3) = qpt_full(1,iqpt)*gprim(3,1) + qpt_full(2,iqpt)*gprim(3,2) + qpt_full(3,iqpt)*gprim(3,3)
6473    do irpt=1,nrpt
6474      ! Product of k and r
6475      kr = kk(1)*rpt(1,irpt)+ kk(2)*rpt(2,irpt)+ kk(3)*rpt(3,irpt)
6476      coskr(iqpt,irpt)=cos(two_pi*kr)
6477      sinkr(iqpt,irpt)=sin(two_pi*kr)
6478    end do
6479  end do
6480 
6481 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
 comm=MPI communicator.

OUTPUT

 atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space.

SOURCE

3508 subroutine ftifc_q2r(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt,comm)
3509 
3510 !Arguments -------------------------------
3511 !scalars
3512  integer,intent(in) :: natom,nqpt,nrpt,comm
3513 !arrays
3514  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
3515  real(dp),intent(out) :: atmfrc(3,natom,3,natom,nrpt)
3516  real(dp),intent(in) :: dynmat(2,3,natom,3,natom,nqpt)
3517 
3518 !Local variables -------------------------
3519 !scalars
3520  integer :: ia,ib,iqpt,irpt,mu,nu,nprocs,my_rank,ierr
3521  real(dp) :: im,kr,re
3522 !arrays
3523  real(dp) :: kk(3)
3524 
3525 ! *********************************************************************
3526 
3527  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3528 
3529  ! Interatomic Forces from Dynamical Matrices
3530  atmfrc = zero
3531  do irpt=1,nrpt
3532    if (mod(irpt, nprocs) /= my_rank) cycle ! mpi-parallelism
3533    do iqpt=1,nqpt
3534 
3535      ! Calculation of the k coordinates in Normalized Reciprocal coordinates
3536      kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
3537      kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
3538      kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
3539 
3540      ! Product of k and r
3541      kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3542 
3543      ! Get the phase factor
3544      re=cos(two_pi*kr)
3545      im=sin(two_pi*kr)
3546 
3547      ! Now, big inner loops on atoms and directions
3548      ! The indices are ordered to give better speed
3549      do ib=1,natom
3550        do nu=1,3
3551          do ia=1,natom
3552            do mu=1,3
3553              ! Real and imaginary part of the interatomic forces
3554              atmfrc(mu,ia,nu,ib,irpt)=atmfrc(mu,ia,nu,ib,irpt) &
3555               +re*dynmat(1,mu,ia,nu,ib,iqpt)&
3556               +im*dynmat(2,mu,ia,nu,ib,iqpt)
3557              !The imaginary part should be equal to zero !!!!!!
3558              !atmfrc(2,mu,ia,nu,ib,irpt)=atmfrc(2,mu,ia,nu,ib,irpt) &
3559              !          +re*dynmat(2,mu,ia,nu,ib,iqpt) &
3560              !          -im*dynmat(1,mu,ia,nu,ib,iqpt)
3561            end do
3562          end do
3563        end do
3564      end do
3565 
3566    end do
3567  end do
3568 
3569  call xmpi_sum(atmfrc, comm, ierr)
3570  !The sumifc has to be weighted by a normalization factor of 1/nqpt
3571  atmfrc = atmfrc/nqpt
3572 
3573 end subroutine ftifc_q2r

m_dynmat/ftifc_r2q [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 ftifc_r2q

FUNCTION

 Generates the Fourier transform of the interatomic forces
 to obtain dynamical matrices in reciprocal space: R --> q.

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
 comm: MPI communicator

OUTPUT

 dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base

SOURCE

3604 subroutine ftifc_r2q(atmfrc, dynmat, gprim, natom, nqpt, nrpt, rpt, spqpt, wghatm, comm)
3605 
3606 !Arguments -------------------------------
3607 !scalars
3608  integer,intent(in) :: natom,nqpt,nrpt,comm
3609 !arrays
3610  real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt)
3611  real(dp),intent(in) :: wghatm(natom,natom,nrpt)
3612  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
3613  real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt)
3614 
3615 !Local variables -------------------------
3616 !scalars
3617  integer :: ia,ib,iqpt,irpt,mu,nu,cnt,my_rank,nprocs, ierr
3618  real(dp) :: facti,factr,im,kr,re
3619  !real(dp) : w(2, natom, natom)
3620 !arrays
3621  real(dp) :: kk(3)
3622 
3623 ! *********************************************************************
3624 
3625  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
3626  dynmat = zero; cnt = 0
3627 
3628  ! MG: This is an hotspot. I don'tknow whether one should rewrite with BLAS1 dot or not.
3629  ! Note, however, that simply removing the check on the weights inside the loop over atoms.
3630  ! leads to a non-negligible speedup with intel (~30% if dipdip -1 is used)
3631  do iqpt=1,nqpt
3632 
3633    ! Calculation of the k coordinates in Normalized Reciprocal coordinates
3634    kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)* gprim(1,2)+spqpt(3,iqpt)*gprim(1,3)
3635    kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)* gprim(2,2)+spqpt(3,iqpt)*gprim(2,3)
3636    kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)* gprim(3,2)+spqpt(3,iqpt)*gprim(3,3)
3637 
3638    do irpt=1,nrpt
3639      cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism.
3640 
3641      ! k.R
3642      kr = kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt)
3643      ! Get phase factor
3644      re = cos(two_pi*kr); im = sin(two_pi*kr)
3645 
3646      ! Inner loop on atoms and directions
3647      do ib=1,natom
3648        do ia=1,natom
3649          !if (abs(wghatm(ia,ib,irpt)) > tol10) then  ! Commented by MG
3650            factr = re * wghatm(ia,ib,irpt)
3651            facti = im * wghatm(ia,ib,irpt)
3652            do nu=1,3
3653              do mu=1,3
3654                ! Real and imaginary part of the dynamical matrices
3655                ! Atmfrc should be real
3656                dynmat(1,mu,ia,nu,ib,iqpt) = dynmat(1,mu,ia,nu,ib,iqpt) + factr * atmfrc(mu,ia,nu,ib,irpt)
3657                dynmat(2,mu,ia,nu,ib,iqpt) = dynmat(2,mu,ia,nu,ib,iqpt) + facti * atmfrc(mu,ia,nu,ib,irpt)
3658              end do
3659            end do
3660          !end if
3661        end do
3662      end do
3663 
3664      ! MG: New version: I don't know if it's faster.
3665      !w(1,:,:) = re * wghatm(:,:,irpt)
3666      !w(2,:,:) = im * wghatm(:,:,irpt)
3667      !do ib=1,natom
3668      !  do nu=1,3
3669      !    do ia=1,natom
3670      !      do mu=1,3
3671      !        ! Real and imaginary part of the dynamical matrices
3672      !        ! Atmfrc should be real
3673      !        dynmat(1,mu,ia,nu,ib,iqpt) = dynmat(1,mu,ia,nu,ib,iqpt) + w(1,ia,ib) * atmfrc(mu,ia,nu,ib,irpt)
3674      !        dynmat(2,mu,ia,nu,ib,iqpt) = dynmat(2,mu,ia,nu,ib,iqpt) + w(2,ia,ib) * atmfrc(mu,ia,nu,ib,irpt)
3675      !     end do
3676      !    end do
3677      !  end do
3678      !end do
3679 
3680    end do
3681  end do
3682 
3683  if (nprocs > 1) call xmpi_sum(dynmat, comm, ierr)
3684 
3685 end subroutine ftifc_r2q

m_dynmat/get_bigbox_and_weights [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 get_bigbox_and_weights

FUNCTION

 Compute the Big Box containing the R points in the cartesian real space needed to Fourier Transform
 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.)
 natom= Number of atoms
 nqbz= Number of q-points in BZ.
 ngqpt(3)= Numbers used to generate the q points to sample the Brillouin zone using an homogeneous grid
 nqshift= number of shifts in q-mesh
 qshift(3, nqshift) = Q-mesh shifts
 rprim(3,3)= Normalized coordinates in real space.
 rprimd, gprimd
 rcan(3,natom)  = Atomic position in canonical coordinates
 cutmode=Define the cutoff used to filter the output R-points according to their weights.
   0 --> No cutoff (mainly for debugging)
   1 --> Include only those R-points for which sum(abs(wg(:,:,irpt)) < tol20
         This is the approach used for the dynamical matrix.
   2 --> Include only those R-points for which the trace over iatom of abs(wg(iat,iat,irpt)) < tol20
         This option is used for objects that depend on a single atomic index.
 comm= MPI communicator

OUTPUT

 nrpt= Total Number of R points in the Big Box
 cell(3,nrpt) Give the index of the the cell and irpt
 rpt(3,nrpt)= Canonical coordinates of the R points in the unit cell. These coordinates are normalized (=> * acell(3)!!)
 r_inscribed_sphere
 wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector

SOURCE

2672 subroutine get_bigbox_and_weights(brav, natom, nqbz, ngqpt, nqshift, qshift, rprim, rprimd, gprim, rcan, &
2673                                   cutmode, nrpt, rpt, cell, wghatm, r_inscribed_sphere, comm)
2674 
2675 !Arguments -------------------------------
2676 !scalars
2677  integer,intent(in) :: brav, natom, nqbz, nqshift, cutmode, comm
2678  integer,intent(out) :: nrpt
2679  real(dp),intent(out) :: r_inscribed_sphere
2680 !arrays
2681  integer,intent(in) :: ngqpt(3)
2682  real(dp),intent(in) :: gprim(3,3),rprim(3,3),rprimd(3,3), rcan(3, natom)
2683  real(dp),intent(in) :: qshift(3, nqshift)
2684  integer,allocatable,intent(out) :: cell(:,:)
2685  real(dp),allocatable,intent(out) :: rpt(:,:), wghatm(:,:,:)
2686 
2687 !Local variables -------------------------
2688 !scalars
2689  integer :: my_ierr, ierr, ii, irpt, all_nrpt
2690  real(dp) :: toldist
2691  integer :: ngqpt9(9)
2692  character(len=500*4) :: msg
2693 !arrays
2694  integer,allocatable :: all_cell(:,:)
2695  real(dp),allocatable :: all_rpt(:,:), all_wghatm(:,:,:)
2696 
2697 ! *********************************************************************
2698 
2699  ABI_CHECK(any(cutmode == [0, 1, 2]), "cutmode should be in [0, 1, 2]")
2700 
2701  ! Create the Big Box of R vectors in real space and compute the number of points (cells) in real space
2702  call make_bigbox(brav, all_cell, ngqpt, nqshift, rprim, all_nrpt, all_rpt)
2703 
2704  ! Weights associated to these R points and to atomic pairs
2705  ABI_MALLOC(all_wghatm, (natom, natom, all_nrpt))
2706 
2707  ! HM: this tolerance is highly dependent on the compilation/architecture
2708  !     numeric errors in the DDB text file. Try a few tolerances and check whether all the weights are found.
2709  ngqpt9 = 0; ngqpt9(1:3) = ngqpt(1:3)
2710  toldist = tol8
2711  do while (toldist <= tol6)
2712    ! Note ngqpt(9) with intent(inout)!
2713    call wght9(brav, gprim, natom, ngqpt9, nqbz, nqshift, all_nrpt, qshift, rcan, &
2714               all_rpt, rprimd, toldist, r_inscribed_sphere, all_wghatm, my_ierr)
2715    call xmpi_max(my_ierr, ierr, comm, ii)
2716    if (ierr > 0) toldist = toldist * 10
2717    if (ierr == 0) exit
2718  end do
2719 
2720  if (ierr > 0) then
2721    write(msg, '(3a,es14.4,2a,i0, 14a)' ) &
2722     'The sum of the weight is not equal to nqpt.',ch10,&
2723     'The sum of the weights is: ',sum(all_wghatm),ch10,&
2724     'The number of q points is: ',nqbz, ch10, &
2725     'This might have several sources.',ch10,&
2726     'If toldist is larger than 1.0e-8, the atom positions might be loose.',ch10,&
2727     'and the q point weights not computed properly.',ch10,&
2728     'Action: make input atomic positions more symmetric.',ch10,&
2729     'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine and recompile.',ch10,&
2730     'Actually, this can also happen when ngqpt is 0 0 0,',ch10,&
2731     'if abs(brav) /= 1, in this case you should change brav to 1. If brav is already set to 1 (default) try -1.'
2732    ABI_ERROR(msg)
2733  end if
2734 
2735  ! Only conserve the necessary points in rpt.
2736  nrpt = 0
2737  do irpt=1,all_nrpt
2738    if (filterw(all_wghatm(:,:,irpt))) cycle
2739    nrpt = nrpt + 1
2740  end do
2741 
2742  ! Allocate output arrays and transfer data.
2743  ABI_MALLOC(rpt, (3, nrpt))
2744  ABI_MALLOC(cell, (3, nrpt))
2745  ABI_MALLOC(wghatm, (natom, natom, nrpt))
2746 
2747  ii = 0
2748  do irpt=1,all_nrpt
2749    if (filterw(all_wghatm(:,:,irpt))) cycle
2750    ii = ii + 1
2751    rpt(:, ii) = all_rpt(:,irpt)
2752    wghatm(:,:,ii) = all_wghatm(:,:,irpt)
2753    cell(:,ii) = all_cell(:,irpt)
2754  end do
2755 
2756  ABI_FREE(all_rpt)
2757  ABI_FREE(all_wghatm)
2758  ABI_FREE(all_cell)
2759 
2760 contains
2761 
2762 logical pure function filterw(wg)
2763 
2764  real(dp),intent(in) :: wg(natom,natom)
2765  integer :: iat
2766  real(dp) :: trace
2767 
2768  select case (cutmode)
2769  case (1)
2770    filterw = sum(abs(wg)) < tol20
2771  case (2)
2772    trace = zero
2773    do iat=1,natom
2774      trace = trace + abs(wg(iat,iat))
2775    end do
2776    filterw = trace < tol20
2777  case default
2778    filterw = .False.
2779  end select
2780 
2781 end function filterw
2782 
2783 end subroutine get_bigbox_and_weights

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
 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
 comm=MPI communicator.
 [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part
 [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part

OUTPUT

 d2cart(2,3,mpert,3,mpert)=dynamical matrix obtained for the wavevector qpt (normalized using qphnrm)

SOURCE

5510 subroutine gtdyn9(acell,atmfrc,dielt,dipdip,dyewq0,d2cart,gmet,gprim,mpert,natom,&
5511                   nrpt,qphnrm,qpt,rmet,rprim,rpt,trans,ucvol,wghatm,xred,zeff,qdrp_cart,ewald_option,comm,&
5512                   dipquad,quadquad)  ! optional
5513 
5514 !Arguments -------------------------------
5515 !scalars
5516  integer,intent(in) :: dipdip,mpert,natom,nrpt,ewald_option,comm
5517  real(dp),intent(in) :: qphnrm,ucvol
5518  integer,optional,intent(in) :: dipquad, quadquad
5519 !arrays
5520  real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qpt(3)
5521  real(dp),intent(in) :: rmet(3,3),rprim(3,3),rpt(3,nrpt)
5522  real(dp),intent(in) :: trans(3,natom),wghatm(natom,natom,nrpt),xred(3,natom)
5523  real(dp),intent(in) :: zeff(3,3,natom)
5524  real(dp),intent(in) :: qdrp_cart(3,3,3,natom)
5525  real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt)
5526  real(dp),intent(in) :: dyewq0(3,3,natom)
5527  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
5528 
5529 !Local variables -------------------------
5530 !scalars
5531  integer,parameter :: nqpt1 = 1, option2 = 2, sumg0 = 0, plus1 = 1, iqpt1 = 1
5532  integer :: i1, i2, ib, nsize, dipquad_, quadquad_
5533 !arrays
5534  real(dp) :: qphon(3) !, tsec(2)
5535  real(dp),allocatable :: dq(:,:,:,:,:),dyew(:,:,:,:,:)
5536 
5537 ! *********************************************************************
5538 
5539  ! Keep track of time spent in gtdyn9
5540  !call timab(1750, 1, tsec)
5541 
5542  ABI_MALLOC(dq,(2,3,natom,3,natom))
5543 
5544  ! Define quadrupolar options
5545  dipquad_=0; if(present(dipquad)) dipquad_=dipquad
5546  quadquad_=0; if(present(quadquad)) quadquad_=quadquad
5547 
5548  ! Get the normalized wavevector
5549  if(abs(qphnrm)<1.0d-7)then
5550    qphon(1:3)=zero
5551  else
5552    qphon(1:3)=qpt(1:3)/qphnrm
5553  end if
5554 
5555  ! Generate the analytical part from the interatomic forces
5556  call ftifc_r2q(atmfrc, dq, gprim, natom, nqpt1, nrpt, rpt, qphon, wghatm, comm)
5557 
5558  ! The analytical dynamical matrix dq has been generated
5559  ! in the normalized canonical coordinate system.
5560  ! Now, the phase is modified, in order to recover the usual (xred) coordinate of atoms.
5561  call dymfz9(dq,natom,nqpt1,gprim,option2,qphon,trans)
5562 
5563  if (dipdip==1.or.dipquad_==1.or.quadquad_==1) then
5564    ! Add the non-analytical part
5565    ! Compute dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
5566    ! second energy derivative wrt xred(3,natom) in Hartrees (Denoted A-bar in the notes)
5567    ABI_MALLOC(dyew,(2,3,natom,3,natom))
5568 
5569    call ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff,&
5570       qdrp_cart,option=ewald_option,dipquad=dipquad_,quadquad=quadquad_)
5571 
5572    call q0dy3_apply(natom,dyewq0,dyew)
5573    call nanal9(dyew,dq,iqpt1,natom,nqpt1,plus1)
5574 
5575    ABI_FREE(dyew)
5576  end if
5577 
5578  ! Copy the dynamical matrix in the proper location
5579  ! First zero all the elements
5580  nsize=2*(3*mpert)**2
5581  d2cart = zero
5582 
5583  ! Copy the elements from dq to d2cart
5584  d2cart(:,:,1:natom,:,1:natom)=dq(:,:,1:natom,:,1:natom)
5585 
5586  ! In case we have the gamma point,
5587  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)then
5588    ! Copy the effective charge and dielectric constant in the final array
5589    do i1=1,3
5590      do i2=1,3
5591        d2cart(1,i1,natom+2,i2,natom+2)=dielt(i1,i2)
5592        do ib=1,natom
5593          d2cart(1,i1,natom+2,i2,ib)=zeff(i1,i2,ib)
5594          d2cart(1,i2,ib,i1,natom+2)=zeff(i1,i2,ib)
5595        end do
5596      end do
5597    end do
5598  end if
5599 
5600  ABI_FREE(dq)
5601 
5602  !call timab(1750, 2, tsec)
5603 
5604 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

SOURCE

3797 subroutine ifclo9(ifccar,ifcloc,vect1,vect2,vect3)
3798 
3799 !Arguments -------------------------------
3800 !arrays
3801  real(dp),intent(in) :: ifccar(3,3),vect1(3),vect2(3),vect3(3)
3802  real(dp),intent(out) :: ifcloc(3,3)
3803 
3804 !Local variables -------------------------
3805 !scalars
3806  integer :: ii,jj
3807 !arrays
3808  real(dp) :: work(3,3)
3809 
3810 ! *********************************************************************
3811 
3812  do jj=1,3
3813    do ii=1,3
3814      work(jj,ii)=zero
3815    end do
3816    do ii=1,3
3817      work(jj,1)=work(jj,1)+ifccar(jj,ii)*vect1(ii)
3818      work(jj,2)=work(jj,2)+ifccar(jj,ii)*vect2(ii)
3819      work(jj,3)=work(jj,3)+ifccar(jj,ii)*vect3(ii)
3820    end do
3821  end do
3822 
3823  do jj=1,3
3824    do ii=1,3
3825      ifcloc(ii,jj)=zero
3826    end do
3827    do ii=1,3
3828      ifcloc(1,jj)=ifcloc(1,jj)+vect1(ii)*work(ii,jj)
3829      ifcloc(2,jj)=ifcloc(2,jj)+vect2(ii)*work(ii,jj)
3830      ifcloc(3,jj)=ifcloc(3,jj)+vect3(ii)*work(ii,jj)
3831    end do
3832  end do
3833 
3834 end subroutine ifclo9

m_dynmat/make_bigbox [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 make_bigbox

FUNCTION

 Helper functions to faciliate 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(3,nrpt)= integer coordinates of the cells (R points) in the rprim basis
 nprt= Number of cells (R points) in the Big Box
 rpt(3,mrpt)= canonical coordinates of the cells (R points)
  These coordinates are normalized (=> * acell(3)!!)
  The array is allocated here with the proper dimension. Client code is responsible
  for the deallocation.

SOURCE

2817 subroutine make_bigbox(brav, cell, ngqpt, nqshft, rprim, nrpt, rpt)
2818 
2819 !Arguments -------------------------------
2820 !scalars
2821  integer,intent(in) :: brav,nqshft
2822  integer,intent(out) :: nrpt
2823 !arrays
2824  integer,intent(in) :: ngqpt(3)
2825  real(dp),intent(in) :: rprim(3,3)
2826  real(dp),allocatable,intent(out) :: rpt(:,:)
2827  integer,allocatable,intent(out) :: cell(:,:)
2828 
2829 !Local variables -------------------------
2830 !scalars
2831  integer :: choice,mrpt
2832 !arrays
2833  real(dp) :: dummy_rpt(3,1)
2834  integer:: dummy_cell(1,3)
2835 
2836 ! *********************************************************************
2837 
2838  ! Compute the number of points (cells) in real space
2839  choice=0
2840  call bigbx9(brav,dummy_cell,choice,1,ngqpt,nqshft,mrpt,rprim,dummy_rpt)
2841 
2842  ! Now we can allocate and calculate the points and the weights.
2843  nrpt = mrpt
2844  ABI_MALLOC(rpt,(3,nrpt))
2845  ABI_MALLOC(cell,(3,nrpt))
2846 
2847  choice=1
2848  call bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt)
2849 
2850 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,...)
  [herm_opt]= 1 to hermitianize mat (default)
          0 if no symmetrization should be performed

SIDE EFFECTS

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

SOURCE

6218 subroutine massmult_and_breaksym(natom, ntypat, typat, amu, mat, &
6219                                  herm_opt) ! optional
6220 
6221 !Arguments -------------------------------
6222 !scalars
6223  integer,intent(in) :: natom,ntypat
6224  integer,optional,intent(in) :: herm_opt
6225 !arrays
6226  integer,intent(in) :: typat(natom)
6227  real(dp),intent(in) :: amu(ntypat)
6228  real(dp),intent(inout) :: mat(2*3*natom*3*natom)
6229 
6230 !Local variables -------------------------
6231 !scalars
6232  integer :: i1,i2,idir1,idir2,index,ipert1,ipert2, herm_opt__
6233  real(dp),parameter :: break_symm=1.0d-12
6234  !real(dp),parameter :: break_symm=zero
6235  real(dp) :: fac
6236 !arrays
6237  real(dp) :: nearidentity(3,3)
6238 
6239 ! *********************************************************************
6240 
6241  herm_opt__ = 1; if (present(herm_opt)) herm_opt__ = herm_opt
6242 
6243  ! This slight breaking of the symmetry allows the results to be more portable between machines
6244  nearidentity(:,:)=one
6245  nearidentity(1,1)=one+break_symm
6246  nearidentity(3,3)=one-break_symm
6247 
6248  ! Include the masses in the dynamical matrix
6249  do ipert1=1,natom
6250    do ipert2=1,natom
6251      fac=1.0_dp/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass
6252      do idir1=1,3
6253        do idir2=1,3
6254          i1=idir1+(ipert1-1)*3
6255          i2=idir2+(ipert2-1)*3
6256          index=i1+3*natom*(i2-1)
6257          mat(2*index-1)=mat(2*index-1)*fac*nearidentity(idir1,idir2)
6258          mat(2*index  )=mat(2*index  )*fac*nearidentity(idir1,idir2)
6259          ! This is to break slightly the translation invariance, and make the automatic tests more portable
6260          if(ipert1==ipert2 .and. idir1==idir2)then
6261            mat(2*index-1)=mat(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp
6262          end if
6263        end do
6264      end do
6265    end do
6266  end do
6267 
6268  ! Make the dynamical matrix hermitian
6269  if (herm_opt__ == 1) call mkherm(mat,3*natom)
6270 
6271 end subroutine massmult_and_breaksym

m_dynmat/massmult_and_breaksym_cplx [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

  mult_masses_and_break_symms_cplx

FUNCTION

  Similar to massmult_and_breaksym, the only difference is that it receives complex array.

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.

 For plus=0, see Eq.(76) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]],
 get the left hand side.

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

SOURCE

5415 subroutine nanal9(dyew,dynmat,iqpt,natom,nqpt,plus)
5416 
5417 !Arguments -------------------------------
5418 !scalars
5419  integer,intent(in) :: iqpt,natom,nqpt,plus
5420 !arrays
5421  real(dp),intent(in) :: dyew(2,3,natom,3,natom)
5422  real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt)
5423 
5424 !Local variables -------------------------
5425 !scalars
5426  integer :: ia,ib,mu,nu
5427  character(len=500) :: msg
5428 
5429 ! *********************************************************************
5430 
5431  if (plus==0) then
5432 
5433    do ia=1,natom
5434      do ib=1,natom
5435        do mu=1,3
5436          do nu=1,3
5437            ! The following four lines are OK
5438            dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) - dyew(1,mu,ia,nu,ib)
5439            dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) - dyew(2,mu,ia,nu,ib)
5440          end do
5441        end do
5442      end do
5443    end do
5444 
5445  else if (plus==1) then
5446    do ia=1,natom
5447      do ib=1,natom
5448        do mu=1,3
5449          do nu=1,3
5450            dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) + dyew(1,mu,ia,nu,ib)
5451            dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) + dyew(2,mu,ia,nu,ib)
5452          end do
5453        end do
5454      end do
5455    end do
5456 
5457  else
5458    write(msg,'(3a,i0,a)' )&
5459     'The argument "plus" must be equal to 0 or 1.',ch10,&
5460     'The value ',plus,' is not available.'
5461    ABI_BUG(msg)
5462  end if
5463 
5464 end subroutine nanal9

m_dynmat/phdispl_from_eigvec [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 phdispl_from_eigvec

FUNCTION

  Phonon displacements in cart coords from eigenvectors

INPUTS

  natom: number of atoms in unit cell
  ntypat=number of atom types
  typat(natom)=integer label of each type of atom (1,2,...)
  amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms)
  eigvec(2*3*natom*3*natom)= eigenvectors of the dynamical matrix in cartesian coordinates.

OUTPUT

  displ(2*3*natom*3*natom)=displacements of atoms in cartesian coordinates.

SOURCE

5946 pure subroutine phdispl_from_eigvec(natom, ntypat, typat, amu, eigvec, displ)
5947 
5948 !Arguments -------------------------------
5949 !scalars
5950  integer,intent(in) :: natom, ntypat
5951 !arrays
5952  integer,intent(in) :: typat(natom)
5953  real(dp),intent(in) :: amu(ntypat)
5954  real(dp),intent(in) :: eigvec(2*3*natom*3*natom)
5955  real(dp),intent(out) :: displ(2*3*natom*3*natom)
5956 
5957 !Local variables -------------------------
5958 !scalars
5959  integer :: i1,idir1,imode,ipert1, index
5960 
5961 ! *********************************************************************
5962 
5963  do imode=1,3*natom
5964 
5965    do idir1=1,3
5966      do ipert1=1,natom
5967        i1=idir1+(ipert1-1)*3
5968        index=i1+3*natom*(imode-1)
5969        displ(2*index-1)=eigvec(2*index-1) / sqrt(amu(typat(ipert1))*amu_emass)
5970        displ(2*index  )=eigvec(2*index  ) / sqrt(amu(typat(ipert1))*amu_emass)
5971      end do
5972    end do
5973 
5974  end do
5975 
5976 end subroutine phdispl_from_eigvec

m_dynmat/pheigvec_normalize [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 pheigvec_normalize

FUNCTION

  Normalize input eigenvectors in cartesian coordinates

INPUTS

  natom: number of atoms in unit cell

SIDE EFFECTS

  eigvec(2*3*natom*3*natom)=in output the normalized eigenvectors in cartesian coordinates.

SOURCE

5884 pure subroutine pheigvec_normalize(natom, eigvec)
5885 
5886 !Arguments -------------------------------
5887 !scalars
5888  integer,intent(in) :: natom
5889 !arrays
5890  real(dp),intent(inout) :: eigvec(2*3*natom*3*natom)
5891 
5892 !Local variables -------------------------
5893 !scalars
5894  integer :: i1,idir1,imode,ipert1,index
5895  real(dp) :: norm
5896 
5897 ! *********************************************************************
5898 
5899  do imode=1,3*natom
5900 
5901    norm=zero
5902    do idir1=1,3
5903      do ipert1=1,natom
5904        i1=idir1+(ipert1-1)*3
5905        index=i1+3*natom*(imode-1)
5906        norm=norm+eigvec(2*index-1)**2+eigvec(2*index)**2
5907      end do
5908    end do
5909    norm=sqrt(norm)
5910 
5911    do idir1=1,3
5912      do ipert1=1,natom
5913        i1=idir1+(ipert1-1)*3
5914        index=i1+3*natom*(imode-1)
5915        eigvec(2*index-1)=eigvec(2*index-1)/norm
5916        eigvec(2*index)=eigvec(2*index)/norm
5917      end do
5918    end do
5919 
5920  end do
5921 
5922 end subroutine pheigvec_normalize

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
 See Eq.(71) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]],
 get the left hand side.

INPUTS

  dyewq0(3,3,natom) = part needed to correct the dynamical matrix for atom self-interaction.
  natom= number of atom in the unit cell

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

SOURCE

1879 subroutine q0dy3_apply(natom,dyewq0,dyew)
1880 
1881 !Arguments -------------------------------
1882 !scalars
1883  integer,intent(in) :: natom
1884 !arrays
1885  real(dp),intent(in) :: dyewq0(3,3,natom)
1886  real(dp),intent(inout) :: dyew(2,3,natom,3,natom)
1887 
1888 !Local variables -------------------------
1889 !scalars
1890  integer :: ia,mu,nu
1891 
1892 ! *********************************************************************
1893 
1894  do mu=1,3
1895    do nu=1,3
1896      do ia=1,natom
1897        dyew(1,mu,ia,nu,ia)=dyew(1,mu,ia,nu,ia)-dyewq0(mu,nu,ia)
1898      end do
1899    end do
1900  end do
1901 
1902 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
 See Eq.(71) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], the sum over \kappa"

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

SOURCE

1949 subroutine q0dy3_calc(natom,dyewq0,dyew,option)
1950 
1951 !Arguments -------------------------------
1952 !scalars
1953  integer,intent(in) :: natom,option
1954 !arrays
1955  real(dp),intent(in) :: dyew(2,3,natom,3,natom)
1956  real(dp),intent(out) :: dyewq0(3,3,natom)
1957 
1958 !Local variables -------------------------
1959 !scalars
1960  integer :: ia,ib,mu,nu
1961  character(len=500) :: msg
1962 
1963 ! *********************************************************************
1964 
1965  if(option==1.or.option==2)then
1966    do mu=1,3
1967      do nu=1,3
1968        do ia=1,natom
1969          dyewq0(mu,nu,ia)=zero
1970          do ib=1,natom
1971            dyewq0(mu,nu,ia)=dyewq0(mu,nu,ia)+dyew(1,mu,ia,nu,ib)
1972          end do
1973        end do
1974      end do
1975    end do
1976  else
1977    write (msg, '(3a)')&
1978 &   'option should be 1 or 2.',ch10,&
1979 &   'action: correct calling routine'
1980    ABI_BUG(msg)
1981  end if
1982 
1983  if(option==2)then
1984    do ia=1,natom
1985      do mu=1,3
1986        do nu=mu,3
1987          dyewq0(mu,nu,ia)=(dyewq0(mu,nu,ia)+dyewq0(nu,mu,ia))/2
1988          dyewq0(nu,mu,ia)=dyewq0(mu,nu,ia)
1989        end do
1990      end do
1991    end do
1992  end if
1993 
1994 end subroutine q0dy3_calc

m_dynmat/sylwtens [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 sylwtens

FUNCTION

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

INPUTS

  has_strain = if .true. i2pert includes strain perturbation
  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 reduced space)
  symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space)
  symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian 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,m_nonlinear,m_respfn_driver

CHILDREN

SOURCE

4980 !subroutine sylwtens(indsym,mpert,natom,nsym,rfpert,symrec,symrel,symrel_cart)
4981 subroutine sylwtens(indsym,mpert,natom,nsym,rfpert,symrec,symrel)
4982 
4983 !Arguments -------------------------------
4984 !scalars
4985  integer,intent(in) :: mpert,natom,nsym
4986 !arrays
4987  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4988  integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert)
4989 ! real(dp),intent(in) :: symrel_cart(3,3,nsym)
4990 
4991 !Local variables -------------------------
4992 !scalars
4993  integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_
4994 ! integer :: i2dir_a,i2dir_b
4995  integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2
4996  integer :: ipesy3,isym
4997 ! integer :: istr,idisy2_a,idisy2_b
4998  logical :: is_strain
4999 ! real(dp) :: flag_dp
5000 !arrays
5001 ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
5002  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
5003  integer,allocatable :: pertsy(:,:,:,:,:,:)
5004 
5005 !***********************************************************************
5006 
5007  ABI_MALLOC(pertsy,(3,mpert,3,mpert,3,mpert))
5008  pertsy(:,:,:,:,:,:) = 0
5009 
5010 !Loop over perturbations
5011 
5012  do i1pert_ = 1, mpert
5013    do i2pert_ = 1, mpert
5014      is_strain=.false.    
5015      do i3pert_ = 1, mpert
5016 
5017        do i1dir_ = 1, 3
5018          do i2dir_ = 1, 3
5019            do i3dir_ = 1, 3
5020 
5021              i1pert = (mpert - i1pert_ + 1)
5022              if (i1pert <= natom) i1pert = natom + 1 - i1pert
5023              i2pert = (mpert - i2pert_ + 1)
5024              if (i2pert <= natom) i2pert = natom + 1 - i2pert
5025              i3pert = (mpert - i3pert_ + 1)
5026              if (i3pert <= natom) i3pert = natom + 1 - i3pert
5027 
5028              if (i1pert <= natom) then
5029                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
5030              else if (i2pert <= natom) then
5031                i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_
5032              else if (i3pert <= natom) then
5033                i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_
5034              else
5035                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
5036              end if
5037 
5038              if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then
5039 
5040 !              Loop over all symmetries
5041 
5042                flag = 0
5043                do isym = 1, nsym
5044 
5045                  found = 1
5046 
5047 !                Select the symmetric element of i1pert,i2pert,i3pert
5048 
5049                  if (i1pert <= natom) then
5050                    ipesy1 = indsym(4,isym,i1pert)
5051                    sym1(:,:) = symrec(:,:,isym)
5052                  else if (i1pert == natom + 2) then
5053                    ipesy1 = i1pert
5054                    sym1(:,:) = symrel(:,:,isym)
5055                  else
5056                    found = 0
5057                  end if
5058 
5059                  if (i2pert <= natom) then
5060                    ipesy2 = indsym(4,isym,i2pert)
5061                    sym2(:,:) = symrec(:,:,isym)
5062                  else if (i2pert == natom + 2) then
5063                    ipesy2 = i2pert
5064                    sym2(:,:) = symrel(:,:,isym)
5065                  else if (i2pert == natom + 3.or. i2pert == natom + 4) then
5066 !                  !TODO: Symmetries on strain perturbation do not work yet.
5067                    found = 0
5068                    is_strain=.true.
5069 !
5070 !                   ipesy2 = i2pert
5071 !                   sym2(:,:) = NINT(symrel_cart(:,:,isym))
5072 !                   if (i2pert==natom+3) istr=i2dir
5073 !                   if (i2pert==natom+4) istr=3+i2dir
5074 !                   i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
5075                  else
5076                    found = 0
5077                  end if
5078 
5079                  if (i3pert == natom + 8) then
5080                    ipesy3 = i3pert
5081                    sym3(:,:) = symrel(:,:,isym)
5082                  else
5083                    found = 0
5084                  end if
5085 
5086 
5087 !                See if the symmetric element is available and check if some
5088 !                of the elements may be zero. In the latter case, they do not need
5089 !                to be computed.
5090 
5091                  if (.not.is_strain) then
5092                    if ((flag /= -1).and.&
5093 &                   (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
5094                      flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir)
5095                    end if
5096 
5097                    do idisy1 = 1, 3
5098                      do idisy2 = 1, 3
5099                        do idisy3 = 1, 3
5100 
5101                          if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.&
5102 &                         (sym3(i3dir,idisy3) /= 0)) then
5103                            if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
5104                              found = 0
5105 !                            exit      ! exit loop over symmetries
5106                            end if
5107                          end if
5108 
5109 
5110                          if ((flag == -1).and.&
5111 &                         ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then
5112                            if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.&
5113 &                           (sym3(i3dir,idisy3)/=0)) then
5114                              flag = 0
5115                            end if
5116                          end if
5117 
5118                        end do
5119                      end do
5120                    end do
5121 !                 else 
5122 !                   if ((flag_dp /= -1).and.&
5123 !&                   (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
5124 !                     flag = sym1(i1dir,i1dir)*sym2(i2dir_a,i2dir_a)* &
5125 !                   & sym2(i2dir_b,i2dir_b)*sym3(i3dir,i3dir)
5126 !                   end if
5127 !
5128 !                   do idisy1 = 1, 3
5129 !                     do idisy2 = 1, 3
5130 !                       if (ipesy2==natom+3) istr=idisy2
5131 !                       if (ipesy2==natom+4) istr=3+idisy2
5132 !                       idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr)
5133 !                       do idisy3 = 1, 3
5134 !
5135 !                         if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir_a,idisy2_a) /= 0).and.&
5136 !&                          (sym2(i2dir_b,idisy2_b) /= 0).and.(sym3(i3dir,idisy3) /= 0)) then
5137 !                           if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
5138 !                             found = 0
5139 !!                            exit      ! exit loop over symmetries
5140 !                           end if
5141 !                         end if
5142 !
5143 !
5144 !                         if ((flag == -1).and.&
5145 !&                         ((idisy1/=i1dir).or.(idisy2_a/=i2dir_a).or.(idisy2_b/=i2dir_b).or.(idisy3/=i3dir))) then
5146 !                           if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir_a,idisy2_a)/=0).and.&
5147 !&                           (sym2(i2dir_b,idisy2_b)/=0).and.(sym3(i3dir,idisy3)/=0)) then
5148 !                             flag = 0
5149 !                           end if
5150 !                         end if
5151 !
5152 !                       end do
5153 !                     end do
5154 !                   end do
5155                  end if
5156 
5157                  if (found == 1) then
5158                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1
5159                  end if
5160 
5161 !                In case a symmetry operation only changes the sign of an
5162 !                element, this element has to be equal to zero
5163 
5164                  if (flag == -1) then
5165                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2
5166                    exit
5167                  end if
5168 
5169                end do    ! close loop on symmetries
5170 
5171 !              If the elemetn i1pert,i2pert,i3pert is not symmetric
5172 !              to a basis element, it is a basis element
5173 
5174                if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then
5175                  pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
5176                end if
5177 
5178              end if ! rfpert /= 0
5179 
5180            end do        ! close loop over perturbations
5181          end do
5182        end do
5183      end do
5184    end do
5185  end do
5186 
5187 !Now, take into account the permutation of (i1pert,i1dir)
5188 !and (i2pert,i2dir)
5189 
5190  do i1pert = 1, mpert
5191    do i2pert = 1, mpert
5192      do i3pert = 1, mpert
5193 
5194        do i1dir = 1, 3
5195          do i2dir = 1, 3
5196            do i3dir = 1, 3
5197 
5198              if ((i1pert /= i2pert).or.(i1dir /= i2dir)) then
5199 
5200                if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.&
5201                 (pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) == 1)) then
5202                  pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = -1
5203                end if
5204 
5205              end if
5206 
5207            end do
5208          end do
5209        end do
5210 
5211      end do
5212    end do
5213  end do
5214 
5215  rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:)
5216 
5217  ABI_FREE(pertsy)
5218 
5219 end subroutine sylwtens

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

SOURCE

2033 subroutine symdyma(dmati,indsym,natom,nsym,qptn,rprimd,symrel,symafm)
2034 
2035 !Arguments -------------------------------
2036 !scalars
2037  integer,intent(in) :: natom,nsym
2038 !arrays
2039  integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
2040  integer,intent(in) :: symafm(nsym)
2041  real(dp),intent(in) :: qptn(3),rprimd(3,3)
2042  real(dp),intent(inout) :: dmati(2*3*natom*3*natom)
2043 
2044 !Local variables -------------------------
2045 !scalars
2046  integer :: i1,i2,iat,idir,ii,index,isgn,isym,itirev,jat,jdir,jj,kk,ll
2047  integer :: niat,njat,timrev
2048  real(dp) :: arg1,arg2,dmint,im,re,sumi,sumr
2049 !arrays
2050  integer :: indij(natom,natom),symq(4,2,nsym),symrec(3,3,nsym)
2051  real(dp) :: TqR(3,3),TqS_(3,3),dynmat(2,3,natom,3,natom)
2052  real(dp) :: dynmatint(2*nsym,2,3,natom,3,natom),gprimd(3,3)
2053  real(dp) :: symcart(3,3,nsym)
2054 
2055 ! *********************************************************************
2056 
2057  ! 0) initializations
2058  call matr3inv(rprimd,gprimd)
2059  do isym=1,nsym
2060    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2061  end do
2062 
2063  TqR=zero
2064  TqS_=zero
2065  dynmat=zero
2066 
2067 !Note: dynmat is used as work space here
2068  i1=0
2069  do iat=1,natom
2070    do idir=1,3
2071      i1=i1+1
2072      i2=0
2073      do jat=1,natom
2074        do jdir=1,3
2075          i2=i2+1
2076          index=i1+3*natom*(i2-1)
2077          dynmat(1,idir,iat,jdir,jat)=dmati(2*index-1)
2078          dynmat(2,idir,iat,jdir,jat)=dmati(2*index  )
2079        end do
2080      end do
2081    end do
2082  end do
2083 
2084 !Transform symrel to cartesian coordinates (RC coding)
2085 !do isym=1,nsym
2086 !symcart(:,:,isym)=matmul(rprimd,matmul(dble(symrel(:,:,isym)),gprimd))
2087 !end do
2088 
2089 !Coding from symdm9
2090  do isym=1,nsym
2091    do jj=1,3
2092      symcart(:,jj,isym)=zero
2093      do kk=1,3
2094        do ll=1,3
2095          symcart(:,jj,isym)=symcart(:,jj,isym)+rprimd(:,kk)*gprimd(jj,ll)*symrel(kk,ll,isym)
2096        end do
2097      end do
2098    end do
2099  end do
2100 
2101  ! Get the symq of the CURRENT Q POINT
2102  ! mjv: set prtvol=0 for production runs.
2103  call littlegroup_q(nsym,qptn,symq,symrec,symafm,timrev,prtvol=0)
2104 
2105  indij(:,:)=0
2106  dynmatint=zero
2107 
2108  do isym=1,nsym  ! loop over all the symmetries
2109    ! write(std_out,*) 'current symmetry',isym
2110    do itirev=1,2  ! loop over the time-reversal symmetry
2111      isgn=3-2*itirev
2112      ! write(std_out,*) 'timereversal',isgn
2113 
2114      if (symq(4,itirev,isym)==1) then ! isym belongs to the wave vector point group
2115        ! write(std_out,*) 'isym belongs to the wave vector point group'
2116        do iat=1,natom
2117          do jat=1,natom
2118            niat=indsym(4,isym,iat)  ! niat={R|t}iat
2119            njat=indsym(4,isym,jat)  ! njat={R|t}jat
2120            indij(niat,njat)=indij(niat,njat)+1
2121            ! write(std_out,'(a,5i5)') 'current status:',iat,jat,niat,njat,indij(niat,njat)
2122            ! phase calculation, arg1 and arg2 because of two-atom derivative
2123            arg1=two_pi*( qptn(1)*indsym(1,isym,iat)+&
2124              qptn(2)*indsym(2,isym,iat)+&
2125              qptn(3)*indsym(3,isym,iat) )
2126            arg2=two_pi*( qptn(1)*indsym(1,isym,jat)+&
2127              qptn(2)*indsym(2,isym,jat)+&
2128              qptn(3)*indsym(3,isym,jat) )
2129 
2130            re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2)
2131            im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2))
2132 
2133            do idir=1,3     ! loop over displacements
2134              do jdir=1,3   ! loop over displacements
2135                ! we pick the (iat,jat) (3x3) block of the dyn.mat.
2136                sumr=zero
2137                sumi=zero
2138                do ii=1,3
2139                  do jj=1,3
2140                    sumr=sumr+symcart(idir,ii,isym)*dynmat(1,ii,niat,jj,njat)*symcart(jdir,jj,isym)
2141                    sumi=sumi+symcart(idir,ii,isym)*dynmat(2,ii,niat,jj,njat)*symcart(jdir,jj,isym)
2142                  end do
2143                end do
2144                sumi=isgn*sumi
2145 
2146                dynmatint(nsym*(itirev-1)+isym,1,idir,iat,jdir,jat)=re*sumr-im*sumi
2147                dynmatint(nsym*(itirev-1)+isym,2,idir,iat,jdir,jat)=re*sumi+im*sumr
2148              end do
2149            end do
2150          end do
2151        end do ! end treatment of the (iat,jat) (3x3) block of dynmat
2152      end if ! symmetry check
2153    end do ! time-reversal
2154  end do ! symmetries
2155 
2156  !4) make the average, get the final symmetric dynamical matrix
2157  do iat=1,natom
2158    do jat=1,natom
2159      do idir=1,3
2160        do jdir=1,3
2161          dmint=zero
2162          do isym=1,2*nsym
2163            dmint=dmint+dynmatint(isym,1,idir,iat,jdir,jat)
2164          end do
2165          dynmat(1,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
2166          dmint=zero
2167          do isym=1,2*nsym
2168            dmint=dmint+dynmatint(isym,2,idir,iat,jdir,jat)
2169          end do
2170          dynmat(2,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat))
2171        end do
2172      end do
2173    end do
2174  end do
2175 
2176  i1=0
2177  do iat=1,natom
2178    do idir=1,3
2179      i1=i1+1
2180      i2=0
2181      do jat=1,natom
2182        do jdir=1,3
2183          i2=i2+1
2184          index=i1+3*natom*(i2-1)
2185          dmati(2*index-1)=dynmat(1,idir,iat,jdir,jat)
2186          dmati(2*index  )=dynmat(2,idir,iat,jdir,jat)
2187        end do
2188      end do
2189    end do
2190  end do
2191 
2192 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

SOURCE

4745 subroutine sytens(indsym,mpert,natom,nsym,rfpert,symrec,symrel)
4746 
4747 !Arguments -------------------------------
4748 !scalars
4749  integer,intent(in) :: mpert,natom,nsym
4750 !arrays
4751  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
4752  integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert)
4753 
4754 !Local variables -------------------------
4755 !scalars
4756  integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_
4757  integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2
4758  integer :: ipesy3,isym
4759 !arrays
4760  integer :: sym1(3,3),sym2(3,3),sym3(3,3)
4761  integer,allocatable :: pertsy(:,:,:,:,:,:)
4762 
4763 !***********************************************************************
4764 
4765  ABI_MALLOC(pertsy,(3,mpert,3,mpert,3,mpert))
4766  pertsy(:,:,:,:,:,:) = 0
4767 
4768 !Loop over perturbations
4769 
4770  do i1pert_ = 1, mpert
4771    do i2pert_ = 1, mpert
4772      do i3pert_ = 1, mpert
4773 
4774        do i1dir_ = 1, 3
4775          do i2dir_ = 1, 3
4776            do i3dir_ = 1, 3
4777 
4778              i1pert = (mpert - i1pert_ + 1)
4779              if (i1pert <= natom) i1pert = natom + 1 - i1pert
4780              i2pert = (mpert - i2pert_ + 1)
4781              if (i2pert <= natom) i2pert = natom + 1 - i2pert
4782              i3pert = (mpert - i3pert_ + 1)
4783              if (i3pert <= natom) i3pert = natom + 1 - i3pert
4784 
4785              if (i1pert <= natom) then
4786                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
4787              else if (i2pert <= natom) then
4788                i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_
4789              else if (i3pert <= natom) then
4790                i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_
4791              else
4792                i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
4793              end if
4794 
4795              if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then
4796 
4797 !              Loop over all symmetries
4798 
4799                flag = 0
4800                do isym = 1, nsym
4801 
4802                  found = 1
4803 
4804 !                Select the symmetric element of i1pert,i2pert,i3pert
4805 
4806                  if (i1pert <= natom) then
4807                    ipesy1 = indsym(4,isym,i1pert)
4808                    sym1(:,:) = symrec(:,:,isym)
4809                  else if (i1pert == natom + 2) then
4810                    ipesy1 = i1pert
4811                    sym1(:,:) = symrel(:,:,isym)
4812                  else
4813                    found = 0
4814                  end if
4815 
4816                  if (i2pert <= natom) then
4817                    ipesy2 = indsym(4,isym,i2pert)
4818                    sym2(:,:) = symrec(:,:,isym)
4819                  else if (i2pert == natom + 2) then
4820                    ipesy2 = i2pert
4821                    sym2(:,:) = symrel(:,:,isym)
4822                  else
4823                    found = 0
4824                  end if
4825 
4826                  if (i3pert <= natom) then
4827                    ipesy3 = indsym(4,isym,i3pert)
4828                    sym3(:,:) = symrec(:,:,isym)
4829                  else if (i3pert == natom + 2) then
4830                    ipesy3 = i3pert
4831                    sym3(:,:) = symrel(:,:,isym)
4832                  else
4833                    found = 0
4834                  end if
4835 
4836 !                See if the symmetric element is available and check if some
4837 !                of the elements may be zeor. In the latter case, they do not need
4838 !                to be computed.
4839 
4840 
4841                  if ((flag /= -1).and.&
4842 &                 (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
4843                    flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir)
4844                  end if
4845 
4846 
4847                  do idisy1 = 1, 3
4848                    do idisy2 = 1, 3
4849                      do idisy3 = 1, 3
4850 
4851                        if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.&
4852 &                       (sym3(i3dir,idisy3) /= 0)) then
4853                          if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
4854                            found = 0
4855 !                          exit      ! exit loop over symmetries
4856                          end if
4857                        end if
4858 
4859 
4860                        if ((flag == -1).and.&
4861 &                       ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then
4862                          if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.&
4863 &                         (sym3(i3dir,idisy3)/=0)) then
4864                            flag = 0
4865                          end if
4866                        end if
4867 
4868                      end do
4869                    end do
4870                  end do
4871 
4872                  if (found == 1) then
4873                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1
4874                  end if
4875 
4876 !                In case a symmetry operation only changes the sign of an
4877 !                element, this element has to be equal to zero
4878 
4879                  if (flag == -1) then
4880                    pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2
4881                    exit
4882                  end if
4883 
4884                end do    ! close loop on symmetries
4885 
4886 !              If the elemetn i1pert,i2pert,i3pert is not symmetric
4887 !              to a basis element, it is a basis element
4888 
4889                if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then
4890                  pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4891                end if
4892 
4893              end if ! rfpert /= 0
4894 
4895            end do        ! close loop over perturbations
4896          end do
4897        end do
4898      end do
4899    end do
4900  end do
4901 
4902 !Now, take into account the permutation of (i1pert,i1dir)
4903 !and (i3pert,i3dir)
4904 
4905  do i1pert = 1, mpert
4906    do i2pert = 1, mpert
4907      do i3pert = 1, mpert
4908 
4909        do i1dir = 1, 3
4910          do i2dir = 1, 3
4911            do i3dir = 1, 3
4912 
4913              if ((i1pert /= i3pert).or.(i1dir /= i3dir)) then
4914 
4915                if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.&
4916 &               (pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) == 1)) then
4917                  pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = -1
4918                end if
4919 
4920              end if
4921 
4922            end do
4923          end do
4924        end do
4925 
4926      end do
4927    end do
4928  end do
4929 
4930  rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:)
4931 
4932  ABI_FREE(pertsy)
4933 
4934 end subroutine sytens

m_dynmat/wght9 [ Functions ]

[ Top ] [ m_dynmat ] [ Functions ]

NAME

 wght9

FUNCTION

 Generates a weight to each R point 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)
 toldist= Tolerance on the distance between two R points.

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

SOURCE

3874 subroutine wght9(brav,gprim,natom,ngqpt,nqpt,nqshft,nrpt,qshft,rcan,rpt,rprimd,toldist,r_inscribed_sphere,wghatm,ierr)
3875 
3876 !Arguments -------------------------------
3877 !scalars
3878  integer,intent(in) :: brav,natom,nqpt,nqshft,nrpt
3879  integer,intent(out) :: ierr
3880  real(dp),intent(out) :: r_inscribed_sphere
3881  real(dp),intent(in) :: toldist
3882 !arrays
3883  integer,intent(inout) :: ngqpt(9)
3884  real(dp),intent(in) :: gprim(3,3),qshft(3,4),rcan(3,natom),rpt(3,nrpt),rprimd(3,3)
3885  real(dp),intent(out) :: wghatm(natom,natom,nrpt)
3886 
3887 !Local variables -------------------------
3888 !scalars
3889  integer :: ia,ib,ii,jj,kk,iqshft,irpt,jqshft,nbordh,tok,nptws,nreq,idir
3890  real(dp) :: factor,sumwght,normsq,proj
3891  character(len=500) :: msg
3892 !arrays
3893  integer :: nbord(9)
3894  real(dp) :: rdiff(9),red(3,3),ptws(4, 729),pp(3),rdiff_tmp(3)
3895 
3896 ! *********************************************************************
3897 
3898  ierr = 0
3899 
3900  ! First analyze the vectors qshft
3901  if (nqshft /= 1) then
3902 
3903    if (brav == 4) then
3904      write(msg,'(3a,i0,3a)' )&
3905      'For the time being, only nqshft=1',ch10,&
3906      'is allowed with brav=4, while it is nqshft=',nqshft,'.',ch10,&
3907      'Action: in the input file, correct either brav or nqshft.'
3908      ABI_ERROR(msg)
3909    end if
3910 
3911    if (nqshft == 2) then
3912      ! Make sure that the q vectors form a BCC lattice
3913      do ii=1,3
3914        if(abs(abs(qshft(ii,1)-qshft(ii,2))-.5_dp)>1.d-10)then
3915          write(msg, '(a,a,a,a,a,a,a)' )&
3916          'The test of the q1shft vectors shows that they',ch10,&
3917          'do not generate a body-centered lattice, which',ch10,&
3918          'is mandatory for nqshft=2.',ch10,&
3919          'Action: change the q1shft vectors in your input file.'
3920          ABI_ERROR(msg)
3921        end if
3922      end do
3923    else if (nqshft == 4) then
3924      ! Make sure that the q vectors form a FCC lattice
3925      do iqshft=1,3
3926        do jqshft=iqshft+1,4
3927          tok=0
3928          do ii=1,3
3929            ! Test on the presence of a +-0.5 difference
3930            if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-.5_dp) <1.d-10) tok=tok+1
3931            ! Test on the presence of a 0 or +-1.0 difference
3932            if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-1._dp) <1.d-10  .or.&
3933               abs(qshft(ii,iqshft)-qshft(ii,jqshft)) < 1.d-10) tok=tok+4
3934          end do
3935          ! Test 1 should be satisfied twice, and test 2 once
3936          if(tok/=6)then
3937            write(msg, '(7a)' )&
3938            'The test of the q1shft vectors shows that they',ch10,&
3939            'do not generate a face-centered lattice, which',ch10,&
3940            'is mandatory for nqshft=4.',ch10,&
3941            'Action: change the q1shft vectors in your input file.'
3942            ABI_ERROR(msg)
3943          end if
3944        end do
3945      end do
3946    else
3947      write(msg, '(a,i4,3a)' )&
3948      'nqshft must be 1, 2 or 4. It is nqshft=',nqshft,'.',ch10,&
3949      'Action: change nqshft in your input file.'
3950      ABI_ERROR(msg)
3951    end if
3952  end if
3953 
3954  factor=0.5_dp
3955  if(brav==2 .or. brav==3) factor=0.25_dp
3956  if(nqshft/=1)factor=factor*2
3957 
3958  if (brav==1) then
3959    ! Does not support multiple shifts
3960    if (nqshft/=1) then
3961      ABI_ERROR('This version of the weights does not support nqshft/=1.')
3962    end if
3963 
3964    ! Find the points of the lattice given by ngqpt*acell. These are used to define
3965    ! a Wigner-Seitz cell around the origin. The origin is excluded from the list.
3966    ! TODO : in principle this should be only -1 to +1 for ii jj kk!
3967    nptws=0
3968    do ii=-2,2
3969      do jj=-2,2
3970        do kk=-2,2
3971          do idir=1,3
3972            pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3)
3973          end do
3974          normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3)
3975          if (normsq > tol6) then
3976            nptws = nptws + 1
3977            ptws(:3,nptws) = pp(:)
3978            ptws(4,nptws) = half*normsq
3979          end if
3980        end do
3981      end do
3982    end do
3983  end if ! end new_wght
3984  !write(std_out,*)'factor,ngqpt',factor,ngqpt(1:3)
3985 
3986  r_inscribed_sphere = sum((matmul(rprimd(:,:),ngqpt(1:3)))**2)
3987  do ii=-1,1
3988    do jj=-1,1
3989      do kk=-1,1
3990        if (ii==0 .and. jj==0 .and. kk==0) cycle
3991        do idir=1,3
3992          pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3)
3993        end do
3994        normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3)
3995        r_inscribed_sphere = min(r_inscribed_sphere, normsq)
3996      end do
3997    end do
3998  end do
3999  r_inscribed_sphere = sqrt(r_inscribed_sphere)
4000 
4001 
4002 !Begin the big loop on ia and ib
4003  do ia=1,natom
4004    do ib=1,natom
4005 
4006      ! Simple Lattice
4007      if (abs(brav)==1) then
4008        ! In this case, it is better to work in reduced coordinates
4009        ! As rcan is in canonical coordinates, => multiplication by gprim
4010        do ii=1,3
4011          red(1,ii)=  rcan(1,ia)*gprim(1,ii) +rcan(2,ia)*gprim(2,ii) +rcan(3,ia)*gprim(3,ii)
4012          red(2,ii)=  rcan(1,ib)*gprim(1,ii) +rcan(2,ib)*gprim(2,ii) +rcan(3,ib)*gprim(3,ii)
4013        end do
4014      end if
4015 
4016      do irpt=1,nrpt
4017 
4018        ! Initialization of the weights to 1.0
4019        wghatm(ia,ib,irpt)=1.0_dp
4020 
4021        ! Compute the difference vector
4022 
4023        ! Simple Cubic Lattice
4024        if (abs(brav)==1) then
4025          ! Change of rpt to reduced coordinates
4026          do ii=1,3
4027            red(3,ii)=  rpt(1,irpt)*gprim(1,ii) +rpt(2,irpt)*gprim(2,ii) +rpt(3,irpt)*gprim(3,ii)
4028            rdiff(ii)=red(2,ii)-red(1,ii)+red(3,ii)
4029          end do
4030          if (brav==1) then
4031            ! rdiff in cartesian coordinates
4032            do ii=1,3
4033              rdiff_tmp(ii)=rdiff(1)*rprimd(ii,1)+rdiff(2)*rprimd(ii,2)+rdiff(3)*rprimd(ii,3)
4034            end do
4035            rdiff(1:3)=rdiff_tmp(1:3)
4036          end if
4037 
4038        else
4039          ! Other lattices
4040          do ii=1,3
4041            rdiff(ii)=rcan(ii,ib)-rcan(ii,ia)+rpt(ii,irpt)
4042          end do
4043        end if
4044 
4045        ! Assignement of weights
4046 
4047        if(nqshft==1 .and. brav/=4)then
4048 
4049          if (brav/=1) then
4050            do ii=1,3
4051              ! If the rpt vector is greater than the allowed space => weight = 0.0
4052              if (abs(rdiff(ii))-tol10>factor*ngqpt(ii)) then
4053                wghatm(ia,ib,irpt)=zero
4054              else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4055                ! If the point is in a boundary position => weight/2
4056                wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4057              end if
4058            end do
4059          else
4060            ! new weights
4061            wghatm(ia,ib,irpt)=zero
4062            nreq = 1
4063            do ii=1,nptws
4064              proj = rdiff(1)*ptws(1,ii)+rdiff(2)*ptws(2,ii)+rdiff(3)*ptws(3,ii)
4065              ! if rdiff closer to ptws than the origin the weight is zero
4066              ! if rdiff close to the origin with respect to all the other ptws the weight is 1
4067              ! if rdiff is equidistant from the origin and N other ptws the weight is 1/(N+1)
4068              if (proj - ptws(4,ii) > toldist) then
4069                nreq = 0
4070                EXIT
4071              else if(abs(proj-ptws(4,ii)) <= toldist) then
4072                nreq=nreq+1
4073              end if
4074            end do
4075            if (nreq>0) then
4076              wghatm(ia,ib,irpt)=one/DBLE(nreq)
4077            end if
4078          end if
4079 
4080        else if(brav==4)then
4081          ! Hexagonal
4082          ! Examination of the X and Y boundaries in order to form an hexagon
4083          ! First generate the relevant boundaries
4084          rdiff(4)=0.5_dp*( rdiff(1)+sqrt(3.0_dp)*rdiff(2) )
4085          ngqpt(4)=ngqpt(1)
4086          rdiff(5)=0.5_dp*( rdiff(1)-sqrt(3.0_dp)*rdiff(2) )
4087          ngqpt(5)=ngqpt(1)
4088 
4089          ! Test the four inequalities
4090          do ii=1,5
4091            if(ii/=2)then
4092 
4093              nbord(ii)=0
4094              ! If the rpt vector is greater than the allowed space => weight = 0.0
4095              if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4096                wghatm(ia,ib,irpt)=zero
4097              else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4098                ! If the point is in a boundary position increment nbord(ii)
4099                nbord(ii)=1
4100              end if
4101 
4102            end if
4103          end do
4104 
4105          ! Computation of weights
4106          nbordh=nbord(1)+nbord(4)+nbord(5)
4107          if (nbordh==1) then
4108            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4109          else if (nbordh==2) then
4110            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4111          else if (nbordh/=0) then
4112            ABI_BUG('There is a problem of borders and weights (hex).')
4113          end if
4114          if (nbord(3)==1)then
4115            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4116          end if
4117 
4118        else if(nqshft==2 .and. brav/=4)then
4119 
4120          ! BCC packing of k-points
4121          ! First, generate the relevant boundaries
4122          rdiff(4)= rdiff(1)+rdiff(2)
4123          rdiff(5)= rdiff(1)-rdiff(2)
4124          rdiff(6)= rdiff(1)+rdiff(3)
4125          rdiff(7)= rdiff(1)-rdiff(3)
4126          rdiff(8)= rdiff(3)+rdiff(2)
4127          rdiff(9)= rdiff(3)-rdiff(2)
4128          if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
4129            write(msg, '(a,a,a,3i6,a,a,a,a)' )&
4130            'In the BCC case, the three ngqpt numbers ',ch10,&
4131            '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
4132            'should be equal.',ch10,&
4133            'Action: use identical ngqpt(1:3) in your input file.'
4134            ABI_ERROR(msg)
4135          end if
4136          do ii=4,9
4137            ngqpt(ii)=ngqpt(1)
4138          end do
4139 
4140          ! Test the relevant inequalities
4141          nbord(1)=0
4142          do ii=4,9
4143            ! If the rpt vector is greater than the allowed space => weight = 0.0
4144            if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4145              wghatm(ia,ib,irpt)=zero
4146            else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4147              ! If the point is in a boundary position increment nbord(1)
4148              nbord(1)=nbord(1)+1
4149            end if
4150          end do
4151 
4152          ! Computation of weights
4153          if (nbord(1)==1) then
4154            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4155          else if (nbord(1)==2) then
4156            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4157          else if (nbord(1)==3) then
4158            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
4159          else if (nbord(1)==4) then
4160            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/6
4161          else if (nbord(1)/=0) then
4162            ABI_ERROR(' There is a problem of borders and weights (BCC).')
4163          end if
4164 
4165        else if(nqshft==4 .and. brav/=4)then
4166 
4167          ! FCC packing of k-points
4168          ! First, generate the relevant boundaries
4169          rdiff(4)= (rdiff(1)+rdiff(2)+rdiff(3))*2._dp/3._dp
4170          rdiff(5)= (rdiff(1)-rdiff(2)+rdiff(3))*2._dp/3._dp
4171          rdiff(6)= (rdiff(1)+rdiff(2)-rdiff(3))*2._dp/3._dp
4172          rdiff(7)= (rdiff(1)-rdiff(2)-rdiff(3))*2._dp/3._dp
4173          if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then
4174            write(msg, '(a,a,a,3i6,a,a,a,a)' )&
4175            'In the FCC case, the three ngqpt numbers ',ch10,&
4176            '    ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,&
4177            'should be equal.',ch10,&
4178            'Action: use identical ngqpt(1:3) in your input file.'
4179            ABI_ERROR(msg)
4180          end if
4181          do ii=4,7
4182            ngqpt(ii)=ngqpt(1)
4183          end do
4184 
4185          ! Test the relevant inequalities
4186          nbord(1)=0
4187          do ii=1,7
4188            ! If the rpt vector is greater than the allowed space => weight = 0.0
4189            if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then
4190              wghatm(ia,ib,irpt)=zero
4191              ! If the point is in a boundary position increment nbord(1)
4192            else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then
4193              nbord(1)=nbord(1)+1
4194            end if
4195          end do
4196 
4197          ! Computation of weights
4198          if (nbord(1)==1) then
4199            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2
4200          else if (nbord(1)==2) then
4201            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3
4202          else if (nbord(1)==3) then
4203            wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4
4204          else if (nbord(1)/=0 .and. wghatm(ia,ib,irpt)>1.d-10) then
4205            ! Interestingly nbord(1)==4 happens for some points outside of the volume
4206            ABI_BUG(' There is a problem of borders and weights (FCC).')
4207          end if
4208 
4209        else
4210          write(msg, '(3a,i0,a)' )&
4211          'One should not arrive here ... ',ch10,&
4212          'The value nqshft ',nqshft,' is not available'
4213          ABI_BUG(msg)
4214        end if
4215      end do ! Assignement of weights is done
4216    end do ! End of the double loop on ia and ib
4217  end do
4218 
4219  ! Check the results
4220  do ia=1,natom
4221    do ib=1,natom
4222      sumwght=zero
4223      do irpt=1,nrpt
4224        ! Check if the sum of the weights is equal to the number of q points
4225        sumwght=sumwght+wghatm(ia,ib,irpt)
4226        !write(std_out,'(a,3(i0,1x))' )' atom1, atom2, irpt ; rpt ; wghatm ',ia,ib,irpt
4227        !write(std_out,'(3es16.6,es18.6)' )rpt(1,irpt),rpt(2,irpt),rpt(3,irpt),wghatm(ia,ib,irpt)
4228      end do
4229      if (abs(sumwght-nqpt)>tol10) ierr = 1
4230    end do
4231  end do
4232 
4233 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

SOURCE

2506 subroutine wings3(carflg,d2cart,mpert)
2507 
2508 !Arguments -------------------------------
2509 !scalars
2510  integer,intent(in) :: mpert
2511 !arrays
2512  integer,intent(inout) :: carflg(3,mpert,3,mpert)
2513  real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert)
2514 
2515 !Local variables -------------------------
2516 !scalars
2517  integer :: idir,idir1,ipert,ipert1
2518 
2519 ! *********************************************************************
2520 
2521  do ipert=1,mpert
2522    do idir=1,3
2523      if(carflg(idir,ipert,idir,ipert)==0)then
2524        do ipert1=1,mpert
2525          do idir1=1,3
2526            carflg(idir,ipert,idir1,ipert1)=0
2527            carflg(idir1,ipert1,idir,ipert)=0
2528            d2cart(1,idir,ipert,idir1,ipert1)=zero
2529            d2cart(2,idir,ipert,idir1,ipert1)=zero
2530            d2cart(1,idir1,ipert1,idir,ipert)=zero
2531            d2cart(2,idir1,ipert1,idir,ipert)=zero
2532          end do
2533        end do
2534      end if
2535    end do
2536  end do
2537 
2538 end subroutine wings3