## 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 (C) 2008-2018 ABINIT group (MG, MT, FJ, TRangel, DCA, XG, AHR, DJA, DRH)
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
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)
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


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


## 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)
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
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)))
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)))
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)))
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)))
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
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
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
3056          if (determinant(is)==1) then
3058          end if
3059        else if (ALL(tnons(:,is2)<tol8)) then
3061        else
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 (C) 2007 Jonathan Yates, Arash Mostofi,
Young-Su Lee, Nicola Marzari, Ivo Souza, David Vanderbilt.
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