TABLE OF CONTENTS


ABINIT/irreducible_set_pert [ Functions ]

[ Top ] [ Functions ]

NAME

 irreducible_set_pert

FUNCTION

 Determines a set of perturbations that form a basis
 in that, using symmetry, they can be used to generate
 all other perturbations that are asked to be calculated (target).

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 iper
  natom= number of atoms
  nsym=number of space group symmetries
  rfdir(3)=direction for the perturbations
  rfpert(mpert)=information on the perturbations
  symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
  symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
  symq(4,2,nsym)= (integer) three first numbers define the G vector;
   fourth number is 0 if the q-vector is not preserved, is 1 otherwise
   second index is one without time-reversal symmetry, two with time-reversal symmetry

OUTPUT

   pertsy(3,mpert)= the target perturbation is described by the two last indices (idir, and ipert),
                    the value is 0, 1 or -1, see notes.

NOTES

 Output will be in the pertsy array,
   0 for non-target perturbations
   1 for basis perturbations
  -1 for perturbations that can be found from basis perturbations

PARENTS

      get_npert_rbz,m_dvdb,respfn

CHILDREN

SOURCE

3712 subroutine irreducible_set_pert(indsym,mpert,natom,nsym,pertsy,rfdir,rfpert,symq,symrec,symrel)
3713 
3714 
3715 !This section has been created automatically by the script Abilint (TD).
3716 !Do not modify the following lines by hand.
3717 #undef ABI_FUNC
3718 #define ABI_FUNC 'irreducible_set_pert'
3719 !End of the abilint section
3720 
3721  implicit none
3722 
3723 !Arguments -------------------------------
3724 !scalars
3725  integer,intent(in) :: mpert,natom,nsym
3726 !arrays
3727  integer,intent(in) :: indsym(4,nsym,natom),rfdir(3),rfpert(mpert)
3728  integer,intent(in) :: symq(4,2,nsym),symrec(3,3,nsym),symrel(3,3,nsym)
3729  integer,intent(out) :: pertsy(3,mpert)
3730 
3731 !Local variables -------------------------
3732 !scalars
3733  integer :: found,idir1,idisy1,ii,ipert1,ipesy1,isign,isym,itirev,jj
3734 !arrays
3735  integer :: sym1(3,3)
3736 
3737 ! *********************************************************************
3738 
3739 !Zero pertsy
3740  pertsy(:,:)=0
3741 
3742  do ipert1=1,mpert
3743    do idir1=1,3
3744      if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
3745 !      write(std_out,*)' for candidate idir =',idir1,' ipert = ',ipert1
3746 
3747 !      Loop on all symmetries, including time-reversal
3748        do isym=1,nsym
3749          do itirev=1,2
3750            isign=3-2*itirev
3751 
3752            if(symq(4,itirev,isym)/=0)then
3753 
3754              found=1
3755 
3756 !            Here select the symmetric of ipert1
3757              if(ipert1<=natom)then
3758                ipesy1=indsym(4,isym,ipert1)
3759                do ii=1,3
3760                  do jj=1,3
3761                    sym1(ii,jj)=symrec(ii,jj,isym)
3762                  end do
3763                end do
3764              else if(ipert1==(natom+2) .or. ipert1==(natom+6))then
3765                ipesy1=ipert1
3766                do ii=1,3
3767                  do jj=1,3
3768                    sym1(ii,jj)=symrel(ii,jj,isym)
3769                  end do
3770                end do
3771              else
3772                found=0
3773              end if
3774 
3775 !            Now that a symmetric perturbation has been obtained,
3776 !            including the expression of the symmetry matrix, see
3777 !            if the symmetric perturbations are available
3778              if( found==1 ) then
3779 
3780                do idisy1=1,3
3781                  if(sym1(idir1,idisy1)/=0)then
3782                    if(pertsy(idisy1,ipesy1)==0)then
3783                      found=0
3784                      exit
3785                    end if
3786                  end if
3787                end do
3788              end if
3789 
3790 !            Now, if still found, then it is a symmetric
3791 !            of some linear combination of existing perturbations
3792              if(found==1)then
3793 
3794 !              DEBUG
3795 !              write(std_out,*)' all found !  isym, isign= ',isym,isign
3796 !              write(std_out,1010)((sym1(ii,jj),ii=1,3),jj=1,3)
3797 !              write(std_out,1010)((sym2(ii,jj),ii=1,3),jj=1,3)
3798 !              write(std_out,*)sumr,sumi
3799 !              1010    format(9i4)
3800 !              ENDDEBUG
3801 
3802                pertsy(idir1,ipert1)=-1
3803                exit ! Exit loop on symmetry operations
3804 
3805              end if
3806 
3807            end if !  End loop on all symmetries + time-reversal
3808          end do
3809        end do
3810 
3811 !      Now that all symmetries have been examined,
3812 !      if still not symmetric of a linear combination
3813 !      of basis perturbations, then it is a basis perturbation
3814        if(pertsy(idir1,ipert1)/=-1) pertsy(idir1,ipert1)=1
3815 !      write(std_out,'(a,3i5)' ) ' irreducible_set_pert :',idir1,ipert1,pertsy(idir1,ipert1)
3816 
3817      end if ! End big loop on all elements
3818    end do
3819  end do
3820 
3821 end subroutine irreducible_set_pert

ABINIT/m_geometry [ Modules ]

[ Top ] [ Modules ]

NAME

  m_geometry

FUNCTION

  This module contains basic tools to operate on vectors expressed in reduced coordinates.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG, MT, FJ, TRangel, DCA, XG, AHR, DJA, DRH)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_geometry
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_atomdata
29  use m_sort
30 
31  use m_io_tools,       only : open_file
32  use m_numeric_tools,  only : uniformrandom, isinteger, set2unit
33  use m_symtk,          only : mati3inv, mati3det, matr3inv, symdet
34  use m_hide_lapack,    only : matr3eigval
35  use m_pptools,        only : prmat
36  use m_numeric_tools,  only : wrap2_pmhalf
37  use m_hide_lapack,    only : matrginv
38 
39  implicit none
40 
41  private
42 
43  public :: normv              ! Norm of vector(s) in reduced coordinates either in real or reciprocal space.
44  public :: vdotw              ! Scalar product between two reduced vectors either in real or reciprocal space.
45  public :: acrossb            ! Cross product of two 3-vectors.
46  public :: wigner_seitz       ! Find the grid of points falling inside the Wigner-Seitz cell.
47  public :: phdispl_cart2red   ! Calculate the displacement vectors for all branches in reduced coordinates.
48  public :: getspinrot         ! Compute the components of the spinor rotation matrix
49  public :: spinrot_cmat       ! Construct 2x2 complex matrix representing rotation operator in spin-space.
50  public :: rotmat             ! Finds the rotation matrix.
51  public :: fixsym             ! Check that iatfix does not break symmetry.
52  public :: metric             ! Compute metric matrices.
53  public :: mkradim            ! Make rprim and acell from rprimd
54  public :: mkrdim             ! Make rprimd from acell from rprim
55  public :: chkrprimd          ! Test if {rprim,acell,rprimd} are consistent
56  public :: chkdilatmx         ! check if dilatation of unit cell is consistent with initial G-sphere
57  public :: xcart2xred         ! From cart coords to reduced
58  public :: xred2xcart         ! From reduced coords to cart.
59  public :: fred2fcart         ! Convert reduced forces into cartesian forces
60  public :: fcart2fred         ! Convert cartesian forces into reduced forces
61  public :: bonds_lgth_angles  ! Write GEO file
62  public :: randomcellpos      ! Creates unit cell with random atomic positions.
63  public :: ioniondist         ! Compute ion-ion distances
64  public :: dist2              ! Calculates the distance of v1 and v2 in a crystal by epeating the unit cell
65  public :: shellstruct        ! Calculates shell structure (multiplicities, radii)
66  public :: remove_inversion   ! Remove the inversion symmetry and improper rotations
67  public :: symredcart         ! Convert a symmetry operation from reduced coordinates (integers) to cart coords (reals)
68  public :: strainsym          ! Symmetrize the strain tensor.
69  public :: stresssym          ! Symmetrize the stress tensor.
70  public :: strconv            ! Convert from symmetric storage mode in reduced coords to cart coords.
71  public :: littlegroup_pert   ! Determines the set of symmetries that leaves a perturbation invariant.
72  public :: irreducible_set_pert  ! Determines a set of perturbations that form a basis
73 
74  interface normv
75   module procedure normv_rdp_vector
76   module procedure normv_int_vector
77   !module procedure normv_int_vector_array  ! WARNING for the time being, do not use these 2 procedures,
78   !module procedure normv_rdp_vector_array  ! sunstudio12 is not able to resolve which sub should be called.
79  end interface normv
80 
81  interface vdotw
82   module procedure vdotw_rr_vector
83   module procedure vdotw_rc_vector
84  end interface vdotw
85 
86 CONTAINS  !===========================================================

m_geometry/acrossb [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 acrossb

FUNCTION

 Calculates the cross product of two 3-vectors

INPUTS

   a(3): real(dp) vector
   b(3): real(dp) vector

OUTPUT

   c(3): real(dp) vector = a X b

PARENTS

      calc_b_matrix,m_abimover,simple_j_dia

CHILDREN

SOURCE

467 subroutine acrossb(a,b,c)
468 
469 
470 !This section has been created automatically by the script Abilint (TD).
471 !Do not modify the following lines by hand.
472 #undef ABI_FUNC
473 #define ABI_FUNC 'acrossb'
474 !End of the abilint section
475 
476  implicit none
477 
478 !Arguments ---------------------------------------------
479 !arrays
480  real(dp),intent(in) :: a(3),b(3)
481  real(dp),intent(out) :: c(3)
482 
483 ! *********************************************************************
484 
485  c(1) =  a(2)*b(3) - a(3)*b(2)
486  c(2) = -a(1)*b(3) + a(3)*b(1)
487  c(3) =  a(1)*b(2) - b(1)*a(2)
488 
489 end subroutine acrossb

m_geometry/bonds_lgth_angles [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 bonds_lgth_angles

FUNCTION

 From list of coordinates and primitive translations, output
 a list of bonds lengths and bond angles.

INPUTS

  coordn = maximum coordination number to be taken into account
  fnameabo_app_geo=name of file for _GEO data
  natom  = number of atoms in unit cell
  ntypat = number of types of atoms in unit cell.
  rprimd(3,3)  = real space dimensional primitive translations (bohr)
  typat(natom) = type integer for each atom in cell
  znucl(ntypat)= real(dp), atomic number of atom type
  xred(3,natom)= reduced coordinates of atoms

OUTPUT

 data written in file fnameabo_app_geo

NOTES

  The tolerance tol8 aims at giving a machine-independent ordering.
  (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)

PARENTS

      outscfcv

CHILDREN

      atomdata_from_znucl,wrtout,xred2xcart

SOURCE

2001 subroutine bonds_lgth_angles(coordn,fnameabo_app_geo,natom,ntypat,rprimd,typat,xred,znucl)
2002 
2003 
2004 !This section has been created automatically by the script Abilint (TD).
2005 !Do not modify the following lines by hand.
2006 #undef ABI_FUNC
2007 #define ABI_FUNC 'bonds_lgth_angles'
2008 !End of the abilint section
2009 
2010  implicit none
2011 
2012 !Arguments ------------------------------------
2013 !scalars
2014  integer,intent(in) :: coordn,natom,ntypat
2015  character(len=*),intent(in) :: fnameabo_app_geo
2016 !arrays
2017  integer,intent(in) :: typat(natom)
2018  real(dp),intent(in) :: rprimd(3,3),znucl(ntypat)
2019  real(dp),intent(inout) :: xred(3,natom)
2020 
2021 !Local variables-------------------------------
2022 !scalars
2023  integer :: done,ia,ib,ic,ii,ineighb,jneighb,mneighb,mu,ndig,nu,t1,t2,t3,tmax,temp_unit
2024  real(dp) :: adotb,asq,bsq,co,length,sq,thdeg
2025 !real(dp)u1,u2,u3,v1,v2,v3
2026  character(len=500) :: message
2027  type(atomdata_t) :: atom
2028 !arrays
2029  integer,allocatable :: list_neighb(:,:,:)
2030  real(dp) :: bab(3),bac(3),dif(3),rmet(3,3)
2031  real(dp),allocatable :: sqrlength(:),xangst(:,:),xcart(:,:)
2032  character(len=8),allocatable :: iden(:)
2033 
2034 ! *************************************************************************
2035 
2036 !Initialize the file
2037  write(message, '(a,a,a)' )' bonds_lgth_angles : about to open file ',trim(fnameabo_app_geo),ch10
2038  call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
2039 
2040  if (open_file(fnameabo_app_geo,message,newunit=temp_unit,status='unknown',form='formatted') /= 0) then
2041    MSG_ERROR(message)
2042  end if
2043  rewind(temp_unit)
2044 
2045  write(message, '(a,a)' ) ch10,' ABINIT package : GEO file '
2046  call wrtout(temp_unit,message,'COLL')
2047 
2048 !Compute maximum number of neighbors is the neighbor list,
2049 !from the indicative coordination number
2050 !Note : the following formula includes next nearest neighbors, but not others
2051  mneighb=1+coordn+coordn*(coordn-1)
2052 
2053  write(message, '(a,a,i2,a,a,i4,a,a,a,i4,a)' ) ch10,&
2054 & ' Maximal coordination number, as estimated by the user : ',coordn,ch10,&
2055 & '  giving a maximum of ',coordn*coordn,&
2056 & ' nearest neighbors and next nearest neighbors, ',ch10,&
2057 & '                  and ',(coordn*(coordn-1))/2,&
2058 & ' distinct angles between nearest neighbors'
2059  call wrtout(temp_unit,message,'COLL')
2060 
2061 !Compute metric tensor in real space rmet
2062  do nu=1,3
2063    do mu=1,3
2064      rmet(mu,nu)=rprimd(1,mu)*rprimd(1,nu)+&
2065 &     rprimd(2,mu)*rprimd(2,nu)+&
2066 &     rprimd(3,mu)*rprimd(3,nu)
2067    end do
2068  end do
2069 
2070  write(message, '(a,a)' )ch10,' Primitive vectors of the periodic cell (bohr)'
2071  call wrtout(temp_unit,message,'COLL')
2072  do nu=1,3
2073    write(message, '(1x,a,i1,a,3f10.5)' ) '  R(',nu,')=',rprimd(:,nu)
2074    call wrtout(temp_unit,message,'COLL')
2075  end do
2076 
2077  write(message, '(a,a)' ) ch10,&
2078 & ' Atom list        Reduced coordinates          Cartesian coordinates (bohr)'
2079  call wrtout(temp_unit,message,'COLL')
2080 
2081 !Set up a list of character identifiers for all atoms : iden(ia)
2082  ABI_ALLOCATE(iden,(natom))
2083  iden(:)='        '
2084  do ia=1,natom
2085    ndig=int(log10(dble(ia)+0.5d0))+1
2086    call atomdata_from_znucl(atom,znucl(typat(ia)))
2087    if(ndig==1) write(iden(ia), '(a,a,i1,a)' )  atom%symbol,'(',ia,')   '
2088    if(ndig==2) write(iden(ia), '(a,a,i2,a)' )  atom%symbol,'(',ia,')  '
2089    if(ndig==3) write(iden(ia), '(a,a,i3,a)' )  atom%symbol,'(',ia,') '
2090    if(ndig==4) write(iden(ia), '(a,a,i4,a)' )  atom%symbol,'(',ia,')'
2091    if(ndig>4)then
2092      close(temp_unit)
2093      write(message, '(a,i8,a,a)' )&
2094 &     'bonds_lgth_angles cannot handle more than 9999 atoms, while natom=',natom,ch10,&
2095 &     'Action: decrease natom, or contact ABINIT group.'
2096      MSG_BUG(message)
2097    end if
2098  end do
2099 
2100 !Compute cartesian coordinates, and print reduced and cartesian coordinates
2101 !then print coordinates in angstrom, with the format neede for xmol
2102  ABI_ALLOCATE(xangst,(3,natom))
2103  ABI_ALLOCATE(xcart,(3,natom))
2104  call xred2xcart(natom,rprimd,xcart,xred)
2105  xangst(:,:)=xcart(:,:)*Bohr_Ang
2106 
2107  do ia=1,natom
2108    write(message, '(a,a,3f10.5,a,3f10.5)' ) &
2109 &   '   ',iden(ia),(xred(ii,ia)+tol10,ii=1,3),&
2110 &   '    ',(xcart(ii,ia)+tol10,ii=1,3)
2111    call wrtout(temp_unit,message,'COLL')
2112  end do
2113 
2114  write(message, '(a,a,a,a,i4,a)' )ch10,&
2115 & ' XMOL data : natom, followed by cartesian coordinates in Angstrom',&
2116 & ch10,ch10,natom,ch10
2117  call wrtout(temp_unit,message,'COLL')
2118 
2119  do ia=1,natom
2120    call atomdata_from_znucl(atom,znucl(typat(ia)))
2121    write(message, '(a,a,3f10.5)' )'   ',atom%symbol,xangst(1:3,ia)
2122    call wrtout(temp_unit,message,'COLL')
2123  end do
2124 
2125  ABI_DEALLOCATE(xangst)
2126  ABI_DEALLOCATE(xcart)
2127 
2128  ABI_ALLOCATE(list_neighb,(0:mneighb+1,4,2))
2129  ABI_ALLOCATE(sqrlength,(0:mneighb+1))
2130 
2131 !Compute list of neighbors
2132  do ia=1,natom
2133 
2134    write(message, '(a,a,a,a,a,a,a,a,a)' ) ch10,'===========',&
2135 &   '=====================================================================',&
2136 &   ch10,' ',iden(ia),ch10,ch10,' Bond lengths '
2137    call wrtout(temp_unit,message,'COLL')
2138 
2139 !  Search other atoms for bonds, but must proceed
2140 !  in such a way to consider a search box sufficiently large,
2141 !  so increase the size of the search box until the
2142 !  final bond length list do not change
2143    do tmax=0,5
2144 
2145 !    Set initial list of neighbors to zero,
2146 !    and initial square of bond lengths to a very large number.
2147 !    Note that the dimension is larger than neighb to ease
2148 !    the later sorting : neighbors 0 and neighb+1 are non-existent, while
2149 !    neighbor 1 will be the atom itself ...
2150      list_neighb(0:mneighb+1,1:4,1)=0
2151      sqrlength(1:mneighb+1)=huge(0.0d0)
2152      sqrlength(0)=-1.0d0
2153 
2154 !    Here search on all atoms inside the box defined by tmax
2155      do ib=1,natom
2156        do t3=-tmax,tmax
2157          do t2=-tmax,tmax
2158            do t1=-tmax,tmax
2159              dif(1)=xred(1,ia)-(xred(1,ib)+dble(t1))
2160              dif(2)=xred(2,ia)-(xred(2,ib)+dble(t2))
2161              dif(3)=xred(3,ia)-(xred(3,ib)+dble(t3))
2162              sq=rsdot(dif(1),dif(2),dif(3),dif(1),dif(2),dif(3),rmet)
2163 
2164 !            Insert the atom at the proper place in the neighbor list.
2165              do ineighb=mneighb,0,-1
2166 !              Note the tolerance
2167                if(sq+tol8>sqrlength(ineighb))then
2168                  sqrlength(ineighb+1)=sq
2169                  list_neighb(ineighb+1,1,1)=ib
2170                  list_neighb(ineighb+1,2,1)=t1
2171                  list_neighb(ineighb+1,3,1)=t2
2172                  list_neighb(ineighb+1,4,1)=t3
2173 !                DEBUG
2174 !                if(ineighb/=mneighb)then
2175 !                write(std_out,*)' '
2176 !                do ii=1,mneighb
2177 !                write(std_out,*)ii,sqrlength(ii)
2178 !                end do
2179 !                end if
2180 !                ENDDEBUG
2181                  exit
2182                else
2183                  sqrlength(ineighb+1)=sqrlength(ineighb)
2184                  list_neighb(ineighb+1,1:4,1)=list_neighb(ineighb,1:4,1)
2185                end if
2186              end do
2187 
2188            end do
2189          end do
2190        end do
2191 !      end ib loop:
2192      end do
2193 
2194 !    Now, check that the box defined by tmax was large enough :
2195 !    require the present and old lists to be the same
2196      done=0
2197 
2198      if(tmax>0)then
2199        done=1
2200        do ineighb=1,mneighb
2201 !        DEBUG
2202 !        write(std_out,'(5i5,f12.5)' )ineighb,list_neighb(ineighb,1:4,1),&
2203 !        &                                    sqrlength(ineighb)
2204 !        write(std_out,'(5i5)' )ineighb,list_neighb(ineighb,1:4,2)
2205 !        ENDDEBUG
2206          if( list_neighb(ineighb,1,1)/=list_neighb(ineighb,1,2) .or. &
2207 &         list_neighb(ineighb,2,1)/=list_neighb(ineighb,2,2) .or. &
2208 &         list_neighb(ineighb,3,1)/=list_neighb(ineighb,3,2) .or. &
2209 &         list_neighb(ineighb,4,1)/=list_neighb(ineighb,4,2)       )then
2210            done=0
2211          end if
2212        end do
2213      end if
2214 
2215 !    If done==1, then one can exit the loop : the correct list of
2216 !    neighbors is contained in list_neighb(1:neighb,1:4,1),
2217 !    with the first neighbor being the atom itself
2218      if(done==1)exit
2219 
2220 !    If the work is not done, while tmax==5, then there is a problem .
2221      if(tmax==5)then
2222        close(temp_unit)
2223        write(message, '(2a)' )&
2224 &       'Did not succeed to generate a reliable list of bonds ',&
2225 &       'since tmax is exceeded.'
2226        MSG_BUG(message)
2227      end if
2228 
2229 !    Copy the new list into the old list.
2230      list_neighb(1:mneighb,1:4,2)=list_neighb(1:mneighb,1:4,1)
2231 
2232 !    Loop on tmax (note that there are exit instruction inside the loop)
2233    end do
2234 
2235 
2236 
2237 !  Output the bond list
2238    do ineighb=2,mneighb
2239      ib=list_neighb(ineighb,1,1)
2240      length=sqrt(sqrlength(ineighb))
2241      write(message, '(a,a,a,a,3i2,t27,a,f10.5,a,f9.5,a)' )&
2242 &     '  ',trim(iden(ia)),' - ',trim(iden(ib)),&
2243 &     list_neighb(ineighb,2:4,1),'bond length is ',&
2244 &     length,' bohr  ( or ',Bohr_Ang*length,' Angst.)'
2245      call wrtout(temp_unit,message,'COLL')
2246    end do
2247 
2248 !  Output the angle list
2249    if(coordn>1)then
2250 
2251      write(message, '(a,a)' ) ch10,' Bond angles '
2252      call wrtout(temp_unit,message,'COLL')
2253 
2254      do ineighb=2,coordn
2255        do jneighb=ineighb+1,coordn+1
2256 
2257          ib=list_neighb(ineighb,1,1)
2258          ic=list_neighb(jneighb,1,1)
2259          do mu=1,3
2260            bab(mu)=xred(mu,ib)+dble(list_neighb(ineighb,1+mu,1))-xred(mu,ia)
2261            bac(mu)=xred(mu,ic)+dble(list_neighb(jneighb,1+mu,1))-xred(mu,ia)
2262          end do
2263          asq=rsdot(bab(1),bab(2),bab(3),bab(1),bab(2),bab(3),rmet)
2264          bsq=rsdot(bac(1),bac(2),bac(3),bac(1),bac(2),bac(3),rmet)
2265          adotb=rsdot(bab(1),bab(2),bab(3),bac(1),bac(2),bac(3),rmet)
2266          co=adotb/sqrt(asq*bsq)
2267          if( abs(co)-1.0d0 >= 0.0d0 )then
2268            if( abs(co)-1.0d0 <= 1.0d-12 )then
2269 !            Allows for a small numerical inaccuracy
2270              thdeg=0.0d0
2271              if(co < 0.0d0) thdeg=180.0d0
2272            else
2273              MSG_BUG('the evaluation of the angle is wrong.')
2274            end if
2275          else
2276            thdeg=acos(co)*180.d0*piinv
2277          end if
2278 
2279          write(message, '(a,a,3i2,a,a,a,a,3i2,t44,a,f13.5,a)' )&
2280 &         '  ',trim(iden(ib)),list_neighb(ineighb,2:4,1),' - ',&
2281 &         trim(iden(ia)),' - ',trim(iden(ic)),&
2282 &         list_neighb(jneighb,2:4,1),'bond angle is ',thdeg,' degrees '
2283          call wrtout(temp_unit,message,'COLL')
2284        end do
2285      end do
2286 
2287    end if
2288  end do !  End big ia loop:
2289 
2290  ABI_DEALLOCATE(iden)
2291  ABI_DEALLOCATE(list_neighb)
2292  ABI_DEALLOCATE(sqrlength)
2293 
2294  close(temp_unit)
2295 
2296  contains
2297 
2298    function rsdot(u1,u2,u3,v1,v2,v3,rmet)
2299 
2300 
2301 !This section has been created automatically by the script Abilint (TD).
2302 !Do not modify the following lines by hand.
2303 #undef ABI_FUNC
2304 #define ABI_FUNC 'rsdot'
2305 !End of the abilint section
2306 
2307    real(dp) :: rsdot
2308    real(dp),intent(in) :: u1,u2,u3,v1,v2,v3
2309    real(dp),intent(in) :: rmet(3,3)
2310    rsdot=rmet(1,1)*u1*v1+rmet(2,1)*u2*v1+&
2311 &   rmet(3,1)*u3*v1+rmet(1,2)*u1*v2+rmet(2,2)*u2*v2+&
2312 &   rmet(3,2)*u3*v2+rmet(1,3)*u1*v3+rmet(2,3)*u2*v3+rmet(3,3)*u3*v3
2313  end function rsdot
2314 
2315 end subroutine bonds_lgth_angles

m_geometry/chkdilatmx [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 chkdilatmx

FUNCTION

 Check whether the new rprimd does not give a too large number
 of plane waves, compared to the one booked for rprimd, taking
 into account the maximal dilatation dilatmx. Actually check whether
 the new Fermi sphere is inside the old one, dilated.

INPUTS

  chkdilatmx_ = if 1, will prevent to have any vector outside the Fermi sphere, possibly
       by rescaling (three times at most), and then stopping the execution
                if 0, simply send a warning, but continues execution
  dilatmx     = maximal dilatation factor (usually the input variable)
  rprimd      = new primitive vectors
  rprimd_orig = original primitive vectors (usually the input variable)

OUTPUT

  dilatmx_errmsg=Emptry string if calculation can continue.
            If the calculation cannot continue, dilatmx_errmsg will contain
            the message that should be reported in the output file.

            Client code should handle a possible problem with the following test:

              if (LEN_TRIM(dilatmx_errmsg) then
                dump dilatmx_errmsg to the main output file.
                handle_error
              end if

PARENTS

      driver,mover

CHILDREN

      matr3eigval,matr3inv

SOURCE

1519 subroutine chkdilatmx(chkdilatmx_,dilatmx,rprimd,rprimd_orig,dilatmx_errmsg)
1520 
1521 
1522 !This section has been created automatically by the script Abilint (TD).
1523 !Do not modify the following lines by hand.
1524 #undef ABI_FUNC
1525 #define ABI_FUNC 'chkdilatmx'
1526 !End of the abilint section
1527 
1528  implicit none
1529 
1530 !Arguments ------------------------------------
1531 !scalars
1532  integer,intent(in) :: chkdilatmx_
1533  real(dp),intent(in) :: dilatmx
1534  character(len=500),intent(out) :: dilatmx_errmsg
1535 !arrays
1536  real(dp),intent(inout) :: rprimd(3,3)
1537  real(dp),intent(in) :: rprimd_orig(3,3)
1538 
1539 !Local variables-------------------------------
1540 !scalars
1541  integer :: ii,jj,mu
1542  real(dp) :: alpha,dilatmx_new
1543 !arrays
1544  real(dp) :: eigval(3),gprimd_orig(3,3),met(3,3),old_to_new(3,3)
1545  character(len=500) :: message
1546 
1547 ! *************************************************************************
1548 
1549 !Generates gprimd
1550  call matr3inv(rprimd_orig,gprimd_orig)
1551 
1552 !Find the matrix that transform an original xcart to xred, then to the new xcart
1553  do mu=1,3
1554    old_to_new(mu,:)=rprimd(mu,1)*gprimd_orig(:,1)+&
1555 &   rprimd(mu,2)*gprimd_orig(:,2)+&
1556 &   rprimd(mu,3)*gprimd_orig(:,3)
1557  end do
1558 
1559 !The largest increase in length will be obtained thanks
1560 !to the diagonalization of the corresponding metric matrix :
1561 !it is the square root of its largest eigenvalue.
1562  do ii=1,3
1563    do jj=1,3
1564      met(ii,jj)=old_to_new(1,ii)*old_to_new(1,jj)+&
1565 &     old_to_new(2,ii)*old_to_new(2,jj)+&
1566 &     old_to_new(3,ii)*old_to_new(3,jj)
1567    end do
1568  end do
1569 
1570  call matr3eigval(eigval,met)
1571 
1572  dilatmx_new=sqrt(maxval(eigval(:)))
1573 
1574  dilatmx_errmsg = ""
1575  if(dilatmx_new>dilatmx+tol6)then
1576 
1577 ! MJV 2014 07 22: correct rprim to maximum jump allowed by dilatmx
1578 ! XG 20171011 : eigenvalues of "old_to_old" tensor are of course the unity !
1579 
1580    if(chkdilatmx_/=0)then
1581      alpha = (dilatmx - one) / (dilatmx_new - one)
1582 !    for safety, only 90 percent of max jump
1583      alpha = 0.9_dp * alpha
1584 
1585      rprimd = alpha * rprimd + (one - alpha) * rprimd_orig
1586 
1587      write(dilatmx_errmsg,'(3a,es16.6,4a,es16.6,2a,es16.6,a)')&
1588 &     'The new primitive vectors rprimd (an evolving quantity)',ch10,&
1589 &     'are too large with respect to the old rprimd and the accompanying dilatmx:',dilatmx,ch10,&
1590 &     'This large change of unit cell parameters is not allowed by the present value of dilatmx.',ch10,&
1591 &     'An adequate value would have been dilatmx_new=',dilatmx_new,ch10,&
1592 &     'Calculation continues with limited jump, by rescaling the projected move by the factor',alpha,'.'
1593    else
1594      write(message, '(3a,es16.6,2a,es16.6,2a)' )&
1595 &     'The new primitive vectors rprimd (an evolving quantity)',ch10,&
1596 &     'are too large, given the initial rprimd and the accompanying dilatmx:',dilatmx,ch10,&
1597 &     'An adequate value would have been dilatmx_new=',dilatmx_new,ch10,&
1598 &     'As chkdilatmx=0, assume experienced user. Execution will continue.'
1599      MSG_WARNING(message)
1600    end if
1601 
1602  end if
1603 
1604 end subroutine chkdilatmx

m_geometry/chkrprimd [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 chkrprimd

FUNCTION

 Test if {rprim,acell,rprimd} are consistent
 It means that rprimd can be reconstructed from the rprim and acell
 Output a message if is not the case

INPUTS

OUTPUT

  (only writing)

PARENTS

CHILDREN

SOURCE

1422 subroutine chkrprimd(acell,rprim,rprimd,iout)
1423 
1424 
1425 !This section has been created automatically by the script Abilint (TD).
1426 !Do not modify the following lines by hand.
1427 #undef ABI_FUNC
1428 #define ABI_FUNC 'chkrprimd'
1429 !End of the abilint section
1430 
1431 implicit none
1432 
1433 !Arguments ------------------------------------
1434 !scalars
1435 integer,intent(in) :: iout
1436 !arrays
1437 real(dp),intent(in) :: rprim(3,3)
1438 real(dp),intent(in) :: rprimd(3,3)
1439 real(dp),intent(in) :: acell(3)
1440 
1441 !Local variables-------------------------------
1442 !scalars
1443 integer :: ii,jj
1444 !arrays
1445 real(dp) :: rprimd_test(3,3)
1446 logical :: equal
1447 
1448 ! ***********************************************************
1449 
1450 !###########################################################
1451 !### 1. Compute rprimd from rprim and acell
1452  do ii=1,3
1453    do jj=1,3
1454      rprimd_test(ii,jj)=rprim(ii,jj)*acell(jj)
1455    end do
1456  end do
1457 
1458 
1459 !###########################################################
1460 !### 2. Compare rprimd and rprimd_test
1461 
1462  equal=.TRUE.
1463  do ii=1,3
1464    do jj=1,3
1465      if (abs(rprimd_test(ii,jj)-rprimd(ii,jj))>1.E-12) then
1466        equal=.FALSE.
1467      end if
1468    end do
1469  end do
1470 
1471  if (equal)then
1472    write(iout,*) 'chkrprimd: rprimd is consistent'
1473  else
1474    write(iout,*) 'chkrprimd: rprimd is NOT consistent ERROR'
1475  end if
1476 
1477 end subroutine chkrprimd

m_geometry/dist2 [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  dist2

FUNCTION

  Calculates the distance of v1 and v2 in a crystal by epeating the unit cell

INPUTS

  v1,v2
  rprimd: dimensions of the unit cell. if not given 1,0,0/0,1,0/0,0,1 is assumed
  option: 0 v1, v2 given in cartesian coordinates (default) / 1 v1,v2 given in reduced coordinates

OUTPUT

  dist2

PARENTS

  ioniondist

CHILDREN

SOURCE

2871 function dist2(v1,v2,rprimd,option)
2872 
2873 
2874 !This section has been created automatically by the script Abilint (TD).
2875 !Do not modify the following lines by hand.
2876 #undef ABI_FUNC
2877 #define ABI_FUNC 'dist2'
2878 !End of the abilint section
2879 
2880  implicit none
2881 
2882 !Arguments ------------------------------------
2883 !scalars
2884  integer,intent(in),optional :: option
2885  real(dp) :: dist2
2886 !arrays
2887  real(dp),intent(in),optional :: rprimd(3,3)
2888  real(dp),intent(in) :: v1(3),v2(3)
2889 
2890 !Local variables-------------------------------
2891 !scalars
2892  integer :: i1,i2,i3,opt,s1,s2,s3
2893  real(dp):: min2,norm2,ucvol
2894 !arrays
2895  integer :: limits(3)
2896  real(dp) :: corner(3),dred(3),dtot(3),dv(3),dwrap(3),sh(3)
2897  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2898  real(dp) :: vprimd(3,3)
2899 
2900 ! *************************************************************************
2901 
2902  if (.not.PRESENT(rprimd)) then
2903    vprimd=reshape((/1,0,0,  0,1,0,  0,0,1/),(/3,3/))
2904  else
2905    vprimd=rprimd
2906  end if
2907 
2908  call metric(gmet,gprimd,-1,rmet,vprimd,ucvol)
2909 
2910  dv(:)=v2(:)-v1(:)
2911 
2912 !If in cartesian coordinates, need to be transformed to reduced coordinates.
2913  opt=0
2914  if(present(option))then
2915    opt=option
2916  end if
2917  if(opt==0)then
2918    dred(:)=gprimd(1,:)*dv(1)+gprimd(2,:)*dv(2)+gprimd(3,:)*dv(3)
2919  else
2920    dred(:)=dv(:)
2921  end if
2922 
2923 !Wrap in the ]-1/2,1/2] interval
2924  call wrap2_pmhalf(dred(1),dwrap(1),sh(1))
2925  call wrap2_pmhalf(dred(2),dwrap(2),sh(2))
2926  call wrap2_pmhalf(dred(3),dwrap(3),sh(3))
2927 
2928 !Compute the limits of the parallelipipedic box that contains the Wigner-Seitz cell
2929 !The reduced coordinates of the corners of the Wigner-Seitz cell are computed (multiplied by two)
2930 !Then, the maximal values of these reduced coordinates are stored.
2931  limits(:)=0
2932  do s1=-1,1,2
2933    do s2=-1,1,2
2934      do s3=-1,1,2
2935        corner(:)=gmet(:,1)*s1*rmet(1,1)+gmet(:,2)*s2*rmet(2,2)+gmet(:,3)*s3*rmet(3,3)
2936        limits(1)=max(limits(1),ceiling(abs(corner(1))+tol14))
2937        limits(2)=max(limits(2),ceiling(abs(corner(2))+tol14))
2938        limits(3)=max(limits(3),ceiling(abs(corner(3))+tol14))
2939      end do
2940    end do
2941  end do
2942 
2943 !Use all relevant primitive real space lattice vectors to find the minimal difference vector
2944  min2=huge(zero)
2945  do i1=-limits(1),limits(1)
2946    do i2=-limits(2),limits(2)
2947      do i3=-limits(3),limits(3)
2948        dtot(1)=dwrap(1)+i1
2949        dtot(2)=dwrap(2)+i2
2950        dtot(3)=dwrap(3)+i3
2951        norm2=dtot(1)*rmet(1,1)*dtot(1)+dtot(2)*rmet(2,2)*dtot(2)+dtot(3)*rmet(3,3)*dtot(3)+&
2952 &       2*(dtot(1)*rmet(1,2)*dtot(2)+dtot(2)*rmet(2,3)*dtot(3)+dtot(3)*rmet(3,1)*dtot(1))
2953        min2=min(norm2,min2)
2954      end do
2955    end do
2956  end do
2957  dist2=sqrt(min2)
2958 
2959 end function dist2

m_geometry/fcart2fred [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 fcart2fred

FUNCTION

 Convert cartesian forces into reduced forces

INPUTS

  fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
  natom=Number of atoms in the unitary cell
  rprimd(3,3)=dimensional primitive

OUTPUT

  fred(3,natom)=symmetrized grtn = d(etotal)/d(xred)

NOTES

  Unlike fred, fcart has been corrected by enforcing
  the translational symmetry, namely that the sum of force
  on all atoms is zero.

PARENTS

      gstateimg,m_abihist,m_effective_potential,m_mep,mover,prec_simple
      pred_bfgs,pred_delocint,pred_lbfgs,pred_verlet,prtxfase

CHILDREN

SOURCE

1922 subroutine fcart2fred(fcart,fred,rprimd,natom)
1923 
1924 
1925 !This section has been created automatically by the script Abilint (TD).
1926 !Do not modify the following lines by hand.
1927 #undef ABI_FUNC
1928 #define ABI_FUNC 'fcart2fred'
1929 !End of the abilint section
1930 
1931  implicit none
1932 
1933 !Arguments ------------------------------------
1934 !scalars
1935  integer,intent(in) :: natom
1936 !arrays
1937  real(dp),intent(in) :: fcart(3,natom)
1938  real(dp),intent(out) :: fred(3,natom)
1939  real(dp),intent(in) :: rprimd(3,3)
1940 
1941 !Local variables-------------------------------
1942 !scalars
1943  integer :: iatom,mu
1944 
1945 ! *************************************************************************
1946 
1947 !MT, april 2012: the coding was not consistent with fred2fcart
1948  do iatom=1,natom
1949    do mu=1,3
1950      fred(mu,iatom)= - (rprimd(1,mu)*fcart(1,iatom)+&
1951 &     rprimd(2,mu)*fcart(2,iatom)+&
1952 &     rprimd(3,mu)*fcart(3,iatom))
1953    end do
1954  end do
1955 
1956 !Previous version
1957 !do iatom=1,natom
1958 !do mu=1,3
1959 !fred(mu,iatom)= - (rprimd(mu,1)*fcart(1,iatom)+&
1960 !&     rprimd(mu,2)*fcart(2,iatom)+&
1961 !&     rprimd(mu,3)*fcart(3,iatom))
1962 !end do
1963 !end do
1964 
1965 end subroutine fcart2fred

m_geometry/fixsym [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 fixsym

FUNCTION

 Using input indsym which tells which atoms are related by symmetry,
 check that iatfix consistently fixes (freezes) all atoms which are
 related by symmetry, i.e. that iatfix does not break symmetry.

INPUTS

 iatfix(3,natom)=integer array with 1 in every position for which
  the atom is to be kept fixed
  NOTE that this is not the input data structure for iatfix but it is
  the internal data structure used through most of the subroutines
 indsym(4,nsym,natom)=indirect indexing array for symmetrically related
  atoms; 4th element is label of symmetrically related atom
 natom=number of atoms
 nsym=number of symmetries (should be > 1 when this is called)

OUTPUT

  (only checking)

 NOTE
  Stops execution with an error message if iatfix breaks symmetry.

PARENTS

      gstate

CHILDREN

SOURCE

1157 subroutine fixsym(iatfix,indsym,natom,nsym)
1158 
1159 
1160 !This section has been created automatically by the script Abilint (TD).
1161 !Do not modify the following lines by hand.
1162 #undef ABI_FUNC
1163 #define ABI_FUNC 'fixsym'
1164 !End of the abilint section
1165 
1166  implicit none
1167 
1168 !Arguments ------------------------------------
1169 !scalars
1170  integer,intent(in) :: natom,nsym
1171 !arrays
1172  integer,intent(in) :: iatfix(3,natom),indsym(4,nsym,natom)
1173 
1174 !Local variables-------------------------------
1175 !scalars
1176  integer :: iatom,isym,jatom
1177  character(len=500) :: message
1178 
1179 ! *************************************************************************
1180 
1181  if (nsym > 1) then
1182    do iatom=1,natom
1183      do isym=1,nsym
1184 !      jatom is the label of a symmetrically related atom
1185        jatom=indsym(4,isym,iatom)
1186 !      Thus the atoms jatom and iatom must be fixed along the same directions
1187        if ( iatfix(1,jatom) /=  iatfix(1,iatom) .or. &
1188 &       iatfix(2,jatom) /=  iatfix(2,iatom) .or. &
1189 &       iatfix(3,jatom) /=  iatfix(3,iatom)       ) then
1190          write(message, '(a,i0,a,i0,7a)' )&
1191 &         'Atom number ',jatom,' is symmetrically  equivalent to atom number ',iatom,',',ch10,&
1192 &         'but according to iatfix, iatfixx, iatfixy and iatfixz, they',ch10,&
1193 &         'are not fixed along the same directions, which is forbidden.',ch10,&
1194 &         'Action: modify either the symmetry or iatfix(x,y,z) and resubmit.'
1195          MSG_ERROR(message)
1196        end if
1197      end do
1198    end do
1199  end if
1200 
1201 end subroutine fixsym

m_geometry/fred2fcart [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 fred2fcart

FUNCTION

 Convert reduced forces into cartesian forces

INPUTS

  fred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
  natom=Number of atoms in the unitary cell
  Favgz_null=TRUE if the average cartesian force has to be set to zero
             FALSE if it is set to zero only in x,y directions (not z)
  gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)

OUTPUT

  fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)

NOTES

    Unlike fred, fcart has been corrected by enforcing
    the translational symmetry, namely that the sum of force
    on all atoms is zero (except is a slab is used)

PARENTS

      forces,m_mep

CHILDREN

SOURCE

1843 subroutine fred2fcart(favg,Favgz_null,fcart,fred,gprimd,natom)
1844 
1845 
1846 !This section has been created automatically by the script Abilint (TD).
1847 !Do not modify the following lines by hand.
1848 #undef ABI_FUNC
1849 #define ABI_FUNC 'fred2fcart'
1850 !End of the abilint section
1851 
1852  implicit none
1853 
1854 !Arguments ------------------------------------
1855 !scalars
1856  integer,intent(in) :: natom
1857  logical :: Favgz_null
1858 !arrays
1859  real(dp),intent(out) :: fcart(3,natom)
1860  real(dp),intent(in) :: fred(3,natom)
1861  real(dp),intent(in) :: gprimd(3,3)
1862  real(dp),intent(out) :: favg(3)
1863 
1864 !Local variables-------------------------------
1865 !scalars
1866  integer :: iatom,mu
1867 
1868 ! *************************************************************************
1869 
1870 !Note conversion to cartesian coordinates (bohr) AND
1871 !negation to make a force out of a gradient
1872  favg(:)=zero
1873  do iatom=1,natom
1874    do mu=1,3
1875      fcart(mu,iatom)= - (gprimd(mu,1)*fred(1,iatom)+&
1876 &     gprimd(mu,2)*fred(2,iatom)+&
1877 &     gprimd(mu,3)*fred(3,iatom))
1878      favg(mu)=favg(mu)+fcart(mu,iatom)
1879    end do
1880  end do
1881 
1882 !Subtract off average force from each force component
1883 !to avoid spurious drifting of atoms across cell.
1884  favg(:)=favg(:)/dble(natom)
1885  if(.not.Favgz_null) favg(3)=zero
1886  do iatom=1,natom
1887    fcart(:,iatom)=fcart(:,iatom)-favg(:)
1888  end do
1889 
1890 end subroutine fred2fcart

m_geometry/getspinrot [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 getspinrot

FUNCTION

 From the symmetry matrix symrel_conv expressed in the coordinate system rprimd,
 compute the components of the spinor rotation matrix

INPUTS

 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 symrel_conv(3,3)=symmetry operation in real space in terms
  of primitive translations rprimd

OUTPUT

 spinrot(4)=components of the spinor rotation matrix :
  spinrot(1)=$\cos \phi / 2$
  spinrot(2)=$\sin \phi / 2 \times u_x$
  spinrot(3)=$\sin \phi / 2 \times u_y$
  spinrot(4)=$\sin \phi / 2 \times u_z$
  where $\phi$ is the angle of rotation, and
  $(u_x,u_y,u_z)$ is the normalized direction of the rotation axis

NOTES

 Only the proper part of the symmetry operation is taken into account:
 pure rotations, while the inversion part is taken away, if present.

 The whole collection of symmetry matrices is call symrel(3,3,nsym)
 symrel1 contains just one of those matrices symrel1(3,3)

PARENTS

      cg_rotate,m_crystal,wfconv

CHILDREN

      mati3det,matr3inv

SOURCE

796 subroutine getspinrot(rprimd,spinrot,symrel_conv)
797 
798 
799 !This section has been created automatically by the script Abilint (TD).
800 !Do not modify the following lines by hand.
801 #undef ABI_FUNC
802 #define ABI_FUNC 'getspinrot'
803 !End of the abilint section
804 
805  implicit none
806 
807 !Arguments ------------------------------------
808 !arrays
809  integer,intent(in) :: symrel_conv(3,3)
810  real(dp),intent(in) :: rprimd(3,3)
811  real(dp),intent(out) :: spinrot(4)
812 
813 !Local variables-------------------------------
814 !scalars
815  integer :: det,ii
816  real(dp) :: cos_phi,norminv,phi,scprod,sin_phi
817  !character(len=500) :: message
818 !arrays
819  integer :: identity(3,3),symrel1(3,3)
820  real(dp) :: axis(3),coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3)
821  real(dp) :: rprimd_invt(3,3),vecta(3),vectb(3),vectc(3)
822 
823 !**************************************************************************
824 
825  symrel1(:,:)=symrel_conv(:,:)
826 
827 !Compute determinant of the matrix
828  call mati3det(symrel1,det)
829 
830 !Produce a rotation from an improper symmetry
831  if(det==-1)symrel1(:,:)=-symrel1(:,:)
832 
833 !Test the possibility of the unit matrix
834  identity(:,:)=0
835  identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
836 
837  if( sum((symrel1(:,:)-identity(:,:))**2)/=0)then
838 
839 !  Transform symmetry matrix in the system defined by rprimd
840    call matr3inv(rprimd,rprimd_invt)
841    do ii=1,3
842      coord(:,ii)=rprimd_invt(ii,:)
843    end do
844    call matr3inv(coord,coordinvt)
845    do ii=1,3
846      matr1(:,ii)=symrel1(:,1)*coord(1,ii)+&
847 &     symrel1(:,2)*coord(2,ii)+&
848 &     symrel1(:,3)*coord(3,ii)
849    end do
850    do ii=1,3
851      matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+&
852 &     coordinvt(2,:)*matr1(2,ii)+&
853 &     coordinvt(3,:)*matr1(3,ii)
854    end do
855 
856 !  Find the eigenvector with unit eigenvalue of the
857 !  rotation matrix in cartesian coordinate, matr2
858 
859    matr1(:,:)=matr2(:,:)
860    matr1(1,1)=matr1(1,1)-one
861    matr1(2,2)=matr1(2,2)-one
862    matr1(3,3)=matr1(3,3)-one
863 
864 !  Compute the axis of rotation and the cos and sin of rotation angle
865    if(matr1(1,1)**2 + matr1(2,1)**2 + matr1(3,1)**2 < tol8 )then
866 !    The first direction is the axis
867      axis(1)=one ; axis(2)=zero ; axis(3)=zero
868      cos_phi=matr2(2,2)
869      sin_phi=matr2(3,2)
870    else if(matr1(1,2)**2 + matr1(2,2)**2 + matr1(3,2)**2 < tol8 )then
871 !    The second direction is the axis
872      axis(1)=zero ; axis(2)=one ; axis(3)=zero
873      cos_phi=matr2(3,3)
874      sin_phi=matr2(1,3)
875    else
876 !    In this case, try use the first and second vector to build the
877 !    rotation axis : compute their cross product
878      axis(1)=matr1(2,1)*matr1(3,2)-matr1(2,2)*matr1(3,1)
879      axis(2)=matr1(3,1)*matr1(1,2)-matr1(3,2)*matr1(1,1)
880      axis(3)=matr1(1,1)*matr1(2,2)-matr1(1,2)*matr1(2,1)
881 !    Then, try to normalize it
882      scprod=sum(axis(:)**2)
883      if(scprod<tol8)then
884 !      The first and second vectors were linearly dependent
885 !      Thus, use the first and third vectors
886        axis(1)=matr1(2,1)*matr1(3,3)-matr1(2,3)*matr1(3,1)
887        axis(2)=matr1(3,1)*matr1(1,3)-matr1(3,3)*matr1(1,1)
888        axis(3)=matr1(1,1)*matr1(2,3)-matr1(1,3)*matr1(2,1)
889 !      Normalize the vector
890        scprod=sum(axis(:)**2)
891        if(scprod<tol8)then
892          MSG_BUG('Cannot find the rotation axis.')
893        end if
894      end if
895      norminv=one/sqrt(scprod)
896      axis(:)=axis(:)*norminv
897 
898 !    Project the axis vector out of the first unit vector,
899 !    and renormalize the projected vector
900 !    (the first vector cannot be the axis, as tested before)
901      vecta(1)=one-axis(1)**2
902      vecta(2)=-axis(1)*axis(2)
903      vecta(3)=-axis(1)*axis(3)
904      scprod=sum(vecta(:)**2)
905      norminv=one/sqrt(scprod)
906      vecta(:)=vecta(:)*norminv
907 !    Rotate the vector A, to get vector B
908      vectb(:)=matr2(:,1)*vecta(1)+matr2(:,2)*vecta(2)+matr2(:,3)*vecta(3)
909 !    Get dot product of vectors A and B, giving cos of the rotation angle
910      cos_phi=sum(vecta(:)*vectb(:))
911 !    Compute the cross product of the axis and vector A
912      vectc(1)=axis(2)*vecta(3)-axis(3)*vecta(2)
913      vectc(2)=axis(3)*vecta(1)-axis(1)*vecta(3)
914      vectc(3)=axis(1)*vecta(2)-axis(2)*vecta(1)
915 !    Get dot product of vectors B and C, giving sin of the rotation angle
916      sin_phi=sum(vectb(:)*vectc(:))
917    end if
918 
919 !  Get the rotation angle, then the parameters of the spinor rotation
920 !  Here, treat possible inaccurate values of the cosine of phi
921    if(cos_phi>  one-tol8 )cos_phi=  one-tol8
922    if(cos_phi<-(one-tol8))cos_phi=-(one-tol8)
923    phi=acos(cos_phi)
924    if(sin_phi<zero)phi=-phi
925 !  Rectify the angle, such that its absolute values corresponds to
926 !  180, 120, 90, 60, or 0 degrees
927    phi=(nint(six*phi/pi))/six*pi
928 !  Compute components of the spinor matrix
929    spinrot(1)=cos(half*phi)
930    spinrot(2)=axis(1)*sin(half*phi)
931    spinrot(3)=axis(2)*sin(half*phi)
932    spinrot(4)=axis(3)*sin(half*phi)
933 
934  else
935 
936 !  Here, the case of the unit matrix
937    axis(:)=zero
938    phi=zero
939    spinrot(1)=one
940    spinrot(2)=zero
941    spinrot(3)=zero
942    spinrot(4)=zero
943 
944  end if ! the case of the identity matrix
945 
946 !DEBUG
947 !write(std_out,*)' getspinrot :'
948 !write(std_out,*)' symrel_conv =',symrel_conv(:,:)
949 !write(std_out,*)' symrel =',symrel1(:,:)
950 !write(std_out,*)' rprimd =',rprimd(:,:)
951 !write(std_out,*)' matr2 =',matr2(:,:)
952 !write(std_out,*)' matr1 =',matr1(:,:)
953 !write(std_out,*)' phi (degree)=',phi*180._dp/pi
954 !write(std_out,'(a,3d16.6)' )' axis=',axis(:)
955 !write(std_out,*)' vecta=',vecta(:)
956 !stop
957 !ENDDEBUG
958 
959 end subroutine getspinrot

m_geometry/ioniondist [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 ioniondist

FUNCTION

  Compute ion-ion distances

INPUTS

  natom= number of atoms in unit cell
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  xred(3,natom)=dimensionless reduced coordinates of atoms
  inm(natom,natom)=index (m,n) of the atom
  option= 1 output ion-ion distances / 2 output ordering of ion-ion
          distances / 3 output variables in varlist
          according to ion-ion distances * magnetic ordering
          magv magnetic ordering of atoms given als 1 and -1, if not
          given fm is assumed
  varlist=List of variables
  magv(natom)= magnetic ordering of atoms
  atp=atom on which the perturbation was done

OUTPUT

PARENTS

      pawuj_utils,shellstruct

CHILDREN

      prmat,wrtout

SOURCE

2733 subroutine ioniondist(natom,rprimd,xred,inm,option,varlist,magv,atp,prtvol)
2734 
2735 
2736 !This section has been created automatically by the script Abilint (TD).
2737 !Do not modify the following lines by hand.
2738 #undef ABI_FUNC
2739 #define ABI_FUNC 'ioniondist'
2740 !End of the abilint section
2741 
2742  implicit none
2743 
2744 !Arguments ------------------------------------
2745 !scalars
2746  integer,intent(in)              :: natom,option
2747  integer,intent(in),optional     :: atp                   !atom on which the perturbation was done
2748 !arrays
2749  real(dp),intent(in)             :: rprimd(3,3)
2750  real(dp),intent(in)             :: xred(3,natom)
2751  real(dp),intent(out)            :: inm(natom,natom)
2752  integer,intent(in),optional     :: magv(natom)
2753  real(dp),intent(in),optional    :: varlist(natom)
2754  integer,intent(in),optional     :: prtvol
2755 
2756 !Local variables-------------------------------
2757 !scalars
2758  integer                      :: iatom,jatom,katom,kdum,atpp,prtvoll
2759  !character(len=500)           :: message
2760 !arrays
2761  integer                      :: interq(natom)
2762  real(dp)                     :: hxcart(3,natom),distm(natom,natom)
2763  real(dp)                     :: magvv(natom)
2764 
2765 ! *************************************************************************
2766 
2767  hxcart=matmul(rprimd,xred)
2768  interq=(/(iatom,iatom=1,natom)/)
2769  inm=0
2770 
2771  if (present(magv)) then
2772    magvv=magv
2773  else
2774    magvv=(/ (1, iatom=1,natom)  /)
2775  end if
2776 
2777  if (present(atp)) then
2778    atpp=atp
2779  else
2780    atpp=1
2781  end if
2782 
2783  if (present(prtvol)) then
2784    prtvoll=prtvol
2785  else
2786    prtvoll=1
2787  end if
2788 
2789  if (option==3.and.(.not.present(varlist))) then
2790    call  wrtout(std_out,'ioniondist error: option=3 but no variable list provided for symmetrization','COLL')
2791    return
2792  end if
2793 
2794 
2795 !DEBUG
2796 !write(message, '(a,a)' ) ch10,' ioniondist start '
2797 !call wrtout(std_out,message,'COLL')
2798 !END DEBUG
2799 
2800  distm=0
2801  katom=atpp-1
2802  do iatom=1,natom
2803    katom=katom+1
2804    if (katom > natom) katom=1
2805    distm(iatom,iatom)=0
2806    do jatom=iatom,natom
2807      distm(iatom,jatom)=dist2(xred(:,katom),xred(:,jatom),rprimd,1)*magvv(katom)*magvv(jatom)
2808      distm(jatom,iatom)=distm(iatom,jatom)
2809    end do
2810  end do
2811 
2812  if (prtvoll>=3) then
2813    call  wrtout(std_out,'ioniondist: ionic distances:','COLL')
2814    call prmat(distm,natom,natom,natom,std_out)
2815  end if
2816 
2817  distm=anint(distm*10000_dp)/10000_dp           ! rounding needed else distm(iatom,jatom)/= distm(1,kdum) sometimes fails
2818 
2819  do iatom=1,natom
2820    if (option==1) then
2821      inm(iatom,:)=distm(iatom,:)
2822    else
2823      do jatom=iatom,natom
2824        kdum=1
2825        do while ( (kdum <= natom) .and. (distm(iatom,jatom)/= distm(1,kdum)) )
2826          kdum=kdum+1
2827        end do
2828        if (option==2) then
2829          inm(iatom,jatom)=interq(kdum)
2830        else if (option==3) then
2831          inm(iatom,jatom)=varlist(kdum)
2832        end if
2833        inm(jatom,iatom)=inm(iatom,jatom)
2834      end do
2835    end if
2836  end do
2837 
2838  if (prtvoll==2) then
2839    call wrtout(std_out,'ioniondist: symmetrized matrix:','COLL')
2840    call prmat(distm,1,natom,natom,std_out)
2841  else if (prtvoll>=3) then
2842    call wrtout(std_out,'ioniondist: symmetrized matrix:','COLL')
2843    call prmat(distm,natom,natom,natom,std_out)
2844  end if
2845 
2846 end subroutine ioniondist

m_geometry/littlegroup_pert [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 littlegroup_pert

FUNCTION

 If syuse==0 and abs(rfmeth)==2, determines the set of symmetries that leaves a perturbation invariant.
 (Actually, all symmetries that leaves a q-wavevector invariant should be used to reduce the number
 of k-points for all perturbations. Unfortunately, one has to take into account the sign reversal of the
 perturbation under the symmetry operations, which makes GS routines not usable for the respfn code.
 The intermediate choice was to select only those that keep also the perturbation invariant.
 Note that the wavevector of the perturbation must also be invariant,
 a translation vector in real space is NOT allowed ).

INPUTS

 gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
 idir=direction of the perturbation
 indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine symatm for definition (if nsym>1)
 iout=if non-zero, output on unit iout
 ipert=characteristics of the perturbation
 natom= number of atoms
 nsym=number of space group symmetries
 rfmeth =
   1 or -1 if non-stationary block
   2 or -2 if stationary block
   3 or -3 if third order derivatives
   positive if symmetries are used to set elements to zero whenever possible, negative to prevent this to happen.
 symq(4,2,nsym)= Table computed by littlegroup_q.
   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
 symafm(nsym)=(anti)ferromagnetic part of the symmetry operations
 symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
 symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
 syuse= flag to use the symmetries or not. If 0 usei it, if 1 do not use it.
 tnons(3,nsym)=nonsymmorphic translations of space group in terms
  of real space primitive translations (may be 0)
 [unit]=By default the routine writes to std_out and this is very annoying if we are inside a big loop.
   Use unit=dev_null or a negative integer to disable writing.

OUTPUT

 nsym1 =number of space group symmetries that leaves the perturbation invariant
 symaf1(nsym1)=(anti)ferromagnetic part of the corresponding symmetry operations
 symrl1(3,3,nsym1)=corresponding 3x3 matrices of the group symmetries (real space)
 tnons1(3,nsym1)=corresponding nonsymmorphic translations of space group in terms
   of real space primitive translations (may be 0)!!

PARENTS

      dfpt_looppert,get_npert_rbz,m_dvdb,read_gkk

CHILDREN

      stresssym,wrtout

SOURCE

3519 subroutine littlegroup_pert(gprimd,idir,indsym,iout,ipert,natom,nsym,nsym1, &
3520 &    rfmeth,symafm,symaf1,symq,symrec,symrel,symrl1,syuse,tnons,tnons1, &
3521 &    unit) ! Optional
3522 
3523 
3524 !This section has been created automatically by the script Abilint (TD).
3525 !Do not modify the following lines by hand.
3526 #undef ABI_FUNC
3527 #define ABI_FUNC 'littlegroup_pert'
3528 !End of the abilint section
3529 
3530  implicit none
3531 
3532 !Arguments -------------------------------
3533 !scalars
3534  integer,intent(in) :: idir,iout,ipert,natom,nsym,rfmeth,syuse
3535  integer,intent(in),optional :: unit
3536  integer,intent(out) :: nsym1
3537 !arrays
3538  integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),symq(4,2,nsym)
3539  integer,intent(in) :: symrec(3,3,nsym),symrel(3,3,nsym)
3540  integer,intent(out) :: symaf1(nsym),symrl1(3,3,nsym)
3541  real(dp),intent(in) :: gprimd(3,3),tnons(3,nsym)
3542  real(dp),intent(out) :: tnons1(3,nsym)
3543 
3544 !Local variables -------------------------
3545 !scalars
3546  integer :: idir1,ii,istr,isym,jj,nsym_test,tok,ount
3547  character(len=500) :: msg
3548 !arrays
3549  integer :: sym_test(3,3,2)
3550  real(dp) :: str_test(6)
3551 
3552 ! *********************************************************************
3553 
3554  ount = std_out; if (present(unit)) ount = unit
3555 
3556  nsym1=0
3557  if((ipert==natom+3 .or. ipert==natom+4) .and. syuse==0 .and. abs(rfmeth)==2) then
3558 !  Strain perturbation section
3559 !  Use ground state routine which symmetrizes cartesian stress as a quick
3560 !  and dirty test for the invariance of the strain (ipert,idir) under
3561 !  each candidate symmetry
3562 !  I am presently assuming that translations are acceptable because I dont
3563 !  see why not.
3564 
3565    istr=3*(ipert-natom-3)+idir
3566    nsym_test=2
3567 !  Store identity as first element for test
3568    sym_test(:,:,1)=0
3569    sym_test(1,1,1)=1; sym_test(2,2,1)=1; sym_test(3,3,1)=1
3570    do isym=1,nsym
3571      sym_test(:,:,2)=symrec(:,:,isym)
3572      str_test(:)=0.0_dp
3573      str_test(istr)=1.0_dp
3574      call stresssym(gprimd,nsym_test,str_test,sym_test)
3575      if(abs(str_test(istr)-1.0_dp)<tol8)then
3576 !      The test has been successful !
3577        nsym1=nsym1+1
3578        symaf1(nsym1)=symafm(isym)
3579        do ii=1,3
3580          tnons1(ii,nsym1)=tnons(ii,isym)
3581          do jj=1,3
3582            symrl1(ii,jj,nsym1)=symrel(ii,jj,isym)
3583          end do
3584        end do
3585      end if
3586    end do
3587 
3588  else if(ipert>natom .or. syuse/=0 .or. abs(rfmeth)/=2)then
3589 
3590 !  Not yet coded for d/dk or electric field perturbations
3591    nsym1=1
3592    do ii=1,3
3593      tnons1(ii,1)=0._dp
3594      symaf1(1)=1
3595      do jj=1,3
3596        symrl1(ii,jj,1)=0
3597        if(ii==jj)symrl1(ii,jj,1)=1
3598      end do
3599    end do
3600 
3601  else
3602 
3603    do isym=1,nsym
3604 !    Check that the symmetry operation preserves the wavevector
3605 !    (a translation is NOT allowed)
3606      if(symq(4,1,isym)==1 .and.&
3607 &     symq(1,1,isym)==0 .and.&
3608 &     symq(2,1,isym)==0 .and.&
3609 &     symq(3,1,isym)==0          )then
3610 !      Check that the symmetry operation preserves the atom
3611        if(ipert==indsym(4,isym,ipert))then
3612 !        Check if the direction is preserved
3613          tok=1
3614          do idir1=1,3
3615            if((idir1==idir.and.symrec(idir,idir1,isym)/=1) .or.&
3616 &           (idir1/=idir.and.symrec(idir,idir1,isym)/=0))then
3617              tok=0
3618            end if
3619          end do
3620          if(tok==1)then
3621 !          All the tests have been successful !
3622            nsym1=nsym1+1
3623            symaf1(nsym1)=symafm(isym)
3624            do ii=1,3
3625              tnons1(ii,nsym1)=tnons(ii,isym)
3626              do jj=1,3
3627                symrl1(ii,jj,nsym1)=symrel(ii,jj,isym)
3628              end do
3629            end do
3630          end if
3631 
3632        end if
3633      end if
3634    end do
3635  end if
3636 
3637  if (nsym1<1) then
3638    write(msg,'(a,i0,a)')&
3639 &   ' The number of selected symmetries should be > 0, while it is nsym=',nsym1,'.'
3640    MSG_BUG(msg)
3641  end if
3642 
3643  if (nsym1 /= 1) then
3644    if (iout /= ount .and. iout > 0) then
3645      write(msg,'(a,i5,a)')' Found ',nsym1,' symmetries that leave the perturbation invariant.'
3646      call wrtout(iout,msg,'COLL')
3647    end if
3648    write(msg,'(a,i5,a)')' littlegroup_pert : found ',nsym1,' symmetries that leave the perturbation invariant :'
3649    call wrtout(ount,msg,'COLL')
3650  else
3651    if (iout /= ount .and. iout > 0) then
3652      write(msg,'(a,a)')' The set of symmetries contains',' only one element for this perturbation.'
3653      call wrtout(iout,msg,'COLL')
3654    end if
3655    write(msg,'(a)')' littlegroup_pert : only one element in the set of symmetries for this perturbation :'
3656    call wrtout(ount,msg,'COLL')
3657  end if
3658 
3659  if (ount > 0) then
3660    do isym=1,nsym1
3661      write(msg, '(9i4)' )((symrl1(ii,jj,isym),ii=1,3),jj=1,3)
3662      call wrtout(ount,msg,'COLL')
3663    end do
3664  end if
3665 
3666 end subroutine littlegroup_pert

m_geometry/metric [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 metric

FUNCTION

 Compute first dimensional primitive translation vectors in reciprocal space
 gprimd from rprimd, and eventually writes out.
 Then, computes metrics for real and recip space rmet and gmet using length
 dimensional primitive translation vectors in columns of rprimd(3,3) and gprimd(3,3).
  gprimd is the inverse transpose of rprimd.
  i.e. $ rmet_{i,j}= \sum_k ( rprimd_{k,i}*rprimd_{k,j} )  $
       $ gmet_{i,j}= \sum_k ( gprimd_{k,i}*gprimd_{k,j} )  $
 Also computes unit cell volume ucvol in $\textrm{bohr}^3$

INPUTS

  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  iout=unit number of output file.  If iout<0, do not write output.

OUTPUT

  gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
  ucvol=unit cell volume ($\textrm{bohr}^{3}$).

PARENTS

      afterscfloop,bethe_salpeter,chkinp,clnup1,conducti_nc,conducti_paw
      conducti_paw_core,cut3d,d2frnl,dfpt_eltfrhar,dfpt_eltfrkin,dfpt_eltfrxc
      dfpt_looppert,dfpt_newvtr,dfpt_scfcv,dist2,elpolariz,emispec_paw,energy
      extrapwf,fftprof,finddistrproc,forces,forstr,get_npert_rbz,getkgrid
      hartre,hartrestr,ingeo,initaim,initberry,inkpts,inqpt,invacuum,invars2m
      ks_ddiago,linear_optics_paw,m_ab7_symmetry,m_crystal,m_cut3d,m_ddb
      m_dens,m_effective_potential,m_effective_potential_file,m_fft
      m_fft_prof,m_fit_data,m_hamiltonian,m_io_kss,m_ioarr,m_mep,m_pawpwij
      m_screening,m_tdep_latt,m_use_ga,m_vcoul,m_wfk,mag_constr,mag_constr_e
      memory_eval,mkcore_wvl,mlwfovlp_qp,moddiel,mpi_setup,mrgscr,newrho
      newvtr,nres2vres,odamix,optic,pawgrnl,prcref,prcref_PMA,pred_bfgs
      pred_delocint,pred_isothermal,pred_langevin,pred_lbfgs,pred_nose
      pred_srkna14,pred_verlet,prt_cif,prtefield,prtimg,psolver_rhohxc
      rhotoxc,scfcv,screening,setup1,setup_bse,setup_screening,setup_sigma
      sigma,smallprim,stress,strhar,symmetrize_rprimd,testkgrid,thmeig
      vdw_dftd2,vdw_dftd3,wrt_moldyn_netcdf,wvl_initro,xchybrid_ncpp_cc
      xfpack_vin2x,xfpack_x2vin

CHILDREN

      matr3inv,wrtout

SOURCE

1252 subroutine metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
1253 
1254 
1255 !This section has been created automatically by the script Abilint (TD).
1256 !Do not modify the following lines by hand.
1257 #undef ABI_FUNC
1258 #define ABI_FUNC 'metric'
1259 !End of the abilint section
1260 
1261  implicit none
1262 
1263 !Arguments ------------------------------------
1264 !scalars
1265  integer,intent(in) :: iout
1266  real(dp),intent(out) :: ucvol
1267 !arrays
1268  real(dp),intent(in) :: rprimd(3,3)
1269  real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1270 
1271 !Local variables-------------------------------
1272 !scalars
1273  integer :: nu
1274  character(len=500) :: message
1275 !arrays
1276  real(dp) :: angle(3)
1277 
1278 ! *************************************************************************
1279 
1280 !Compute unit cell volume
1281  ucvol=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
1282 & rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
1283 & rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
1284 
1285 !Check that the input primitive translations are not linearly dependent (and none is zero); i.e. ucvol~=0
1286 !Also ask that the mixed product is positive.
1287  if (abs(ucvol)<tol12) then
1288 !  write(std_out,*)"rprimd",rprimd,"ucvol",ucvol
1289    write(message,'(5a)')&
1290 &   'Input rprim and acell gives vanishing unit cell volume.',ch10,&
1291 &   'This indicates linear dependency between primitive lattice vectors',ch10,&
1292 &   'Action: correct either rprim or acell in input file.'
1293    MSG_ERROR(message)
1294  end if
1295  if (ucvol<zero)then
1296    write(message,'(2a,3(a,3es16.6,a),7a)')&
1297 &   'Current rprimd gives negative (R1xR2).R3 . ',ch10,&
1298 &   'Rprimd =',rprimd(:,1),ch10,&
1299 &   '        ',rprimd(:,2),ch10,&
1300 &   '        ',rprimd(:,3),ch10,&
1301 &   'Action: if the cell size and shape are fixed (optcell==0),',ch10,&
1302 &   '        exchange two of the input rprim vectors;',ch10,&
1303 &   '        if you are optimizing the cell size and shape (optcell/=0),',ch10,&
1304 &   '        maybe the move was too large, and you might try to decrease strprecon.'
1305    MSG_ERROR(message)
1306  end if
1307 
1308 !Generates gprimd
1309  call matr3inv(rprimd,gprimd)
1310 
1311 !Write out rprimd, gprimd and ucvol
1312  if (iout>=0) then
1313    write(message,'(2a)')' Real(R)+Recip(G) ','space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):'
1314    call wrtout(iout,message,'COLL')
1315    do nu=1,3
1316      write(message, '(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)' ) &
1317 &     'R(',nu,')=',rprimd(:,nu)+tol10,&
1318 &     'G(',nu,')=',gprimd(:,nu)+tol10
1319      call wrtout(iout,message,'COLL')
1320    end do
1321    write(message,'(a,1p,e15.7,a)') ' Unit cell volume ucvol=',ucvol+tol10,' bohr^3'
1322    call wrtout(iout,message,'COLL')
1323    call wrtout(std_out,message,'COLL')
1324  end if
1325 
1326 !Compute real space metric.
1327  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
1328 
1329 !Compute reciprocal space metric.
1330  gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
1331 
1332 !Write out the angles
1333  if (iout>=0) then
1334    angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0
1335    angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0
1336    angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0
1337    write(message, '(a,3es16.8,a)' )' Angles (23,13,12)=',angle(1:3),' degrees'
1338    call wrtout(iout,message,'COLL')
1339    call wrtout(std_out,message,'COLL')
1340  end if
1341 
1342 end subroutine metric

m_geometry/mkradim [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 mkradim

FUNCTION

  Not so trivial subroutine to make dimensionless real space
  primitive translations rprim(3,3) from dimensional rprimd(3).
  also make acell(3).

INPUTS

  rprimd(3,3)=dimensional real space primitive translations (bohr)
              where: rprimd(i,j)=rprim(i,j)*acell(j)

OUTPUT

  acell(3)=unit cell length scales (bohr)
  rprim(3,3)=dimensionless real space primitive translations

PARENTS

      gstate,gstateimg,ingeo,m_ddk,m_pimd,m_use_ga,pred_steepdesc
      predict_pimd,wvl_memory,xfpack_vin2x

CHILDREN

SOURCE

1370 subroutine mkradim(acell,rprim,rprimd)
1371 
1372 
1373 !This section has been created automatically by the script Abilint (TD).
1374 !Do not modify the following lines by hand.
1375 #undef ABI_FUNC
1376 #define ABI_FUNC 'mkradim'
1377 !End of the abilint section
1378 
1379  implicit none
1380 
1381 !Arguments ------------------------------------
1382 !arrays
1383  real(dp),intent(out) :: acell(3),rprim(3,3)
1384  real(dp),intent(in) :: rprimd(3,3)
1385 
1386 !Local variables-------------------------------
1387 !scalars
1388  integer :: ii
1389 
1390 ! *************************************************************************
1391 
1392 !Use a representation based on normalised rprim vectors
1393  do ii=1,3
1394    acell(ii)=sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2)
1395    rprim(:,ii)=rprimd(:,ii)/acell(ii)
1396  end do
1397 
1398 end subroutine mkradim

m_geometry/mkrdim [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 mkrdim

FUNCTION

  Trivial subroutine to make dimensional real space
  primitive translations from length scales acell(3)
  and dimensionless translations rprim(3,3).

INPUTS

  acell(3)=unit cell length scales (bohr)
  rprim(3,3)=dimensionless real space primitive translations

OUTPUT

  rprimd(3,3)=dimensional real space primitive translations (bohr)
              where: rprimd(i,j)=rprim(i,j)*acell(j)

PARENTS

      bethe_salpeter,dfpt_looppert,dfpt_symph,driver,finddistrproc
      get_npert_rbz,gstateimg,harmonic_thermo,ingeo,invars1,invars2m,m_ddb
      m_ifc,m_results_img,m_use_ga,memory_eval,mpi_setup,outvar_o_z,pred_bfgs
      pred_isothermal,pred_lbfgs,pred_steepdesc,pred_verlet,predict_pimd
      randomcellpos,screening,setup1,setup_bse,setup_screening,setup_sigma
      sigma,thmeig,wvl_setboxgeometry,xfpack_x2vin

CHILDREN

SOURCE

1636 subroutine mkrdim(acell,rprim,rprimd)
1637 
1638 
1639 !This section has been created automatically by the script Abilint (TD).
1640 !Do not modify the following lines by hand.
1641 #undef ABI_FUNC
1642 #define ABI_FUNC 'mkrdim'
1643 !End of the abilint section
1644 
1645  implicit none
1646 
1647 !Arguments ------------------------------------
1648 !arrays
1649  real(dp),intent(in) :: acell(3),rprim(3,3)
1650  real(dp),intent(out) :: rprimd(3,3)
1651 
1652 !Local variables-------------------------------
1653 !scalars
1654  integer :: ii,jj
1655 
1656 ! *************************************************************************
1657 
1658  do ii=1,3
1659    do jj=1,3
1660      rprimd(ii,jj)=rprim(ii,jj)*acell(jj)
1661    end do
1662  end do
1663 
1664 end subroutine mkrdim

m_geometry/normv_int_vector [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  normv_int_vector

FUNCTION

  Returns the norm of an integer 3D vector expressed in reduced coordinates.
  either in real or reciprocal space. In the later case the factor 2pi has
  to be included, due to the conventions used in abinit to define the reciprocal lattice.

INPUTS

OUTPUT

SOURCE

168 function normv_int_vector(xv,met,space) result(res)
169 
170 
171 !This section has been created automatically by the script Abilint (TD).
172 !Do not modify the following lines by hand.
173 #undef ABI_FUNC
174 #define ABI_FUNC 'normv_int_vector'
175 !End of the abilint section
176 
177  implicit none
178 
179 !Arguments ------------------------------------
180 !scalars
181  real(dp) :: res
182  character(len=1),intent(in) :: space
183 !arrays
184  real(dp),intent(in) :: met(3,3)
185  integer,intent(in) :: xv(3)
186 
187 ! *************************************************************************
188 
189  res =  ( xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3)  &
190 &  +two*( xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) )
191 
192  select case (space)
193  case ('r','R')
194    res=SQRT(res)
195  case ('g','G')
196    res=two_pi*SQRT(res)
197  case default
198    MSG_BUG('Wrong value for space')
199  end select
200 
201 end function normv_int_vector

m_geometry/normv_int_vector_array [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  normv_int_vector_array

FUNCTION

  Returns the norm of an array of integer 3D vectors expressed in reduced coordinates.
  either in real or reciprocal space. In the later case the factor 2pi has
  to be included, due to the conventions used in abinit to define the reciprocal lattice.

INPUTS

OUTPUT

SOURCE

221 function normv_int_vector_array(xv,met,space) result(res)
222 
223 
224 !This section has been created automatically by the script Abilint (TD).
225 !Do not modify the following lines by hand.
226 #undef ABI_FUNC
227 #define ABI_FUNC 'normv_int_vector_array'
228 !End of the abilint section
229 
230  implicit none
231 
232 !Arguments ------------------------------------
233 !scalars
234  character(len=1),intent(in) :: space
235 !arrays
236  real(dp),intent(in) :: met(3,3)
237  integer,intent(in) :: xv(:,:)
238  !this awful trick is needed to avoid problems with abilint
239  real(dp) :: res(SIZE(xv(1,:)))
240  !real(dp) :: res(SIZE(xv,DIM=2))
241 
242 ! *************************************************************************
243 
244  res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:)  &
245 &     +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) )
246 
247  select case (space)
248  case ('r','R')
249    res(:)=SQRT(res(:))
250  case ('g','G')
251    res(:)=two_pi*SQRT(res(:))
252  case default
253    MSG_BUG('Wrong value for space')
254  end select
255 
256 end function normv_int_vector_array

m_geometry/normv_rdp_vector [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 normv_rdp_vector

FUNCTION

 Compute the norm of a vector expressed in reduced coordinates using the metric met.
 The result is multiplied by 2pi in case of a vector in reciprocal space
 to take into account the correct normalisation of the reciprocal lattice vectors

INPUTS

  xv(3)=Vector in reduced coordinates
  met(3,3)=Metric tensor
  space=Character defining whether we are working in real (r|R) or reciprocal space (g|G)

OUTPUT

  normv_rdp_vector=norm of xv

NOTES

  The routine is able to deal both with a single vector as well as arrays of vectors.
  Versions for integer and real vectors are provided.

PARENTS

CHILDREN

SOURCE

116 function normv_rdp_vector(xv,met,space) result(res)
117 
118 
119 !This section has been created automatically by the script Abilint (TD).
120 !Do not modify the following lines by hand.
121 #undef ABI_FUNC
122 #define ABI_FUNC 'normv_rdp_vector'
123 !End of the abilint section
124 
125  implicit none
126 
127 !Arguments ------------------------------------
128 !scalars
129  real(dp) :: res
130  character(len=1),intent(in) :: space
131 !arrays
132  real(dp),intent(in) :: met(3,3),xv(3)
133 
134 ! *************************************************************************
135 
136  res =  (xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3)  &
137 &  +two*(xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) )
138 
139  select case (space)
140  case ('r','R')
141    res=SQRT(res)
142  case ('g','G')
143    res=two_pi*SQRT(res)
144  case default
145    MSG_BUG('Wrong value for space')
146  end select
147 
148 end function normv_rdp_vector

m_geometry/normv_rdp_vector_array [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  normv_rdp_vector_array

FUNCTION

  Returns the norm of an array of real 3D vectors expressed in reduced coordinates.
  either in real or reciprocal space. In the later case the factor 2pi has
  to be included, due to the conventions used in abinit to define the reciprocal lattice.

INPUTS

OUTPUT

SOURCE

276 function normv_rdp_vector_array(xv,met,space) result(res)
277 
278 
279 !This section has been created automatically by the script Abilint (TD).
280 !Do not modify the following lines by hand.
281 #undef ABI_FUNC
282 #define ABI_FUNC 'normv_rdp_vector_array'
283 !End of the abilint section
284 
285  implicit none
286 
287 !Arguments ------------------------------------
288 !scalars
289  character(len=1),intent(in) :: space
290 !arrays
291  real(dp),intent(in) :: met(3,3)
292  real(dp),intent(in) :: xv(:,:)
293  !this awful trick is needed to avoid problems with abilint
294  real(dp) :: res(SIZE(xv(1,:)))
295  !real(dp) :: res(SIZE(xv,DIM=2))
296 
297 ! *************************************************************************
298 
299  res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:)  &
300 &     +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) )
301 
302  select case (space)
303  case ('r','R')
304    res(:)=SQRT(res(:))
305  case ('g','G')
306    res(:)=two_pi*SQRT(res)
307  case default
308    MSG_BUG('Wrong value for space')
309  end select
310 
311 end function normv_rdp_vector_array

m_geometry/phdispl_cart2red [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  phdispl_cart2red

FUNCTION

  Calculates the displacement vectors for all branches in reduced coordinates.
  $ displ_red = displ_cart \cdot gprimd $ for each phonon branch.

INPUTS

  natom=Number of atoms.
  gprimd(3,3)Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  displ_cart(2,3*natom,3*natom)=Phonon displacement in Cartesian coordinates.

OUTPUT

  displ_red(2,3*natom,3*natom)=Phonon displacement in reduded coordinates.

PARENTS

      get_tau_k,m_ddb,m_ifc,mka2f,mkph_linwid,read_gkk

CHILDREN

SOURCE

702 subroutine phdispl_cart2red(natom,gprimd,displ_cart,displ_red)
703 
704 
705 !This section has been created automatically by the script Abilint (TD).
706 !Do not modify the following lines by hand.
707 #undef ABI_FUNC
708 #define ABI_FUNC 'phdispl_cart2red'
709 !End of the abilint section
710 
711  implicit none
712 
713 !Arguments ------------------------------------
714 !scalars
715  integer,intent(in) :: natom
716 !arrays
717  real(dp),intent(in) :: gprimd(3,3)
718  real(dp),intent(in) :: displ_cart(2,3*natom,3*natom)
719  real(dp),intent(out) :: displ_red(2,3*natom,3*natom)
720 
721 !Local variables-------------------------
722 !scalars
723  integer :: nbranch,jbranch,iatom,idir,ibranch,kdir,k1
724 
725 ! *************************************************************************
726 
727  displ_red = zero
728 
729  nbranch=3*natom
730 
731  do jbranch=1,nbranch
732    !
733    do iatom=1,natom
734      do idir=1,3
735        ibranch=idir+3*(iatom-1)
736        do kdir=1,3
737          k1 = kdir+3*(iatom-1)
738          ! WARNING: could be non-transpose of rprimd matrix : to be checked.
739          ! 23 june 2004: rprimd becomes gprimd
740          ! could be gprim and then multiply by acell...
741          ! Nope, checked and ok with gprimd 24 jun 2004
742          displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + &
743 &         gprimd(kdir,idir) * displ_cart(1,k1,jbranch)
744 
745          displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + &
746 &         gprimd(kdir,idir) * displ_cart(2,k1,jbranch)
747 
748        end do !kdir
749      end do !idir
750    end do !iatom
751    !
752  end do !jbranch
753 
754 end subroutine phdispl_cart2red

m_geometry/randomcellpos [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  randomcellpos

FUNCTION

  This subroutine creates a unit cell with random atomic positions. It is
  assumed that the cell parameters are given and fixed. Several methods are
  used to generate the cell.

INPUTS

 natom=number of atoms
 npsp=number of pseudopotentials (needed for the dimension of znucl)
 ntypat=number of type of atoms
 random_atpos=input variable
   0 no generation of random atomic potision
   1 completely random atomic potisions
   2 random atomic positions, avoiding too close atoms
     (prevent coming closer than a fraction of the sum of covalent radii)
   3 same than 2 but also generates the rprim and acell randomly
    within some given ranges (angles between 50 and 130)
 ratsph(1:ntypat)=radius of the atomic sphere
 rprimd(3,3)=dimensional primitive translations in real space (bohr)
 typat(1:natom)= input variable giving the type of each atom
 znucl(1:npsp)=nuclear number of atom as specified in psp file

OUTPUT

 xred(3,natom)=reduced dimensionless atomic coordinates

SIDE EFFECTS

NOTES

PARENTS

      ingeo

CHILDREN

      atomdata_from_znucl,mkrdim,xred2xcart

SOURCE

2358 subroutine randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd,typat,xred,znucl,acell)
2359 
2360 
2361 !This section has been created automatically by the script Abilint (TD).
2362 !Do not modify the following lines by hand.
2363 #undef ABI_FUNC
2364 #define ABI_FUNC 'randomcellpos'
2365 !End of the abilint section
2366 
2367  implicit none
2368 
2369 !Arguments ------------------------------------
2370 !scalars
2371  integer,intent(in) :: natom,npsp,ntypat,random_atpos
2372 !arrays
2373  integer, intent(in)   :: typat(natom)
2374  real(dp),intent(in)   :: ratsph(ntypat)
2375  real(dp), intent(inout)  :: rprim(3,3)
2376  real(dp), intent(inout)  :: rprimd(3,3)
2377  real(dp), intent(inout) :: xred(3,natom)
2378  real(dp), intent(in) :: znucl(npsp)
2379  real(dp), intent(inout) :: acell(3)
2380 
2381 !Local variables-------------------------------
2382  integer ::   iatom=0,ii,idum=-20
2383  real(dp) ::  rij(3), rijd(3), radiuscovi, radiuscovj, dist, rati, ratj, angdeg(3)
2384  real(dp) ::  cosang,aa,cc,a2
2385  character(len=500) :: message
2386  type(atomdata_t) :: atom
2387 
2388 ! *************************************************************************
2389 
2390 !DEBUG
2391 !For the time being, print rprimd to keep it as an argument, in spite of abirule checking.
2392 !write (std_out,*) ' randomcellpos : enter'
2393 !write(std_out,*)' rprimd=',rprimd
2394 !write(std_out,*)' znucl=',znucl
2395 !write(std_out,*)' typat=',typat
2396 !write(std_out,*)' random_atpos=',random_atpos
2397 !ENDDEBUG
2398 
2399  if(random_atpos==2 .and. npsp/=ntypat)then
2400    write(message, '(a,i5,2a,i5,a,i5,4a)' )&
2401 &   'Input variable random_atpos= ',random_atpos,ch10,&
2402 &   'However, the number of pseudopotentials ',npsp,', is not equal to the number of type of atoms ',ntypat,ch10,&
2403 &   'The use of alchemical mixing cannot be combined with the constraint based on the mixing of covalent radii.',ch10,&
2404 &   'Action: switch to another value of random_atpos.'
2405    MSG_ERROR(message)
2406  end if
2407 
2408 !random_atpos = 0   Default value, no random initialisation
2409 !random_atpos = 1   Fully random (Is it really useful ???)
2410 !random_atpos = 2   Random, but the sum of the two covalent radii is
2411 !less than the interatomic distance
2412 !random_atpos = 3   Random, but the sum of the two (other type of)
2413 !radii is less than the interatomic distance
2414 !random_atpos = 4   Random, but the sum of the two pseudopotential
2415 !radii is less than the interatomic distance
2416 !random_atpos = 5   Random, but the interatomic distance must be bigger
2417 !than the sum of
2418 !some input variable (well, instead of defining a new variable, why
2419 !not use ratsph ?)
2420 !Right now we are not using a factor for the tested distance.. something to be done, after a new variable has been defined
2421 
2422  if (random_atpos /= 0) then
2423    select case (random_atpos)
2424    case (1)
2425      do ii=1,natom
2426        xred(1,ii)=uniformrandom(idum)
2427        xred(2,ii)=uniformrandom(idum)
2428        xred(3,ii)=uniformrandom(idum)
2429      end do
2430    case (2)
2431      iatom=0
2432      do
2433        iatom=iatom+1
2434        xred(1,iatom)=uniformrandom(idum)
2435        xred(2,iatom)=uniformrandom(idum)
2436        xred(3,iatom)=uniformrandom(idum)
2437        call atomdata_from_znucl(atom,znucl(typat(iatom)))
2438        radiuscovi = atom%rcov
2439        do ii=1,iatom-1
2440          rij=xred(:,iatom)-xred(:,ii)
2441 !          periodic boundary conditions
2442          rij = rij - 0.5
2443          rij = rij - anint (rij)
2444 !          coming back to cube between (0,1)
2445          rij = rij + 0.5
2446 !          convert reduced coordinates to cartesian coordinates
2447          call xred2xcart(1,rprimd,rijd,rij)
2448          dist=dot_product(rijd,rijd)
2449          call atomdata_from_znucl(atom,znucl(typat(ii)))
2450          radiuscovj = atom%rcov
2451          if (dist<(radiuscovj+radiuscovi)) then
2452            iatom = iatom -1
2453            EXIT
2454          end if
2455        end do
2456        if (iatom>=natom) EXIT
2457      end do
2458    case(3)
2459      iatom=0
2460      do
2461        iatom=iatom+1
2462        xred(1,iatom)=uniformrandom(idum)
2463        xred(2,iatom)=uniformrandom(idum)
2464        xred(3,iatom)=uniformrandom(idum)
2465        call atomdata_from_znucl(atom,znucl(typat(iatom)))
2466        radiuscovi = atom%rcov
2467        do ii=1,iatom-1
2468          rij=xred(:,iatom)-xred(:,ii)
2469 !          periodic boundary conditions
2470          rij = rij - 0.5
2471          rij = rij - anint (rij)
2472 !          coming back to cube between (0,1)
2473          rij = rij + 0.5
2474 !          convert reduced coordinates to cartesian coordinates
2475          call xred2xcart(1,rprimd,rijd,rij)
2476          dist=dot_product(rijd,rijd)
2477          call atomdata_from_znucl(atom,znucl(typat(ii)))
2478          radiuscovj = atom%rcov
2479          if (dist<(radiuscovj+radiuscovi)) then
2480            iatom = iatom -1
2481            EXIT
2482          end if
2483        end do
2484        if (iatom>=natom) EXIT
2485      end do
2486      do ii=1,3
2487 !        generates cells with angles between 60 and 120 degrees
2488        angdeg(ii)=60_dp+uniformrandom(idum)*60.0_dp
2489      end do
2490      if (angdeg(1)+angdeg(2)+angdeg(3)>360._dp) then
2491        angdeg(3)=360._dp-angdeg(1)-angdeg(2)
2492      end if
2493 !      check if angles are between the limits and create rprim
2494      if( abs(angdeg(1)-angdeg(2))<tol12 .and. &
2495 &     abs(angdeg(2)-angdeg(3))<tol12 .and. &
2496 &     abs(angdeg(1)-90._dp)+abs(angdeg(2)-90._dp)+abs(angdeg(3)-90._dp)>tol12 )then
2497 !        Treat the case of equal angles (except all right angles) :
2498 !        generates trigonal symmetry wrt third axis
2499        cosang=cos(pi*angdeg(1)/180.0_dp)
2500        a2=2.0_dp/3.0_dp*(1.0_dp-cosang)
2501        aa=sqrt(a2)
2502        cc=sqrt(1.0_dp-a2)
2503        rprim(1,1)=aa        ; rprim(2,1)=0.0_dp                 ; rprim(3,1)=cc
2504        rprim(1,2)=-0.5_dp*aa ; rprim(2,2)= sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,2)=cc
2505        rprim(1,3)=-0.5_dp*aa ; rprim(2,3)=-sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,3)=cc
2506 !        DEBUG
2507 !        write(std_out,*)' ingeo : angdeg=',angdeg(1:3)
2508 !        write(std_out,*)' ingeo : aa,cc=',aa,cc
2509 !        ENDDEBUG
2510      else
2511 !        Treat all the other cases
2512        rprim(:,:)=0.0_dp
2513        rprim(1,1)=1.0_dp
2514        rprim(1,2)=cos(pi*angdeg(3)/180.0_dp)
2515        rprim(2,2)=sin(pi*angdeg(3)/180.0_dp)
2516        rprim(1,3)=cos(pi*angdeg(2)/180.0_dp)
2517        rprim(2,3)=(cos(pi*angdeg(1)/180.0_dp)-rprim(1,2)*rprim(1,3))/rprim(2,2)
2518        rprim(3,3)=sqrt(1.0_dp-rprim(1,3)**2-rprim(2,3)**2)
2519      end if
2520 !      generate acell
2521      aa=zero
2522      do ii=1,npsp
2523        aa=znucl(ii)
2524      end do
2525      do ii=1,3
2526        acell(ii)=aa+uniformrandom(idum)*4.0
2527      end do
2528      call mkrdim(acell,rprim,rprimd)
2529    case(4)
2530      write(std_out,*) 'Not implemented yet'
2531    case(5)
2532      iatom=0
2533      do
2534        iatom=iatom+1
2535        xred(1,iatom)=uniformrandom(idum)
2536        xred(2,iatom)=uniformrandom(idum)
2537        xred(3,iatom)=uniformrandom(idum)
2538        rati=ratsph(typat(iatom))
2539        do ii=1,iatom-1
2540          ratj=ratsph(typat(ii))
2541 !          apply periodic boundary conditions
2542          rij=(xred(:,iatom)-xred(:,ii))-0.5
2543          rij = rij - ANINT ( rij )
2544          rij = rij + 0.5
2545          call xred2xcart(natom,rprimd,rijd,rij)
2546          dist=dot_product(rijd,rijd)
2547          if (dist<(rati+ratj)) EXIT
2548        end do
2549        if (iatom==natom) EXIT
2550        if (ii<(iatom-1)) iatom=iatom-1
2551      end do
2552    end select
2553  end if
2554 
2555 end subroutine randomcellpos

m_geometry/remove_inversion [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 remove_inversion

FUNCTION

  Remove the inversion symmetry from a symmetry set as well
  all the improper rotations (if present)

INPUTS

  nsym=initial number of symmetries
  symrel(3,3,nsym)=Initial set of symmetry operarations in real space
  tnons(3,nsym)=Initial fractional translations

OUTPUT

  nsym_out=Number of symmetries in the set without improper rotation
  symrel_out(:,:) [pointer] = output symmetries without improper rotations
  tnons_out(:) [pointer] = fractional translations associated to symrel_out
  pinv=-1 if the inversion has been removed, 1 otherwise

NOTES

  Note the use of pointers, memory is allocated inside the procedure and passed back
  to the caller. Thus memory deallocation is relegated to the caller. To be on the safe side
  the pointers should be nullified before entering.

PARENTS

      m_crystal,m_io_kss,outkss

CHILDREN

      set2unit,symdet,wrtout

SOURCE

2994 subroutine remove_inversion(nsym,symrel,tnons,nsym_out,symrel_out,tnons_out,pinv)
2995 
2996 
2997 !This section has been created automatically by the script Abilint (TD).
2998 !Do not modify the following lines by hand.
2999 #undef ABI_FUNC
3000 #define ABI_FUNC 'remove_inversion'
3001 !End of the abilint section
3002 
3003  implicit none
3004 
3005 !Arguments ------------------------------------
3006 !scalars
3007  integer,intent(in) :: nsym
3008  integer,intent(out) :: nsym_out,pinv
3009 !arrays
3010  integer,intent(in) :: symrel(3,3,nsym)
3011  integer,pointer :: symrel_out(:,:,:)
3012  real(dp),intent(in) :: tnons(3,nsym)
3013  real(dp),pointer :: tnons_out(:,:)
3014 
3015 !Local variables-------------------------------
3016 !scalars
3017  integer :: is,is2,is_discarded,is_inv,is_retained,nsym2
3018  logical :: found
3019  character(len=500) :: msg
3020 !arrays
3021  integer :: determinant(nsym),inversion(3,3),symrel2(3,3,nsym)
3022  real(dp) :: dtnons(3),tnons2(3,nsym)
3023 
3024 ! *********************************************************************
3025 
3026  MSG_WARNING('Removing inversion related symmetrie from initial set')
3027 
3028  ! Find the occurence of the inversion symmetry.
3029  call set2unit(inversion) ; inversion=-inversion
3030 
3031  is_inv=0; found=.FALSE.
3032  do while (is_inv<nsym .and. .not.found)
3033    is_inv=is_inv+1; found=ALL(symrel(:,:,is_inv)==inversion)
3034  end do
3035  if (found) then
3036    write(msg,'(a,i3)')' The inversion is symmetry operation no. ',is_inv
3037  else
3038    write(msg,'(a)')' The inversion was not found in the symmetries list.'
3039  end if
3040  call wrtout(std_out,msg,'COLL')
3041 
3042  ! Find the symmetries that are related through the inversion symmetry
3043  call symdet(determinant,nsym,symrel)
3044  nsym2=0
3045  do is=1,nsym-1
3046    do is2=is+1,nsym
3047 
3048      dtnons(:)=tnons(:,is2)-tnons(:,is)-tnons(:,is_inv)
3049      found=ALL(symrel(:,:,is)==-symrel(:,:,is2)).and.isinteger(dtnons,tol8)
3050 
3051      if (found) then
3052        nsym2=nsym2+1
3053        ! Retain symmetries with positive determinant
3054        if (ALL(tnons(:,is2)<tol8).and.ALL(tnons(:,is)<tol8)) then
3055          is_retained=is2 ; is_discarded=is
3056          if (determinant(is)==1) then
3057            is_retained=is  ; is_discarded=is2
3058          end if
3059        else if (ALL(tnons(:,is2)<tol8)) then
3060          is_retained=is2 ; is_discarded=is
3061        else
3062          is_retained=is ;  is_discarded=is2
3063        end if
3064 
3065        symrel2(:,:,nsym2)=symrel(:,:,is_retained)
3066        tnons2   (:,nsym2)=tnons   (:,is_retained)
3067        write(msg,'(a,i3,a,i3,3a,i3,a)')&
3068 &       ' Symmetry operations no. ',is,' and no. ',is2,&
3069 &       ' are related through the inversion.',ch10,&
3070 &       ' Symmetry operation no. ',is_discarded,' will be suppressed.'
3071        call wrtout(std_out,msg,'COLL')
3072      end if ! found
3073 
3074    end do !is2
3075  end do !is
3076 
3077  if (nsym2/=(nsym/2).or.nsym==1) then
3078    call wrtout(std_out, ' Program uses the original set of symmetries ', 'COLL')
3079    nsym_out=nsym
3080    ABI_ALLOCATE(symrel_out,(3,3,nsym))
3081    ABI_ALLOCATE(tnons_out,(3,nsym))
3082    symrel_out(:,:,:)=symrel(:,:,1:nsym)
3083    tnons_out(:,:)=tnons(:,1:nsym)
3084    pinv=1
3085  else
3086    write(msg,'(a)')' Inversion related operations have been suppressed from symmetries list.'
3087    call wrtout(std_out,msg,'COLL')
3088    nsym_out=nsym2
3089    ABI_ALLOCATE(symrel_out,(3,3,nsym2))
3090    ABI_ALLOCATE(tnons_out,(3,nsym2))
3091    symrel_out(:,:,:)=symrel2(:,:,1:nsym2)
3092    tnons_out(:,:)=tnons(:,1:nsym2)
3093    pinv=-1
3094  end if
3095 
3096 end subroutine remove_inversion

m_geometry/rotmat [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 rotmat

FUNCTION

 Finds the rotation matrix.

INPUTS

  xaxis(3)= vectors defining the x axis
  zaxis(3)= vectors defining the z axis

OUTPUT

  inversion_flag = flag that indicates that an inversion operation
   on the coordinate system should be done
  umat(3,3)= matrix that rotates the x=(1 0 0) and z=(0 0 1) to the new
   values defined in xaxis and zaxis

NOTES

 Here I set that the axe x is originally at the 1 0 0 direction and z is originally 0 0 1.
 So calling rotmat(x',z') will find the rotation
 matrix for the case in which we rotate the x and z
 axes from their default values to x' and z'.

PARENTS

      mlwfovlp_ylmfac,mlwfovlp_ylmfar

CHILDREN

      wrtout

SOURCE

1046 subroutine rotmat(xaxis,zaxis,inversion_flag,umat)
1047 
1048 
1049 !This section has been created automatically by the script Abilint (TD).
1050 !Do not modify the following lines by hand.
1051 #undef ABI_FUNC
1052 #define ABI_FUNC 'rotmat'
1053 !End of the abilint section
1054 
1055  implicit none
1056 
1057 !Arguments ------------------------------------
1058 !scalars
1059  integer,intent(out) :: inversion_flag
1060 !arrays
1061  real(dp),intent(in) :: xaxis(3),zaxis(3)
1062  real(dp),intent(out) :: umat(3,3)
1063 
1064 !Local variables-------------------------------
1065 !scalars
1066  real(dp) :: cosine,xmod,zmod
1067  character(len=500) :: message
1068 !arrays
1069  real(dp) :: yaxis(3)
1070 
1071 ! *************************************************************************
1072 
1073  xmod = sqrt(xaxis(1)**2 + xaxis(2)**2 + xaxis(3)**2)
1074  zmod = sqrt(zaxis(1)**2 + zaxis(2)**2 + zaxis(3)**2)
1075 
1076  if(xmod < 1.d-8)then
1077    write(message,'(a,a,a,i6)')&
1078 &   'The module of the xaxis should be greater than 1.d-8,',ch10,&
1079 &   'however, |xaxis|=',xmod
1080    MSG_BUG(message)
1081  end if
1082 
1083  if(zmod < 1.d-8)then
1084    write(message,'(a,a,a,i6)')&
1085 &   'The module of the zaxis should be greater than 1.d-8,',ch10,&
1086 &   'however, |zaxis|=',zmod
1087    MSG_ERROR(message)
1088  end if
1089 
1090 !verify that both axis are perpendicular
1091  cosine = (xaxis(1)*zaxis(1) + xaxis(2)*zaxis(2) &
1092 & + xaxis(3)*zaxis(3))/(xmod*zmod)
1093 
1094  if(abs(cosine) > 1.d-8)then
1095    write(message,'(a,a,a,i6)')&
1096 &   'xaxis and zaxis should be perpendicular,',ch10,&
1097 &   'however, cosine=',cosine
1098    MSG_BUG(message)
1099  end if
1100 
1101 !new y axis as cross product
1102  yaxis(1) = (zaxis(2)*xaxis(3) - xaxis(2)*zaxis(3))/(xmod*zmod)
1103  yaxis(2) = (zaxis(3)*xaxis(1) - xaxis(3)*zaxis(1))/(xmod*zmod)
1104  yaxis(3) = (zaxis(1)*xaxis(2) - xaxis(1)*zaxis(2))/(xmod*zmod)
1105 
1106 !hack to allow inversion operation on coordinate transformation
1107 !uses unlikely large but legal values of proj_x and/or proj_z
1108 !to flag inversion
1109  inversion_flag=0
1110  if(xmod>10._dp .or. zmod>10._dp) then
1111    inversion_flag=1
1112    write(message, '(4a)' )&
1113 &   'inversion operation will be appended to axis transformation',ch10,&
1114 &   'Action: If you did not intend this, make |z|<10 and |x|<10 ',ch10
1115    call wrtout(std_out,message,'COLL')
1116  end if
1117 
1118  umat(1,:) = xaxis(:)/xmod
1119  umat(2,:) = yaxis(:)
1120  umat(3,:) = zaxis(:)/zmod
1121 
1122 end subroutine rotmat

m_geometry/shellstruct [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  shellstruct

FUNCTION

  Calculates shell structure (multiplicities, radii) of an atomic configuration

INPUTS

  natom=number of atoms in unit cell
  xred=reduced coordinates of atoms
  rprimd=unit cell vectors
  magv = magnetic ordering of atoms given as 1 and -1, if not given fm is assumed
  atp = atom on which the perturbation was done

OUTPUT

  sdisv(nat)= distance of each shell to central atom (only the first nsh entries are relevant)
  nsh= number of shells
  mult(nat) = number of atoms on shell (only the first nsh entries are relevant)

PARENTS

      pawuj_det

CHILDREN

      ioniondist,prmat,sort_dp,sort_int,wrtout

SOURCE

2585 subroutine shellstruct(xred,rprimd,natom,magv,distv,smult,sdisv,nsh,atp,prtvol)
2586 
2587 
2588 !This section has been created automatically by the script Abilint (TD).
2589 !Do not modify the following lines by hand.
2590 #undef ABI_FUNC
2591 #define ABI_FUNC 'shellstruct'
2592 !End of the abilint section
2593 
2594  implicit none
2595 
2596 !Arguments ------------------------------------
2597 !scalars
2598  integer,intent(in)              :: natom
2599  integer,intent(in),optional     :: atp
2600  integer,intent(in),optional     :: prtvol
2601  integer,intent(out)             :: nsh
2602 !arrays
2603  real(dp),intent(in)             :: rprimd(3,3)
2604  real(dp),intent(in)             :: xred(3,natom)
2605  integer,intent(out)             :: smult(natom)
2606  integer,intent(in),optional     :: magv(natom)
2607  real(dp),intent(out)            :: sdisv(natom)
2608  real(dp),intent(out)            :: distv(natom)
2609 
2610 !Local variables-------------------------------
2611 !scalars
2612  integer                      :: iatom,atpp,ish,prtvoll
2613  character(len=500)           :: message
2614  real(dp),parameter           :: rndfact=10000_dp
2615 !arrays
2616  integer                      :: iperm(natom),jperm(natom)
2617  real(dp)                     :: distvh(natom,natom)
2618  real(dp)                     :: magvv(natom)
2619 
2620 ! *************************************************************************
2621 
2622  if (present(magv)) then
2623    magvv=magv
2624  else
2625    magvv=(/ (1, iatom=1,natom)  /)
2626  end if
2627 
2628  if (present(atp)) then
2629    atpp=atp
2630  else
2631    atpp=1
2632  end if
2633 
2634  if (present(prtvol)) then
2635    prtvoll=prtvol
2636  else
2637    prtvoll=1
2638  end if
2639 
2640 !DEBUB
2641  write(std_out,*)'shellstruct start'
2642 !END DEBUG
2643 
2644 !Calculate ionic distances
2645  call ioniondist(natom,rprimd,xred,distvh,1,magv=int(magvv),atp=atpp)
2646  distv=distvh(1,:)
2647 
2648  if (prtvol>2) then
2649    write(std_out,'(a)')' shellstruct ionic distances in cell (distv) : '
2650    call prmat(distv(1:natom),1,natom,1,std_out)
2651  end if
2652 
2653  iperm=(/ (iatom, iatom=1,natom ) /)
2654  jperm=iperm
2655  distv=anint(distv*rndfact)/rndfact
2656 !Sort distances
2657  call sort_dp(natom,distv,iperm,10d-5)
2658  call sort_int(natom,iperm,jperm)
2659 
2660  smult=0
2661  sdisv=dot_product(rprimd(1,:),rprimd(1,:))+dot_product(rprimd(2,:),rprimd(2,:))+dot_product(rprimd(3,:),rprimd(3,:))
2662 
2663  nsh=1
2664  smult(1)=1
2665  sdisv(1)=distv(1)
2666 
2667  do iatom=2,natom
2668    do ish=1,natom
2669      if (distv(iatom)>sdisv(ish)) then
2670        cycle
2671      else if (distv(iatom)==sdisv(ish)) then
2672        smult(ish)=smult(ish)+1
2673        exit
2674      else if (distv(iatom)<sdisv(ish)) then
2675        smult(ish+1:natom)=smult(ish:natom-1)
2676        sdisv(ish+1:natom)=sdisv(ish:natom-1)
2677        smult(ish)=1
2678        sdisv(ish)=distv(iatom)
2679        nsh=nsh+1
2680        exit
2681      end if
2682    end do
2683  end do
2684 
2685  distv=(/ ( distv(jperm(iatom)),iatom=1,natom ) /)
2686 
2687  if (prtvoll>2) then
2688    write(message,'(a,i4,a)')' shellstruct found ',nsh,' shells at distances (sdisv) '
2689    call wrtout(std_out,message,'COLL')
2690    call prmat(sdisv(1:nsh),1,nsh,1,std_out)
2691    write(message,fmt='(a,150i4)')' and multiplicities (smult) ', smult(1:nsh)
2692    call wrtout(std_out,message,'COLL')
2693  end if
2694 
2695 !DEBUB
2696  write(std_out,*)'shellstruct leave'
2697 !END DEBUG
2698 
2699 end subroutine shellstruct

m_geometry/spinrot_cmat [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

  spinrot_cmat

FUNCTION

  Construct 2x2 complex matrix representing rotation operator in spin-space.

INPUTS

  spinrot(4)=components of the spinor rotation matrix computed by getspinrot

OUTPUT

  spinrot(2,2)=Rotation matrix (complex array)

PARENTS

CHILDREN

SOURCE

 981 pure function spinrot_cmat(spinrot)
 982 
 983 
 984 !This section has been created automatically by the script Abilint (TD).
 985 !Do not modify the following lines by hand.
 986 #undef ABI_FUNC
 987 #define ABI_FUNC 'spinrot_cmat'
 988 !End of the abilint section
 989 
 990  implicit none
 991 
 992 !Arguments ------------------------------------
 993  real(dp),intent(in) :: spinrot(4)
 994  complex(dpc) :: spinrot_cmat(2,2)
 995 
 996 ! *************************************************************************
 997 
 998  ! Rotation in spinor space (same equations as in wfconv)
 999  spinrot_cmat(1,1) = spinrot(1) + j_dpc*spinrot(4)
1000  spinrot_cmat(1,2) = spinrot(3) + j_dpc*spinrot(2)
1001  spinrot_cmat(2,1) =-spinrot(3) + j_dpc*spinrot(2)
1002  spinrot_cmat(2,2) = spinrot(1) - j_dpc*spinrot(4)
1003 
1004  ! My equation
1005  !spinrot_cmat(1,1) = spinrot(1) - j_dpc*spinrot(4)
1006  !spinrot_cmat(1,2) =-spinrot(3) - j_dpc*spinrot(2)
1007  !spinrot_cmat(2,1) = spinrot(3) - j_dpc*spinrot(2)
1008  !spinrot_cmat(2,2) = spinrot(1) + j_dpc*spinrot(4)
1009 
1010 end function spinrot_cmat

m_geometry/strainsym [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 strainsym

FUNCTION

 For given order of point group, symmetrizes the strain tensor,
 then produce primitive vectors based on the symmetrized strain.

INPUTS

 nsym=order of group.
 rprimd(3,3)= primitive vectors, to be symmetrized
 rprimd0(3,3)= reference primitive vectors, already symmetrized
 symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations

OUTPUT

 rprimd_symm(3,3)= symmetrized primitive vectors

PARENTS

      xfpack_vin2x,xfpack_x2vin

CHILDREN

      dgemm,mati3inv,matrginv

SOURCE

3200 subroutine strainsym(nsym,rprimd0,rprimd,rprimd_symm,symrel)
3201 
3202  use m_linalg_interfaces
3203 
3204 !This section has been created automatically by the script Abilint (TD).
3205 !Do not modify the following lines by hand.
3206 #undef ABI_FUNC
3207 #define ABI_FUNC 'strainsym'
3208 !End of the abilint section
3209 
3210  implicit none
3211 
3212 !Arguments ------------------------------------
3213 !scalars
3214  integer,intent(in) :: nsym
3215 !arrays
3216  integer,intent(in) :: symrel(3,3,nsym)
3217  real(dp),intent(in) :: rprimd(3,3),rprimd0(3,3)
3218  real(dp),intent(out) :: rprimd_symm(3,3)
3219 
3220 !Local variables-------------------------------
3221 !scalars
3222  integer :: isym
3223 !arrays
3224  integer :: symrel_it(3,3)
3225  real(dp) :: rprimd0_inv(3,3),strain(3,3),strain_symm(3,3),tmp_mat(3,3)
3226 
3227 !**************************************************************************
3228 
3229 !copy initial rprimd input and construct inverse
3230  rprimd0_inv = rprimd0
3231  call matrginv(rprimd0_inv,3,3)
3232 
3233 !define strain as rprimd = strain * rprimd0 (in cartesian frame)
3234 !so strain = rprimd * rprimd0^{-1}
3235 !transform to triclinic frame with rprimd0^{-1} * strain * rprimd0
3236 !giving strain as rprimd0^{-1} * rprimd
3237  call dgemm('N','N',3,3,3,one,rprimd0_inv,3,rprimd,3,zero,strain,3)
3238 
3239 !loop over symmetry elements to obtain symmetrized strain matrix
3240  strain_symm = zero
3241  do isym = 1, nsym
3242 
3243 !  this loop accumulates symrel^{-1}*strain*symrel into strain_symm
3244 
3245 !  mati3inv gives the inverse transpose of symrel
3246    call mati3inv(symrel(:,:,isym),symrel_it)
3247    call dgemm('N','N',3,3,3,one,strain,3,dble(symrel(:,:,isym)),3,zero,tmp_mat,3)
3248    call dgemm('T','N',3,3,3,one,dble(symrel_it),3,tmp_mat,3,one,strain_symm,3)
3249 
3250  end do
3251 
3252 !normalize by number of symmetry operations
3253  strain_symm = strain_symm/dble(nsym)
3254 
3255 !this step is equivalent to r_new = r_old * strain * r_old^{-1} * r_old,
3256 !that is, convert strain back to cartesian frame and then multipy by r_old,
3257 !to get the r_new primitive vectors
3258 
3259  call dgemm('N','N',3,3,3,one,rprimd0,3,strain_symm,3,zero,rprimd_symm,3)
3260 
3261 end subroutine strainsym

m_geometry/strconv [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 strconv

FUNCTION

 If original gprimd is input, convert from symmetric storage mode
 3x3 tensor in reduced coordinates "frac" to symmetric storage mode
 symmetric tensor in cartesian coordinates "cart".

INPUTS

  frac(6)=3x3 tensor in symmetric storage mode, reduced coordinates
  gprimd(3,3)=reciprocal space dimensional primitive translations (bohr^-1)

OUTPUT

  cart(6)=symmetric storage mode for symmetric 3x3 tensor in cartesian coords.

NOTES

 $cart(i,j)=G(i,a) G(j,b) frac(a,b)$
 "Symmetric" storage mode for 3x3 tensor is 6 element array with
 elements 11, 22, 33, 32, 31, and 21.
 "cart" may be same array as "frac".
 If rprimd transpose is input instead of gprimd, then convert tensor
 in cartesian coordinates to reduced coordinates

PARENTS

      ctocprj,d2frnl,mkcore,mkcore_paw,mkcore_wvl,nonlop_pl,nonlop_ylm
      stresssym

CHILDREN

SOURCE

3409 subroutine strconv(frac,gprimd,cart)
3410 
3411 
3412 !This section has been created automatically by the script Abilint (TD).
3413 !Do not modify the following lines by hand.
3414 #undef ABI_FUNC
3415 #define ABI_FUNC 'strconv'
3416 !End of the abilint section
3417 
3418  implicit none
3419 
3420 !Arguments ------------------------------------
3421 !arrays
3422  real(dp),intent(in) :: frac(6),gprimd(3,3)
3423  real(dp),intent(inout) :: cart(6) ! alias of frac   !vz_i
3424 
3425 !Local variables-------------------------------
3426 !scalars
3427  integer :: ii,jj
3428 !arrays
3429  real(dp) :: work1(3,3),work2(3,3)
3430 
3431 ! *************************************************************************
3432 
3433  work1(1,1)=frac(1)
3434  work1(2,2)=frac(2)
3435  work1(3,3)=frac(3)
3436  work1(3,2)=frac(4) ; work1(2,3)=frac(4)
3437  work1(3,1)=frac(5) ; work1(1,3)=frac(5)
3438  work1(2,1)=frac(6) ; work1(1,2)=frac(6)
3439 
3440  do ii=1,3
3441    work2(:,ii)=zero
3442    do jj=1,3
3443      work2(:,ii)=work2(:,ii)+gprimd(ii,jj)*work1(:,jj)
3444    end do
3445  end do
3446 
3447  do ii=1,3
3448    work1(ii,:)=zero
3449    do jj=1,3
3450      work1(ii,:)=work1(ii,:)+gprimd(ii,jj)*work2(jj,:)
3451    end do
3452  end do
3453 
3454  cart(1)=work1(1,1)
3455  cart(2)=work1(2,2)
3456  cart(3)=work1(3,3)
3457  cart(4)=work1(2,3)
3458  cart(5)=work1(1,3)
3459  cart(6)=work1(1,2)
3460 
3461 end subroutine strconv

m_geometry/stresssym [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 stresssym

FUNCTION

 For given order of point group, symmetrizes the stress tensor,
 in symmetrized storage mode and cartesian coordinates, using input
 3x3 symmetry operators in reduced coordinates.
 symmetrized tensor replaces input tensor.

INPUTS

 gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
 nsym=order of group.
 sym(3,3,nsym)=symmetry operators (usually symrec=expressed in terms
               of action on reciprocal lattice primitive translations);
               integers.

OUTPUT

 stress(6)=stress tensor, in cartesian coordinates, in symmetric storage mode

SIDE EFFECTS

PARENTS

      dfpt_nselt,dfpt_nstpaw,forstrnps,littlegroup_pert,pawgrnl,stress

CHILDREN

      matr3inv,strconv

SOURCE

3294 subroutine stresssym(gprimd,nsym,stress,sym)
3295 
3296 
3297 !This section has been created automatically by the script Abilint (TD).
3298 !Do not modify the following lines by hand.
3299 #undef ABI_FUNC
3300 #define ABI_FUNC 'stresssym'
3301 !End of the abilint section
3302 
3303  implicit none
3304 
3305 !Arguments ------------------------------------
3306 !scalars
3307  integer,intent(in) :: nsym
3308 !arrays
3309  integer,intent(in) :: sym(3,3,nsym)
3310  real(dp),intent(in) :: gprimd(3,3)
3311  real(dp),intent(inout) :: stress(6)
3312 
3313 !Local variables-------------------------------
3314 !scalars
3315  integer :: ii,isym,mu,nu
3316  real(dp) :: summ,tmp
3317 !arrays
3318  real(dp) :: rprimd(3,3),rprimdt(3,3),strfrac(6),tensor(3,3),tt(3,3)
3319 
3320 !*************************************************************************
3321 
3322 !Obtain matrix of real space dimensional primitive translations
3323 !(inverse tranpose of gprimd), and its transpose
3324  call matr3inv(gprimd,rprimd)
3325  rprimdt=transpose(rprimd)
3326 
3327 !Compute stress tensor in reduced coordinates
3328  call strconv(stress,rprimdt,strfrac)
3329 
3330 !Switch to full storage mode
3331  tensor(1,1)=strfrac(1)
3332  tensor(2,2)=strfrac(2)
3333  tensor(3,3)=strfrac(3)
3334  tensor(3,2)=strfrac(4)
3335  tensor(3,1)=strfrac(5)
3336  tensor(2,1)=strfrac(6)
3337  tensor(2,3)=tensor(3,2)
3338  tensor(1,3)=tensor(3,1)
3339  tensor(1,2)=tensor(2,1)
3340 
3341  do nu=1,3
3342    do mu=1,3
3343      tt(mu,nu)=tensor(mu,nu)/dble(nsym)
3344      tensor(mu,nu)=0.0_dp
3345    end do
3346  end do
3347 
3348 !loop over all symmetry operations:
3349  do isym=1,nsym
3350    do mu=1,3
3351      do nu=1,3
3352        summ=0._dp
3353        do ii=1,3
3354          tmp=tt(ii,1)*sym(nu,1,isym)+tt(ii,2)*sym(nu,2,isym)+&
3355 &         tt(ii,3)*sym(nu,3,isym)
3356          summ=summ+sym(mu,ii,isym)*tmp
3357        end do
3358        tensor(mu,nu)=tensor(mu,nu)+summ
3359      end do
3360    end do
3361  end do
3362 
3363 !Switch back to symmetric storage mode
3364  strfrac(1)=tensor(1,1)
3365  strfrac(2)=tensor(2,2)
3366  strfrac(3)=tensor(3,3)
3367  strfrac(4)=tensor(3,2)
3368  strfrac(5)=tensor(3,1)
3369  strfrac(6)=tensor(2,1)
3370 
3371 !Convert back stress tensor (symmetrized) in cartesian coordinates
3372  call strconv(strfrac,gprimd,stress)
3373 
3374 end subroutine stresssym

m_geometry/symredcart [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 symredcart

FUNCTION

 Convert a symmetry operation from reduced coordinates (integers)
 to cartesian coordinates (reals). Can operate in real or reciprocal space

INPUTS

 symred(3,3)=symmetry matrice in reduced coordinates (integers) (real or reciprocal space)
 aprim(3,3)=real or reciprocal space dimensional primitive translations (see below)
 bprim(3,3)=real or reciprocal space dimensional primitive translations (see below)

OUTPUT

 symcart(3,3)=symmetry matrice in cartesian coordinates (reals)

NOTES

 When aprim=rprimd and bprim=gprimd, the routine operates in real space (on a real space symmetry)
 When aprim=gprimd and bprim=rprimd, the routine operates in reciprocal space (on a real space symmetry)

PARENTS

      m_matlu,m_phonons,symrhg

CHILDREN

SOURCE

3126 subroutine symredcart(aprim,bprim,symcart,symred)
3127 
3128 
3129 !This section has been created automatically by the script Abilint (TD).
3130 !Do not modify the following lines by hand.
3131 #undef ABI_FUNC
3132 #define ABI_FUNC 'symredcart'
3133 !End of the abilint section
3134 
3135  implicit none
3136 
3137 !Arguments ------------------------------------
3138 !arrays
3139  integer,intent(in) :: symred(3,3)
3140  real(dp),intent(in) :: aprim(3,3),bprim(3,3)
3141  real(dp),intent(out) :: symcart(3,3)
3142 
3143 !Local variables-------------------------------
3144 !scalars
3145  integer :: ii,jj,kk
3146  real(dp) :: symtmp
3147 !arrays
3148  real(dp) :: work(3,3)
3149 
3150 ! *************************************************************************
3151 
3152  work=zero
3153  do kk=1,3
3154    do jj=1,3
3155      symtmp=dble(symred(jj,kk))
3156      do ii=1,3
3157        work(ii,jj)=work(ii,jj)+bprim(ii,kk)*symtmp
3158      end do
3159    end do
3160  end do
3161 
3162  symcart=zero
3163  do kk=1,3
3164    do jj=1,3
3165      symtmp=work(jj,kk)
3166      do ii=1,3
3167        symcart(ii,jj)=symcart(ii,jj)+aprim(ii,kk)*symtmp
3168      end do
3169    end do
3170  end do
3171 
3172 end subroutine symredcart

m_geometry/vdotw_rc_vector [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 vdotw_rc_vector

FUNCTION

 Compute the scalar product between two vectors expressed in reduced coordinates
 First vector is real, the second one is complex.
 The result is multiplied by (2pi)**2 in case of vectors in reciprocal space
 to take into account the correct normalisation of the reciprocal lattice vectors

INPUTS

  xv(3),xw(3)=Vectors in reduced coordinates
  met(3,3)=Metric tensor
  space=Character defining whether we are working in real (r) or reciprocal space (g)

OUTPUT

  res=complex scalar product of xv and xw

PARENTS

CHILDREN

SOURCE

404 function vdotw_rc_vector(xv,xw,met,space) result(res)
405 
406 
407 !This section has been created automatically by the script Abilint (TD).
408 !Do not modify the following lines by hand.
409 #undef ABI_FUNC
410 #define ABI_FUNC 'vdotw_rc_vector'
411 !End of the abilint section
412 
413  implicit none
414 
415 !Arguments ------------------------------------
416 !scalars
417  complex(dpc) :: res
418  character(len=1),intent(in) :: space
419 !arrays
420  real(dp),intent(in) :: met(3,3),xv(3)
421  complex(dpc),intent(in) :: xw(3)
422 
423 ! *************************************************************************
424 
425  res = (  met(1,1)* xv(1)*xw(1)                &
426 &        +met(2,2)* xv(2)*xw(2)                &
427 &        +met(3,3)* xv(3)*xw(3)                &
428 &        +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) &
429 &        +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) &
430 &        +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) )
431 
432  select case (space)
433  case ('r','R')
434    return
435  case ('g','G')
436    res= res * (two_pi**2)
437  case default
438    MSG_BUG('Wrong value for space')
439  end select
440 
441 end function vdotw_rc_vector

m_geometry/vdotw_rr_vector [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 vdotw_rr_vector

FUNCTION

 Compute the scalar product between two vectors expressed in reduced coordinates
 The result is multiplied by (2pi)**2 in case of vectors in reciprocal space
 to take into account the correct normalisation of the reciprocal lattice vectors

INPUTS

  xv(3),xw(3)=Vectors in reduced coordinates
  met(3,3)=Metric tensor
  space=Character defining whether we are working in real (r) or reciprocal space (g)

OUTPUT

  res=scalar product of xv and xw

PARENTS

CHILDREN

SOURCE

339 function vdotw_rr_vector(xv,xw,met,space) result(res)
340 
341 
342 !This section has been created automatically by the script Abilint (TD).
343 !Do not modify the following lines by hand.
344 #undef ABI_FUNC
345 #define ABI_FUNC 'vdotw_rr_vector'
346 !End of the abilint section
347 
348  implicit none
349 
350 !Arguments ------------------------------------
351 !scalars
352  real(dp) :: res
353  character(len=1),intent(in) :: space
354 !arrays
355  real(dp),intent(in) :: met(3,3),xv(3),xw(3)
356 
357 ! *************************************************************************
358 
359  res = (  met(1,1)* xv(1)*xw(1)                &
360 &        +met(2,2)* xv(2)*xw(2)                &
361 &        +met(3,3)* xv(3)*xw(3)                &
362 &        +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) &
363 &        +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) &
364 &        +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) )
365 
366  select case (space)
367  case ('r','R')
368    return
369  case ('g','G')
370    res= res * (two_pi**2)
371  case default
372    MSG_BUG('Wrong value for space')
373  end select
374 
375 end function vdotw_rr_vector

m_geometry/wigner_seitz [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 wigner_seitz

FUNCTION

 Calculates a grid of points that falls inside of (and eventually on the surface of)
 the Wigner-Seitz supercell centered on the origin of the B lattice with primitive
 translations nmonkh(1)*a_1+nmonkh(2)*a_2+nmonkh(3)*a_3.
 Subroutine taken from the Wannier90 code. Modified by MG to fulfil abinit coding rules.
 API slightly changed wrt the wannier90 version.

COPYRIGHT

 Copyright (C) 2007 Jonathan Yates, Arash Mostofi,
 Young-Su Lee, Nicola Marzari, Ivo Souza, David Vanderbilt.
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  center(3)=The Wigner-Seitz cell is centered on this point in reduced coordinates.
  rmet(3,3)=Real space metric ($\textrm{bohr}^{2}$).
  kptrlatt(3)=Values defining the supercell.
  prtvol=If different from 0 print out the points falling inside the W-S cell and the correponding weights.
  lmax(3)=see Notes below.

OUTPUT

  npts=number of points falling inside the Wigner-Seitz cell
  irvec(3,npts)=Reduced coordinated of the points inside the W-S cell
  ndegen(npts)=Weigths associated to each point.

SIDE EFFECTS

  irvec and ndegen are are allocated with the correct
  size inside the routine and returned to the caller.

NOTES

 The Wannier functions live in a supercell of the real space unit cell.
 This supercell is mp_grid unit cells long in each direction
 The algorithm loops over grid points r on a unit cell that is 8 times larger than this
 primitive supercell.
 One of these points is in the W-S cell if it is closer to center(:)
 than any of the other points R where R are the translation vectors of the supercell.
 In the end npts contains the total number of grid points that have been found in the Wigner-Seitz cell
 The number of lattice vectors R along each direction of the supercell is defined by lmax.

PARENTS

CHILDREN

SOURCE

542 subroutine wigner_seitz(center, lmax, kptrlatt, rmet, npts, irvec, ndegen, prtvol)
543 
544 
545 !This section has been created automatically by the script Abilint (TD).
546 !Do not modify the following lines by hand.
547 #undef ABI_FUNC
548 #define ABI_FUNC 'wigner_seitz'
549 !End of the abilint section
550 
551  implicit none
552 
553 !Arguments ------------------------------------
554 !scalars
555  integer,optional,intent(in) :: prtvol
556  integer,intent(out) :: npts
557 !arrays
558  integer,intent(in) :: kptrlatt(3,3),lmax(3)
559  integer,allocatable,intent(out) :: irvec(:,:),ndegen(:)
560  real(dp),intent(in) :: center(3),rmet(3,3)
561 
562 !Local variables-------------------------------
563 !scalars
564  integer :: in1,in2,in3,l1,l2,l3,ii,icount,n1,n2,n3
565  integer :: l0,l1_max,l2_max,l3_max,nl,verbose,mm1,mm2,mm3
566  real(dp) :: tot,dist_min
567  real(dp),parameter :: TOL_DIST=tol7
568  character(len=500) :: msg
569 !arrays
570  real(dp) :: diff(3)
571  real(dp),allocatable :: dist(:),swap2(:,:),swap1(:)
572 
573 ! *************************************************************************
574 
575  verbose = 0; if (PRESENT(prtvol)) verbose = prtvol
576 
577  if (kptrlatt(1,2)/=0 .or. kptrlatt(2,1)/=0 .or. &
578 &    kptrlatt(1,3)/=0 .or. kptrlatt(3,1)/=0 .or. &
579 &    kptrlatt(2,3)/=0 .or. kptrlatt(3,2)/=0 ) then
580    MSG_ERROR('Off-diagonal elements of kptrlatt must be zero')
581  end if
582 
583  n1 = kptrlatt(1,1)
584  n2 = kptrlatt(2,2)
585  n3 = kptrlatt(3,3)
586 
587  l1_max = lmax(1)
588  l2_max = lmax(2)
589  l3_max = lmax(3)
590 
591  nl=(2*l1_max+1)*(2*l2_max+1)*(2*l3_max+1)
592  l0=1+l1_max*(1+(2*l2_max+1)**2+(2*l3_max+1)) ! Index of the origin.
593  ABI_MALLOC(dist,(nl))
594 
595  ! Allocate with maximum size
596  mm1 = 2 * n1 + 1
597  mm2 = 2 * n2 + 1
598  mm3 = 2 * n3 + 1
599  ABI_MALLOC(irvec, (3, mm1*mm2*mm3))
600  ABI_MALLOC(ndegen, (mm1*mm2*mm3))
601 
602  npts = 0
603  do in1=-n1,n1
604    do in2=-n2,n2
605      do in3=-n3,n3
606 
607       ! Loop over the nl points R. R=0 corresponds to l1=l2=l3=1, or icount=l0
608       icount = 0
609       do l1=-l1_max,l1_max
610         do l2=-l2_max,l2_max
611           do l3=-l3_max,l3_max
612             ! Calculate |r - R -r0|^2.
613             diff(1) = in1 - l1 * n1 - center(1)
614             diff(2) = in2 - l2 * n2 - center(2)
615             diff(3) = in3 - l3 * n3 - center(3)
616             icount = icount+1
617             dist(icount) = DOT_PRODUCT(diff, MATMUL(rmet, diff))
618           end do
619         end do
620       end do
621 
622       dist_min = MINVAL(dist)
623 
624       if (ABS(dist(l0) - dist_min) < TOL_DIST) then
625         npts = npts + 1
626         ndegen (npts) = 0
627         do ii=1,nl
628           if (ABS(dist(ii) - dist_min) < TOL_DIST) ndegen(npts) = ndegen(npts) + 1
629         end do
630         irvec(1, npts) = in1
631         irvec(2, npts) = in2
632         irvec(3, npts) = in3
633       end if
634      end do !in3
635    end do !in2
636  end do !in1
637 
638  if (verbose>=1) then
639    write(msg,'(a,i0)')' lattice points in Wigner-Seitz supercell: ',npts
640    call wrtout(std_out,msg,'COLL')
641    do ii=1,npts
642      write(msg,'(a,3(i3),a,i4)')'  vector ', irvec(:,ii),' degeneracy: ', ndegen(ii)
643      call wrtout(std_out,msg,'COLL')
644    end do
645  end if
646 
647  ! Check the "sum rule"
648  tot = zero
649  do ii=1,npts
650    tot = tot + one/ndegen(ii)
651  end do
652  if (ABS(tot-(n1*n2*n3))>tol8) then
653    write(msg,'(a,es16.8,a,i0)')'Something wrong in the generation of WS mesh: tot ',tot,' /= ',n1*n2*n3
654    MSG_ERROR(msg)
655  end if
656 
657  ABI_FREE(dist)
658 
659  ! Reallocate the output with correct size.
660  ABI_MALLOC(swap1,(npts))
661  swap1(:)=ndegen(1:npts)
662  ABI_FREE(ndegen)
663  ABI_MALLOC(ndegen,(npts))
664  ndegen=swap1
665  ABI_FREE(swap1)
666 
667  ABI_MALLOC(swap2,(3,npts))
668  swap2(:,:)=irvec(1:3,1:npts)
669  ABI_FREE(irvec)
670  ABI_MALLOC(irvec,(3,npts))
671  irvec=swap2
672  ABI_FREE(swap2)
673 
674 end subroutine wigner_seitz

m_geometry/xcart2xred [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 xcart2xred

FUNCTION

 Convert from cartesian coordinates xcart(3,natom) in bohr to
 dimensionless reduced coordinates xred(3,natom) by using
 xred(mu,ia)=gprimd(1,mu)*xcart(1,ia)
            +gprimd(2,mu)*xcart(2,ia)
            +gprimd(3,mu)*xcart(3,ia)
 where gprimd is the inverse of rprimd
 Note that the reverse operation is deon by xred2xcart

INPUTS

  natom=number of atoms in unit cell
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  xcart(3,natom)=cartesian coordinates of atoms (bohr)

OUTPUT

  xred(3,natom)=dimensionless reduced coordinates of atoms

PARENTS

      driver,evdw_wannier,ingeo,m_cut3d,m_dens,m_effective_potential
      m_effective_potential_file,m_mep,m_paw_pwaves_lmn,m_pred_lotf
      mkcore_paw,mkcore_wvl,mover_effpot,pawmkaewf,pimd_langevin_npt
      pimd_langevin_nvt,pimd_nosehoover_npt,pimd_nosehoover_nvt,prcref
      prcref_PMA,pred_delocint,pred_diisrelax,pred_isokinetic,pred_isothermal
      pred_langevin,pred_moldyn,pred_nose,pred_srkna14,pred_steepdesc
      pred_velverlet,pred_verlet,relaxpol,wrt_moldyn_netcdf
      wvl_setboxgeometry

CHILDREN

      matr3inv

SOURCE

1703 subroutine xcart2xred(natom,rprimd,xcart,xred)
1704 
1705 
1706 !This section has been created automatically by the script Abilint (TD).
1707 !Do not modify the following lines by hand.
1708 #undef ABI_FUNC
1709 #define ABI_FUNC 'xcart2xred'
1710 !End of the abilint section
1711 
1712  implicit none
1713 
1714 !Arguments ------------------------------------
1715 !scalars
1716  integer,intent(in) :: natom
1717 !arrays
1718  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1719  real(dp),intent(out) :: xred(3,natom)
1720 
1721 !Local variables-------------------------------
1722 !scalars
1723  integer :: iatom,mu
1724 !arrays
1725  real(dp) :: gprimd(3,3)
1726 
1727 ! *************************************************************************
1728 
1729  call matr3inv(rprimd,gprimd)
1730  do iatom=1,natom
1731    do mu=1,3
1732      xred(mu,iatom)= gprimd(1,mu)*xcart(1,iatom)+gprimd(2,mu)*xcart(2,iatom)+&
1733 &     gprimd(3,mu)*xcart(3,iatom)
1734    end do
1735  end do
1736 
1737 end subroutine xcart2xred

m_geometry/xred2xcart [ Functions ]

[ Top ] [ m_geometry ] [ Functions ]

NAME

 xred2xcart

FUNCTION

 Convert from dimensionless reduced coordinates xred(3,natom)
 to cartesian coordinates xcart(3,natom) in bohr by using
 xcart(mu,ia)=rprimd(mu,1)*xred(1,ia)
             +rprimd(mu,2)*xred(2,ia)
             +rprimd(mu,3)*xred(3,ia)
 Note that the reverse operation is done by xcart2xred.F90

INPUTS

  natom=number of atoms in unit cell
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  xred(3,natom)=dimensionless reduced coordinates of atoms

OUTPUT

  xcart(3,natom)=cartesian coordinates of atoms (bohr)

PARENTS

      afterscfloop,berryphase,berryphase_new,bonds_lgth_angles,constrf,cut3d
      denfgr,driver,evdw_wannier,forstr,ingeo,ionion_realspace,ionion_surface
      m_abihist,m_crystal,m_ddb,m_effective_potential,m_fit_polynomial_coeff
      m_mep,m_pred_lotf,m_results_img,m_tdep_abitypes,make_efg_el
      make_efg_ion,mkcore_paw,mkcore_wvl,mkgrid_fft,mklocl,mklocl_realspace
      mlwfovlp_projpaw,mover_effpot,out1dm,outqmc,outvar_o_z,outxml
      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_npt
      pimd_nosehoover_nvt,prec_simple,pred_delocint,pred_diisrelax,pred_hmc
      pred_isokinetic,pred_isothermal,pred_langevin,pred_moldyn,pred_nose
      pred_srkna14,pred_steepdesc,pred_velverlet,pred_verlet,prtimg
      prtspgroup,prtxfase,randomcellpos,rhotov,setvtr,spin_current,symspgr
      thmeig,vso_realspace_local,vtorho,wrt_moldyn_netcdf,wvl_denspot_set
      wvl_initro,wvl_memory,wvl_nhatgrid,wvl_projectors_set,wvl_rwwf
      wvl_setboxgeometry,wvl_wfs_set,wvl_wfsinp_reformat,wvl_wfsinp_scratch
      xfh_recover_deloc

CHILDREN

SOURCE

1781 subroutine xred2xcart(natom,rprimd,xcart,xred)
1782 
1783 
1784 !This section has been created automatically by the script Abilint (TD).
1785 !Do not modify the following lines by hand.
1786 #undef ABI_FUNC
1787 #define ABI_FUNC 'xred2xcart'
1788 !End of the abilint section
1789 
1790  implicit none
1791 
1792 !Arguments ------------------------------------
1793 !scalars
1794  integer,intent(in) :: natom
1795 !arrays
1796  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1797  real(dp),intent(out) :: xcart(3,natom)
1798 
1799 !Local variables-------------------------------
1800 !scalars
1801  integer :: iatom,mu
1802 
1803 ! *************************************************************************
1804 
1805  do iatom=1,natom
1806    do mu=1,3
1807      xcart(mu,iatom)=rprimd(mu,1)*xred(1,iatom)+rprimd(mu,2)*xred(2,iatom)+rprimd(mu,3)*xred(3,iatom)
1808    end do
1809  end do
1810 
1811 end subroutine xred2xcart